30 f_y, f_z, c_Xh, gs_Xh, bc_prs_surface,bc_sym_surface, Ax, bd, dt, mu, rho)
31 type(
field_t),
intent(inout) :: p, u, v, w
32 type(
field_t),
intent(inout) :: u_e, v_e, w_e
33 type(
field_t),
intent(inout) :: p_res
34 type(
field_t),
intent(inout) :: f_x, f_y, f_z
35 type(
coef_t),
intent(inout) :: c_Xh
36 type(
gs_t),
intent(inout) :: gs_Xh
39 class(
ax_t),
intent(inout) :: Ax
40 real(kind=
rp),
intent(inout) :: bd
41 real(kind=
rp),
intent(in) :: dt
42 real(kind=
rp),
intent(in) :: mu
43 real(kind=
rp),
intent(in) :: rho
47 type(
field_t),
pointer :: ta1, ta2, ta3, wa1, wa2, wa3, work1, work2
48 integer :: temp_indices(8)
62 c_xh%h1(i,1,1,1) = 1.0_rp / rho
63 c_xh%h2(i,1,1,1) = 0.0_rp
67 call curl(ta1, ta2, ta3, u_e, v_e, w_e, work1, work2, c_xh)
68 call curl(wa1, wa2, wa3, ta1, ta2, ta3, work1, work2, c_xh)
70 do concurrent(i = 1:n)
71 ta1%x(i,1,1,1) = f_x%x(i,1,1,1) / rho &
72 - ((wa1%x(i,1,1,1) * (mu / rho)) * c_xh%B(i,1,1,1))
73 ta2%x(i,1,1,1) = f_y%x(i,1,1,1) / rho &
74 - ((wa2%x(i,1,1,1) * (mu / rho)) * c_xh%B(i,1,1,1))
75 ta3%x(i,1,1,1) = f_z%x(i,1,1,1) / rho &
76 - ((wa3%x(i,1,1,1) * (mu / rho)) * c_xh%B(i,1,1,1))
79 call gs_xh%op(ta1, gs_op_add)
80 call gs_xh%op(ta2, gs_op_add)
81 call gs_xh%op(ta3, gs_op_add)
83 do concurrent(i = 1:n)
84 ta1%x(i,1,1,1) = ta1%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
85 ta2%x(i,1,1,1) = ta2%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
86 ta3%x(i,1,1,1) = ta3%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
89 call cdtp(wa1%x, ta1%x, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
90 call cdtp(wa2%x, ta2%x, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
91 call cdtp(wa3%x, ta3%x, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
93 call ax%compute(p_res%x,p%x,c_xh,p%msh,p%Xh)
95 do concurrent(i = 1:n)
96 p_res%x(i,1,1,1) = (-p_res%x(i,1,1,1)) &
97 + wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1)
103 do concurrent(i = 1:n)
104 wa1%x(i,1,1,1) = 0.0_rp
105 wa2%x(i,1,1,1) = 0.0_rp
106 wa3%x(i,1,1,1) = 0.0_rp
109 call bc_sym_surface%apply_surfvec(wa1%x,wa2%x,wa3%x,ta1%x, ta2%x, ta3%x, n)
112 do concurrent(i = 1:n)
113 ta1%x(i,1,1,1) = 0.0_rp
114 ta2%x(i,1,1,1) = 0.0_rp
115 ta3%x(i,1,1,1) = 0.0_rp
118 call bc_prs_surface%apply_surfvec(ta1%x, ta2%x, ta3%x, u%x, v%x, w%x, n)
120 do concurrent(i = 1:n)
121 p_res%x(i,1,1,1) = p_res%x(i,1,1,1) &
122 - (dtbd * (ta1%x(i,1,1,1) + ta2%x(i,1,1,1) + ta3%x(i,1,1,1)))&
123 - (wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1))
131 p, f_x, f_y, f_z, c_Xh, msh, Xh, mu, rho, bd, dt, n)
132 class(
ax_t),
intent(in) :: Ax
133 type(
mesh_t),
intent(inout) :: msh
134 type(
space_t),
intent(inout) :: Xh
135 type(
field_t),
intent(inout) :: p, u, v, w
136 type(
field_t),
intent(inout) :: u_res, v_res, w_res
137 type(
field_t),
intent(inout) :: f_x, f_y, f_z
138 type(
coef_t),
intent(inout) :: c_Xh
139 real(kind=
rp),
intent(in) :: mu
140 real(kind=
rp),
intent(in) :: rho
141 real(kind=
rp),
intent(in) :: bd
142 real(kind=
rp),
intent(in) :: dt
143 integer :: temp_indices(3)
144 type(
field_t),
pointer :: ta1, ta2, ta3
145 integer,
intent(in) :: n
149 c_xh%h1(i,1,1,1) = mu
150 c_xh%h2(i,1,1,1) = rho * (bd / dt)
154 call ax%compute(u_res%x, u%x, c_xh, msh, xh)
155 call ax%compute(v_res%x, v%x, c_xh, msh, xh)
156 call ax%compute(w_res%x, w%x, c_xh, msh, xh)
162 call opgrad(ta1%x, ta2%x, ta3%x, p%x, c_xh)
164 do concurrent(i = 1:n)
165 u_res%x(i,1,1,1) = (-u_res%x(i,1,1,1)) - ta1%x(i,1,1,1) + f_x%x(i,1,1,1)
166 v_res%x(i,1,1,1) = (-v_res%x(i,1,1,1)) - ta2%x(i,1,1,1) + f_y%x(i,1,1,1)
167 w_res%x(i,1,1,1) = (-w_res%x(i,1,1,1)) - ta3%x(i,1,1,1) + f_z%x(i,1,1,1)
Defines a Matrix-vector product.
Dirichlet condition applied in the facet normal direction.
integer, parameter, public rp
Global precision used in computations.
subroutine, public cdtp(dtx, x, dr, ds, dt, coef)
Apply D^T to a scalar field, where D is the derivative matrix.
subroutine, public opgrad(ux, uy, uz, u, coef, es, ee)
Compute the gradient of a scalar field.
subroutine, public curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
Residuals in the Pn-Pn formulation (CPU version)
subroutine pnpn_vel_res_cpu_compute(Ax, u, v, w, u_res, v_res, w_res, p, f_x, f_y, f_z, c_Xh, msh, Xh, mu, rho, bd, dt, n)
subroutine pnpn_prs_res_cpu_compute(p, p_res, u, v, w, u_e, v_e, w_e, f_x, f_y, f_z, c_Xh, gs_Xh, bc_prs_surface, bc_sym_surface, Ax, bd, dt, mu, rho)
Defines Pressure and velocity residuals in the Pn-Pn formulation.
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
type(scratch_registry_t), target, public neko_scratch_registry
Global scratch registry.
Defines a function space.
Base type for a matrix-vector product providing .
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Dirichlet condition in facet normal direction.
Abstract type to compute pressure residual.
Abstract type to compute velocity residual.
The function space for the SEM solution fields.