Neko  0.8.99
A portable framework for high-order spectral element flow simulations
dynamic_smagorinsky_cpu.f90
Go to the documentation of this file.
1 ! Copyright (c) 2024, The Neko Authors
2 ! All rights reserved.
3 !
4 ! Redistribution and use in source and binary forms, with or without
5 ! modification, are permitted provided that the following conditions
6 ! are met:
7 !
8 ! * Redistributions of source code must retain the above copyright
9 ! notice, this list of conditions and the following disclaimer.
10 !
11 ! * Redistributions in binary form must reproduce the above
12 ! copyright notice, this list of conditions and the following
13 ! disclaimer in the documentation and/or other materials provided
14 ! with the distribution.
15 !
16 ! * Neither the name of the authors nor the names of its
17 ! contributors may be used to endorse or promote products derived
18 ! from this software without specific prior written permission.
19 !
20 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24 ! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25 ! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26 ! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27 ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 ! POSSIBILITY OF SUCH DAMAGE.
32 !
34 module dynamic_smagorinsky_cpu
35  use num_types, only : rp
36  use field_list, only : field_list_t
37  use math, only : cadd, neko_eps, col2, sub2, col3, cmult
38  use scratch_registry, only : neko_scratch_registry
39  use field_registry, only : neko_field_registry
40  use field, only : field_t
41  use operators, only : strain_rate
42  use coefs, only : coef_t
43  use elementwise_filter, only : elementwise_filter_t
44  use gs_ops, only : gs_op_add
45  implicit none
46  private
47 
48  public :: dynamic_smagorinsky_compute_cpu
49 
50 contains
51 
64  subroutine dynamic_smagorinsky_compute_cpu(t, tstep, coef, nut, delta, &
65  c_dyn, test_filter, mij, lij, num, den)
66  real(kind=rp), intent(in) :: t
67  integer, intent(in) :: tstep
68  type(coef_t), intent(in) :: coef
69  type(field_t), intent(inout) :: nut
70  type(field_t), intent(in) :: delta
71  type(field_t), intent(inout) :: c_dyn
72  type(elementwise_filter_t), intent(inout) :: test_filter
73  type(field_t), intent(inout) :: mij(6), lij(6)
74  type(field_t), intent(inout) :: num, den
75 
76  type(field_t), pointer :: u, v, w
77  type(field_t) :: c_dyn_curr
78  ! the strain rate tensor
79  type(field_t), pointer :: s11, s22, s33, s12, s13, s23, s_abs
80  real(kind=rp) :: alpha ! running averaging coefficient
81  integer :: temp_indices(7)
82  integer :: i
83 
84  if (tstep .eq. 1) then
85  alpha = 1.0_rp
86  else
87  alpha = 0.9_rp
88  end if
89 
90  u => neko_field_registry%get_field_by_name("u")
91  v => neko_field_registry%get_field_by_name("v")
92  w => neko_field_registry%get_field_by_name("w")
93 
94  call neko_scratch_registry%request_field(s11, temp_indices(1))
95  call neko_scratch_registry%request_field(s22, temp_indices(2))
96  call neko_scratch_registry%request_field(s33, temp_indices(3))
97  call neko_scratch_registry%request_field(s12, temp_indices(4))
98  call neko_scratch_registry%request_field(s13, temp_indices(5))
99  call neko_scratch_registry%request_field(s23, temp_indices(6))
100  call neko_scratch_registry%request_field(s_abs, temp_indices(7))
101 
102  ! Compute the strain rate tensor
103  call strain_rate(s11%x, s22%x, s33%x, s12%x, s13%x, s23%x, u, v, w, coef)
104 
105  call coef%gs_h%op(s11%x, s11%dof%size(), gs_op_add)
106  call coef%gs_h%op(s22%x, s11%dof%size(), gs_op_add)
107  call coef%gs_h%op(s33%x, s11%dof%size(), gs_op_add)
108  call coef%gs_h%op(s12%x, s11%dof%size(), gs_op_add)
109  call coef%gs_h%op(s13%x, s11%dof%size(), gs_op_add)
110  call coef%gs_h%op(s23%x, s11%dof%size(), gs_op_add)
111 
112  do concurrent(i = 1:s11%dof%size())
113  s11%x(i,1,1,1) = s11%x(i,1,1,1) * coef%mult(i,1,1,1)
114  s22%x(i,1,1,1) = s22%x(i,1,1,1) * coef%mult(i,1,1,1)
115  s33%x(i,1,1,1) = s33%x(i,1,1,1) * coef%mult(i,1,1,1)
116  s12%x(i,1,1,1) = s12%x(i,1,1,1) * coef%mult(i,1,1,1)
117  s13%x(i,1,1,1) = s13%x(i,1,1,1) * coef%mult(i,1,1,1)
118  s23%x(i,1,1,1) = s23%x(i,1,1,1) * coef%mult(i,1,1,1)
119  end do
120 
121  do concurrent(i = 1:u%dof%size())
122  s_abs%x(i,1,1,1) = sqrt(2.0_rp * (s11%x(i,1,1,1)*s11%x(i,1,1,1) + &
123  s22%x(i,1,1,1)*s22%x(i,1,1,1) + &
124  s33%x(i,1,1,1)*s33%x(i,1,1,1)) + &
125  4.0_rp * (s12%x(i,1,1,1)*s12%x(i,1,1,1) + &
126  s13%x(i,1,1,1)*s13%x(i,1,1,1) + &
127  s23%x(i,1,1,1)*s23%x(i,1,1,1)))
128  end do
129 
130  call compute_lij_cpu(lij, u, v, w, test_filter, u%dof%size(), u%msh%nelv)
131  call compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
132  s_abs, test_filter, delta, u%dof%size(), u%msh%nelv)
133  call compute_num_den_cpu(num, den, lij, mij, alpha, u%dof%size())
134 
135  do concurrent(i =1:u%dof%size())
136  if (den%x(i,1,1,1) .gt. 0.0_rp) then
137  c_dyn%x(i,1,1,1) = 0.5_rp * (num%x(i,1,1,1)/den%x(i,1,1,1))
138  else
139  c_dyn%x(i,1,1,1) = 0.0_rp
140  end if
141  c_dyn%x(i,1,1,1) = max(c_dyn%x(i,1,1,1),0.0_rp)
142  nut%x(i,1,1,1) = c_dyn%x(i,1,1,1) * delta%x(i,1,1,1)**2 * s_abs%x(i,1,1,1)
143  end do
144 
145 
146  call neko_scratch_registry%relinquish_field(temp_indices)
147  end subroutine dynamic_smagorinsky_compute_cpu
148 
157  subroutine compute_lij_cpu(lij, u, v, w, test_filter, n, nelv)
158  type(field_t), intent(inout) :: lij(6)
159  type(field_t), pointer, intent(in) :: u, v, w
160  type(elementwise_filter_t), intent(inout) :: test_filter
161  integer, intent(in) :: n
162  integer, intent(inout) :: nelv
163  integer :: i
165  real(kind=rp), dimension(u%dof%size()) :: fu, fv, fw
166 
167  ! Use test filter for the velocity fields
168  call test_filter%filter_3d(fu, u%x, nelv)
169  call test_filter%filter_3d(fv, v%x, nelv)
170  call test_filter%filter_3d(fw, w%x, nelv)
171 
172  !! The first term
173  do concurrent(i = 1:n)
174  lij(1)%x(i,1,1,1) = fu(i) * fu(i)
175  lij(2)%x(i,1,1,1) = fv(i) * fv(i)
176  lij(3)%x(i,1,1,1) = fw(i) * fw(i)
177  lij(4)%x(i,1,1,1) = fu(i) * fv(i)
178  lij(5)%x(i,1,1,1) = fu(i) * fw(i)
179  lij(6)%x(i,1,1,1) = fv(i) * fw(i)
180  end do
181 
182  !! Subtract the second term:
183  !! use test filter for the cross terms
184  !! fu and fv are used as work array
185  call col3(fu, u%x, u%x, n)
186  call test_filter%filter_3d(fv, fu, nelv)
187  call sub2(lij(1)%x, fv, n)
188 
189  call col3(fu, v%x, v%x, n)
190  call test_filter%filter_3d(fv, fu, nelv)
191  call sub2(lij(2)%x, fv, n)
192 
193  call col3(fu, w%x, w%x, n)
194  call test_filter%filter_3d(fv, fu, nelv)
195  call sub2(lij(3)%x, fv, n)
196 
197  call col3(fu, u%x, v%x, n)
198  call test_filter%filter_3d(fv, fu, nelv)
199  call sub2(lij(4)%x, fv, n)
200 
201  call col3(fu, u%x, w%x, n)
202  call test_filter%filter_3d(fv, fu, nelv)
203  call sub2(lij(5)%x, fv, n)
204 
205  call col3(fu, v%x, w%x, n)
206  call test_filter%filter_3d(fv, fu, nelv)
207  call sub2(lij(6)%x, fv, n)
208 
209  end subroutine compute_lij_cpu
210 
219  subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
220  s_abs, test_filter, delta, n, nelv)
221  type(field_t), intent(inout) :: mij(6)
222  type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23, s_abs
223  type(elementwise_filter_t), intent(inout) :: test_filter
224  type(field_t), intent(in) :: delta
225  integer, intent(in) :: n
226  integer, intent(inout) :: nelv
227 
228  real(kind=rp), dimension(n) :: fs11, fs22, fs33, fs12, fs13, fs23, fs_abs
229  real(kind=rp) :: delta_ratio2 !test- to grid- filter ratio, squared
230  integer :: i
231  real(kind=rp) :: delta2
232 
233  delta_ratio2 = ((test_filter%nx-1)/(test_filter%nt-1))**2
234 
235  !! The first term:
236  !! _____ ____
237  !! (delta_test/delta)^2 s_abs*s_ij
238  call test_filter%filter_3d(fs_abs, s_abs%x, nelv)
239 
240  call test_filter%filter_3d(fs11, s11%x, nelv)
241  call col3(mij(1)%x, fs_abs, fs11, n)
242  call cmult(mij(1)%x, delta_ratio2, n)
243 
244  call test_filter%filter_3d(fs22, s22%x, nelv)
245  call col3(mij(2)%x, fs_abs, fs11, n)
246  call cmult(mij(2)%x, delta_ratio2, n)
247 
248  call test_filter%filter_3d(fs33, s33%x, nelv)
249  call col3(mij(3)%x, fs_abs, fs11, n)
250  call cmult(mij(3)%x, delta_ratio2, n)
251 
252  call test_filter%filter_3d(fs12, s12%x, nelv)
253  call col3(mij(4)%x, fs_abs, fs11, n)
254  call cmult(mij(4)%x, delta_ratio2, n)
255 
256  call test_filter%filter_3d(fs13, s13%x, nelv)
257  call col3(mij(5)%x, fs_abs, fs11, n)
258  call cmult(mij(5)%x, delta_ratio2, n)
259 
260  call test_filter%filter_3d(fs23, s23%x, nelv)
261  call col3(mij(6)%x, fs_abs, fs11, n)
262  call cmult(mij(6)%x, delta_ratio2, n)
263 
264  !! Substract the second term:
265  !! _____ ____ __________
266  !! (delta_test/delta)^2 s_abs*s_ij - s_abs*s_ij
267  !! fs11 and fs22 are used as work array
268  call col3(fs11, s_abs%x, s11%x, n)
269  call test_filter%filter_3d(fs22, fs11, nelv)
270  call sub2(mij(1)%x, fs22, n)
271 
272  call col3(fs11, s_abs%x, s22%x, n)
273  call test_filter%filter_3d(fs22, fs11, nelv)
274  call sub2(mij(2)%x, fs22, n)
275 
276  call col3(fs11, s_abs%x, s33%x, n)
277  call test_filter%filter_3d(fs22, fs11, nelv)
278  call sub2(mij(3)%x, fs22, n)
279 
280  call col3(fs11, s_abs%x, s12%x, n)
281  call test_filter%filter_3d(fs22, fs11, nelv)
282  call sub2(mij(4)%x, fs22, n)
283 
284  call col3(fs11, s_abs%x, s13%x, n)
285  call test_filter%filter_3d(fs22, fs11, nelv)
286  call sub2(mij(5)%x, fs22, n)
287 
288  call col3(fs11, s_abs%x, s23%x, n)
289  call test_filter%filter_3d(fs22, fs11, nelv)
290  call sub2(mij(6)%x, fs22, n)
291 
292  !! Lastly multiplied by delta^2
293  do concurrent(i = 1:n)
294  delta2 = delta%x(i,1,1,1)**2
295  mij(1)%x(i,1,1,1) = mij(1)%x(i,1,1,1) * delta2
296  mij(2)%x(i,1,1,1) = mij(2)%x(i,1,1,1) * delta2
297  mij(3)%x(i,1,1,1) = mij(3)%x(i,1,1,1) * delta2
298  mij(4)%x(i,1,1,1) = mij(4)%x(i,1,1,1) * delta2
299  mij(5)%x(i,1,1,1) = mij(5)%x(i,1,1,1) * delta2
300  mij(6)%x(i,1,1,1) = mij(6)%x(i,1,1,1) * delta2
301  end do
302 
303  end subroutine compute_mij_cpu
304 
311  subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
312  type(field_t), intent(inout) :: num, den
313  type(field_t), intent(in) :: lij(6), mij(6)
314  real(kind=rp), intent(in) :: alpha
315  integer, intent(in) :: n
316 
317  real(kind=rp), dimension(n) :: num_curr, den_curr
318  integer :: i
319 
320  do concurrent(i = 1:n)
321  num_curr(i) = mij(1)%x(i,1,1,1)*lij(1)%x(i,1,1,1) + &
322  mij(2)%x(i,1,1,1)*lij(2)%x(i,1,1,1) + &
323  mij(3)%x(i,1,1,1)*lij(3)%x(i,1,1,1) + &
324  2.0_rp*(mij(4)%x(i,1,1,1)*lij(4)%x(i,1,1,1) + &
325  mij(5)%x(i,1,1,1)*lij(5)%x(i,1,1,1) + &
326  mij(6)%x(i,1,1,1)*lij(6)%x(i,1,1,1))
327  den_curr(i) = mij(1)%x(i,1,1,1)*mij(1)%x(i,1,1,1) + &
328  mij(2)%x(i,1,1,1)*mij(2)%x(i,1,1,1) + &
329  mij(3)%x(i,1,1,1)*mij(3)%x(i,1,1,1) + &
330  2.0_rp*(mij(4)%x(i,1,1,1)*mij(4)%x(i,1,1,1) + &
331  mij(5)%x(i,1,1,1)*mij(5)%x(i,1,1,1) + &
332  mij(6)%x(i,1,1,1)*mij(6)%x(i,1,1,1))
333  end do
334 
335  ! running average over time
336  do concurrent(i = 1:n)
337  num%x(i,1,1,1) = alpha * num%x(i,1,1,1) + (1.0_rp - alpha) * num_curr(i)
338  den%x(i,1,1,1) = alpha * den%x(i,1,1,1) + (1.0_rp - alpha) * den_curr(i)
339  end do
340 
341  end subroutine compute_num_den_cpu
342 
343 end module dynamic_smagorinsky_cpu
344 
coefs
Coefficients.
Definition: coef.f90:34
dynamic_smagorinsky_cpu
Implements the CPU kernel for the smagorinsky_t type.
Definition: dynamic_smagorinsky_cpu.f90:34
dynamic_smagorinsky_cpu::compute_num_den_cpu
subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
Compute numerator and denominator for c_dyn on the CPU.
Definition: dynamic_smagorinsky_cpu.f90:312
dynamic_smagorinsky_cpu::compute_mij_cpu
subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, s_abs, test_filter, delta, n, nelv)
Compute M_ij on the CPU.
Definition: dynamic_smagorinsky_cpu.f90:221
dynamic_smagorinsky_cpu::compute_lij_cpu
subroutine compute_lij_cpu(lij, u, v, w, test_filter, n, nelv)
Compute Germano Identity on the CPU.
Definition: dynamic_smagorinsky_cpu.f90:158
dynamic_smagorinsky_cpu::dynamic_smagorinsky_compute_cpu
subroutine, public dynamic_smagorinsky_compute_cpu(t, tstep, coef, nut, delta, c_dyn, test_filter, mij, lij, num, den)
Compute eddy viscosity on the CPU.
Definition: dynamic_smagorinsky_cpu.f90:66
elementwise_filter
Implements explicit_filter_t.
Definition: elementwise_filter.f90:35
field_list
Definition: field_list.f90:1
field_registry
Defines a registry for storing solution fields.
Definition: field_registry.f90:35
field_registry::neko_field_registry
type(field_registry_t), target, public neko_field_registry
Global field registry.
Definition: field_registry.f90:75
field
Defines a field.
Definition: field.f90:34
gs_ops
Defines Gather-scatter operations.
Definition: gs_ops.f90:34
gs_ops::gs_op_add
integer, parameter, public gs_op_add
Definition: gs_ops.f90:36
math
Definition: math.f90:60
math::cmult
subroutine, public cmult(a, c, n)
Multiplication by constant c .
Definition: math.f90:277
math::cadd
subroutine, public cadd(a, s, n)
Add a scalar to vector .
Definition: math.f90:289
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition: math.f90:689
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition: math.f90:702
math::neko_eps
real(kind=rp), parameter, public neko_eps
Machine epsilon .
Definition: math.f90:67
math::sub2
subroutine, public sub2(a, b, n)
Vector substraction .
Definition: math.f90:589
num_types
Definition: num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
operators
Operators.
Definition: operators.f90:34
operators::strain_rate
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:428
scratch_registry
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
Definition: scratch_registry.f90:37
scratch_registry::neko_scratch_registry
type(scratch_registry_t), target, public neko_scratch_registry
Global scratch registry.
Definition: scratch_registry.f90:85
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition: coef.f90:55
elementwise_filter::elementwise_filter_t
Implements the explicit filter for SEM.
Definition: elementwise_filter.f90:47
field::field_t
Definition: field.f90:46
field_list::field_list_t
field_list_t, To be able to group fields together
Definition: field_list.f90:13
max
#define max(a, b)
Definition: tensor.cu:40