Neko 1.99.1
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
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_field_registry%get_field_by_name("u_e")
94 v => neko_field_registry%get_field_by_name("v_e")
95 w => neko_field_registry%get_field_by_name("w_e")
96 else
97 u => neko_field_registry%get_field_by_name("u")
98 v => neko_field_registry%get_field_by_name("v")
99 w => neko_field_registry%get_field_by_name("w")
100 end if
101
102
103 call neko_scratch_registry%request_field(s11, temp_indices(1))
104 call neko_scratch_registry%request_field(s22, temp_indices(2))
105 call neko_scratch_registry%request_field(s33, temp_indices(3))
106 call neko_scratch_registry%request_field(s12, temp_indices(4))
107 call neko_scratch_registry%request_field(s13, temp_indices(5))
108 call neko_scratch_registry%request_field(s23, temp_indices(6))
109 call neko_scratch_registry%request_field(s_abs, temp_indices(7))
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
169 subroutine compute_lij_cpu(lij, u, v, w, test_filter, n)
170 type(field_t), intent(inout) :: lij(6)
171 type(field_t), pointer, intent(in) :: u, v, w
172 type(elementwise_filter_t), intent(inout) :: test_filter
173 integer, intent(in) :: n
174 integer :: i
176 integer :: temp_indices(3)
177 type(field_t), pointer :: fu, fv, fw
178
179 ! Use test filter for the velocity fields
180 call neko_scratch_registry%request_field(fu, temp_indices(1))
181 call neko_scratch_registry%request_field(fv, temp_indices(2))
182 call neko_scratch_registry%request_field(fw, temp_indices(3))
183 call test_filter%apply(fu, u)
184 call test_filter%apply(fv, v)
185 call test_filter%apply(fw, w)
186
187 !! The first term
188 do concurrent(i = 1:n)
189 lij(1)%x(i,1,1,1) = fu%x(i,1,1,1) * fu%x(i,1,1,1)
190 lij(2)%x(i,1,1,1) = fv%x(i,1,1,1) * fv%x(i,1,1,1)
191 lij(3)%x(i,1,1,1) = fw%x(i,1,1,1) * fw%x(i,1,1,1)
192 lij(4)%x(i,1,1,1) = fu%x(i,1,1,1) * fv%x(i,1,1,1)
193 lij(5)%x(i,1,1,1) = fu%x(i,1,1,1) * fw%x(i,1,1,1)
194 lij(6)%x(i,1,1,1) = fv%x(i,1,1,1) * fw%x(i,1,1,1)
195 end do
196
197 !! Subtract the second term:
198 !! use test filter for the cross terms
199 !! fu and fv are used as work array
200 call col3(fu%x, u%x, u%x, n)
201 call test_filter%apply(fv, fu)
202 call sub2(lij(1)%x, fv%x, n)
203
204 call col3(fu%x, v%x, v%x, n)
205 call test_filter%apply(fv, fu)
206 call sub2(lij(2)%x, fv%x, n)
207
208 call col3(fu%x, w%x, w%x, n)
209 call test_filter%apply(fv, fu)
210 call sub2(lij(3)%x, fv%x, n)
211
212 call col3(fu%x, u%x, v%x, n)
213 call test_filter%apply(fv, fu)
214 call sub2(lij(4)%x, fv%x, n)
215
216 call col3(fu%x, u%x, w%x, n)
217 call test_filter%apply(fv, fu)
218 call sub2(lij(5)%x, fv%x, n)
219
220 call col3(fu%x, v%x, w%x, n)
221 call test_filter%apply(fv, fu)
222 call sub2(lij(6)%x, fv%x, n)
223
224 end subroutine compute_lij_cpu
225
234 subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
235 s_abs, test_filter, delta, n)
236 type(field_t), intent(inout) :: mij(6)
237 type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23, s_abs
238 type(elementwise_filter_t), intent(inout) :: test_filter
239 type(field_t), intent(in) :: delta
240 integer, intent(in) :: n
241
242 integer :: temp_indices(7)
243 type(field_t), pointer :: fs11, fs22, fs33, fs12, fs13, fs23, fs_abs
244 real(kind=rp) :: delta_ratio2 !test- to grid- filter ratio, squared
245 integer :: i
246 real(kind=rp) :: delta2
247
248 delta_ratio2 = ((test_filter%nx-1.0_rp)/(test_filter%nt-1.0_rp))**2
249
250 call neko_scratch_registry%request_field(fs11, temp_indices(1))
251 call neko_scratch_registry%request_field(fs22, temp_indices(2))
252 call neko_scratch_registry%request_field(fs33, temp_indices(3))
253 call neko_scratch_registry%request_field(fs12, temp_indices(4))
254 call neko_scratch_registry%request_field(fs13, temp_indices(5))
255 call neko_scratch_registry%request_field(fs23, temp_indices(6))
256 call neko_scratch_registry%request_field(fs_abs, temp_indices(7))
257 !! The first term:
258 !! _____ ____
259 !! (delta_test/delta)^2 s_abs*s_ij
260 call test_filter%apply(fs_abs, s_abs)
261
262 call test_filter%apply(fs11, s11)
263 call col3(mij(1)%x, fs_abs%x, fs11%x, n)
264 call cmult(mij(1)%x, delta_ratio2, n)
265
266 call test_filter%apply(fs22, s22)
267 call col3(mij(2)%x, fs_abs%x, fs22%x, n)
268 call cmult(mij(2)%x, delta_ratio2, n)
269
270 call test_filter%apply(fs33, s33)
271 call col3(mij(3)%x, fs_abs%x, fs33%x, n)
272 call cmult(mij(3)%x, delta_ratio2, n)
273
274 call test_filter%apply(fs12, s12)
275 call col3(mij(4)%x, fs_abs%x, fs12%x, n)
276 call cmult(mij(4)%x, delta_ratio2, n)
277
278 call test_filter%apply(fs13, s13)
279 call col3(mij(5)%x, fs_abs%x, fs13%x, n)
280 call cmult(mij(5)%x, delta_ratio2, n)
281
282 call test_filter%apply(fs23, s23)
283 call col3(mij(6)%x, fs_abs%x, fs23%x, n)
284 call cmult(mij(6)%x, delta_ratio2, n)
285
286 !! Substract the second term:
287 !! _____ ____ __________
288 !! (delta_test/delta)^2 s_abs*s_ij - s_abs*s_ij
289 !! fs11 and fs22 are used as work array
290 call col3(fs11%x, s_abs%x, s11%x, n)
291 call test_filter%apply(fs22, fs11)
292 call sub2(mij(1)%x, fs22%x, n)
293
294 call col3(fs11%x, s_abs%x, s22%x, n)
295 call test_filter%apply(fs22, fs11)
296 call sub2(mij(2)%x, fs22%x, n)
297
298 call col3(fs11%x, s_abs%x, s33%x, n)
299 call test_filter%apply(fs22, fs11)
300 call sub2(mij(3)%x, fs22%x, n)
301
302 call col3(fs11%x, s_abs%x, s12%x, n)
303 call test_filter%apply(fs22, fs11)
304 call sub2(mij(4)%x, fs22%x, n)
305
306 call col3(fs11%x, s_abs%x, s13%x, n)
307 call test_filter%apply(fs22, fs11)
308 call sub2(mij(5)%x, fs22%x, n)
309
310 call col3(fs11%x, s_abs%x, s23%x, n)
311 call test_filter%apply(fs22, fs11)
312 call sub2(mij(6)%x, fs22%x, n)
313
314 !! Lastly multiplied by delta^2
315 do concurrent(i = 1:n)
316 delta2 = delta%x(i,1,1,1)**2
317 mij(1)%x(i,1,1,1) = mij(1)%x(i,1,1,1) * delta2
318 mij(2)%x(i,1,1,1) = mij(2)%x(i,1,1,1) * delta2
319 mij(3)%x(i,1,1,1) = mij(3)%x(i,1,1,1) * delta2
320 mij(4)%x(i,1,1,1) = mij(4)%x(i,1,1,1) * delta2
321 mij(5)%x(i,1,1,1) = mij(5)%x(i,1,1,1) * delta2
322 mij(6)%x(i,1,1,1) = mij(6)%x(i,1,1,1) * delta2
323 end do
324
325 end subroutine compute_mij_cpu
326
333 subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
334 type(field_t), intent(inout) :: num, den
335 type(field_t), intent(in) :: lij(6), mij(6)
336 real(kind=rp), intent(in) :: alpha
337 integer, intent(in) :: n
338
339 real(kind=rp), dimension(n) :: num_curr, den_curr
340 integer :: i
341
342 do concurrent(i = 1:n)
343 num_curr(i) = mij(1)%x(i,1,1,1)*lij(1)%x(i,1,1,1) + &
344 mij(2)%x(i,1,1,1)*lij(2)%x(i,1,1,1) + &
345 mij(3)%x(i,1,1,1)*lij(3)%x(i,1,1,1) + &
346 2.0_rp*(mij(4)%x(i,1,1,1)*lij(4)%x(i,1,1,1) + &
347 mij(5)%x(i,1,1,1)*lij(5)%x(i,1,1,1) + &
348 mij(6)%x(i,1,1,1)*lij(6)%x(i,1,1,1))
349 den_curr(i) = mij(1)%x(i,1,1,1)*mij(1)%x(i,1,1,1) + &
350 mij(2)%x(i,1,1,1)*mij(2)%x(i,1,1,1) + &
351 mij(3)%x(i,1,1,1)*mij(3)%x(i,1,1,1) + &
352 2.0_rp*(mij(4)%x(i,1,1,1)*mij(4)%x(i,1,1,1) + &
353 mij(5)%x(i,1,1,1)*mij(5)%x(i,1,1,1) + &
354 mij(6)%x(i,1,1,1)*mij(6)%x(i,1,1,1))
355 end do
356
357 ! running average over time
358 do concurrent(i = 1:n)
359 num%x(i,1,1,1) = alpha * num%x(i,1,1,1) + (1.0_rp - alpha) * num_curr(i)
360 den%x(i,1,1,1) = alpha * den%x(i,1,1,1) + (1.0_rp - alpha) * den_curr(i)
361 end do
362
363 end subroutine compute_num_den_cpu
364
365end module dynamic_smagorinsky_cpu
366
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:170
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:236
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:334
elementwise_filter
Implements elementwise_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:94
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:489
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:84
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 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