Neko 1.99.3
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, &
113 u%x, v%x, w%x, coef)
114
115 call coef%gs_h%op(s11%x, s11%dof%size(), gs_op_add)
116 call coef%gs_h%op(s22%x, s11%dof%size(), gs_op_add)
117 call coef%gs_h%op(s33%x, s11%dof%size(), gs_op_add)
118 call coef%gs_h%op(s12%x, s11%dof%size(), gs_op_add)
119 call coef%gs_h%op(s13%x, s11%dof%size(), gs_op_add)
120 call coef%gs_h%op(s23%x, s11%dof%size(), gs_op_add)
121
122 do concurrent(i = 1:u%dof%size())
123 s11%x(i,1,1,1) = s11%x(i,1,1,1) * coef%mult(i,1,1,1)
124 s22%x(i,1,1,1) = s22%x(i,1,1,1) * coef%mult(i,1,1,1)
125 s33%x(i,1,1,1) = s33%x(i,1,1,1) * coef%mult(i,1,1,1)
126 s12%x(i,1,1,1) = s12%x(i,1,1,1) * coef%mult(i,1,1,1)
127 s13%x(i,1,1,1) = s13%x(i,1,1,1) * coef%mult(i,1,1,1)
128 s23%x(i,1,1,1) = s23%x(i,1,1,1) * coef%mult(i,1,1,1)
129 end do
130
131 do concurrent(i = 1:u%dof%size())
132 s_abs%x(i,1,1,1) = sqrt(2.0_rp * (s11%x(i,1,1,1)*s11%x(i,1,1,1) + &
133 s22%x(i,1,1,1)*s22%x(i,1,1,1) + &
134 s33%x(i,1,1,1)*s33%x(i,1,1,1)) + &
135 4.0_rp * (s12%x(i,1,1,1)*s12%x(i,1,1,1) + &
136 s13%x(i,1,1,1)*s13%x(i,1,1,1) + &
137 s23%x(i,1,1,1)*s23%x(i,1,1,1)))
138 end do
139
140 call compute_lij_cpu(lij, u, v, w, test_filter, u%dof%size())
141 call compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
142 s_abs, test_filter, delta, u%dof%size())
143 call compute_num_den_cpu(num, den, lij, mij, alpha, u%dof%size())
144
145 do concurrent(i =1:u%dof%size())
146 if (den%x(i,1,1,1) .gt. 0.0_rp) then
147 c_dyn%x(i,1,1,1) = 0.5_rp * (num%x(i,1,1,1)/den%x(i,1,1,1))
148 else
149 c_dyn%x(i,1,1,1) = 0.0_rp
150 end if
151 c_dyn%x(i,1,1,1) = max(c_dyn%x(i,1,1,1),0.0_rp)
152 nut%x(i,1,1,1) = c_dyn%x(i,1,1,1) * delta%x(i,1,1,1)**2 &
153 * s_abs%x(i,1,1,1)
154 end do
155
156 call coef%gs_h%op(nut, gs_op_add)
157 call col2(nut%x, coef%mult, nut%dof%size())
158
159 call neko_scratch_registry%relinquish_field(temp_indices)
160 end subroutine dynamic_smagorinsky_compute_cpu
161
171 subroutine compute_lij_cpu(lij, u, v, w, test_filter, n)
172 type(field_t), intent(inout) :: lij(6)
173 type(field_t), pointer, intent(in) :: u, v, w
174 type(elementwise_filter_t), intent(inout) :: test_filter
175 integer, intent(in) :: n
176 integer :: i
178 integer :: temp_indices(3)
179 type(field_t), pointer :: fu, fv, fw
180
181 ! Use test filter for the velocity fields
182 call neko_scratch_registry%request_field(fu, temp_indices(1), .false.)
183 call neko_scratch_registry%request_field(fv, temp_indices(2), .false.)
184 call neko_scratch_registry%request_field(fw, temp_indices(3), .false.)
185 call test_filter%apply(fu, u)
186 call test_filter%apply(fv, v)
187 call test_filter%apply(fw, w)
188
189 !! The first term
190 do concurrent(i = 1:n)
191 lij(1)%x(i,1,1,1) = fu%x(i,1,1,1) * fu%x(i,1,1,1)
192 lij(2)%x(i,1,1,1) = fv%x(i,1,1,1) * fv%x(i,1,1,1)
193 lij(3)%x(i,1,1,1) = fw%x(i,1,1,1) * fw%x(i,1,1,1)
194 lij(4)%x(i,1,1,1) = fu%x(i,1,1,1) * fv%x(i,1,1,1)
195 lij(5)%x(i,1,1,1) = fu%x(i,1,1,1) * fw%x(i,1,1,1)
196 lij(6)%x(i,1,1,1) = fv%x(i,1,1,1) * fw%x(i,1,1,1)
197 end do
198
199 !! Subtract the second term:
200 !! use test filter for the cross terms
201 !! fu and fv are used as work array
202 call col3(fu%x, u%x, u%x, n)
203 call test_filter%apply(fv, fu)
204 call sub2(lij(1)%x, fv%x, n)
205
206 call col3(fu%x, v%x, v%x, n)
207 call test_filter%apply(fv, fu)
208 call sub2(lij(2)%x, fv%x, n)
209
210 call col3(fu%x, w%x, w%x, n)
211 call test_filter%apply(fv, fu)
212 call sub2(lij(3)%x, fv%x, n)
213
214 call col3(fu%x, u%x, v%x, n)
215 call test_filter%apply(fv, fu)
216 call sub2(lij(4)%x, fv%x, n)
217
218 call col3(fu%x, u%x, w%x, n)
219 call test_filter%apply(fv, fu)
220 call sub2(lij(5)%x, fv%x, n)
221
222 call col3(fu%x, v%x, w%x, n)
223 call test_filter%apply(fv, fu)
224 call sub2(lij(6)%x, fv%x, n)
225
226 end subroutine compute_lij_cpu
227
241 subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
242 s_abs, test_filter, delta, n)
243 type(field_t), intent(inout) :: mij(6)
244 type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23, s_abs
245 type(elementwise_filter_t), intent(inout) :: test_filter
246 type(field_t), intent(in) :: delta
247 integer, intent(in) :: n
248
249 integer :: temp_indices(7)
250 type(field_t), pointer :: fs11, fs22, fs33, fs12, fs13, fs23, fs_abs
251 real(kind=rp) :: delta_ratio2 !test- to grid- filter ratio, squared
252 integer :: i
253 real(kind=rp) :: delta2
254
255 delta_ratio2 = ((test_filter%nx-1.0_rp)/(test_filter%nt-1.0_rp))**2
256
257 call neko_scratch_registry%request_field(fs11, temp_indices(1), .false.)
258 call neko_scratch_registry%request_field(fs22, temp_indices(2), .false.)
259 call neko_scratch_registry%request_field(fs33, temp_indices(3), .false.)
260 call neko_scratch_registry%request_field(fs12, temp_indices(4), .false.)
261 call neko_scratch_registry%request_field(fs13, temp_indices(5), .false.)
262 call neko_scratch_registry%request_field(fs23, temp_indices(6), .false.)
263 call neko_scratch_registry%request_field(fs_abs, temp_indices(7), .false.)
264 !! The first term:
265 !! _____ ____
266 !! (delta_test/delta)^2 s_abs*s_ij
267 call test_filter%apply(fs_abs, s_abs)
268
269 call test_filter%apply(fs11, s11)
270 call col3(mij(1)%x, fs_abs%x, fs11%x, n)
271 call cmult(mij(1)%x, delta_ratio2, n)
272
273 call test_filter%apply(fs22, s22)
274 call col3(mij(2)%x, fs_abs%x, fs22%x, n)
275 call cmult(mij(2)%x, delta_ratio2, n)
276
277 call test_filter%apply(fs33, s33)
278 call col3(mij(3)%x, fs_abs%x, fs33%x, n)
279 call cmult(mij(3)%x, delta_ratio2, n)
280
281 call test_filter%apply(fs12, s12)
282 call col3(mij(4)%x, fs_abs%x, fs12%x, n)
283 call cmult(mij(4)%x, delta_ratio2, n)
284
285 call test_filter%apply(fs13, s13)
286 call col3(mij(5)%x, fs_abs%x, fs13%x, n)
287 call cmult(mij(5)%x, delta_ratio2, n)
288
289 call test_filter%apply(fs23, s23)
290 call col3(mij(6)%x, fs_abs%x, fs23%x, n)
291 call cmult(mij(6)%x, delta_ratio2, n)
292
293 !! Substract the second term:
294 !! _____ ____ __________
295 !! (delta_test/delta)^2 s_abs*s_ij - s_abs*s_ij
296 !! fs11 and fs22 are used as work array
297 call col3(fs11%x, s_abs%x, s11%x, n)
298 call test_filter%apply(fs22, fs11)
299 call sub2(mij(1)%x, fs22%x, n)
300
301 call col3(fs11%x, s_abs%x, s22%x, n)
302 call test_filter%apply(fs22, fs11)
303 call sub2(mij(2)%x, fs22%x, n)
304
305 call col3(fs11%x, s_abs%x, s33%x, n)
306 call test_filter%apply(fs22, fs11)
307 call sub2(mij(3)%x, fs22%x, n)
308
309 call col3(fs11%x, s_abs%x, s12%x, n)
310 call test_filter%apply(fs22, fs11)
311 call sub2(mij(4)%x, fs22%x, n)
312
313 call col3(fs11%x, s_abs%x, s13%x, n)
314 call test_filter%apply(fs22, fs11)
315 call sub2(mij(5)%x, fs22%x, n)
316
317 call col3(fs11%x, s_abs%x, s23%x, n)
318 call test_filter%apply(fs22, fs11)
319 call sub2(mij(6)%x, fs22%x, n)
320
321 !! Lastly multiplied by delta^2
322 do concurrent(i = 1:n)
323 delta2 = delta%x(i,1,1,1)**2
324 mij(1)%x(i,1,1,1) = mij(1)%x(i,1,1,1) * delta2
325 mij(2)%x(i,1,1,1) = mij(2)%x(i,1,1,1) * delta2
326 mij(3)%x(i,1,1,1) = mij(3)%x(i,1,1,1) * delta2
327 mij(4)%x(i,1,1,1) = mij(4)%x(i,1,1,1) * delta2
328 mij(5)%x(i,1,1,1) = mij(5)%x(i,1,1,1) * delta2
329 mij(6)%x(i,1,1,1) = mij(6)%x(i,1,1,1) * delta2
330 end do
331
332 end subroutine compute_mij_cpu
333
341 subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
342 type(field_t), intent(inout) :: num, den
343 type(field_t), intent(in) :: lij(6), mij(6)
344 real(kind=rp), intent(in) :: alpha
345 integer, intent(in) :: n
346
347 real(kind=rp), dimension(n) :: num_curr, den_curr
348 integer :: i
349
350 do concurrent(i = 1:n)
351 num_curr(i) = mij(1)%x(i,1,1,1)*lij(1)%x(i,1,1,1) + &
352 mij(2)%x(i,1,1,1)*lij(2)%x(i,1,1,1) + &
353 mij(3)%x(i,1,1,1)*lij(3)%x(i,1,1,1) + &
354 2.0_rp*(mij(4)%x(i,1,1,1)*lij(4)%x(i,1,1,1) + &
355 mij(5)%x(i,1,1,1)*lij(5)%x(i,1,1,1) + &
356 mij(6)%x(i,1,1,1)*lij(6)%x(i,1,1,1))
357 den_curr(i) = mij(1)%x(i,1,1,1)*mij(1)%x(i,1,1,1) + &
358 mij(2)%x(i,1,1,1)*mij(2)%x(i,1,1,1) + &
359 mij(3)%x(i,1,1,1)*mij(3)%x(i,1,1,1) + &
360 2.0_rp*(mij(4)%x(i,1,1,1)*mij(4)%x(i,1,1,1) + &
361 mij(5)%x(i,1,1,1)*mij(5)%x(i,1,1,1) + &
362 mij(6)%x(i,1,1,1)*mij(6)%x(i,1,1,1))
363 end do
364
365 ! running average over time
366 do concurrent(i = 1:n)
367 num%x(i,1,1,1) = alpha * num%x(i,1,1,1) + (1.0_rp - alpha) * num_curr(i)
368 den%x(i,1,1,1) = alpha * den%x(i,1,1,1) + (1.0_rp - alpha) * den_curr(i)
369 end do
370
371 end subroutine compute_num_den_cpu
372
373end module dynamic_smagorinsky_cpu
374
operators::strain_rate
Compute the strain rate tensor of a vector field.
Definition operators.f90:119
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:172
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:243
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:342
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:502
math::cadd
subroutine, public cadd(a, s, n)
Add a scalar to vector .
Definition math.f90:564
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:1044
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.f90:1059
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:946
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
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:144
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:138
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:63
elementwise_filter::elementwise_filter_t
Implements the elementwise filter for SEM.
Definition elementwise_filter.f90:56
field::field_t
Definition field.f90:48
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