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())
131 call compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
132 s_abs, test_filter, delta, u%dof%size())
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 coef%gs_h%op(nut, gs_op_add)
146 call col2(nut%x, coef%mult, nut%dof%size())
147
148 call neko_scratch_registry%relinquish_field(temp_indices)
149 end subroutine dynamic_smagorinsky_compute_cpu
150
159 subroutine compute_lij_cpu(lij, u, v, w, test_filter, n)
160 type(field_t), intent(inout) :: lij(6)
161 type(field_t), pointer, intent(in) :: u, v, w
162 type(elementwise_filter_t), intent(inout) :: test_filter
163 integer, intent(in) :: n
164 integer :: i
166 integer :: temp_indices(3)
167 type(field_t), pointer :: fu, fv, fw
168
169 ! Use test filter for the velocity fields
170 call neko_scratch_registry%request_field(fu, temp_indices(1))
171 call neko_scratch_registry%request_field(fv, temp_indices(2))
172 call neko_scratch_registry%request_field(fw, temp_indices(3))
173 call test_filter%apply(fu, u)
174 call test_filter%apply(fv, v)
175 call test_filter%apply(fw, w)
176
177 !! The first term
178 do concurrent(i = 1:n)
179 lij(1)%x(i,1,1,1) = fu%x(i,1,1,1) * fu%x(i,1,1,1)
180 lij(2)%x(i,1,1,1) = fv%x(i,1,1,1) * fv%x(i,1,1,1)
181 lij(3)%x(i,1,1,1) = fw%x(i,1,1,1) * fw%x(i,1,1,1)
182 lij(4)%x(i,1,1,1) = fu%x(i,1,1,1) * fv%x(i,1,1,1)
183 lij(5)%x(i,1,1,1) = fu%x(i,1,1,1) * fw%x(i,1,1,1)
184 lij(6)%x(i,1,1,1) = fv%x(i,1,1,1) * fw%x(i,1,1,1)
185 end do
186
187 !! Subtract the second term:
188 !! use test filter for the cross terms
189 !! fu and fv are used as work array
190 call col3(fu%x, u%x, u%x, n)
191 call test_filter%apply(fv, fu)
192 call sub2(lij(1)%x, fv%x, n)
193
194 call col3(fu%x, v%x, v%x, n)
195 call test_filter%apply(fv, fu)
196 call sub2(lij(2)%x, fv%x, n)
197
198 call col3(fu%x, w%x, w%x, n)
199 call test_filter%apply(fv, fu)
200 call sub2(lij(3)%x, fv%x, n)
201
202 call col3(fu%x, u%x, v%x, n)
203 call test_filter%apply(fv, fu)
204 call sub2(lij(4)%x, fv%x, n)
205
206 call col3(fu%x, u%x, w%x, n)
207 call test_filter%apply(fv, fu)
208 call sub2(lij(5)%x, fv%x, n)
209
210 call col3(fu%x, v%x, w%x, n)
211 call test_filter%apply(fv, fu)
212 call sub2(lij(6)%x, fv%x, n)
213
214 end subroutine compute_lij_cpu
215
224 subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
225 s_abs, test_filter, delta, n)
226 type(field_t), intent(inout) :: mij(6)
227 type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23, s_abs
228 type(elementwise_filter_t), intent(inout) :: test_filter
229 type(field_t), intent(in) :: delta
230 integer, intent(in) :: n
231
232 integer :: temp_indices(7)
233 type(field_t), pointer :: fs11, fs22, fs33, fs12, fs13, fs23, fs_abs
234 real(kind=rp) :: delta_ratio2 !test- to grid- filter ratio, squared
235 integer :: i
236 real(kind=rp) :: delta2
237
238 delta_ratio2 = ((test_filter%nx-1.0_rp)/(test_filter%nt-1.0_rp))**2
239
240 call neko_scratch_registry%request_field(fs11, temp_indices(1))
241 call neko_scratch_registry%request_field(fs22, temp_indices(2))
242 call neko_scratch_registry%request_field(fs33, temp_indices(3))
243 call neko_scratch_registry%request_field(fs12, temp_indices(4))
244 call neko_scratch_registry%request_field(fs13, temp_indices(5))
245 call neko_scratch_registry%request_field(fs23, temp_indices(6))
246 call neko_scratch_registry%request_field(fs_abs, temp_indices(7))
247 !! The first term:
248 !! _____ ____
249 !! (delta_test/delta)^2 s_abs*s_ij
250 call test_filter%apply(fs_abs, s_abs)
251
252 call test_filter%apply(fs11, s11)
253 call col3(mij(1)%x, fs_abs%x, fs11%x, n)
254 call cmult(mij(1)%x, delta_ratio2, n)
255
256 call test_filter%apply(fs22, s22)
257 call col3(mij(2)%x, fs_abs%x, fs22%x, n)
258 call cmult(mij(2)%x, delta_ratio2, n)
259
260 call test_filter%apply(fs33, s33)
261 call col3(mij(3)%x, fs_abs%x, fs33%x, n)
262 call cmult(mij(3)%x, delta_ratio2, n)
263
264 call test_filter%apply(fs12, s12)
265 call col3(mij(4)%x, fs_abs%x, fs12%x, n)
266 call cmult(mij(4)%x, delta_ratio2, n)
267
268 call test_filter%apply(fs13, s13)
269 call col3(mij(5)%x, fs_abs%x, fs13%x, n)
270 call cmult(mij(5)%x, delta_ratio2, n)
271
272 call test_filter%apply(fs23, s23)
273 call col3(mij(6)%x, fs_abs%x, fs23%x, n)
274 call cmult(mij(6)%x, delta_ratio2, n)
275
276 !! Substract the second term:
277 !! _____ ____ __________
278 !! (delta_test/delta)^2 s_abs*s_ij - s_abs*s_ij
279 !! fs11 and fs22 are used as work array
280 call col3(fs11%x, s_abs%x, s11%x, n)
281 call test_filter%apply(fs22, fs11)
282 call sub2(mij(1)%x, fs22%x, n)
283
284 call col3(fs11%x, s_abs%x, s22%x, n)
285 call test_filter%apply(fs22, fs11)
286 call sub2(mij(2)%x, fs22%x, n)
287
288 call col3(fs11%x, s_abs%x, s33%x, n)
289 call test_filter%apply(fs22, fs11)
290 call sub2(mij(3)%x, fs22%x, n)
291
292 call col3(fs11%x, s_abs%x, s12%x, n)
293 call test_filter%apply(fs22, fs11)
294 call sub2(mij(4)%x, fs22%x, n)
295
296 call col3(fs11%x, s_abs%x, s13%x, n)
297 call test_filter%apply(fs22, fs11)
298 call sub2(mij(5)%x, fs22%x, n)
299
300 call col3(fs11%x, s_abs%x, s23%x, n)
301 call test_filter%apply(fs22, fs11)
302 call sub2(mij(6)%x, fs22%x, n)
303
304 !! Lastly multiplied by delta^2
305 do concurrent(i = 1:n)
306 delta2 = delta%x(i,1,1,1)**2
307 mij(1)%x(i,1,1,1) = mij(1)%x(i,1,1,1) * delta2
308 mij(2)%x(i,1,1,1) = mij(2)%x(i,1,1,1) * delta2
309 mij(3)%x(i,1,1,1) = mij(3)%x(i,1,1,1) * delta2
310 mij(4)%x(i,1,1,1) = mij(4)%x(i,1,1,1) * delta2
311 mij(5)%x(i,1,1,1) = mij(5)%x(i,1,1,1) * delta2
312 mij(6)%x(i,1,1,1) = mij(6)%x(i,1,1,1) * delta2
313 end do
314
315 end subroutine compute_mij_cpu
316
323 subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
324 type(field_t), intent(inout) :: num, den
325 type(field_t), intent(in) :: lij(6), mij(6)
326 real(kind=rp), intent(in) :: alpha
327 integer, intent(in) :: n
328
329 real(kind=rp), dimension(n) :: num_curr, den_curr
330 integer :: i
331
332 do concurrent(i = 1:n)
333 num_curr(i) = mij(1)%x(i,1,1,1)*lij(1)%x(i,1,1,1) + &
334 mij(2)%x(i,1,1,1)*lij(2)%x(i,1,1,1) + &
335 mij(3)%x(i,1,1,1)*lij(3)%x(i,1,1,1) + &
336 2.0_rp*(mij(4)%x(i,1,1,1)*lij(4)%x(i,1,1,1) + &
337 mij(5)%x(i,1,1,1)*lij(5)%x(i,1,1,1) + &
338 mij(6)%x(i,1,1,1)*lij(6)%x(i,1,1,1))
339 den_curr(i) = mij(1)%x(i,1,1,1)*mij(1)%x(i,1,1,1) + &
340 mij(2)%x(i,1,1,1)*mij(2)%x(i,1,1,1) + &
341 mij(3)%x(i,1,1,1)*mij(3)%x(i,1,1,1) + &
342 2.0_rp*(mij(4)%x(i,1,1,1)*mij(4)%x(i,1,1,1) + &
343 mij(5)%x(i,1,1,1)*mij(5)%x(i,1,1,1) + &
344 mij(6)%x(i,1,1,1)*mij(6)%x(i,1,1,1))
345 end do
346
347 ! running average over time
348 do concurrent(i = 1:n)
349 num%x(i,1,1,1) = alpha * num%x(i,1,1,1) + (1.0_rp - alpha) * num_curr(i)
350 den%x(i,1,1,1) = alpha * den%x(i,1,1,1) + (1.0_rp - alpha) * den_curr(i)
351 end do
352
353 end subroutine compute_num_den_cpu
354
355end module dynamic_smagorinsky_cpu
356
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_lij_cpu
subroutine compute_lij_cpu(lij, u, v, w, test_filter, n)
Compute Germano Identity on the CPU.
Definition dynamic_smagorinsky_cpu.f90:160
dynamic_smagorinsky_cpu::compute_mij_cpu
subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, s_abs, test_filter, delta, n)
Compute M_ij on the CPU.
Definition dynamic_smagorinsky_cpu.f90:226
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:324
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:438
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:56
field::field_t
Definition field.f90:47
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