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