63 real(kind=
rp),
intent(in) :: t
64 integer,
intent(in) :: tstep
65 type(
coef_t),
intent(in) :: coef
66 type(
field_t),
intent(inout) :: nut
67 type(
field_t),
intent(in) :: delta
68 real(kind=
rp),
intent(in) :: c
70 type(
field_t),
pointer :: g11, g12, g13, g21, g22, g23, g31, g32, g33
71 type(
field_t),
pointer :: u, v, w
73 real(kind=
rp) :: sigg11, sigg12, sigg13, sigg22, sigg23, sigg33
74 real(kind=
rp) :: sigma1, sigma2, sigma3
75 real(kind=
rp) :: invariant1, invariant2, invariant3
76 real(kind=
rp) :: alpha1, alpha2, alpha3
77 real(kind=
rp) :: dsigma
78 real(kind=
rp) :: pi_3 = 4.0_rp/3.0_rp*atan(1.0_rp)
82 integer :: temp_indices(9)
106 call dudxyz (g11%x, u%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
107 call dudxyz (g12%x, u%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
108 call dudxyz (g13%x, u%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
110 call dudxyz (g21%x, v%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
111 call dudxyz (g22%x, v%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
112 call dudxyz (g23%x, v%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
114 call dudxyz (g31%x, w%x, coef%drdx, coef%dsdx, coef%dtdx, coef)
115 call dudxyz (g32%x, w%x, coef%drdy, coef%dsdy, coef%dtdy, coef)
116 call dudxyz (g33%x, w%x, coef%drdz, coef%dsdz, coef%dtdz, coef)
128 do concurrent(i = 1:g11%dof%size())
129 g11%x(i,1,1,1) = g11%x(i,1,1,1) * coef%mult(i,1,1,1)
130 g12%x(i,1,1,1) = g12%x(i,1,1,1) * coef%mult(i,1,1,1)
131 g13%x(i,1,1,1) = g13%x(i,1,1,1) * coef%mult(i,1,1,1)
132 g21%x(i,1,1,1) = g21%x(i,1,1,1) * coef%mult(i,1,1,1)
133 g22%x(i,1,1,1) = g22%x(i,1,1,1) * coef%mult(i,1,1,1)
134 g23%x(i,1,1,1) = g23%x(i,1,1,1) * coef%mult(i,1,1,1)
135 g31%x(i,1,1,1) = g31%x(i,1,1,1) * coef%mult(i,1,1,1)
136 g32%x(i,1,1,1) = g32%x(i,1,1,1) * coef%mult(i,1,1,1)
137 g33%x(i,1,1,1) = g33%x(i,1,1,1) * coef%mult(i,1,1,1)
140 do concurrent(e = 1:coef%msh%nelv)
141 do concurrent(i = 1:coef%Xh%lxyz)
143 sigg11 = g11%x(i,1,1,e)**2 + g21%x(i,1,1,e)**2 + g31%x(i,1,1,e)**2
144 sigg22 = g12%x(i,1,1,e)**2 + g22%x(i,1,1,e)**2 + g32%x(i,1,1,e)**2
145 sigg33 = g13%x(i,1,1,e)**2 + g23%x(i,1,1,e)**2 + g33%x(i,1,1,e)**2
146 sigg12 = g11%x(i,1,1,e)*g12%x(i,1,1,e) + &
147 g21%x(i,1,1,e)*g22%x(i,1,1,e) + &
148 g31%x(i,1,1,e)*g32%x(i,1,1,e)
149 sigg13 = g11%x(i,1,1,e)*g13%x(i,1,1,e) + &
150 g21%x(i,1,1,e)*g23%x(i,1,1,e) + &
151 g31%x(i,1,1,e)*g33%x(i,1,1,e)
152 sigg23 = g12%x(i,1,1,e)*g13%x(i,1,1,e) + &
153 g22%x(i,1,1,e)*g23%x(i,1,1,e) + &
154 g32%x(i,1,1,e)*g33%x(i,1,1,e)
162 if (abs(sigg11) .lt. eps)
then
165 if (abs(sigg12) .lt. eps)
then
168 if (abs(sigg13) .lt. eps)
then
171 if (abs(sigg22) .lt. eps)
then
174 if (abs(sigg23) .lt. eps)
then
177 if (abs(sigg33) .lt. eps)
then
181 if (abs(sigg12*sigg12 + &
182 sigg13*sigg13 + sigg23*sigg23) .lt. eps)
then
185 sigma1 = sqrt(
max(
max(
max(sigg11, sigg22), sigg33), 0.0_rp))
186 sigma3 = sqrt(
max(min(min(sigg11, sigg22), sigg33), 0.0_rp))
187 invariant1 = sigg11 + sigg22 + sigg33
188 sigma2 = sqrt(abs(invariant1 - sigma1*sigma1 - sigma3*sigma3))
192 invariant1 = sigg11 + sigg22 + sigg33
193 invariant2 = sigg11*sigg22 + sigg11*sigg33 + sigg22*sigg33 - &
194 (sigg12*sigg12 + sigg13*sigg13 + sigg23*sigg23)
195 invariant3 = sigg11*sigg22*sigg33 + &
196 2.0_rp*sigg12*sigg13*sigg23 - &
197 (sigg11*sigg23*sigg23 + sigg22*sigg13*sigg13 + &
198 sigg33*sigg12*sigg12)
203 invariant1 =
max(invariant1, 0.0_rp)
204 invariant2 =
max(invariant2, 0.0_rp)
205 invariant3 =
max(invariant3, 0.0_rp)
208 alpha1 = invariant1*invariant1/9.0_rp - invariant2/3.0_rp
212 alpha1 =
max(alpha1, 0.0_rp)
214 alpha2 = invariant1*invariant1*invariant1/27.0_rp - &
215 invariant1*invariant2/6.0_rp + invariant3/2.0_rp
220 tmp1 = alpha2/sqrt(alpha1 * alpha1 * alpha1)
222 if (tmp1 .le. -1.0_rp)
then
225 sigma1 = sqrt(
max(invariant1/3.0_rp + sqrt(alpha1), 0.0_rp))
227 sigma3 = sqrt(invariant1/3.0_rp - 2.0_rp*sqrt(alpha1))
229 elseif (tmp1 .ge. 1.0_rp)
then
231 sigma1 = sqrt(
max(invariant1/3.0_rp + 2.0_rp*sqrt(alpha1), &
233 sigma2 = sqrt(invariant1/3.0_rp - sqrt(alpha1))
236 alpha3 = acos(tmp1)/3.0_rp
238 if (abs(invariant3) .lt. eps)
then
241 sigma1 = sqrt(
max(invariant1/3.0_rp + &
242 2.0_rp*sqrt(alpha1)*cos(alpha3), 0.0_rp))
243 sigma2 = sqrt(abs(invariant1 - sigma1*sigma1))
246 sigma1 = sqrt(
max(invariant1/3.0_rp + &
247 2.0_rp*sqrt(alpha1)*cos(alpha3), 0.0_rp))
248 sigma2 = sqrt(invariant1/3.0_rp - &
249 2.0_rp*sqrt(alpha1)*cos(pi_3 + alpha3))
250 sigma3 = sqrt(abs(invariant1 - &
251 sigma1*sigma1-sigma2*sigma2))
257 if (sigma1 .gt. 0.0_rp)
then
259 sigma3*(sigma1 - sigma2)*(sigma2 - sigma3)/(sigma1*sigma1)
265 dsigma =
max(dsigma, 0.0_rp)
269 nut%x(i,1,1,e) = (c*delta%x(i,1,1,e))**2 * dsigma