31 f_y, f_z, c_Xh, gs_Xh, bc_prs_surface, bc_sym_surface, Ax, bd, dt, mu, &
33 type(
field_t),
intent(inout) :: p, u, v, w
34 type(
field_t),
intent(inout) :: u_e, v_e, w_e
35 type(
field_t),
intent(inout) :: p_res
36 type(
field_t),
intent(inout) :: f_x, f_y, f_z
37 type(
coef_t),
intent(inout) :: c_Xh
38 type(
gs_t),
intent(inout) :: gs_Xh
41 class(
ax_t),
intent(inout) :: Ax
42 real(kind=
rp),
intent(inout) :: bd
43 real(kind=
rp),
intent(in) :: dt
45 type(
field_t),
intent(in) :: rho
46 real(kind=
rp) :: dtbd, rho_val, mu_val
49 type(
field_t),
pointer :: ta1, ta2, ta3, wa1, wa2, wa3, work1, work2
50 integer :: temp_indices(8)
64 rho_val = rho%x(1,1,1,1)
65 mu_val = mu%x(1,1,1,1)
67 c_xh%h1(i,1,1,1) = 1.0_rp / rho_val
68 c_xh%h2(i,1,1,1) = 0.0_rp
72 call curl(ta1, ta2, ta3, u_e, v_e, w_e, work1, work2, c_xh)
73 call curl(wa1, wa2, wa3, ta1, ta2, ta3, work1, work2, c_xh)
76 do concurrent(i = 1:n)
77 ta1%x(i,1,1,1) = f_x%x(i,1,1,1) / rho_val &
78 - ((wa1%x(i,1,1,1) * (mu_val / rho_val)) * c_xh%B(i,1,1,1))
79 ta2%x(i,1,1,1) = f_y%x(i,1,1,1) / rho_val &
80 - ((wa2%x(i,1,1,1) * (mu_val / rho_val)) * c_xh%B(i,1,1,1))
81 ta3%x(i,1,1,1) = f_z%x(i,1,1,1) / rho_val &
82 - ((wa3%x(i,1,1,1) * (mu_val / rho_val)) * c_xh%B(i,1,1,1))
85 call gs_xh%op(ta1, gs_op_add)
86 call gs_xh%op(ta2, gs_op_add)
87 call gs_xh%op(ta3, gs_op_add)
89 do concurrent(i = 1:n)
90 ta1%x(i,1,1,1) = ta1%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
91 ta2%x(i,1,1,1) = ta2%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
92 ta3%x(i,1,1,1) = ta3%x(i,1,1,1) * c_xh%Binv(i,1,1,1)
95 call cdtp(wa1%x, ta1%x, c_xh%drdx, c_xh%dsdx, c_xh%dtdx, c_xh)
96 call cdtp(wa2%x, ta2%x, c_xh%drdy, c_xh%dsdy, c_xh%dtdy, c_xh)
97 call cdtp(wa3%x, ta3%x, c_xh%drdz, c_xh%dsdz, c_xh%dtdz, c_xh)
99 call ax%compute(p_res%x, p%x, c_xh, p%msh, p%Xh)
101 do concurrent(i = 1:n)
102 p_res%x(i,1,1,1) = (-p_res%x(i,1,1,1)) &
103 + wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1)
109 do concurrent(i = 1:n)
110 wa1%x(i,1,1,1) = 0.0_rp
111 wa2%x(i,1,1,1) = 0.0_rp
112 wa3%x(i,1,1,1) = 0.0_rp
115 call bc_sym_surface%apply_surfvec(wa1%x, wa2%x, wa3%x, ta1%x, ta2%x, ta3%x,&
119 do concurrent(i = 1:n)
120 ta1%x(i,1,1,1) = 0.0_rp
121 ta2%x(i,1,1,1) = 0.0_rp
122 ta3%x(i,1,1,1) = 0.0_rp
125 call bc_prs_surface%apply_surfvec(ta1%x, ta2%x, ta3%x, u%x, v%x, w%x, n)
127 do concurrent(i = 1:n)
128 p_res%x(i,1,1,1) = p_res%x(i,1,1,1) &
129 - (dtbd * (ta1%x(i,1,1,1) + ta2%x(i,1,1,1) + ta3%x(i,1,1,1)))&
130 - (wa1%x(i,1,1,1) + wa2%x(i,1,1,1) + wa3%x(i,1,1,1))
138 p, f_x, f_y, f_z, c_Xh, msh, Xh, mu, rho, bd, dt, n)
139 class(
ax_t),
intent(in) :: Ax
140 type(
mesh_t),
intent(inout) :: msh
141 type(
space_t),
intent(inout) :: Xh
142 type(
field_t),
intent(inout) :: p, u, v, w
143 type(
field_t),
intent(inout) :: u_res, v_res, w_res
144 type(
field_t),
intent(inout) :: f_x, f_y, f_z
145 type(
coef_t),
intent(inout) :: c_Xh
146 type(
field_t),
intent(in) :: mu
147 type(
field_t),
intent(in) :: rho
148 real(kind=
rp),
intent(in) :: bd
149 real(kind=
rp),
intent(in) :: dt
150 real(kind=
rp) :: rho_val, mu_val
151 integer :: temp_indices(3)
152 type(
field_t),
pointer :: ta1, ta2, ta3
153 integer,
intent(in) :: n
157 rho_val = rho%x(1,1,1,1)
158 mu_val = mu%x(1,1,1,1)
161 c_xh%h1(i,1,1,1) = mu_val
162 c_xh%h2(i,1,1,1) = rho_val * bd / dt
166 call ax%compute(u_res%x, u%x, c_xh, msh, xh)
167 call ax%compute(v_res%x, v%x, c_xh, msh, xh)
168 call ax%compute(w_res%x, w%x, c_xh, msh, xh)
173 call opgrad(ta1%x, ta2%x, ta3%x, p%x, c_xh)
175 do concurrent(i = 1:n)
176 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)
177 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)
178 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 cmult2(a, b, c, n)
Multiplication by constant c .
subroutine, public invers2(a, b, n)
Compute inverted vector .
subroutine, public copy(a, b, n)
Copy a vector .
subroutine, public rzero(a, n)
Zero a real vector.
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 cdtp(dtx, x, dr, ds, dt, coef, es, ee)
Apply D^T to a scalar field, where D is the derivative matrix.
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.