35 f_x, f_y, f_z, c_Xh, gs_Xh, bc_prs_surface, bc_sym_surface, Ax, bd, dt,&
37 type(
field_t),
intent(inout) :: p, u, v, w
38 type(
field_t),
intent(inout) :: u_e, v_e, w_e
39 type(
field_t),
intent(inout) :: p_res
40 type(
field_t),
intent(inout) :: f_x, f_y, f_z
41 type(
coef_t),
intent(inout) :: c_Xh
42 type(
gs_t),
intent(inout) :: gs_Xh
45 class(
ax_t),
intent(inout) :: Ax
46 real(kind=
rp),
intent(inout) :: bd
47 real(kind=
rp),
intent(in) :: dt
49 type(
field_t),
intent(in) :: rho
51 integer :: n, nelv, lxyz
54 type(
field_t),
pointer :: ta1, ta2, ta3, wa1, wa2, wa3, work1, work2, work3
55 type(
field_t),
pointer :: s11, s22, s33, s12, s13, s23
56 integer :: temp_indices(15)
82 call rzero(c_xh%h2, n)
86 call curl(ta1, ta2, ta3, u_e, v_e, w_e, work1, work2, c_xh)
87 call curl(wa1, wa2, wa3, ta1, ta2, ta3, work1, work2, c_xh)
89 call col2(wa1%x, mu%x, n)
90 call col2(wa2%x, mu%x, n)
91 call col2(wa3%x, mu%x, n)
95 call strain_rate(s11%x, s22%x, s33%x, s12%x, s13%x, s23%x, &
100 call dudxyz(ta1%x, mu%x, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
101 call dudxyz(ta2%x, mu%x, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
102 call dudxyz(ta3%x, mu%x, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
104 call cmult(ta1%x, 2.0_rp, n)
105 call cmult(ta2%x, 2.0_rp, n)
106 call cmult(ta3%x, 2.0_rp, n)
110 call vdot3(work1%x(:, :, :, e), &
111 ta1%x(:, :, :, e), ta2%x(:, :, :, e), ta3%x(:, :, :, e), &
112 s11%x(:, :, :, e), s12%x(:, :, :, e), s13%x(:, :, :, e), &
115 call vdot3 (work2%x(:, :, :, e), &
116 ta1%x(:, :, :, e), ta2%x(:, :, :, e), ta3%x(:, :, :, e), &
117 s12%x(:, :, :, e), s22%x(:, :, :, e), s23%x(:, :, :, e), &
120 call vdot3 (work3%x(:, :, :, e), &
121 ta1%x(:, :, :, e), ta2%x(:, :, :, e), ta3%x(:, :, :, e), &
122 s13%x(:, :, :, e), s23%x(:, :, :, e), s33%x(:, :, :, e), &
130 call sub2(wa1%x, work1%x, n)
131 call sub2(wa2%x, work2%x, n)
132 call sub2(wa3%x, work3%x, n)
134 do concurrent(i = 1:n)
135 ta1%x(i,1,1,1) = f_x%x(i,1,1,1) / rho%x(i,1,1,1) &
136 - ((wa1%x(i,1,1,1) / rho%x(i,1,1,1)) * c_xh%B(i,1,1,1))
137 ta2%x(i,1,1,1) = f_y%x(i,1,1,1) / rho%x(i,1,1,1) &
138 - ((wa2%x(i,1,1,1) / rho%x(i,1,1,1)) * c_xh%B(i,1,1,1))
139 ta3%x(i,1,1,1) = f_z%x(i,1,1,1) / rho%x(i,1,1,1) &
140 - ((wa3%x(i,1,1,1) / rho%x(i,1,1,1)) * c_xh%B(i,1,1,1))
143 call gs_xh%op(ta1, gs_op_add)
144 call gs_xh%op(ta2, gs_op_add)
145 call gs_xh%op(ta3, gs_op_add)
148 ta1%x(i,1,1,1) = ta1%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
149 ta2%x(i,1,1,1) = ta2%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
150 ta3%x(i,1,1,1) = ta3%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
154 call cdtp(wa1%x, ta1%x, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
155 call cdtp(wa2%x, ta2%x, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
156 call cdtp(wa3%x, ta3%x, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
159 call ax%compute(p_res%x, p%x, c_xh, p%msh, p%Xh)
162 p_res%x(i,1,1,1) = (-p_res%x(i,1,1,1)) &
163 + wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1)
170 wa1%x(i,1,1,1) = 0.0_rp
171 wa2%x(i,1,1,1) = 0.0_rp
172 wa3%x(i,1,1,1) = 0.0_rp
175 call bc_sym_surface%apply_surfvec(wa1%x, wa2%x, wa3%x, ta1%x, ta2%x, ta3%x,&
180 ta1%x(i,1,1,1) = 0.0_rp
181 ta2%x(i,1,1,1) = 0.0_rp
182 ta3%x(i,1,1,1) = 0.0_rp
185 call bc_prs_surface%apply_surfvec(ta1%x, ta2%x, ta3%x, u%x, v%x, w%x, n)
188 p_res%x(i,1,1,1) = p_res%x(i,1,1,1) &
189 - (dtbd * (ta1%x(i,1,1,1) + ta2%x(i,1,1,1) + ta3%x(i,1,1,1)))&
190 - (wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1))
198 p, f_x, f_y, f_z, c_Xh, msh, Xh, mu, rho, bd, dt, n)
199 class(
ax_t),
intent(in) :: Ax
200 type(
mesh_t),
intent(inout) :: msh
201 type(
space_t),
intent(inout) :: Xh
202 type(
field_t),
intent(inout) :: p, u, v, w
203 type(
field_t),
intent(inout) :: u_res, v_res, w_res
204 type(
field_t),
intent(inout) :: f_x, f_y, f_z
205 type(
coef_t),
intent(inout) :: c_Xh
206 type(
field_t),
intent(in) :: mu
207 type(
field_t),
intent(in) :: rho
208 real(kind=
rp),
intent(in) :: bd
209 real(kind=
rp),
intent(in) :: dt
210 integer :: temp_indices(3)
211 type(
field_t),
pointer :: ta1, ta2, ta3
212 integer,
intent(in) :: n
215 call copy(c_xh%h1, mu%x, n)
216 call cmult2(c_xh%h2, rho%x, bd / dt, n)
220 call ax%compute_vector(u_res%x, v_res%x, w_res%x, u%x, v%x, w%x, c_xh,&
228 call opgrad(ta1%x, ta2%x, ta3%x, p%x, c_xh)
232 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)
233 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)
234 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.
subroutine, public cmult(a, c, n)
Multiplication by constant c .
subroutine, public cmult2(a, b, c, n)
Multiplication by constant c .
subroutine, public invers2(a, b, n)
Compute inverted vector .
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
subroutine, public col2(a, b, n)
Vector multiplication .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public vdot3(dot, u1, u2, u3, v1, v2, v3, n)
Compute a dot product (3-d version) assuming vector components etc.
subroutine, public rzero(a, n)
Zero a real vector.
subroutine, public sub2(a, b, n)
Vector substraction .
integer, parameter, public rp
Global precision used in computations.
subroutine, public opgrad(ux, uy, uz, u, coef, es, ee)
Compute the weak gradient of a scalar field, i.e. the gradient multiplied by the mass matrix.
subroutine, public curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
subroutine, public strain_rate(s11, s22, s33, s12, s13, s23, u, v, w, coef)
Compute the strain rate tensor, i.e 0.5 * du_i/dx_j + du_j/dx_i.
subroutine, public cdtp(dtx, x, dr, ds, dt, coef, es, ee)
Apply D^T to a scalar field, where D is the derivative matrix.
subroutine, public dudxyz(du, u, dr, ds, dt, coef)
Compute derivative of a scalar field along a single direction.
Residuals in the Pn-Pn formulation (CPU version)
subroutine pnpn_prs_res_stress_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)
subroutine pnpn_vel_res_stress_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)
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.
CPU implementation of the pressure residual for the PnPn fluid with full viscous stress formulation.
CPU implementation of the velocity residual for the PnPn fluid with full viscous stress formulation.
Abstract type to compute pressure residual.
Abstract type to compute velocity residual.
The function space for the SEM solution fields.