Neko  0.8.99
A portable framework for high-order spectral element flow simulations
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
8  use facet_normal, only : facet_normal_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 
32 contains
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(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
43  type(facet_normal_t), intent(inout) :: bc_prs_surface
44  type(facet_normal_t), intent(inout) :: bc_sym_surface
45  class(ax_t), intent(inout) :: Ax
46  real(kind=rp), intent(inout) :: 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 
195  end subroutine pnpn_prs_res_stress_cpu_compute
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(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
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 
239  end subroutine pnpn_vel_res_stress_cpu_compute
240 
241 end 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:311
subroutine, public cmult2(a, b, c, n)
Multiplication by constant c .
Definition: math.f90:701
subroutine, public invers2(a, b, n)
Compute inverted vector .
Definition: math.f90:500
subroutine, public cfill(a, c, n)
Set all elements to a constant c .
Definition: math.f90:348
subroutine, public col2(a, b, n)
Vector multiplication .
Definition: math.f90:729
subroutine, public copy(a, b, n)
Copy a vector .
Definition: math.f90:239
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:545
subroutine, public rzero(a, n)
Zero a real vector.
Definition: math.f90:195
subroutine, public sub2(a, b, n)
Vector substraction .
Definition: math.f90:629
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.
Definition: operators.f90:171
subroutine, public curl(w1, w2, w3, u1, u2, u3, work1, work2, coef)
Definition: operators.f90:362
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.
Definition: operators.f90:430
subroutine, public cdtp(dtx, x, dr, ds, dt, coef, es, ee)
Apply D^T to a scalar field, where D is the derivative matrix.
Definition: operators.f90:230
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_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.
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