Neko 1.99.9
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, &
66 delta, 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 = 0.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 call neko_scratch_registry%relinquish_field(temp_indices)
251
252 end subroutine compute_lij_cpu
253
268 subroutine compute_mij_cpu(mij, s11, s22, s33, s12, s13, s23, &
269 s_abs, test_filter, delta, n)
270 type(field_t), intent(inout) :: mij(6)
271 type(field_t), intent(inout) :: s11, s22, s33, s12, s13, s23, s_abs
272 type(elementwise_filter_t), intent(inout) :: test_filter
273 type(field_t), intent(in) :: delta
274 integer, intent(in) :: n
275
276 integer :: temp_indices(6)
277 type(field_t), pointer :: fs11, fs22, fs33, fs12, fs13, fs23
278 real(kind=rp) :: fs_abs, delta2
279 real(kind=rp) :: delta_ratio2 !test- to grid- filter ratio, squared
280 integer :: i
281
282 delta_ratio2 = ((test_filter%nx-1.0_rp)/(test_filter%nt-1.0_rp))**2
283
284 call neko_scratch_registry%request_field(fs11, temp_indices(1), .false.)
285 call neko_scratch_registry%request_field(fs22, temp_indices(2), .false.)
286 call neko_scratch_registry%request_field(fs33, temp_indices(3), .false.)
287 call neko_scratch_registry%request_field(fs12, temp_indices(4), .false.)
288 call neko_scratch_registry%request_field(fs13, temp_indices(5), .false.)
289 call neko_scratch_registry%request_field(fs23, temp_indices(6), .false.)
290
291 !! The first term:
292 !! _____ ____
293 !! (delta_test/delta)^2 s_abs*s_ij
294 call test_filter%apply(fs11, s11)
295 call test_filter%apply(fs22, s22)
296 call test_filter%apply(fs33, s33)
297 call test_filter%apply(fs12, s12)
298 call test_filter%apply(fs13, s13)
299 call test_filter%apply(fs23, s23)
300
301 ! Compute the test-filtered contribution and include delta before the
302 ! test filter to support both pointwise and elementwise-constant scales.
303 !OCL NORECURRENCE, NOVREC, NOALIAS
304 !DIR$ CONCURRENT
305 !DIR$ IVDEP
306 !GCC$ ivdep
307 !$omp parallel do private(fs_abs, delta2)
308 do i = 1, n
309 delta2 = delta%x(i,1,1,1)**2
310 fs_abs = sqrt(2.0_rp * &
311 (fs11%x(i,1,1,1)*fs11%x(i,1,1,1) + &
312 fs22%x(i,1,1,1)*fs22%x(i,1,1,1) + &
313 fs33%x(i,1,1,1)*fs33%x(i,1,1,1)) + &
314 4.0_rp * (fs12%x(i,1,1,1)*fs12%x(i,1,1,1) + &
315 fs13%x(i,1,1,1)*fs13%x(i,1,1,1) + &
316 fs23%x(i,1,1,1)*fs23%x(i,1,1,1)))
317 mij(1)%x(i,1,1,1) = delta_ratio2 * delta2 * fs_abs * fs11%x(i,1,1,1)
318 mij(2)%x(i,1,1,1) = delta_ratio2 * delta2 * fs_abs * fs22%x(i,1,1,1)
319 mij(3)%x(i,1,1,1) = delta_ratio2 * delta2 * fs_abs * fs33%x(i,1,1,1)
320 mij(4)%x(i,1,1,1) = delta_ratio2 * delta2 * fs_abs * fs12%x(i,1,1,1)
321 mij(5)%x(i,1,1,1) = delta_ratio2 * delta2 * fs_abs * fs13%x(i,1,1,1)
322 mij(6)%x(i,1,1,1) = delta_ratio2 * delta2 * fs_abs * fs23%x(i,1,1,1)
323 fs33%x(i,1,1,1) = s_abs%x(i,1,1,1) * delta2
324 end do
325 !$omp end parallel do
326
327 ! fs11, fs22, fs33 are reused as work arrays for the second term of M_ij
328 call col3(fs11%x, fs33%x, s11%x, n)
329 call test_filter%apply(fs22, fs11)
330 call sub2(mij(1)%x, fs22%x, n)
331
332 call col3(fs11%x, fs33%x, s22%x, n)
333 call test_filter%apply(fs22, fs11)
334 call sub2(mij(2)%x, fs22%x, n)
335
336 call col3(fs11%x, fs33%x, s33%x, n)
337 call test_filter%apply(fs22, fs11)
338 call sub2(mij(3)%x, fs22%x, n)
339
340 call col3(fs11%x, fs33%x, s12%x, n)
341 call test_filter%apply(fs22, fs11)
342 call sub2(mij(4)%x, fs22%x, n)
343
344 call col3(fs11%x, fs33%x, s13%x, n)
345 call test_filter%apply(fs22, fs11)
346 call sub2(mij(5)%x, fs22%x, n)
347
348 call col3(fs11%x, fs33%x, s23%x, n)
349 call test_filter%apply(fs22, fs11)
350 call sub2(mij(6)%x, fs22%x, n)
351
352 call neko_scratch_registry%relinquish_field(temp_indices)
353
354 end subroutine compute_mij_cpu
355
363 subroutine compute_num_den_cpu(num, den, lij, mij, alpha, n)
364 type(field_t), intent(inout) :: num, den
365 type(field_t), intent(in) :: lij(6), mij(6)
366 real(kind=rp), intent(in) :: alpha
367 integer, intent(in) :: n
368
369 real(kind=rp) :: num_curr, den_curr
370 integer :: i
371
372 !OCL NORECURRENCE, NOVREC, NOALIAS
373 !DIR$ CONCURRENT
374 !DIR$ IVDEP
375 !GCC$ ivdep
376 !$omp parallel do private(i, num_curr, den_curr)
377 do i = 1, n
378 num_curr = mij(1)%x(i,1,1,1)*lij(1)%x(i,1,1,1) + &
379 mij(2)%x(i,1,1,1)*lij(2)%x(i,1,1,1) + &
380 mij(3)%x(i,1,1,1)*lij(3)%x(i,1,1,1) + &
381 2.0_rp*(mij(4)%x(i,1,1,1)*lij(4)%x(i,1,1,1) + &
382 mij(5)%x(i,1,1,1)*lij(5)%x(i,1,1,1) + &
383 mij(6)%x(i,1,1,1)*lij(6)%x(i,1,1,1))
384 den_curr = mij(1)%x(i,1,1,1)*mij(1)%x(i,1,1,1) + &
385 mij(2)%x(i,1,1,1)*mij(2)%x(i,1,1,1) + &
386 mij(3)%x(i,1,1,1)*mij(3)%x(i,1,1,1) + &
387 2.0_rp*(mij(4)%x(i,1,1,1)*mij(4)%x(i,1,1,1) + &
388 mij(5)%x(i,1,1,1)*mij(5)%x(i,1,1,1) + &
389 mij(6)%x(i,1,1,1)*mij(6)%x(i,1,1,1))
390
391 ! running average over time
392 num%x(i,1,1,1) = alpha * num%x(i,1,1,1) + (1.0_rp - alpha) * num_curr
393 den%x(i,1,1,1) = alpha * den%x(i,1,1,1) + (1.0_rp - alpha) * den_curr
394 end do
395 !$omp end parallel do
396
397 end subroutine compute_num_den_cpu
398
399end module dynamic_smagorinsky_cpu
400
operators::strain_rate
Compute the strain rate tensor of a vector field.
Definition operators.f90:120
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. M_ij = delta_test^2 |S_test| S_ij_test.
Definition dynamic_smagorinsky_cpu.f90:270
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:364
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:518
math::cadd
subroutine, public cadd(a, s, n)
Add a scalar to vector .
Definition math.F90:580
math::col2
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.F90:1095
math::col3
subroutine, public col3(a, b, c, n)
Vector multiplication with 3 vectors .
Definition math.F90:1110
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:997
num_types
Definition num_types.f90:1
num_types::rp
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:14
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:158
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:201
coefs::coef_t
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:136
elementwise_filter::elementwise_filter_t
Implements the elementwise filter for SEM.
Definition elementwise_filter.f90:57
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