Neko 0.9.99
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
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 call neko_scratch_registry%relinquish_field(temp_indices)
146 end subroutine dynamic_smagorinsky_compute_cpu
147
156 subroutine compute_lij_cpu(lij, u, v, w, test_filter, n, nelv)
157 type(field_t), intent(inout) :: lij(6)
158 type(field_t), pointer, intent(in) :: u, v, w
159 type(elementwise_filter_t), intent(inout) :: test_filter
160 integer, intent(in) :: n
161 integer, intent(inout) :: nelv
162 integer :: i
164 real(kind=rp), dimension(u%dof%size()) :: fu, fv, fw
165
166 ! Use test filter for the velocity fields
167 call test_filter%filter_3d(fu, u%x, nelv)
168 call test_filter%filter_3d(fv, v%x, nelv)
169 call test_filter%filter_3d(fw, w%x, nelv)
170
171 !! The first term
172 do concurrent(i = 1:n)
173 lij(1)%x(i,1,1,1) = fu(i) * fu(i)
174 lij(2)%x(i,1,1,1) = fv(i) * fv(i)
175 lij(3)%x(i,1,1,1) = fw(i) * fw(i)
176 lij(4)%x(i,1,1,1) = fu(i) * fv(i)
177 lij(5)%x(i,1,1,1) = fu(i) * fw(i)
178 lij(6)%x(i,1,1,1) = fv(i) * fw(i)
179 end do
180
181 !! Subtract the second term:
182 !! use test filter for the cross terms
183 !! fu and fv are used as work array
184 call col3(fu, u%x, u%x, n)
185 call test_filter%filter_3d(fv, fu, nelv)
186 call sub2(lij(1)%x, fv, n)
187
188 call col3(fu, v%x, v%x, n)
189 call test_filter%filter_3d(fv, fu, nelv)
190 call sub2(lij(2)%x, fv, n)
191
192 call col3(fu, w%x, w%x, n)
193 call test_filter%filter_3d(fv, fu, nelv)
194 call sub2(lij(3)%x, fv, n)
195
196 call col3(fu, u%x, v%x, n)
197 call test_filter%filter_3d(fv, fu, nelv)
198 call sub2(lij(4)%x, fv, n)
199
200 call col3(fu, u%x, w%x, n)
201 call test_filter%filter_3d(fv, fu, nelv)
202 call sub2(lij(5)%x, fv, n)
203
204 call col3(fu, v%x, w%x, n)
205 call test_filter%filter_3d(fv, fu, nelv)
206 call sub2(lij(6)%x, fv, n)
207
208 end subroutine compute_lij_cpu
209
218 subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
219 s_abs, test_filter, delta, n, nelv)
220 type(field_t), intent(inout) :: mij(6)
221 type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23, s_abs
222 type(elementwise_filter_t), intent(inout) :: test_filter
223 type(field_t), intent(in) :: delta
224 integer, intent(in) :: n
225 integer, intent(inout) :: nelv
226
227 real(kind=rp), dimension(n) :: fs11, fs22, fs33, fs12, fs13, fs23, fs_abs
228 real(kind=rp) :: delta_ratio2 !test- to grid- filter ratio, squared
229 integer :: i
230 real(kind=rp) :: delta2
231
232 delta_ratio2 = ((test_filter%nx-1.0_rp)/(test_filter%nt-1.0_rp))**2
233
234 !! The first term:
235 !! _____ ____
236 !! (delta_test/delta)^2 s_abs*s_ij
237 call test_filter%filter_3d(fs_abs, s_abs%x, nelv)
238
239 call test_filter%filter_3d(fs11, s11%x, nelv)
240 call col3(mij(1)%x, fs_abs, fs11, n)
241 call cmult(mij(1)%x, delta_ratio2, n)
242
243 call test_filter%filter_3d(fs22, s22%x, nelv)
244 call col3(mij(2)%x, fs_abs, fs22, n)
245 call cmult(mij(2)%x, delta_ratio2, n)
246
247 call test_filter%filter_3d(fs33, s33%x, nelv)
248 call col3(mij(3)%x, fs_abs, fs33, n)
249 call cmult(mij(3)%x, delta_ratio2, n)
250
251 call test_filter%filter_3d(fs12, s12%x, nelv)
252 call col3(mij(4)%x, fs_abs, fs12, n)
253 call cmult(mij(4)%x, delta_ratio2, n)
254
255 call test_filter%filter_3d(fs13, s13%x, nelv)
256 call col3(mij(5)%x, fs_abs, fs13, n)
257 call cmult(mij(5)%x, delta_ratio2, n)
258
259 call test_filter%filter_3d(fs23, s23%x, nelv)
260 call col3(mij(6)%x, fs_abs, fs23, n)
261 call cmult(mij(6)%x, delta_ratio2, n)
262
263 !! Substract the second term:
264 !! _____ ____ __________
265 !! (delta_test/delta)^2 s_abs*s_ij - s_abs*s_ij
266 !! fs11 and fs22 are used as work array
267 call col3(fs11, s_abs%x, s11%x, n)
268 call test_filter%filter_3d(fs22, fs11, nelv)
269 call sub2(mij(1)%x, fs22, n)
270
271 call col3(fs11, s_abs%x, s22%x, n)
272 call test_filter%filter_3d(fs22, fs11, nelv)
273 call sub2(mij(2)%x, fs22, n)
274
275 call col3(fs11, s_abs%x, s33%x, n)
276 call test_filter%filter_3d(fs22, fs11, nelv)
277 call sub2(mij(3)%x, fs22, n)
278
279 call col3(fs11, s_abs%x, s12%x, n)
280 call test_filter%filter_3d(fs22, fs11, nelv)
281 call sub2(mij(4)%x, fs22, n)
282
283 call col3(fs11, s_abs%x, s13%x, n)
284 call test_filter%filter_3d(fs22, fs11, nelv)
285 call sub2(mij(5)%x, fs22, n)
286
287 call col3(fs11, s_abs%x, s23%x, n)
288 call test_filter%filter_3d(fs22, fs11, nelv)
289 call sub2(mij(6)%x, fs22, n)
290
291 !! Lastly multiplied by delta^2
292 do concurrent(i = 1:n)
293 delta2 = delta%x(i,1,1,1)**2
294 mij(1)%x(i,1,1,1) = mij(1)%x(i,1,1,1) * delta2
295 mij(2)%x(i,1,1,1) = mij(2)%x(i,1,1,1) * delta2
296 mij(3)%x(i,1,1,1) = mij(3)%x(i,1,1,1) * delta2
297 mij(4)%x(i,1,1,1) = mij(4)%x(i,1,1,1) * delta2
298 mij(5)%x(i,1,1,1) = mij(5)%x(i,1,1,1) * delta2
299 mij(6)%x(i,1,1,1) = mij(6)%x(i,1,1,1) * delta2
300 end do
301
302 end subroutine compute_mij_cpu
303
310 subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
311 type(field_t), intent(inout) :: num, den
312 type(field_t), intent(in) :: lij(6), mij(6)
313 real(kind=rp), intent(in) :: alpha
314 integer, intent(in) :: n
315
316 real(kind=rp), dimension(n) :: num_curr, den_curr
317 integer :: i
318
319 do concurrent(i = 1:n)
320 num_curr(i) = mij(1)%x(i,1,1,1)*lij(1)%x(i,1,1,1) + &
321 mij(2)%x(i,1,1,1)*lij(2)%x(i,1,1,1) + &
322 mij(3)%x(i,1,1,1)*lij(3)%x(i,1,1,1) + &
323 2.0_rp*(mij(4)%x(i,1,1,1)*lij(4)%x(i,1,1,1) + &
324 mij(5)%x(i,1,1,1)*lij(5)%x(i,1,1,1) + &
325 mij(6)%x(i,1,1,1)*lij(6)%x(i,1,1,1))
326 den_curr(i) = mij(1)%x(i,1,1,1)*mij(1)%x(i,1,1,1) + &
327 mij(2)%x(i,1,1,1)*mij(2)%x(i,1,1,1) + &
328 mij(3)%x(i,1,1,1)*mij(3)%x(i,1,1,1) + &
329 2.0_rp*(mij(4)%x(i,1,1,1)*mij(4)%x(i,1,1,1) + &
330 mij(5)%x(i,1,1,1)*mij(5)%x(i,1,1,1) + &
331 mij(6)%x(i,1,1,1)*mij(6)%x(i,1,1,1))
332 end do
333
334 ! running average over time
335 do concurrent(i = 1:n)
336 num%x(i,1,1,1) = alpha * num%x(i,1,1,1) + (1.0_rp - alpha) * num_curr(i)
337 den%x(i,1,1,1) = alpha * den%x(i,1,1,1) + (1.0_rp - alpha) * den_curr(i)
338 end do
339
340 end subroutine compute_num_den_cpu
341
342end module dynamic_smagorinsky_cpu
343
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:311
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:220
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:157
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:310
math::cadd
subroutine, public cadd(a, s, n)
Add a scalar to vector .
Definition math.f90:322
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:728
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:741
math::neko_eps
real(kind=rp), parameter, public neko_eps
Machine epsilon .
Definition math.f90:69
math::sub2
subroutine, public sub2(a, b, n)
Vector substraction .
Definition math.f90:628
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:436
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:51
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