Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
sigma_cpu.f90
Go to the documentation of this file.
1! Copyright (c) 2023-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!
37
39 use num_types, only : rp
40 use field_list, only : field_list_t
43 use field, only : field_t
44 use operators, only : dudxyz
45 use coefs, only : coef_t
46 use gs_ops, only : gs_op_add
47 use math, only : neko_eps, col2
48 implicit none
49 private
50
51 public :: sigma_compute_cpu
52
53contains
54
63 subroutine sigma_compute_cpu(if_ext, t, tstep, coef, nut, delta, c)
64 logical, intent(in) :: if_ext
65 real(kind=rp), intent(in) :: t
66 integer, intent(in) :: tstep
67 type(coef_t), intent(in) :: coef
68 type(field_t), intent(inout) :: nut
69 type(field_t), intent(in) :: delta
70 real(kind=rp), intent(in) :: c
71 ! This is the velocity gradient tensor
72 type(field_t), pointer :: g11, g12, g13, g21, g22, g23, g31, g32, g33
73 type(field_t), pointer :: u, v, w
74
75 real(kind=rp) :: sigg11, sigg12, sigg13, sigg22, sigg23, sigg33
76 real(kind=rp) :: sigma1, sigma2, sigma3
77 real(kind=rp) :: invariant1, invariant2, invariant3
78 real(kind=rp) :: alpha1, alpha2, alpha3
79 real(kind=rp) :: dsigma
80 real(kind=rp) :: pi_3 = 4.0_rp/3.0_rp*atan(1.0_rp)
81 real(kind=rp) :: tmp1
82 real(kind=rp) :: eps
83
84 integer :: temp_indices(9)
85 integer :: e, i
86
87 ! some constant
88 eps = neko_eps
89
90
91 ! get fields from registry
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 call neko_scratch_registry%request_field(g11, temp_indices(1))
103 call neko_scratch_registry%request_field(g12, temp_indices(2))
104 call neko_scratch_registry%request_field(g13, temp_indices(3))
105 call neko_scratch_registry%request_field(g21, temp_indices(4))
106 call neko_scratch_registry%request_field(g22, temp_indices(5))
107 call neko_scratch_registry%request_field(g23, temp_indices(6))
108 call neko_scratch_registry%request_field(g31, temp_indices(7))
109 call neko_scratch_registry%request_field(g32, temp_indices(8))
110 call neko_scratch_registry%request_field(g33, temp_indices(9))
111
112
113 ! Compute the derivatives of the velocity (the components of the g tensor)
114 call dudxyz (g11%x, u%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
115 call dudxyz (g12%x, u%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
116 call dudxyz (g13%x, u%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
117
118 call dudxyz (g21%x, v%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
119 call dudxyz (g22%x, v%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
120 call dudxyz (g23%x, v%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
121
122 call dudxyz (g31%x, w%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
123 call dudxyz (g32%x, w%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
124 call dudxyz (g33%x, w%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
125
126 call coef%gs_h%op(g11, gs_op_add)
127 call coef%gs_h%op(g12, gs_op_add)
128 call coef%gs_h%op(g13, gs_op_add)
129 call coef%gs_h%op(g21, gs_op_add)
130 call coef%gs_h%op(g22, gs_op_add)
131 call coef%gs_h%op(g23, gs_op_add)
132 call coef%gs_h%op(g31, gs_op_add)
133 call coef%gs_h%op(g32, gs_op_add)
134 call coef%gs_h%op(g33, gs_op_add)
135
136 do concurrent(e = 1:coef%msh%nelv)
137 do concurrent(i = 1:coef%Xh%lxyz)
138 ! G_ij = g^t g = g_mi g_mj
139 sigg11 = g11%x(i,1,1,e)**2 + g21%x(i,1,1,e)**2 + g31%x(i,1,1,e)**2
140 sigg22 = g12%x(i,1,1,e)**2 + g22%x(i,1,1,e)**2 + g32%x(i,1,1,e)**2
141 sigg33 = g13%x(i,1,1,e)**2 + g23%x(i,1,1,e)**2 + g33%x(i,1,1,e)**2
142 sigg12 = g11%x(i,1,1,e)*g12%x(i,1,1,e) + &
143 g21%x(i,1,1,e)*g22%x(i,1,1,e) + &
144 g31%x(i,1,1,e)*g32%x(i,1,1,e)
145 sigg13 = g11%x(i,1,1,e)*g13%x(i,1,1,e) + &
146 g21%x(i,1,1,e)*g23%x(i,1,1,e) + &
147 g31%x(i,1,1,e)*g33%x(i,1,1,e)
148 sigg23 = g12%x(i,1,1,e)*g13%x(i,1,1,e) + &
149 g22%x(i,1,1,e)*g23%x(i,1,1,e) + &
150 g32%x(i,1,1,e)*g33%x(i,1,1,e)
151
152 ! If LAPACK compute eigenvalues of the semi-definite positive matrix G
153 ! ..........to be done later on......
154 ! ELSE use the analytical method as done in the following
155
156 ! eigenvalues with the analytical method of Hasan et al. (2001)
157 ! doi:10.1006/jmre.2001.2400
158 if (abs(sigg11) .lt. eps) then
159 sigg11 = 0.0_rp
160 end if
161 if (abs(sigg12) .lt. eps) then
162 sigg12 = 0.0_rp
163 end if
164 if (abs(sigg13) .lt. eps) then
165 sigg13 = 0.0_rp
166 end if
167 if (abs(sigg22) .lt. eps) then
168 sigg22 = 0.0_rp
169 end if
170 if (abs(sigg23) .lt. eps) then
171 sigg23 = 0.0_rp
172 end if
173 if (abs(sigg33) .lt. eps) then
174 sigg33 = 0.0_rp
175 end if
176
177 if (abs(sigg12*sigg12 + &
178 sigg13*sigg13 + sigg23*sigg23) .lt. eps) then
179 ! G is diagonal
180 ! estimate the singular values according to:
181 sigma1 = sqrt(max(max(max(sigg11, sigg22), sigg33), 0.0_rp))
182 sigma3 = sqrt(max(min(min(sigg11, sigg22), sigg33), 0.0_rp))
183 invariant1 = sigg11 + sigg22 + sigg33
184 sigma2 = sqrt(abs(invariant1 - sigma1*sigma1 - sigma3*sigma3))
185 else
186
187 ! estimation of invariants
188 invariant1 = sigg11 + sigg22 + sigg33
189 invariant2 = sigg11*sigg22 + sigg11*sigg33 + sigg22*sigg33 - &
190 (sigg12*sigg12 + sigg13*sigg13 + sigg23*sigg23)
191 invariant3 = sigg11*sigg22*sigg33 + &
192 2.0_rp*sigg12*sigg13*sigg23 - &
193 (sigg11*sigg23*sigg23 + sigg22*sigg13*sigg13 + &
194 sigg33*sigg12*sigg12)
195
196 ! G is symmetric semi-definite positive matrix:
197 ! the invariants have to be larger-equal zero
198 ! which is obtained via forcing
199 invariant1 = max(invariant1, 0.0_rp)
200 invariant2 = max(invariant2, 0.0_rp)
201 invariant3 = max(invariant3, 0.0_rp)
202
203 ! compute the following angles from the invariants
204 alpha1 = invariant1*invariant1/9.0_rp - invariant2/3.0_rp
205
206 ! since alpha1 is always positive (see Hasan et al. (2001))
207 ! forcing is applied
208 alpha1 = max(alpha1, 0.0_rp)
209
210 alpha2 = invariant1*invariant1*invariant1/27.0_rp - &
211 invariant1*invariant2/6.0_rp + invariant3/2.0_rp
212
213 ! since acos(alpha2/(alpha1^(3/2)))/3.0_rp only valid for
214 ! alpha2^2 < alpha1^3.0_rp and arccos(x) only valid for -1<=x<=1
215 ! alpha3 is between 0 and pi/3
216 tmp1 = alpha2/sqrt(alpha1 * alpha1 * alpha1)
217
218 if (tmp1 .le. -1.0_rp) then
219 ! alpha3=pi/3 -> cos(alpha3)=0.5
220 ! compute the singular values
221 sigma1 = sqrt(max(invariant1/3.0_rp + sqrt(alpha1), 0.0_rp))
222 sigma2 = sigma1
223 sigma3 = sqrt(invariant1/3.0_rp - 2.0_rp*sqrt(alpha1))
224
225 elseif (tmp1 .ge. 1.0_rp) then
226 ! alpha3=0.0_rp -> cos(alpha3)=1.0
227 sigma1 = sqrt(max(invariant1/3.0_rp + 2.0_rp*sqrt(alpha1), &
228 0.0_rp))
229 sigma2 = sqrt(invariant1/3.0_rp - sqrt(alpha1))
230 sigma3 = sigma2
231 else
232 alpha3 = acos(tmp1)/3.0_rp
233
234 if (abs(invariant3) .lt. eps) then
235 ! In case of Invariant3=0, one or more eigenvalues are equal
236 ! to zero. Therefore force sigma3 to 0 and compute sigma1 and
237 ! sigma2
238 sigma1 = sqrt(max(invariant1/3.0_rp + &
239 2.0_rp*sqrt(alpha1)*cos(alpha3), 0.0_rp))
240 sigma2 = sqrt(abs(invariant1 - sigma1*sigma1))
241 sigma3 = 0.0_rp
242 else
243 sigma1 = sqrt(max(invariant1/3.0_rp + &
244 2.0_rp*sqrt(alpha1)*cos(alpha3), 0.0_rp))
245 sigma2 = sqrt(invariant1/3.0_rp - &
246 2.0_rp*sqrt(alpha1)*cos(pi_3 + alpha3))
247 sigma3 = sqrt(abs(invariant1 - &
248 sigma1*sigma1-sigma2*sigma2))
249 end if ! Invariant3=0 ?
250 end if ! tmp1
251 end if ! G diagonal ?
252
253 ! Estimate Dsigma
254 if (sigma1 .gt. 0.0_rp) then
255 dsigma = &
256 sigma3*(sigma1 - sigma2)*(sigma2 - sigma3)/(sigma1*sigma1)
257 else
258 dsigma = 0.0_rp
259 end if
260
261 !clipping to avoid negative values
262 dsigma = max(dsigma, 0.0_rp)
263
264 ! estimate turbulent viscosity
265
266 nut%x(i,1,1,e) = (c*delta%x(i,1,1,e))**2 * dsigma &
267 * coef%mult(i,1,1,1)
268
269 end do
270 end do
271
272 call coef%gs_h%op(nut, gs_op_add)
273 call col2(nut%x, coef%mult, nut%dof%size())
274
275 call neko_scratch_registry%relinquish_field(temp_indices)
276 end subroutine sigma_compute_cpu
277
278end module sigma_cpu
279
Coefficients.
Definition coef.f90:34
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Defines a field.
Definition field.f90:34
Defines Gather-scatter operations.
Definition gs_ops.f90:34
integer, parameter, public gs_op_add
Definition gs_ops.f90:36
Definition math.f90:60
subroutine, public col2(a, b, n)
Vector multiplication .
Definition math.f90:860
real(kind=rp), parameter, public neko_eps
Machine epsilon .
Definition math.f90:69
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Operators.
Definition operators.f90:34
subroutine, public dudxyz(du, u, dr, ds, dt, coef)
Compute derivative of a scalar field along a single direction.
Definition operators.f90:80
Defines a registry for storing and requesting temporary fields This can be used when you have a funct...
type(scratch_registry_t), target, public neko_scratch_registry
Global scratch registry.
Implements the CPU kernel for the sigma_t type. Following Nicoud et al. "Using singular values to bui...
Definition sigma_cpu.f90:38
subroutine, public sigma_compute_cpu(if_ext, t, tstep, coef, nut, delta, c)
Compute eddy viscosity on the CPU.
Definition sigma_cpu.f90:64
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
field_list_t, To be able to group fields together
#define max(a, b)
Definition tensor.cu:40