Neko 1.99.2
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 registry, only : neko_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
65 subroutine dynamic_smagorinsky_compute_cpu(if_ext, t, tstep, coef, nut, delta, &
66 c_dyn, test_filter, mij, lij, num, den)
67 logical, intent(in) :: if_ext
68 real(kind=rp), intent(in) :: t
69 integer, intent(in) :: tstep
70 type(coef_t), intent(in) :: coef
71 type(field_t), intent(inout) :: nut
72 type(field_t), intent(in) :: delta
73 type(field_t), intent(inout) :: c_dyn
74 type(elementwise_filter_t), intent(inout) :: test_filter
75 type(field_t), intent(inout) :: mij(6), lij(6)
76 type(field_t), intent(inout) :: num, den
77
78 type(field_t), pointer :: u, v, w
79 type(field_t) :: c_dyn_curr
80 ! the strain rate tensor
81 type(field_t), pointer :: s11, s22, s33, s12, s13, s23, s_abs
82 real(kind=rp) :: alpha ! running averaging coefficient
83 integer :: temp_indices(7)
84 integer :: i
85
86 if (tstep .eq. 1) then
87 alpha = 1.0_rp
88 else
89 alpha = 0.9_rp
90 end if
91
92 if (if_ext .eqv. .true.) then
93 u => neko_registry%get_field_by_name("u_e")
94 v => neko_registry%get_field_by_name("v_e")
95 w => neko_registry%get_field_by_name("w_e")
96 else
97 u => neko_registry%get_field_by_name("u")
98 v => neko_registry%get_field_by_name("v")
99 w => neko_registry%get_field_by_name("w")
100 end if
101
102
103 call neko_scratch_registry%request_field(s11, temp_indices(1), .false.)
104 call neko_scratch_registry%request_field(s22, temp_indices(2), .false.)
105 call neko_scratch_registry%request_field(s33, temp_indices(3), .false.)
106 call neko_scratch_registry%request_field(s12, temp_indices(4), .false.)
107 call neko_scratch_registry%request_field(s13, temp_indices(5), .false.)
108 call neko_scratch_registry%request_field(s23, temp_indices(6), .false.)
109 call neko_scratch_registry%request_field(s_abs, temp_indices(7), .false.)
110
111 ! Compute the strain rate tensor
112 call strain_rate(s11%x, s22%x, s33%x, s12%x, s13%x, s23%x, u, v, w, coef)
113
114 call coef%gs_h%op(s11%x, s11%dof%size(), gs_op_add)
115 call coef%gs_h%op(s22%x, s11%dof%size(), gs_op_add)
116 call coef%gs_h%op(s33%x, s11%dof%size(), gs_op_add)
117 call coef%gs_h%op(s12%x, s11%dof%size(), gs_op_add)
118 call coef%gs_h%op(s13%x, s11%dof%size(), gs_op_add)
119 call coef%gs_h%op(s23%x, s11%dof%size(), gs_op_add)
120
121 do concurrent(i = 1:u%dof%size())
122 s11%x(i,1,1,1) = s11%x(i,1,1,1) * coef%mult(i,1,1,1)
123 s22%x(i,1,1,1) = s22%x(i,1,1,1) * coef%mult(i,1,1,1)
124 s33%x(i,1,1,1) = s33%x(i,1,1,1) * coef%mult(i,1,1,1)
125 s12%x(i,1,1,1) = s12%x(i,1,1,1) * coef%mult(i,1,1,1)
126 s13%x(i,1,1,1) = s13%x(i,1,1,1) * coef%mult(i,1,1,1)
127 s23%x(i,1,1,1) = s23%x(i,1,1,1) * coef%mult(i,1,1,1)
128 end do
129
130 do concurrent(i = 1:u%dof%size())
131 s_abs%x(i,1,1,1) = sqrt(2.0_rp * (s11%x(i,1,1,1)*s11%x(i,1,1,1) + &
132 s22%x(i,1,1,1)*s22%x(i,1,1,1) + &
133 s33%x(i,1,1,1)*s33%x(i,1,1,1)) + &
134 4.0_rp * (s12%x(i,1,1,1)*s12%x(i,1,1,1) + &
135 s13%x(i,1,1,1)*s13%x(i,1,1,1) + &
136 s23%x(i,1,1,1)*s23%x(i,1,1,1)))
137 end do
138
139 call compute_lij_cpu(lij, u, v, w, test_filter, u%dof%size())
140 call compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
141 s_abs, test_filter, delta, u%dof%size())
142 call compute_num_den_cpu(num, den, lij, mij, alpha, u%dof%size())
143
144 do concurrent(i =1:u%dof%size())
145 if (den%x(i,1,1,1) .gt. 0.0_rp) then
146 c_dyn%x(i,1,1,1) = 0.5_rp * (num%x(i,1,1,1)/den%x(i,1,1,1))
147 else
148 c_dyn%x(i,1,1,1) = 0.0_rp
149 end if
150 c_dyn%x(i,1,1,1) = max(c_dyn%x(i,1,1,1),0.0_rp)
151 nut%x(i,1,1,1) = c_dyn%x(i,1,1,1) * delta%x(i,1,1,1)**2 &
152 * s_abs%x(i,1,1,1)
153 end do
154
155 call coef%gs_h%op(nut, gs_op_add)
156 call col2(nut%x, coef%mult, nut%dof%size())
157
158 call neko_scratch_registry%relinquish_field(temp_indices)
159 end subroutine dynamic_smagorinsky_compute_cpu
160
170 subroutine compute_lij_cpu(lij, u, v, w, test_filter, n)
171 type(field_t), intent(inout) :: lij(6)
172 type(field_t), pointer, intent(in) :: u, v, w
173 type(elementwise_filter_t), intent(inout) :: test_filter
174 integer, intent(in) :: n
175 integer :: i
177 integer :: temp_indices(3)
178 type(field_t), pointer :: fu, fv, fw
179
180 ! Use test filter for the velocity fields
181 call neko_scratch_registry%request_field(fu, temp_indices(1), .false.)
182 call neko_scratch_registry%request_field(fv, temp_indices(2), .false.)
183 call neko_scratch_registry%request_field(fw, temp_indices(3), .false.)
184 call test_filter%apply(fu, u)
185 call test_filter%apply(fv, v)
186 call test_filter%apply(fw, w)
187
188 !! The first term
189 do concurrent(i = 1:n)
190 lij(1)%x(i,1,1,1) = fu%x(i,1,1,1) * fu%x(i,1,1,1)
191 lij(2)%x(i,1,1,1) = fv%x(i,1,1,1) * fv%x(i,1,1,1)
192 lij(3)%x(i,1,1,1) = fw%x(i,1,1,1) * fw%x(i,1,1,1)
193 lij(4)%x(i,1,1,1) = fu%x(i,1,1,1) * fv%x(i,1,1,1)
194 lij(5)%x(i,1,1,1) = fu%x(i,1,1,1) * fw%x(i,1,1,1)
195 lij(6)%x(i,1,1,1) = fv%x(i,1,1,1) * fw%x(i,1,1,1)
196 end do
197
198 !! Subtract the second term:
199 !! use test filter for the cross terms
200 !! fu and fv are used as work array
201 call col3(fu%x, u%x, u%x, n)
202 call test_filter%apply(fv, fu)
203 call sub2(lij(1)%x, fv%x, n)
204
205 call col3(fu%x, v%x, v%x, n)
206 call test_filter%apply(fv, fu)
207 call sub2(lij(2)%x, fv%x, n)
208
209 call col3(fu%x, w%x, w%x, n)
210 call test_filter%apply(fv, fu)
211 call sub2(lij(3)%x, fv%x, n)
212
213 call col3(fu%x, u%x, v%x, n)
214 call test_filter%apply(fv, fu)
215 call sub2(lij(4)%x, fv%x, n)
216
217 call col3(fu%x, u%x, w%x, n)
218 call test_filter%apply(fv, fu)
219 call sub2(lij(5)%x, fv%x, n)
220
221 call col3(fu%x, v%x, w%x, n)
222 call test_filter%apply(fv, fu)
223 call sub2(lij(6)%x, fv%x, n)
224
225 end subroutine compute_lij_cpu
226
240 subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
241 s_abs, test_filter, delta, n)
242 type(field_t), intent(inout) :: mij(6)
243 type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23, s_abs
244 type(elementwise_filter_t), intent(inout) :: test_filter
245 type(field_t), intent(in) :: delta
246 integer, intent(in) :: n
247
248 integer :: temp_indices(7)
249 type(field_t), pointer :: fs11, fs22, fs33, fs12, fs13, fs23, fs_abs
250 real(kind=rp) :: delta_ratio2 !test- to grid- filter ratio, squared
251 integer :: i
252 real(kind=rp) :: delta2
253
254 delta_ratio2 = ((test_filter%nx-1.0_rp)/(test_filter%nt-1.0_rp))**2
255
256 call neko_scratch_registry%request_field(fs11, temp_indices(1), .false.)
257 call neko_scratch_registry%request_field(fs22, temp_indices(2), .false.)
258 call neko_scratch_registry%request_field(fs33, temp_indices(3), .false.)
259 call neko_scratch_registry%request_field(fs12, temp_indices(4), .false.)
260 call neko_scratch_registry%request_field(fs13, temp_indices(5), .false.)
261 call neko_scratch_registry%request_field(fs23, temp_indices(6), .false.)
262 call neko_scratch_registry%request_field(fs_abs, temp_indices(7), .false.)
263 !! The first term:
264 !! _____ ____
265 !! (delta_test/delta)^2 s_abs*s_ij
266 call test_filter%apply(fs_abs, s_abs)
267
268 call test_filter%apply(fs11, s11)
269 call col3(mij(1)%x, fs_abs%x, fs11%x, n)
270 call cmult(mij(1)%x, delta_ratio2, n)
271
272 call test_filter%apply(fs22, s22)
273 call col3(mij(2)%x, fs_abs%x, fs22%x, n)
274 call cmult(mij(2)%x, delta_ratio2, n)
275
276 call test_filter%apply(fs33, s33)
277 call col3(mij(3)%x, fs_abs%x, fs33%x, n)
278 call cmult(mij(3)%x, delta_ratio2, n)
279
280 call test_filter%apply(fs12, s12)
281 call col3(mij(4)%x, fs_abs%x, fs12%x, n)
282 call cmult(mij(4)%x, delta_ratio2, n)
283
284 call test_filter%apply(fs13, s13)
285 call col3(mij(5)%x, fs_abs%x, fs13%x, n)
286 call cmult(mij(5)%x, delta_ratio2, n)
287
288 call test_filter%apply(fs23, s23)
289 call col3(mij(6)%x, fs_abs%x, fs23%x, n)
290 call cmult(mij(6)%x, delta_ratio2, n)
291
292 !! Substract the second term:
293 !! _____ ____ __________
294 !! (delta_test/delta)^2 s_abs*s_ij - s_abs*s_ij
295 !! fs11 and fs22 are used as work array
296 call col3(fs11%x, s_abs%x, s11%x, n)
297 call test_filter%apply(fs22, fs11)
298 call sub2(mij(1)%x, fs22%x, n)
299
300 call col3(fs11%x, s_abs%x, s22%x, n)
301 call test_filter%apply(fs22, fs11)
302 call sub2(mij(2)%x, fs22%x, n)
303
304 call col3(fs11%x, s_abs%x, s33%x, n)
305 call test_filter%apply(fs22, fs11)
306 call sub2(mij(3)%x, fs22%x, n)
307
308 call col3(fs11%x, s_abs%x, s12%x, n)
309 call test_filter%apply(fs22, fs11)
310 call sub2(mij(4)%x, fs22%x, n)
311
312 call col3(fs11%x, s_abs%x, s13%x, n)
313 call test_filter%apply(fs22, fs11)
314 call sub2(mij(5)%x, fs22%x, n)
315
316 call col3(fs11%x, s_abs%x, s23%x, n)
317 call test_filter%apply(fs22, fs11)
318 call sub2(mij(6)%x, fs22%x, n)
319
320 !! Lastly multiplied by delta^2
321 do concurrent(i = 1:n)
322 delta2 = delta%x(i,1,1,1)**2
323 mij(1)%x(i,1,1,1) = mij(1)%x(i,1,1,1) * delta2
324 mij(2)%x(i,1,1,1) = mij(2)%x(i,1,1,1) * delta2
325 mij(3)%x(i,1,1,1) = mij(3)%x(i,1,1,1) * delta2
326 mij(4)%x(i,1,1,1) = mij(4)%x(i,1,1,1) * delta2
327 mij(5)%x(i,1,1,1) = mij(5)%x(i,1,1,1) * delta2
328 mij(6)%x(i,1,1,1) = mij(6)%x(i,1,1,1) * delta2
329 end do
330
331 end subroutine compute_mij_cpu
332
340 subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
341 type(field_t), intent(inout) :: num, den
342 type(field_t), intent(in) :: lij(6), mij(6)
343 real(kind=rp), intent(in) :: alpha
344 integer, intent(in) :: n
345
346 real(kind=rp), dimension(n) :: num_curr, den_curr
347 integer :: i
348
349 do concurrent(i = 1:n)
350 num_curr(i) = mij(1)%x(i,1,1,1)*lij(1)%x(i,1,1,1) + &
351 mij(2)%x(i,1,1,1)*lij(2)%x(i,1,1,1) + &
352 mij(3)%x(i,1,1,1)*lij(3)%x(i,1,1,1) + &
353 2.0_rp*(mij(4)%x(i,1,1,1)*lij(4)%x(i,1,1,1) + &
354 mij(5)%x(i,1,1,1)*lij(5)%x(i,1,1,1) + &
355 mij(6)%x(i,1,1,1)*lij(6)%x(i,1,1,1))
356 den_curr(i) = mij(1)%x(i,1,1,1)*mij(1)%x(i,1,1,1) + &
357 mij(2)%x(i,1,1,1)*mij(2)%x(i,1,1,1) + &
358 mij(3)%x(i,1,1,1)*mij(3)%x(i,1,1,1) + &
359 2.0_rp*(mij(4)%x(i,1,1,1)*mij(4)%x(i,1,1,1) + &
360 mij(5)%x(i,1,1,1)*mij(5)%x(i,1,1,1) + &
361 mij(6)%x(i,1,1,1)*mij(6)%x(i,1,1,1))
362 end do
363
364 ! running average over time
365 do concurrent(i = 1:n)
366 num%x(i,1,1,1) = alpha * num%x(i,1,1,1) + (1.0_rp - alpha) * num_curr(i)
367 den%x(i,1,1,1) = alpha * den%x(i,1,1,1) + (1.0_rp - alpha) * den_curr(i)
368 end do
369
370 end subroutine compute_num_den_cpu
371
372end module dynamic_smagorinsky_cpu
373
coefs
Coefficients.
Definition coef.f90:34
dynamic_smagorinsky_cpu
Implements the CPU kernel for the dynamic_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:171
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:242
dynamic_smagorinsky_cpu::dynamic_smagorinsky_compute_cpu
subroutine, public dynamic_smagorinsky_compute_cpu(if_ext, t, tstep, coef, nut, delta, c_dyn, test_filter, mij, lij, num, den)
Compute eddy viscosity on the CPU.
Definition dynamic_smagorinsky_cpu.f90:67
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:341
elementwise_filter
Implements elementwise_filter_t.
Definition elementwise_filter.f90:35
field_list
Definition field_list.f90:1
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:411
math::cadd
subroutine, public cadd(a, s, n)
Add a scalar to vector .
Definition math.f90:462
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:854
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:867
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:768
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:497
registry
Defines a registry for storing solution fields.
Definition registry.f90:34
registry::neko_registry
type(registry_t), target, public neko_registry
Global field registry.
Definition registry.f90:128
scratch_registry
Defines a registry for storing and requesting temporary objects This can be used when you have a func...
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:117
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:56
elementwise_filter::elementwise_filter_t
Implements the elementwise 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:15
max
#define max(a, b)
Definition tensor.cu:40