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 !OCL NORECURRENCE, NOVREC, NOALIAS
123 !DIR$ CONCURRENT
124 !GCC$ ivdep
125 !$omp parallel do
126 do i = 1, u%dof%size()
127 s11%x(i,1,1,1) = s11%x(i,1,1,1) * coef%mult(i,1,1,1)
128 s22%x(i,1,1,1) = s22%x(i,1,1,1) * coef%mult(i,1,1,1)
129 s33%x(i,1,1,1) = s33%x(i,1,1,1) * coef%mult(i,1,1,1)
130 s12%x(i,1,1,1) = s12%x(i,1,1,1) * coef%mult(i,1,1,1)
131 s13%x(i,1,1,1) = s13%x(i,1,1,1) * coef%mult(i,1,1,1)
132 s23%x(i,1,1,1) = s23%x(i,1,1,1) * coef%mult(i,1,1,1)
133 end do
134 !$omp end parallel do
135
136 !OCL NORECURRENCE, NOVREC, NOALIAS
137 !DIR$ CONCURRENT
138 !GCC$ ivdep
139 !$omp parallel do
140 do i = 1, u%dof%size()
141 s_abs%x(i,1,1,1) = sqrt(2.0_rp * (s11%x(i,1,1,1)*s11%x(i,1,1,1) + &
142 s22%x(i,1,1,1)*s22%x(i,1,1,1) + &
143 s33%x(i,1,1,1)*s33%x(i,1,1,1)) + &
144 4.0_rp * (s12%x(i,1,1,1)*s12%x(i,1,1,1) + &
145 s13%x(i,1,1,1)*s13%x(i,1,1,1) + &
146 s23%x(i,1,1,1)*s23%x(i,1,1,1)))
147 end do
148 !$omp end parallel do
149
150 call compute_lij_cpu(lij, u, v, w, test_filter, u%dof%size())
151 call compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
152 s_abs, test_filter, delta, u%dof%size())
153 call compute_num_den_cpu(num, den, lij, mij, alpha, u%dof%size())
154
155 !OCL NORECURRENCE, NOVREC, NOALIAS
156 !DIR$ CONCURRENT
157 !GCC$ ivdep
158 !$omp parallel do
159 do i = 1, u%dof%size()
160 if (den%x(i,1,1,1) .gt. 0.0_rp) then
161 c_dyn%x(i,1,1,1) = 0.5_rp * (num%x(i,1,1,1)/den%x(i,1,1,1))
162 else
163 c_dyn%x(i,1,1,1) = 0.0_rp
164 end if
165 c_dyn%x(i,1,1,1) = max(c_dyn%x(i,1,1,1),0.0_rp)
166 nut%x(i,1,1,1) = c_dyn%x(i,1,1,1) * delta%x(i,1,1,1)**2 &
167 * s_abs%x(i,1,1,1)
168 end do
169 !$omp end parallel do
170
171 call coef%gs_h%op(nut, gs_op_add)
172 call col2(nut%x, coef%mult, nut%dof%size())
173
174 call neko_scratch_registry%relinquish_field(temp_indices)
175 end subroutine dynamic_smagorinsky_compute_cpu
176
186 subroutine compute_lij_cpu(lij, u, v, w, test_filter, n)
187 type(field_t), intent(inout) :: lij(6)
188 type(field_t), pointer, intent(in) :: u, v, w
189 type(elementwise_filter_t), intent(inout) :: test_filter
190 integer, intent(in) :: n
191 integer :: i
193 integer :: temp_indices(3)
194 type(field_t), pointer :: fu, fv, fw
195
196 ! Use test filter for the velocity fields
197 call neko_scratch_registry%request_field(fu, temp_indices(1), .false.)
198 call neko_scratch_registry%request_field(fv, temp_indices(2), .false.)
199 call neko_scratch_registry%request_field(fw, temp_indices(3), .false.)
200 call test_filter%apply(fu, u)
201 call test_filter%apply(fv, v)
202 call test_filter%apply(fw, w)
203
204 !! The first term
205 !OCL NORECURRENCE, NOVREC, NOALIAS
206 !DIR$ CONCURRENT
207 !GCC$ ivdep
208 !$omp parallel do
209 do i = 1, n
210 lij(1)%x(i,1,1,1) = fu%x(i,1,1,1) * fu%x(i,1,1,1)
211 lij(2)%x(i,1,1,1) = fv%x(i,1,1,1) * fv%x(i,1,1,1)
212 lij(3)%x(i,1,1,1) = fw%x(i,1,1,1) * fw%x(i,1,1,1)
213 lij(4)%x(i,1,1,1) = fu%x(i,1,1,1) * fv%x(i,1,1,1)
214 lij(5)%x(i,1,1,1) = fu%x(i,1,1,1) * fw%x(i,1,1,1)
215 lij(6)%x(i,1,1,1) = fv%x(i,1,1,1) * fw%x(i,1,1,1)
216 end do
217 !$omp end parallel do
218
219 !! Subtract the second term:
220 !! use test filter for the cross terms
221 !! fu and fv are used as work array
222 call col3(fu%x, u%x, u%x, n)
223 call test_filter%apply(fv, fu)
224 call sub2(lij(1)%x, fv%x, n)
225
226 call col3(fu%x, v%x, v%x, n)
227 call test_filter%apply(fv, fu)
228 call sub2(lij(2)%x, fv%x, n)
229
230 call col3(fu%x, w%x, w%x, n)
231 call test_filter%apply(fv, fu)
232 call sub2(lij(3)%x, fv%x, n)
233
234 call col3(fu%x, u%x, v%x, n)
235 call test_filter%apply(fv, fu)
236 call sub2(lij(4)%x, fv%x, n)
237
238 call col3(fu%x, u%x, w%x, n)
239 call test_filter%apply(fv, fu)
240 call sub2(lij(5)%x, fv%x, n)
241
242 call col3(fu%x, v%x, w%x, n)
243 call test_filter%apply(fv, fu)
244 call sub2(lij(6)%x, fv%x, n)
245
246 end subroutine compute_lij_cpu
247
261 subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
262 s_abs, test_filter, delta, n)
263 type(field_t), intent(inout) :: mij(6)
264 type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23, s_abs
265 type(elementwise_filter_t), intent(inout) :: test_filter
266 type(field_t), intent(in) :: delta
267 integer, intent(in) :: n
268
269 integer :: temp_indices(7)
270 type(field_t), pointer :: fs11, fs22, fs33, fs12, fs13, fs23, fs_abs
271 real(kind=rp) :: delta_ratio2 !test- to grid- filter ratio, squared
272 integer :: i
273 real(kind=rp) :: delta2
274
275 delta_ratio2 = ((test_filter%nx-1.0_rp)/(test_filter%nt-1.0_rp))**2
276
277 call neko_scratch_registry%request_field(fs11, temp_indices(1), .false.)
278 call neko_scratch_registry%request_field(fs22, temp_indices(2), .false.)
279 call neko_scratch_registry%request_field(fs33, temp_indices(3), .false.)
280 call neko_scratch_registry%request_field(fs12, temp_indices(4), .false.)
281 call neko_scratch_registry%request_field(fs13, temp_indices(5), .false.)
282 call neko_scratch_registry%request_field(fs23, temp_indices(6), .false.)
283 call neko_scratch_registry%request_field(fs_abs, temp_indices(7), .false.)
284 !! The first term:
285 !! _____ ____
286 !! (delta_test/delta)^2 s_abs*s_ij
287 call test_filter%apply(fs_abs, s_abs)
288
289 call test_filter%apply(fs11, s11)
290 call col3(mij(1)%x, fs_abs%x, fs11%x, n)
291 call cmult(mij(1)%x, delta_ratio2, n)
292
293 call test_filter%apply(fs22, s22)
294 call col3(mij(2)%x, fs_abs%x, fs22%x, n)
295 call cmult(mij(2)%x, delta_ratio2, n)
296
297 call test_filter%apply(fs33, s33)
298 call col3(mij(3)%x, fs_abs%x, fs33%x, n)
299 call cmult(mij(3)%x, delta_ratio2, n)
300
301 call test_filter%apply(fs12, s12)
302 call col3(mij(4)%x, fs_abs%x, fs12%x, n)
303 call cmult(mij(4)%x, delta_ratio2, n)
304
305 call test_filter%apply(fs13, s13)
306 call col3(mij(5)%x, fs_abs%x, fs13%x, n)
307 call cmult(mij(5)%x, delta_ratio2, n)
308
309 call test_filter%apply(fs23, s23)
310 call col3(mij(6)%x, fs_abs%x, fs23%x, n)
311 call cmult(mij(6)%x, delta_ratio2, n)
312
313 !! Substract the second term:
314 !! _____ ____ __________
315 !! (delta_test/delta)^2 s_abs*s_ij - s_abs*s_ij
316 !! fs11 and fs22 are used as work array
317 call col3(fs11%x, s_abs%x, s11%x, n)
318 call test_filter%apply(fs22, fs11)
319 call sub2(mij(1)%x, fs22%x, n)
320
321 call col3(fs11%x, s_abs%x, s22%x, n)
322 call test_filter%apply(fs22, fs11)
323 call sub2(mij(2)%x, fs22%x, n)
324
325 call col3(fs11%x, s_abs%x, s33%x, n)
326 call test_filter%apply(fs22, fs11)
327 call sub2(mij(3)%x, fs22%x, n)
328
329 call col3(fs11%x, s_abs%x, s12%x, n)
330 call test_filter%apply(fs22, fs11)
331 call sub2(mij(4)%x, fs22%x, n)
332
333 call col3(fs11%x, s_abs%x, s13%x, n)
334 call test_filter%apply(fs22, fs11)
335 call sub2(mij(5)%x, fs22%x, n)
336
337 call col3(fs11%x, s_abs%x, s23%x, n)
338 call test_filter%apply(fs22, fs11)
339 call sub2(mij(6)%x, fs22%x, n)
340
341 !! Lastly multiplied by delta^2
342 !OCL NORECURRENCE, NOVREC, NOALIAS
343 !DIR$ CONCURRENT
344 !GCC$ ivdep
345 !$omp parallel do private(i, delta2)
346 do i = 1, n
347 delta2 = delta%x(i,1,1,1)**2
348 mij(1)%x(i,1,1,1) = mij(1)%x(i,1,1,1) * delta2
349 mij(2)%x(i,1,1,1) = mij(2)%x(i,1,1,1) * delta2
350 mij(3)%x(i,1,1,1) = mij(3)%x(i,1,1,1) * delta2
351 mij(4)%x(i,1,1,1) = mij(4)%x(i,1,1,1) * delta2
352 mij(5)%x(i,1,1,1) = mij(5)%x(i,1,1,1) * delta2
353 mij(6)%x(i,1,1,1) = mij(6)%x(i,1,1,1) * delta2
354 end do
355 !$omp end parallel do
356
357 end subroutine compute_mij_cpu
358
366 subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
367 type(field_t), intent(inout) :: num, den
368 type(field_t), intent(in) :: lij(6), mij(6)
369 real(kind=rp), intent(in) :: alpha
370 integer, intent(in) :: n
371
372 real(kind=rp), dimension(n) :: num_curr, den_curr
373 integer :: i
374
375 !OCL NORECURRENCE, NOVREC, NOALIAS
376 !DIR$ CONCURRENT
377 !GCC$ ivdep
378 !$omp parallel do
379 do i = 1, n
380 num_curr(i) = mij(1)%x(i,1,1,1)*lij(1)%x(i,1,1,1) + &
381 mij(2)%x(i,1,1,1)*lij(2)%x(i,1,1,1) + &
382 mij(3)%x(i,1,1,1)*lij(3)%x(i,1,1,1) + &
383 2.0_rp*(mij(4)%x(i,1,1,1)*lij(4)%x(i,1,1,1) + &
384 mij(5)%x(i,1,1,1)*lij(5)%x(i,1,1,1) + &
385 mij(6)%x(i,1,1,1)*lij(6)%x(i,1,1,1))
386 den_curr(i) = mij(1)%x(i,1,1,1)*mij(1)%x(i,1,1,1) + &
387 mij(2)%x(i,1,1,1)*mij(2)%x(i,1,1,1) + &
388 mij(3)%x(i,1,1,1)*mij(3)%x(i,1,1,1) + &
389 2.0_rp*(mij(4)%x(i,1,1,1)*mij(4)%x(i,1,1,1) + &
390 mij(5)%x(i,1,1,1)*mij(5)%x(i,1,1,1) + &
391 mij(6)%x(i,1,1,1)*mij(6)%x(i,1,1,1))
392 end do
393 !$omp end parallel do
394
395 ! running average over time
396 !OCL NORECURRENCE, NOVREC, NOALIAS
397 !DIR$ CONCURRENT
398 !GCC$ ivdep
399 !$omp parallel do
400 do i = 1, n
401 num%x(i,1,1,1) = alpha * num%x(i,1,1,1) + (1.0_rp - alpha) * num_curr(i)
402 den%x(i,1,1,1) = alpha * den%x(i,1,1,1) + (1.0_rp - alpha) * den_curr(i)
403 end do
404 !$omp end parallel do
405
406 end subroutine compute_num_den_cpu
407
408end module dynamic_smagorinsky_cpu
409
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:187
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:263
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:367
elementwise_filter
Implements elementwise_filter_t.
Definition elementwise_filter.f90:35
field_list
Definition field_list.f90:33
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:47
max
#define max(a, b)
Definition tensor.cu:40