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
72 type(
field_t),
pointer :: g11, g12, g13, g21, g22, g23, g31, g32, g33
73 type(
field_t),
pointer :: u, v, w
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)
84 integer :: temp_indices(9)
92 if (if_ext .eqv. .true.)
then
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)
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)
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)
136 do concurrent(e = 1:coef%msh%nelv)
137 do concurrent(i = 1:coef%Xh%lxyz)
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)
158 if (abs(sigg11) .lt. eps)
then
161 if (abs(sigg12) .lt. eps)
then
164 if (abs(sigg13) .lt. eps)
then
167 if (abs(sigg22) .lt. eps)
then
170 if (abs(sigg23) .lt. eps)
then
173 if (abs(sigg33) .lt. eps)
then
177 if (abs(sigg12*sigg12 + &
178 sigg13*sigg13 + sigg23*sigg23) .lt. eps)
then
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))
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)
199 invariant1 =
max(invariant1, 0.0_rp)
200 invariant2 =
max(invariant2, 0.0_rp)
201 invariant3 =
max(invariant3, 0.0_rp)
204 alpha1 = invariant1*invariant1/9.0_rp - invariant2/3.0_rp
208 alpha1 =
max(alpha1, 0.0_rp)
210 alpha2 = invariant1*invariant1*invariant1/27.0_rp - &
211 invariant1*invariant2/6.0_rp + invariant3/2.0_rp
216 tmp1 = alpha2/sqrt(alpha1 * alpha1 * alpha1)
218 if (tmp1 .le. -1.0_rp)
then
221 sigma1 = sqrt(
max(invariant1/3.0_rp + sqrt(alpha1), 0.0_rp))
223 sigma3 = sqrt(invariant1/3.0_rp - 2.0_rp*sqrt(alpha1))
225 elseif (tmp1 .ge. 1.0_rp)
then
227 sigma1 = sqrt(
max(invariant1/3.0_rp + 2.0_rp*sqrt(alpha1), &
229 sigma2 = sqrt(invariant1/3.0_rp - sqrt(alpha1))
232 alpha3 = acos(tmp1)/3.0_rp
234 if (abs(invariant3) .lt. eps)
then
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))
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))
254 if (sigma1 .gt. 0.0_rp)
then
256 sigma3*(sigma1 - sigma2)*(sigma2 - sigma3)/(sigma1*sigma1)
262 dsigma =
max(dsigma, 0.0_rp)
266 nut%x(i,1,1,e) = (c*delta%x(i,1,1,e))**2 * dsigma &
273 call col2(nut%x, coef%mult, nut%dof%size())