62 real(kind=
rp),
intent(in) :: t
63 integer,
intent(in) :: tstep
64 type(
coef_t),
intent(in) :: coef
65 type(
field_t),
intent(inout) :: nut
66 type(
field_t),
intent(in) :: delta
67 real(kind=
rp),
intent(in) :: c
69 type(
field_t),
pointer :: g11, g12, g13, g21, g22, g23, g31, g32, g33
70 type(
field_t),
pointer :: u, v, w
72 real(kind=
rp) :: sigg11, sigg12, sigg13, sigg22, sigg23, sigg33
73 real(kind=
rp) :: sigma1, sigma2, sigma3
74 real(kind=
rp) :: invariant1, invariant2, invariant3
75 real(kind=
rp) :: alpha1, alpha2, alpha3
76 real(kind=
rp) :: dsigma
77 real(kind=
rp) :: pi_3 = 4.0_rp/3.0_rp*atan(1.0_rp)
81 integer :: temp_indices(9)
105 call dudxyz (g11%x, u%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
106 call dudxyz (g12%x, u%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
107 call dudxyz (g13%x, u%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
109 call dudxyz (g21%x, v%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
110 call dudxyz (g22%x, v%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
111 call dudxyz (g23%x, v%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
113 call dudxyz (g31%x, w%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
114 call dudxyz (g32%x, w%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
115 call dudxyz (g33%x, w%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
127 do concurrent(i = 1:g11%dof%size())
128 g11%x(i,1,1,1) = g11%x(i,1,1,1) * coef%mult(i,1,1,1)
129 g12%x(i,1,1,1) = g12%x(i,1,1,1) * coef%mult(i,1,1,1)
130 g13%x(i,1,1,1) = g13%x(i,1,1,1) * coef%mult(i,1,1,1)
131 g21%x(i,1,1,1) = g21%x(i,1,1,1) * coef%mult(i,1,1,1)
132 g22%x(i,1,1,1) = g22%x(i,1,1,1) * coef%mult(i,1,1,1)
133 g23%x(i,1,1,1) = g23%x(i,1,1,1) * coef%mult(i,1,1,1)
134 g31%x(i,1,1,1) = g31%x(i,1,1,1) * coef%mult(i,1,1,1)
135 g32%x(i,1,1,1) = g32%x(i,1,1,1) * coef%mult(i,1,1,1)
136 g33%x(i,1,1,1) = g33%x(i,1,1,1) * coef%mult(i,1,1,1)
139 do concurrent(e = 1:coef%msh%nelv)
140 do concurrent(i = 1:coef%Xh%lxyz)
142 sigg11 = g11%x(i,1,1,e)**2 + g21%x(i,1,1,e)**2 + g31%x(i,1,1,e)**2
143 sigg22 = g12%x(i,1,1,e)**2 + g22%x(i,1,1,e)**2 + g32%x(i,1,1,e)**2
144 sigg33 = g13%x(i,1,1,e)**2 + g23%x(i,1,1,e)**2 + g33%x(i,1,1,e)**2
145 sigg12 = g11%x(i,1,1,e)*g12%x(i,1,1,e) + &
146 g21%x(i,1,1,e)*g22%x(i,1,1,e) + &
147 g31%x(i,1,1,e)*g32%x(i,1,1,e)
148 sigg13 = g11%x(i,1,1,e)*g13%x(i,1,1,e) + &
149 g21%x(i,1,1,e)*g23%x(i,1,1,e) + &
150 g31%x(i,1,1,e)*g33%x(i,1,1,e)
151 sigg23 = g12%x(i,1,1,e)*g13%x(i,1,1,e) + &
152 g22%x(i,1,1,e)*g23%x(i,1,1,e) + &
153 g32%x(i,1,1,e)*g33%x(i,1,1,e)
161 if (abs(sigg11) .lt. eps)
then
164 if (abs(sigg12) .lt. eps)
then
167 if (abs(sigg13) .lt. eps)
then
170 if (abs(sigg22) .lt. eps)
then
173 if (abs(sigg23) .lt. eps)
then
176 if (abs(sigg33) .lt. eps)
then
180 if (abs(sigg12*sigg12 + &
181 sigg13*sigg13 + sigg23*sigg23) .lt. eps)
then
184 sigma1 = sqrt(
max(
max(
max(sigg11, sigg22), sigg33), 0.0_rp))
185 sigma3 = sqrt(
max(min(min(sigg11, sigg22), sigg33), 0.0_rp))
186 invariant1 = sigg11 + sigg22 + sigg33
187 sigma2 = sqrt(abs(invariant1 - sigma1*sigma1 - sigma3*sigma3))
191 invariant1 = sigg11 + sigg22 + sigg33
192 invariant2 = sigg11*sigg22 + sigg11*sigg33 + sigg22*sigg33 - &
193 (sigg12*sigg12 + sigg13*sigg13 + sigg23*sigg23)
194 invariant3 = sigg11*sigg22*sigg33 + &
195 2.0_rp*sigg12*sigg13*sigg23 - &
196 (sigg11*sigg23*sigg23 + sigg22*sigg13*sigg13 + &
197 sigg33*sigg12*sigg12)
202 invariant1 =
max(invariant1, 0.0_rp)
203 invariant2 =
max(invariant2, 0.0_rp)
204 invariant3 =
max(invariant3, 0.0_rp)
207 alpha1 = invariant1*invariant1/9.0_rp - invariant2/3.0_rp
211 alpha1 =
max(alpha1, 0.0_rp)
213 alpha2 = invariant1*invariant1*invariant1/27.0_rp - &
214 invariant1*invariant2/6.0_rp + invariant3/2.0_rp
219 tmp1 = alpha2/sqrt(alpha1 * alpha1 * alpha1)
221 if (tmp1 .le. -1.0_rp)
then
224 sigma1 = sqrt(
max(invariant1/3.0_rp + sqrt(alpha1), 0.0_rp))
226 sigma3 = sqrt(invariant1/3.0_rp - 2.0_rp*sqrt(alpha1))
228 elseif (tmp1 .ge. 1.0_rp)
then
230 sigma1 = sqrt(
max(invariant1/3.0_rp + 2.0_rp*sqrt(alpha1), &
232 sigma2 = sqrt(invariant1/3.0_rp - sqrt(alpha1))
235 alpha3 = acos(tmp1)/3.0_rp
237 if (abs(invariant3) .lt. eps)
then
240 sigma1 = sqrt(
max(invariant1/3.0_rp + &
241 2.0_rp*sqrt(alpha1)*cos(alpha3), 0.0_rp))
242 sigma2 = sqrt(abs(invariant1 - sigma1*sigma1))
245 sigma1 = sqrt(
max(invariant1/3.0_rp + &
246 2.0_rp*sqrt(alpha1)*cos(alpha3), 0.0_rp))
247 sigma2 = sqrt(invariant1/3.0_rp - &
248 2.0_rp*sqrt(alpha1)*cos(pi_3 + alpha3))
249 sigma3 = sqrt(abs(invariant1 - &
250 sigma1*sigma1-sigma2*sigma2))
256 if (sigma1 .gt. 0.0_rp)
then
258 sigma3*(sigma1 - sigma2)*(sigma2 - sigma3)/(sigma1*sigma1)
264 dsigma =
max(dsigma, 0.0_rp)
268 nut%x(i,1,1,e) = (c*delta%x(i,1,1,e))**2 * dsigma
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Defines Gather-scatter operations.
integer, parameter, public gs_op_add
integer, parameter, public rp
Global precision used in computations.
subroutine, public dudxyz(du, u, dr, ds, dt, coef)
Compute derivative of a scalar field along a single direction.
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...
subroutine, public sigma_compute_cpu(t, tstep, coef, nut, delta, c)
Compute eddy viscosity on the CPU.
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
field_list_t, To be able to group fields together