Loading [MathJax]/extensions/tex2jax.js
Neko 0.9.1
A portable framework for high-order spectral element flow simulations
•All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Pages
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!
34module 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
50contains
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
343end 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:311
math::cadd
subroutine, public cadd(a, s, n)
Add a scalar to vector .
Definition math.f90:323
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:729
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:742
math::neko_eps
real(kind=rp), parameter, public neko_eps
Machine epsilon .
Definition math.f90:70
math::sub2
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:629
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:430
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