Loading [MathJax]/extensions/tex2jax.js
Neko 0.9.99
A portable framework for high-order spectral element flow simulations
All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Pages
pnpn_res_stress_cpu.f90
Go to the documentation of this file.
1
3 use gather_scatter, only : gs_t, gs_op_add
5 use field, only : field_t
6 use ax_product, only : ax_t
7 use coefs, only : coef_t
11 use mesh, only : mesh_t
12 use num_types, only : rp
13 use space, only : space_t
14 use math, only : rzero, vdot3, cmult, sub2, col2, copy, cfill, invers2, cmult2
15 implicit none
16 private
17
20 type, public, extends(pnpn_prs_res_t) :: pnpn_prs_res_stress_cpu_t
21 contains
22 procedure, nopass :: compute => pnpn_prs_res_stress_cpu_compute
24
27 type, public, extends(pnpn_vel_res_t) :: pnpn_vel_res_stress_cpu_t
28 contains
29 procedure, nopass :: compute => pnpn_vel_res_stress_cpu_compute
31
32contains
33
34 subroutine pnpn_prs_res_stress_cpu_compute(p, p_res, u, v, w, u_e, v_e, w_e,&
35 f_x, f_y, f_z, c_Xh, gs_Xh, bc_prs_surface, bc_sym_surface, Ax, bd, dt,&
36 mu, rho)
37 type(field_t), intent(inout) :: p, u, v, w
38 type(field_t), intent(in) :: u_e, v_e, w_e
39 type(field_t), intent(inout) :: p_res
40 type(field_t), intent(in) :: f_x, f_y, f_z
41 type(coef_t), intent(inout) :: c_Xh
42 type(gs_t), intent(inout) :: gs_Xh
43 type(facet_normal_t), intent(in) :: bc_prs_surface
44 type(facet_normal_t), intent(in) :: bc_sym_surface
45 class(ax_t), intent(inout) :: Ax
46 real(kind=rp), intent(in) :: bd
47 real(kind=rp), intent(in) :: dt
48 type(field_t), intent(in) :: mu
49 type(field_t), intent(in) :: rho
50 real(kind=rp) :: dtbd
51 integer :: n, nelv, lxyz
52 integer :: i, e
53 ! Work arrays
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)
57
58 ! Work arrays
59 call neko_scratch_registry%request_field(ta1, temp_indices(1))
60 call neko_scratch_registry%request_field(ta2, temp_indices(2))
61 call neko_scratch_registry%request_field(ta3, temp_indices(3))
62 call neko_scratch_registry%request_field(wa1, temp_indices(4))
63 call neko_scratch_registry%request_field(wa2, temp_indices(5))
64 call neko_scratch_registry%request_field(wa3, temp_indices(6))
65 call neko_scratch_registry%request_field(work1, temp_indices(7))
66 call neko_scratch_registry%request_field(work2, temp_indices(8))
67 call neko_scratch_registry%request_field(work3, temp_indices(9))
68
69 ! Stress tensor
70 call neko_scratch_registry%request_field(s11, temp_indices(10))
71 call neko_scratch_registry%request_field(s22, temp_indices(11))
72 call neko_scratch_registry%request_field(s33, temp_indices(12))
73 call neko_scratch_registry%request_field(s12, temp_indices(13))
74 call neko_scratch_registry%request_field(s13, temp_indices(14))
75 call neko_scratch_registry%request_field(s23, temp_indices(15))
76
77 n = c_xh%dof%size()
78 lxyz = c_xh%Xh%lxyz
79 nelv = c_xh%msh%nelv
80
81 call invers2(c_xh%h1, rho%x, n)
82 call rzero(c_xh%h2, n)
83 c_xh%ifh2 = .false.
84
85 ! mu times the double curl of the velocity
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)
88
89 call col2(wa1%x, mu%x, n)
90 call col2(wa2%x, mu%x, n)
91 call col2(wa3%x, mu%x, n)
92
93
94 ! The strain rate tensor
95 call strain_rate(s11%x, s22%x, s33%x, s12%x, s13%x, s23%x, &
96 u_e, v_e, w_e, c_xh)
97
98
99 ! Gradient of viscosity * 2
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)
103
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)
107
108 ! S^T grad \mu
109 do e = 1, nelv
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), &
113 lxyz)
114
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), &
118 lxyz)
119
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), &
123 lxyz)
124 end do
125
126 ! Subtract the two terms of the viscous stress to get
127 ! \nabla x \nabla u - S^T \nabla \mu
128 ! The sign is consitent with the fact that we subtract the term
129 ! below.
130 call sub2(wa1%x, work1%x, n)
131 call sub2(wa2%x, work2%x, n)
132 call sub2(wa3%x, work3%x, n)
133
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))
141 end do
142
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)
146
147 do i = 1, n
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)
151 end do
152
153 ! Compute the components of the divergence of the rhs
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)
157
158 ! The laplacian of the pressure
159 call ax%compute(p_res%x, p%x, c_xh, p%msh, p%Xh)
160
161 do i = 1, n
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)
164 end do
165
166 !
167 ! Surface velocity terms
168 !
169 do i = 1, n
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
173 end do
174
175 call bc_sym_surface%apply_surfvec(wa1%x, wa2%x, wa3%x, ta1%x, ta2%x, ta3%x,&
176 n)
177
178 dtbd = bd / dt
179 do i = 1, n
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
183 end do
184
185 call bc_prs_surface%apply_surfvec(ta1%x, ta2%x, ta3%x, u%x, v%x, w%x, n)
186
187 do i = 1, 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))
191 end do
192
193 call neko_scratch_registry%relinquish_field(temp_indices)
194
196
197 subroutine pnpn_vel_res_stress_cpu_compute(Ax, u, v, w, u_res, v_res, w_res, &
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(in) :: 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
213 integer :: i
214
215 call copy(c_xh%h1, mu%x, n)
216 call cmult2(c_xh%h2, rho%x, bd / dt, n)
217 c_xh%ifh2 = .true.
218
219 ! Viscous stresses
220 call ax%compute_vector(u_res%x, v_res%x, w_res%x, u%x, v%x, w%x, c_xh,&
221 msh, xh)
222
223 call neko_scratch_registry%request_field(ta1, temp_indices(1))
224 call neko_scratch_registry%request_field(ta2, temp_indices(2))
225 call neko_scratch_registry%request_field(ta3, temp_indices(3))
226
227 ! Pressure gradient
228 call opgrad(ta1%x, ta2%x, ta3%x, p%x, c_xh)
229
230 ! Sum all the terms
231 do i = 1, n
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)
235 end do
236
237 call neko_scratch_registry%relinquish_field(temp_indices)
238
240
241end module pnpn_res_stress_cpu
Defines a Matrix-vector product.
Definition ax.f90:34
Coefficients.
Definition coef.f90:34
Dirichlet condition applied in the facet normal direction.
Defines a field.
Definition field.f90:34
Gather-scatter.
Definition math.f90:60
subroutine, public cmult(a, c, n)
Multiplication by constant c .
Definition math.f90:310
subroutine, public cmult2(a, b, c, n)
Multiplication by constant c .
Definition math.f90:700
subroutine, public invers2(a, b, n)
Compute inverted vector .
Definition math.f90:499
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
Definition math.f90:347
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:728
subroutine, public copy(a, b, n)
Copy a vector .
Definition math.f90:238
subroutine, public vdot3(dot, u1, u2, u3, v1, v2, v3, n)
Compute a dot product (3-d version) assuming vector components etc.
Definition math.f90:544
subroutine, public rzero(a, n)
Zero a real vector.
Definition math.f90:194
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:628
Defines a mesh.
Definition mesh.f90:34
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Operators.
Definition operators.f90:34
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.
Definition operators.f90:76
Residuals in the Pn-Pn formulation (CPU version)
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)
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)
Defines Pressure and velocity residuals in the Pn-Pn formulation.
Definition pnpn_res.f90:34
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.
Definition space.f90:34
Base type for a matrix-vector product providing .
Definition ax.f90:43
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
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.
Definition pnpn_res.f90:47
Abstract type to compute velocity residual.
Definition pnpn_res.f90:53
The function space for the SEM solution fields.
Definition space.f90:62