ADflow  v1.0
ADflow is a finite volume RANS solver tailored for gradient-based aerodynamic design optimization.
sa.F90
Go to the documentation of this file.
1 ! This module contains the source code related to the SA turbulence
2 ! model. It is slightly more modularized than the original which makes
3 ! performing reverse mode AD simplier.
4 
5 module sa
6 
7  use constants
8  real(kind=realtype) :: cv13, kar2inv, cw36, cb3inv
9  real(kind=realtype), dimension(:, :, :), allocatable :: qq
10  real(kind=realtype), dimension(:, :, :), pointer :: ddw, ww, ddvt
11  real(kind=realtype), dimension(:, :), pointer :: rrlv
12  real(kind=realtype), dimension(:, :), pointer :: dd2wall
13 
14 contains
15 #ifndef USE_TAPENADE
16  subroutine sa_block(resOnly)
17  !
18  ! sa solves the transport equation for the Spalart-Allmaras
19  ! turbulence model in a decoupled manner using a diagonal
20  ! dominant ADI-scheme. Note that the scratch and boundary
21  ! matrix values are not strictly, but tapande would like to
22  ! see them becuase it must save them.
23  !
24  use constants
25  use blockpointers, only: ndom, il, jl, kl, scratch, bmtj1, bmtj2, &
26  bmti1, bmti2, bmtk1, bmtk2
27  use inputtimespectral, only: ntimeintervalsspectral
28  use iteration, only: currentlevel
29  use inputphysics, only: turbprod
30  use paramturb
31  use turbutils
32  use turbbcroutines
33  implicit none
34  !
35  ! Subroutine argument.
36  !
37  logical, intent(in) :: resOnly
38  !
39  ! Local variables.
40  !
41  integer(kind=intType) :: nn, sps
42 
43  ! Set the arrays for the boundary condition treatment.
44  call bcturbtreatment
45 
46  ! Alloc central jacobian memory
47  allocate (qq(2:il, 2:jl, 2:kl))
48 
49  ! Source Terms
50  call sasource
51 
52  ! Advection Term
53  nn = itu1 - 1
54  call turbadvection(1_inttype, 1_inttype, nn, qq)
55 
56  ! Unsteady Term
57  call unsteadyturbterm(1_inttype, 1_inttype, nn, qq)
58 
59  ! Viscous Terms
60  call saviscous
61 
62  ! Perform the residual scaling
63  call saresscale
64 
65  ! We need to do an acutal solve. Solve and update the eddy
66  ! viscosity and the boundary conditions
67 
68  if (.not. resonly) then
69 
70  ! Do solve
71  call sasolve
72 
73  ! Compute the corresponding eddy viscosity.
74 
75  call saeddyviscosity(2, il, 2, jl, 2, kl)
76 
77  ! Set the halo values for the turbulent variables.
78  ! We are on the finest mesh, so the second layer of halo
79  ! cells must be computed as well.
80 
81  call applyallturbbcthisblock(.true.)
82  end if
83 
84  deallocate (qq)
85  end subroutine sa_block
86 #endif
87 
88  subroutine sasource
89  !
90  ! Source terms.
91  ! Determine the source term and its derivative w.r.t. nuTilde
92  ! for all internal cells of the block.
93  ! Remember that the SA field variable nuTilde = w(i,j,k,itu1)
94 
95  use blockpointers
96  use constants
97  use paramturb
98  use section
99  use inputphysics
100  use inputdiscretization, only: approxsa
101  use flowvarrefstate
102  implicit none
103 
104  ! Local parameters
105  real(kind=realtype), parameter :: f23 = two * third
106 
107  ! Local variables.
108  integer(kind=intType) :: i, j, k, nn, ii
109  real(kind=realtype) :: fv1, fv2, ft2
110  real(kind=realtype) :: ss, sst, nu, dist2inv, chi, chi2, chi3
111  real(kind=realtype) :: rr, gg, gg6, termfw, fwsa, term1, term2
112  real(kind=realtype) :: dfv1, dfv2, dft2, drr, dgg, dfw
113  real(kind=realtype) :: uux, uuy, uuz, vvx, vvy, vvz, wwx, wwy, wwz
114  real(kind=realtype) :: div2, fact, sxx, syy, szz, sxy, sxz, syz
115  real(kind=realtype) :: vortx, vorty, vortz
116  real(kind=realtype) :: omegax, omegay, omegaz
117  real(kind=realtype) :: strainmag2, strainprod, vortprod
118  real(kind=realtype), parameter :: xminn = 1.e-10_realtype
119 
120  ! Set model constants
121  cv13 = rsacv1**3
122  kar2inv = one / (rsak**2)
123  cw36 = rsacw3**6
124  cb3inv = one / rsacb3
125 
126  ! Determine the non-dimensional wheel speed of this block.
127 
128  omegax = timeref * sections(sectionid)%rotRate(1)
129  omegay = timeref * sections(sectionid)%rotRate(2)
130  omegaz = timeref * sections(sectionid)%rotRate(3)
131 
132  ! Create switches to production term depending on the variable that
133  ! should be used
134  if (turbprod .eq. katolaunder) then
135  print *, 'katoLaunder production term not supported for SA'
136  stop
137  end if
138 
139 #ifdef TAPENADE_REVERSE
140  !$AD II-LOOP
141  do ii = 0, nx * ny * nz - 1
142  i = mod(ii, nx) + 2
143  j = mod(ii / nx, ny) + 2
144  k = ii / (nx * ny) + 2
145 #else
146  do k = 2, kl
147  do j = 2, jl
148  do i = 2, il
149 #endif
150  ! Compute the gradient of u in the cell center. Use is made
151  ! of the fact that the surrounding normals sum up to zero,
152  ! such that the cell i,j,k does not give a contribution.
153  ! The gradient is scaled by the factor 2*vol.
154 
155  uux = w(i + 1, j, k, ivx) * si(i, j, k, 1) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 1) &
156  + w(i, j + 1, k, ivx) * sj(i, j, k, 1) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 1) &
157  + w(i, j, k + 1, ivx) * sk(i, j, k, 1) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 1)
158  uuy = w(i + 1, j, k, ivx) * si(i, j, k, 2) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 2) &
159  + w(i, j + 1, k, ivx) * sj(i, j, k, 2) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 2) &
160  + w(i, j, k + 1, ivx) * sk(i, j, k, 2) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 2)
161  uuz = w(i + 1, j, k, ivx) * si(i, j, k, 3) - w(i - 1, j, k, ivx) * si(i - 1, j, k, 3) &
162  + w(i, j + 1, k, ivx) * sj(i, j, k, 3) - w(i, j - 1, k, ivx) * sj(i, j - 1, k, 3) &
163  + w(i, j, k + 1, ivx) * sk(i, j, k, 3) - w(i, j, k - 1, ivx) * sk(i, j, k - 1, 3)
164 
165  ! Idem for the gradient of v.
166 
167  vvx = w(i + 1, j, k, ivy) * si(i, j, k, 1) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 1) &
168  + w(i, j + 1, k, ivy) * sj(i, j, k, 1) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 1) &
169  + w(i, j, k + 1, ivy) * sk(i, j, k, 1) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 1)
170  vvy = w(i + 1, j, k, ivy) * si(i, j, k, 2) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 2) &
171  + w(i, j + 1, k, ivy) * sj(i, j, k, 2) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 2) &
172  + w(i, j, k + 1, ivy) * sk(i, j, k, 2) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 2)
173  vvz = w(i + 1, j, k, ivy) * si(i, j, k, 3) - w(i - 1, j, k, ivy) * si(i - 1, j, k, 3) &
174  + w(i, j + 1, k, ivy) * sj(i, j, k, 3) - w(i, j - 1, k, ivy) * sj(i, j - 1, k, 3) &
175  + w(i, j, k + 1, ivy) * sk(i, j, k, 3) - w(i, j, k - 1, ivy) * sk(i, j, k - 1, 3)
176 
177  ! And for the gradient of w.
178 
179  wwx = w(i + 1, j, k, ivz) * si(i, j, k, 1) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 1) &
180  + w(i, j + 1, k, ivz) * sj(i, j, k, 1) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 1) &
181  + w(i, j, k + 1, ivz) * sk(i, j, k, 1) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 1)
182  wwy = w(i + 1, j, k, ivz) * si(i, j, k, 2) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 2) &
183  + w(i, j + 1, k, ivz) * sj(i, j, k, 2) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 2) &
184  + w(i, j, k + 1, ivz) * sk(i, j, k, 2) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 2)
185  wwz = w(i + 1, j, k, ivz) * si(i, j, k, 3) - w(i - 1, j, k, ivz) * si(i - 1, j, k, 3) &
186  + w(i, j + 1, k, ivz) * sj(i, j, k, 3) - w(i, j - 1, k, ivz) * sj(i, j - 1, k, 3) &
187  + w(i, j, k + 1, ivz) * sk(i, j, k, 3) - w(i, j, k - 1, ivz) * sk(i, j, k - 1, 3)
188 
189  ! Compute the components of the stress tensor.
190  ! The combination of the current scaling of the velocity
191  ! gradients (2*vol) and the definition of the stress tensor,
192  ! leads to the factor 1/(4*vol).
193 
194  fact = fourth / vol(i, j, k)
195 
196  if (turbprod .eq. strain) then
197 
198  sxx = two * fact * uux
199  syy = two * fact * vvy
200  szz = two * fact * wwz
201 
202  sxy = fact * (uuy + vvx)
203  sxz = fact * (uuz + wwx)
204  syz = fact * (vvz + wwy)
205 
206  ! Compute 2/3 * divergence of velocity squared
207 
208  div2 = f23 * (sxx + syy + szz)**2
209 
210  ! Compute strain production term
211 
212  strainmag2 = two * (sxy**2 + sxz**2 + syz**2) &
213  + sxx**2 + syy**2 + szz**2
214 
215  strainprod = two * strainmag2 - div2
216 
217  ss = sqrt(strainprod)
218 
219  else if (turbprod .eq. vorticity) then
220 
221  ! Compute the three components of the vorticity vector.
222  ! Substract the part coming from the rotating frame.
223 
224  vortx = two * fact * (wwy - vvz) - two * omegax
225  vorty = two * fact * (uuz - wwx) - two * omegay
226  vortz = two * fact * (vvx - uuy) - two * omegaz
227 
228  ! Compute the vorticity production term
229 
230  vortprod = vortx**2 + vorty**2 + vortz**2
231 
232  ! First take the square root of the production term to
233  ! obtain the correct production term for spalart-allmaras.
234  ! We do this to avoid if statements.
235 
236  ss = sqrt(vortprod)
237 
238  end if
239 
240  ! Compute the laminar kinematic viscosity, the inverse of
241  ! wall distance squared, the ratio chi (ratio of nuTilde
242  ! and nu) and the functions fv1 and fv2. The latter corrects
243  ! the production term near a viscous wall.
244 
245  nu = rlv(i, j, k) / w(i, j, k, irho)
246  dist2inv = one / (d2wall(i, j, k)**2)
247  chi = w(i, j, k, itu1) / nu
248  chi2 = chi * chi
249  chi3 = chi * chi2
250  fv1 = chi3 / (chi3 + cv13)
251  fv2 = one - chi / (one + chi * fv1)
252 
253  ! The function ft2, which is designed to keep a laminar
254  ! solution laminar. When running in fully turbulent mode
255  ! this function should be set to 0.0.
256 
257  if (useft2sa) then
258  ft2 = rsact3 * exp(-rsact4 * chi2)
259  else
260  ft2 = zero
261  end if
262 
263  ! Correct the production term to account for the influence
264  ! of the wall.
265 
266  sst = ss + w(i, j, k, itu1) * fv2 * kar2inv * dist2inv
267 
268  ! Add rotation term (useRotationSA defined in inputParams.F90)
269 
270  if (userotationsa) then
271  sst = sst + rsacrot * min(zero, sqrt(two * strainmag2))
272  end if
273 
274  ! Make sure that this term remains positive
275  ! (the function fv2 is negative between chi = 1 and 18.4,
276  ! which can cause sst to go negative, which is undesirable).
277 
278  sst = max(sst, xminn)
279 
280  ! Compute the function fw. The argument rr is cut off at 10
281  ! to avoid numerical problems. This is ok, because the
282  ! asymptotical value of fw is then already reached.
283 
284  rr = w(i, j, k, itu1) * kar2inv * dist2inv / sst
285  rr = min(rr, 10.0_realtype)
286  gg = rr + rsacw2 * (rr**6 - rr)
287  gg6 = gg**6
288  termfw = ((one + cw36) / (gg6 + cw36))**sixth
289  fwsa = gg * termfw
290 
291  ! Compute the source term; some terms are saved for the
292  ! linearization. The source term is stored in dvt.
293 
294  if (approxsa) then
295  term1 = zero
296  else
297  term1 = rsacb1 * (one - ft2) * ss
298  end if
299  term2 = dist2inv * (kar2inv * rsacb1 * ((one - ft2) * fv2 + ft2) &
300  - rsacw1 * fwsa)
301 
302  scratch(i, j, k, idvt) = (term1 + term2 * w(i, j, k, itu1)) * w(i, j, k, itu1)
303 
304 #ifndef USE_TAPENADE
305  ! Compute some derivatives w.r.t. nuTilde. These will occur
306  ! in the left hand side, i.e. the matrix for the implicit
307  ! treatment.
308 
309  dfv1 = three * chi2 * cv13 / ((chi3 + cv13)**2)
310  dfv2 = (chi2 * dfv1 - one) / (nu * ((one + chi * fv1)**2))
311  dft2 = -two * rsact4 * chi * ft2 / nu
312 
313  drr = (one - rr * (fv2 + w(i, j, k, itu1) * dfv2)) &
314  * kar2inv * dist2inv / sst
315  dgg = (one - rsacw2 + six * rsacw2 * (rr**5)) * drr
316  dfw = (cw36 / (gg6 + cw36)) * termfw * dgg
317 
318  ! Compute the source term jacobian. Note that the part
319  ! containing term1 is treated explicitly. The reason is that
320  ! implicit treatment of this part leads to a decrease of the
321  ! diagonal dominance of the jacobian and it thus decreases
322  ! the stability. You may want to play around and try to
323  ! take this term into account in the jacobian.
324  ! Note that -dsource/dnu is stored.
325  qq(i, j, k) = -two * term2 * w(i, j, k, itu1) &
326  - dist2inv * w(i, j, k, itu1) * w(i, j, k, itu1) &
327  * (rsacb1 * kar2inv * (dfv2 - ft2 * dfv2 - fv2 * dft2 + dft2) &
328  - rsacw1 * dfw)
329 
330  ! A couple of terms in qq may lead to a negative
331  ! contribution. Clip qq to zero, if the total is negative.
332 
333  qq(i, j, k) = max(qq(i, j, k), zero)
334 #endif
335 #ifdef TAPENADE_REVERSE
336  end do
337 #else
338  end do
339  end do
340  end do
341 #endif
342  end subroutine sasource
343 
344  subroutine saviscous
345  !
346  ! Viscous term.
347  ! Determine the viscous contribution to the residual
348  ! for all internal cells of the block.
349 
350  use blockpointers
351  use paramturb
352  implicit none
353  ! Local variables.
354  integer(kind=intType) :: i, j, k, nn, ii
355  real(kind=realtype) :: nu
356  real(kind=realtype) :: fv1, fv2, ft2
357  real(kind=realtype) :: voli, volmi, volpi, xm, ym, zm, xp, yp, zp
358  real(kind=realtype) :: xa, ya, za, ttm, ttp, cnud, cam, cap
359  real(kind=realtype) :: nutm, nutp, num, nup, cdm, cdp
360  real(kind=realtype) :: c1m, c1p, c10, b1, c1, d1, qs
361 
362  ! Set model constants
363  cv13 = rsacv1**3
364  kar2inv = one / (rsak**2)
365  cw36 = rsacw3**6
366  cb3inv = one / rsacb3
367 
368  !
369  ! Viscous terms in k-direction.
370  !
371 #ifdef TAPENADE_REVERSE
372  !$AD II-LOOP
373  do ii = 0, nx * ny * nz - 1
374  i = mod(ii, nx) + 2
375  j = mod(ii / nx, ny) + 2
376  k = ii / (nx * ny) + 2
377 #else
378  do k = 2, kl
379  do j = 2, jl
380  do i = 2, il
381 #endif
382  ! Compute the metrics in zeta-direction, i.e. along the
383  ! line k = constant.
384 
385  voli = one / vol(i, j, k)
386  volmi = two / (vol(i, j, k) + vol(i, j, k - 1))
387  volpi = two / (vol(i, j, k) + vol(i, j, k + 1))
388 
389  xm = sk(i, j, k - 1, 1) * volmi
390  ym = sk(i, j, k - 1, 2) * volmi
391  zm = sk(i, j, k - 1, 3) * volmi
392  xp = sk(i, j, k, 1) * volpi
393  yp = sk(i, j, k, 2) * volpi
394  zp = sk(i, j, k, 3) * volpi
395 
396  xa = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1)) * voli
397  ya = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2)) * voli
398  za = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3)) * voli
399  ttm = xm * xa + ym * ya + zm * za
400  ttp = xp * xa + yp * ya + zp * za
401 
402  ! ttm and ttp ~ 1/deltaX^2
403 
404  ! Computation of the viscous terms in zeta-direction; note
405  ! that cross-derivatives are neglected, i.e. the mesh is
406  ! assumed to be orthogonal.
407  ! Furthermore, the grad(nu)**2 has been rewritten as
408  ! div(nu grad(nu)) - nu div(grad nu) to enhance stability.
409  ! The second derivative in zeta-direction is constructed as
410  ! the central difference of the first order derivatives, i.e.
411  ! d^2/dzeta^2 = d/dzeta (d/dzeta k+1/2 - d/dzeta k-1/2).
412  ! In this way the metric can be taken into account.
413 
414  ! Compute the diffusion coefficients multiplying the nodes
415  ! k+1, k and k-1 in the second derivative. Make sure that
416  ! these coefficients are nonnegative.
417 
418  cnud = -rsacb2 * w(i, j, k, itu1) * cb3inv
419  cam = ttm * cnud
420  cap = ttp * cnud
421 
422  ! Compute nuTilde at the faces
423 
424  nutm = half * (w(i, j, k - 1, itu1) + w(i, j, k, itu1))
425  nutp = half * (w(i, j, k + 1, itu1) + w(i, j, k, itu1))
426 
427  ! Compute nu at the faces
428 
429  nu = rlv(i, j, k) / w(i, j, k, irho)
430  num = half * (rlv(i, j, k - 1) / w(i, j, k - 1, irho) + nu)
431  nup = half * (rlv(i, j, k + 1) / w(i, j, k + 1, irho) + nu)
432 
433  cdm = (num + (one + rsacb2) * nutm) * ttm * cb3inv
434  cdp = (nup + (one + rsacb2) * nutp) * ttp * cb3inv
435 
436  c1m = max(cdm + cam, zero)
437  c1p = max(cdp + cap, zero)
438  c10 = c1m + c1p
439 
440  ! Update the residual for this cell and store the possible
441  ! coefficients for the matrix in b1, c1 and d1.
442 
443  scratch(i, j, k, idvt) = scratch(i, j, k, idvt) + c1m * w(i, j, k - 1, itu1) &
444  - c10 * w(i, j, k, itu1) + c1p * w(i, j, k + 1, itu1)
445 #ifndef USE_TAPENADE
446  b1 = -c1m
447  c1 = c10
448  d1 = -c1p
449 
450  ! Update the central jacobian. For nonboundary cells this
451  ! is simply c1. For boundary cells this is slightly more
452  ! complicated, because the boundary conditions are treated
453  ! implicitly and the off-diagonal terms b1 and d1 must be
454  ! taken into account.
455  ! The boundary conditions are only treated implicitly if
456  ! the diagonal dominance of the matrix is increased.
457 
458  if (k == 2) then
459  qq(i, j, k) = qq(i, j, k) + c1 &
460  - b1 * max(bmtk1(i, j, itu1, itu1), zero)
461  else if (k == kl) then
462  qq(i, j, k) = qq(i, j, k) + c1 &
463  - d1 * max(bmtk2(i, j, itu1, itu1), zero)
464  else
465  qq(i, j, k) = qq(i, j, k) + c1
466  end if
467 #endif
468 #ifdef TAPENADE_REVERSE
469  end do
470 #else
471  end do
472  end do
473  end do
474 #endif
475  !
476  ! Viscous terms in j-direction.
477  !
478 #ifdef TAPENADE_REVERSE
479  !$AD II-LOOP
480  do ii = 0, nx * ny * nz - 1
481  i = mod(ii, nx) + 2
482  j = mod(ii / nx, ny) + 2
483  k = ii / (nx * ny) + 2
484 #else
485  do k = 2, kl
486  do j = 2, jl
487  do i = 2, il
488 #endif
489  ! Compute the metrics in eta-direction, i.e. along the
490  ! line j = constant.
491 
492  voli = one / vol(i, j, k)
493  volmi = two / (vol(i, j, k) + vol(i, j - 1, k))
494  volpi = two / (vol(i, j, k) + vol(i, j + 1, k))
495 
496  xm = sj(i, j - 1, k, 1) * volmi
497  ym = sj(i, j - 1, k, 2) * volmi
498  zm = sj(i, j - 1, k, 3) * volmi
499  xp = sj(i, j, k, 1) * volpi
500  yp = sj(i, j, k, 2) * volpi
501  zp = sj(i, j, k, 3) * volpi
502 
503  xa = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1)) * voli
504  ya = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2)) * voli
505  za = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3)) * voli
506  ttm = xm * xa + ym * ya + zm * za
507  ttp = xp * xa + yp * ya + zp * za
508 
509  ! Computation of the viscous terms in eta-direction; note
510  ! that cross-derivatives are neglected, i.e. the mesh is
511  ! assumed to be orthogonal.
512  ! Furthermore, the grad(nu)**2 has been rewritten as
513  ! div(nu grad(nu)) - nu div(grad nu) to enhance stability.
514  ! The second derivative in eta-direction is constructed as
515  ! the central difference of the first order derivatives, i.e.
516  ! d^2/deta^2 = d/deta (d/deta j+1/2 - d/deta j-1/2).
517  ! In this way the metric can be taken into account.
518 
519  ! Compute the diffusion coefficients multiplying the nodes
520  ! j+1, j and j-1 in the second derivative. Make sure that
521  ! these coefficients are nonnegative.
522 
523  cnud = -rsacb2 * w(i, j, k, itu1) * cb3inv
524  cam = ttm * cnud
525  cap = ttp * cnud
526 
527  nutm = half * (w(i, j - 1, k, itu1) + w(i, j, k, itu1))
528  nutp = half * (w(i, j + 1, k, itu1) + w(i, j, k, itu1))
529  nu = rlv(i, j, k) / w(i, j, k, irho)
530  num = half * (rlv(i, j - 1, k) / w(i, j - 1, k, irho) + nu)
531  nup = half * (rlv(i, j + 1, k) / w(i, j + 1, k, irho) + nu)
532  cdm = (num + (one + rsacb2) * nutm) * ttm * cb3inv
533  cdp = (nup + (one + rsacb2) * nutp) * ttp * cb3inv
534 
535  c1m = max(cdm + cam, zero)
536  c1p = max(cdp + cap, zero)
537  c10 = c1m + c1p
538 
539  ! Update the residual for this cell and store the possible
540  ! coefficients for the matrix in b1, c1 and d1.
541 
542  scratch(i, j, k, idvt) = scratch(i, j, k, idvt) + c1m * w(i, j - 1, k, itu1) &
543  - c10 * w(i, j, k, itu1) + c1p * w(i, j + 1, k, itu1)
544 #ifndef USE_TAPENADE
545  b1 = -c1m
546  c1 = c10
547  d1 = -c1p
548 
549  ! Update the central jacobian. For nonboundary cells this
550  ! is simply c1. For boundary cells this is slightly more
551  ! complicated, because the boundary conditions are treated
552  ! implicitly and the off-diagonal terms b1 and d1 must be
553  ! taken into account.
554  ! The boundary conditions are only treated implicitly if
555  ! the diagonal dominance of the matrix is increased.
556 
557  if (j == 2) then
558  qq(i, j, k) = qq(i, j, k) + c1 &
559  - b1 * max(bmtj1(i, k, itu1, itu1), zero)
560  else if (j == jl) then
561  qq(i, j, k) = qq(i, j, k) + c1 &
562  - d1 * max(bmtj2(i, k, itu1, itu1), zero)
563  else
564  qq(i, j, k) = qq(i, j, k) + c1
565  end if
566 #endif
567 #ifdef TAPENADE_REVERSE
568  end do
569 #else
570  end do
571  end do
572  end do
573 #endif
574  !
575  ! Viscous terms in i-direction.
576  !
577 #ifdef TAPENADE_REVERSE
578  !$AD II-LOOP
579  do ii = 0, nx * ny * nz - 1
580  i = mod(ii, nx) + 2
581  j = mod(ii / nx, ny) + 2
582  k = ii / (nx * ny) + 2
583 #else
584  do k = 2, kl
585  do j = 2, jl
586  do i = 2, il
587 #endif
588  ! Compute the metrics in xi-direction, i.e. along the
589  ! line i = constant.
590 
591  voli = one / vol(i, j, k)
592  volmi = two / (vol(i, j, k) + vol(i - 1, j, k))
593  volpi = two / (vol(i, j, k) + vol(i + 1, j, k))
594 
595  xm = si(i - 1, j, k, 1) * volmi
596  ym = si(i - 1, j, k, 2) * volmi
597  zm = si(i - 1, j, k, 3) * volmi
598  xp = si(i, j, k, 1) * volpi
599  yp = si(i, j, k, 2) * volpi
600  zp = si(i, j, k, 3) * volpi
601 
602  xa = half * (si(i, j, k, 1) + si(i - 1, j, k, 1)) * voli
603  ya = half * (si(i, j, k, 2) + si(i - 1, j, k, 2)) * voli
604  za = half * (si(i, j, k, 3) + si(i - 1, j, k, 3)) * voli
605  ttm = xm * xa + ym * ya + zm * za
606  ttp = xp * xa + yp * ya + zp * za
607 
608  ! Computation of the viscous terms in xi-direction; note
609  ! that cross-derivatives are neglected, i.e. the mesh is
610  ! assumed to be orthogonal.
611  ! Furthermore, the grad(nu)**2 has been rewritten as
612  ! div(nu grad(nu)) - nu div(grad nu) to enhance stability.
613  ! The second derivative in xi-direction is constructed as
614  ! the central difference of the first order derivatives, i.e.
615  ! d^2/dxi^2 = d/dxi (d/dxi i+1/2 - d/dxi i-1/2).
616  ! In this way the metric can be taken into account.
617 
618  ! Compute the diffusion coefficients multiplying the nodes
619  ! i+1, i and i-1 in the second derivative. Make sure that
620  ! these coefficients are nonnegative.
621 
622  cnud = -rsacb2 * w(i, j, k, itu1) * cb3inv
623  cam = ttm * cnud
624  cap = ttp * cnud
625 
626  nutm = half * (w(i - 1, j, k, itu1) + w(i, j, k, itu1))
627  nutp = half * (w(i + 1, j, k, itu1) + w(i, j, k, itu1))
628  nu = rlv(i, j, k) / w(i, j, k, irho)
629  num = half * (rlv(i - 1, j, k) / w(i - 1, j, k, irho) + nu)
630  nup = half * (rlv(i + 1, j, k) / w(i + 1, j, k, irho) + nu)
631  cdm = (num + (one + rsacb2) * nutm) * ttm * cb3inv
632  cdp = (nup + (one + rsacb2) * nutp) * ttp * cb3inv
633 
634  c1m = max(cdm + cam, zero)
635  c1p = max(cdp + cap, zero)
636  c10 = c1m + c1p
637 
638  ! Update the residual for this cell and store the possible
639  ! coefficients for the matrix in b1, c1 and d1.
640 
641  scratch(i, j, k, idvt) = scratch(i, j, k, idvt) + c1m * w(i - 1, j, k, itu1) &
642  - c10 * w(i, j, k, itu1) + c1p * w(i + 1, j, k, itu1)
643 #ifndef USE_TAPENADE
644  b1 = -c1m
645  c1 = c10
646  d1 = -c1p
647 
648  ! Update the central jacobian. For nonboundary cells this
649  ! is simply c1. For boundary cells this is slightly more
650  ! complicated, because the boundary conditions are treated
651  ! implicitly and the off-diagonal terms b1 and d1 must be
652  ! taken into account.
653  ! The boundary conditions are only treated implicitly if
654  ! the diagonal dominance of the matrix is increased.
655 
656  if (i == 2) then
657  qq(i, j, k) = qq(i, j, k) + c1 &
658  - b1 * max(bmti1(j, k, itu1, itu1), zero)
659  else if (i == il) then
660  qq(i, j, k) = qq(i, j, k) + c1 &
661  - d1 * max(bmti2(j, k, itu1, itu1), zero)
662  else
663  qq(i, j, k) = qq(i, j, k) + c1
664  end if
665 #endif
666 #ifdef TAPENADE_REVERSE
667  end do
668 #else
669  end do
670  end do
671  end do
672 #endif
673  end subroutine saviscous
674 
675  subroutine saresscale
676 
677  !
678  ! Multiply the residual by the volume and store this in dw; this
679  ! * is done for monitoring reasons only. The multiplication with the
680  ! * volume is present to be consistent with the flow residuals; also
681  ! the negative value is taken, again to be consistent with the
682  ! * flow equations. Also multiply by iblank so that no updates occur
683  ! in holes or the overset boundary.
684  use blockpointers
685  implicit none
686 
687  ! Local variables
688  integer(kind=intType) :: i, j, k, ii
689  real(kind=realtype) :: rblank
690 
691 #ifdef TAPENADE_REVERSE
692  !$AD II-LOOP
693  do ii = 0, nx * ny * nz - 1
694  i = mod(ii, nx) + 2
695  j = mod(ii / nx, ny) + 2
696  k = ii / (nx * ny) + 2
697 #else
698  do k = 2, kl
699  do j = 2, jl
700  do i = 2, il
701 #endif
702  rblank = max(real(iblank(i, j, k), realtype), zero)
703  dw(i, j, k, itu1) = -volref(i, j, k) * scratch(i, j, k, idvt) * rblank
704 #ifdef TAPENADE_REVERSE
705  end do
706 #else
707  end do
708  end do
709  end do
710 #endif
711  end subroutine saresscale
712 
713 #ifndef USE_TAPENADE
714  subroutine sasolve
715  !
716  ! saSolve solves the turbulent transport equation for the
717  ! original Spalart-Allmaras model in a decoupled manner using
718  ! a diagonal dominant ADI-scheme.
719  use blockpointers
720  use inputiteration
721  use inputphysics
722  use paramturb
723  use turbutils
724  use turbcurvefits, only: curvetupyp
725  implicit none
726 
727  integer(kind=intType) :: i, j, k, nn, ii
728  real(kind=realtype), dimension(2:max(kl, il, jl)) :: bb, cc, dd, ff
729  real(kind=realtype) :: voli, volmi, volpi, xm, ym, zm, xp, yp, zp
730  real(kind=realtype) :: xa, ya, za, ttm, ttp, cnud, cam, cap
731  real(kind=realtype) :: nutm, nutp, num, nup, cdm, cdp
732  real(kind=realtype) :: c1m, c1p, c10, b1, c1, d1, qs
733  real(kind=realtype) :: uu, um, up, factor, f, tu1p(1), nu, rblank
734 
735  logical, dimension(2:jl, 2:kl), target :: flagI2, flagIl
736  logical, dimension(2:il, 2:kl), target :: flagJ2, flagJl
737  logical, dimension(2:il, 2:jl), target :: flagK2, flagKl
738  logical, dimension(:, :), pointer :: flag
739 
740  ! Initialize the wall function flags to .false.
741 
742  flagi2 = .false.
743  flagil = .false.
744  flagj2 = .false.
745  flagjl = .false.
746  flagk2 = .false.
747  flagkl = .false.
748 
749  ! Modify the rhs of the 1st internal cell, if wall functions
750  ! are used; their value is determined by the table.
751 
752  testwallfunctions: if (wallfunctions) then
753 
754  bocos: do nn = 1, nviscbocos
755 
756  ! Determine the block face on which the subface is located
757  ! and set some variables. As flag points to the entire array
758  ! flagI2, etc., its starting indices are the starting indices
759  ! of its target and not 1.
760 
761  select case (bcfaceid(nn))
762  case (imin)
763  flag => flagi2
764  ddw => dw(2, 1:, 1:, 1:); ddvt => scratch(2, 1:, 1:, idvt:)
765  ww => w(2, 1:, 1:, 1:); rrlv => rlv(2, 1:, 1:)
766  dd2wall => d2wall(2, :, :)
767 
768  case (imax)
769  flag => flagil
770  ddw => dw(il, 1:, 1:, 1:); ddvt => scratch(il, 1:, 1:, idvt:)
771  ww => w(il, 1:, 1:, 1:); rrlv => rlv(il, 1:, 1:)
772  dd2wall => d2wall(il, :, :)
773 
774  case (jmin)
775  flag => flagj2
776  ddw => dw(1:, 2, 1:, 1:); ddvt => scratch(1:, 2, 1:, idvt:)
777  ww => w(1:, 2, 1:, 1:); rrlv => rlv(1:, 2, 1:)
778  dd2wall => d2wall(:, 2, :)
779 
780  case (jmax)
781  flag => flagjl
782  ddw => dw(1:, jl, 1:, 1:); ddvt => scratch(1:, jl, 1:, idvt:)
783  ww => w(1:, jl, 1:, 1:); rrlv => rlv(1:, jl, 1:)
784  dd2wall => d2wall(:, jl, :)
785 
786  case (kmin)
787  flag => flagk2
788  ddw => dw(1:, 1:, 2, 1:); ddvt => scratch(1:, 1:, 2, idvt:)
789  ww => w(1:, 1:, 2, 1:); rrlv => rlv(1:, 1:, 2)
790  dd2wall => d2wall(:, :, 2)
791 
792  case (kmax)
793  flag => flagkl
794  ddw => dw(1:, 1:, kl, :); ddvt => scratch(1:, 1:, kl, idvt:)
795  ww => w(1:, 1:, kl, 1:); rrlv => rlv(1:, 1:, kl)
796  dd2wall => d2wall(:, :, kl)
797 
798  end select
799 
800  ! Loop over the owned faces of this subface. Therefore the
801  ! nodal range of BCData must be used. The offset of +1 is
802  ! present, because the starting index of the cell range is
803  ! 1 larger than the starting index of the nodal range.
804 
805  do j = (bcdata(nn)%jnBeg + 1), bcdata(nn)%jnEnd
806  do i = (bcdata(nn)%inBeg + 1), bcdata(nn)%inEnd
807 
808  ! Set ddw to zero.
809 
810  ddw(i, j, itu1) = zero
811 
812  ! Enforce nu tilde in the 1st internal cell from the
813  ! wall function table. There is an offset of -1 in the
814  ! wall distance. Note that the offset compared to the
815  ! current value must be stored, because dvt contains
816  ! the update. Also note that the curve fits contain the
817  ! non-dimensional value.
818 
819  yp = ww(i, j, irho) * dd2wall(i - 1, j - 1) &
820  * viscsubface(nn)%utau(i, j) / rrlv(i, j)
821 
822  call curvetupyp(tu1p, yp, itu1, itu1)
823  ddvt(i, j, 1) = tu1p(1) * rrlv(i, j) / ww(i, j, irho) - ww(i, j, itu1)
824 
825  ! Set the wall flag to .true.
826 
827  flag(i, j) = .true.
828 
829  end do
830  end do
831 
832  end do bocos
833  end if testwallfunctions
834 
835  ! For implicit relaxation take the local time step into account,
836  ! where dt is the inverse of the central jacobian times the cfl
837  ! number. The following system is solved:
838  ! (I/dt + cc + bb + dd)*dw = rhs, in which I/dt = cc/cfl. As in
839  ! the rest of the algorithm only the modified central jacobian is
840  ! used, stored it now.
841 
842  ! Compute the factor multiplying the central jacobian, which
843  ! is 1 + 1/cfl (implicit relaxation only).
844 
845  factor = one
846  if (turbrelax == turbrelaximplicit) &
847  factor = one + (one - alfaturb) / alfaturb
848 
849  do k = 2, kl
850  do j = 2, jl
851  do i = 2, il
852 
853  qq(i, j, k) = factor * qq(i, j, k)
854 
855  ! Set qq to 1 if the value is determined by the
856  ! wall function table.
857 
858  if ((i == 2 .and. flagi2(j, k)) .or. &
859  (i == il .and. flagil(j, k)) .or. &
860  (j == 2 .and. flagj2(i, k)) .or. &
861  (j == jl .and. flagjl(i, k)) .or. &
862  (k == 2 .and. flagk2(i, j)) .or. &
863  (k == kl .and. flagkl(i, j))) qq(i, j, k) = one
864 
865  end do
866  end do
867  end do
868 
869  ! Initialize the grid velocity to zero. This value will be used
870  ! if the block is not moving.
871 
872  qs = zero
873  !
874  ! dd-ADI step in j-direction. There is no particular reason to
875  ! start in j-direction, it just happened to be so. As we solve
876  ! in j-direction, the j-loop is the innermost loop.
877  !
878  do k = 2, kl
879  do i = 2, il
880  do j = 2, jl
881 
882  ! More or less the same code is executed here as above when
883  ! the residual was built. However, now the off-diagonal
884  ! terms for the dd-ADI must be built and stored. This could
885  ! have been done earlier, but then all the coefficients had
886  ! to be stored. To save memory, they are recomputed.
887  ! Consequently, see the j-loop to build the residual for
888  ! the comments.
889 
890  voli = one / vol(i, j, k)
891  volmi = two / (vol(i, j, k) + vol(i, j - 1, k))
892  volpi = two / (vol(i, j, k) + vol(i, j + 1, k))
893 
894  xm = sj(i, j - 1, k, 1) * volmi
895  ym = sj(i, j - 1, k, 2) * volmi
896  zm = sj(i, j - 1, k, 3) * volmi
897  xp = sj(i, j, k, 1) * volpi
898  yp = sj(i, j, k, 2) * volpi
899  zp = sj(i, j, k, 3) * volpi
900 
901  xa = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1)) * voli
902  ya = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2)) * voli
903  za = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3)) * voli
904  ttm = xm * xa + ym * ya + zm * za
905  ttp = xp * xa + yp * ya + zp * za
906 
907  cnud = -rsacb2 * w(i, j, k, itu1) * cb3inv
908  cam = ttm * cnud
909  cap = ttp * cnud
910 
911  ! Off-diagonal terms due to the diffusion terms
912  ! in j-direction.
913 
914  nutm = half * (w(i, j - 1, k, itu1) + w(i, j, k, itu1))
915  nutp = half * (w(i, j + 1, k, itu1) + w(i, j, k, itu1))
916  nu = rlv(i, j, k) / w(i, j, k, irho)
917  num = half * (rlv(i, j - 1, k) / w(i, j - 1, k, irho) + nu)
918  nup = half * (rlv(i, j + 1, k) / w(i, j + 1, k, irho) + nu)
919  cdm = (num + (one + rsacb2) * nutm) * ttm * cb3inv
920  cdp = (nup + (one + rsacb2) * nutp) * ttp * cb3inv
921 
922  c1m = max(cdm + cam, zero)
923  c1p = max(cdp + cap, zero)
924 
925  bb(j) = -c1m
926  dd(j) = -c1p
927 
928  ! Compute the grid velocity if present.
929  ! It is taken as the average of j and j-1,
930 
931  if (addgridvelocities) &
932  qs = half * (sfacej(i, j, k) + sfacej(i, j - 1, k)) * voli
933 
934  ! Off-diagonal terms due to the advection term in
935  ! j-direction. First order approximation.
936 
937  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
938  um = zero
939  up = zero
940  if (uu < zero) um = uu
941  if (uu > zero) up = uu
942 
943  bb(j) = bb(j) - up
944  dd(j) = dd(j) + um
945 
946  ! Store the central jacobian and rhs in cc and ff.
947  ! Multiply the off-diagonal terms and rhs by the iblank
948  ! value so the update determined for iblank = 0 is zero.
949 
950  rblank = max(real(iblank(i, j, k), realtype), zero)
951 
952  cc(j) = qq(i, j, k)
953  ff(j) = scratch(i, j, k, idvt) * rblank
954 
955  bb(j) = bb(j) * rblank
956  dd(j) = dd(j) * rblank
957 
958  ! Set the off diagonal terms to zero if the wall is flagged.
959 
960  if ((i == 2 .and. flagi2(j, k)) .or. &
961  (i == il .and. flagil(j, k)) .or. &
962  (j == 2 .and. flagj2(i, k)) .or. &
963  (j == jl .and. flagjl(i, k)) .or. &
964  (k == 2 .and. flagk2(i, j)) .or. &
965  (k == kl .and. flagkl(i, j))) then
966  bb(j) = zero
967  dd(j) = zero
968  end if
969 
970  end do
971 
972  ! Solve the tri-diagonal system in j-direction.
973  ! First the backward sweep to eliMinate the upper diagonal dd.
974 
975  do j = ny, 2, -1
976  f = dd(j) / cc(j + 1)
977  cc(j) = cc(j) - f * bb(j + 1)
978  ff(j) = ff(j) - f * ff(j + 1)
979  end do
980 
981  ! The matrix is now in lower block bi-diagonal form.
982  ! Perform a forward sweep to compute the solution.
983 
984  ff(2) = ff(2) / cc(2)
985  do j = 3, jl
986  ff(j) = ff(j) - bb(j) * ff(j - 1)
987  ff(j) = ff(j) / cc(j)
988  end do
989 
990  ! Determine the new rhs for the next direction.
991 
992  do j = 2, jl
993  scratch(i, j, k, idvt) = ff(j) * qq(i, j, k)
994  end do
995 
996  end do
997  end do
998  !
999  ! dd-ADI step in i-direction. As we solve in i-direction, the
1000  ! i-loop is the innermost loop.
1001  !
1002  do k = 2, kl
1003  do j = 2, jl
1004  do i = 2, il
1005 
1006  ! More or less the same code is executed here as above when
1007  ! the residual was built. However, now the off-diagonal
1008  ! terms for the dd-ADI must be built and stored. This could
1009  ! have been done earlier, but then all the coefficients had
1010  ! to be stored. To save memory, they are recomputed.
1011  ! Consequently, see the i-loop to build the residual for
1012  ! the comments.
1013 
1014  voli = one / vol(i, j, k)
1015  volmi = two / (vol(i, j, k) + vol(i - 1, j, k))
1016  volpi = two / (vol(i, j, k) + vol(i + 1, j, k))
1017 
1018  xm = si(i - 1, j, k, 1) * volmi
1019  ym = si(i - 1, j, k, 2) * volmi
1020  zm = si(i - 1, j, k, 3) * volmi
1021  xp = si(i, j, k, 1) * volpi
1022  yp = si(i, j, k, 2) * volpi
1023  zp = si(i, j, k, 3) * volpi
1024 
1025  xa = half * (si(i, j, k, 1) + si(i - 1, j, k, 1)) * voli
1026  ya = half * (si(i, j, k, 2) + si(i - 1, j, k, 2)) * voli
1027  za = half * (si(i, j, k, 3) + si(i - 1, j, k, 3)) * voli
1028  ttm = xm * xa + ym * ya + zm * za
1029  ttp = xp * xa + yp * ya + zp * za
1030 
1031  cnud = -rsacb2 * w(i, j, k, itu1) * cb3inv
1032  cam = ttm * cnud
1033  cap = ttp * cnud
1034 
1035  ! Off-diagonal terms due to the diffusion terms
1036  ! in i-direction.
1037 
1038  nutm = half * (w(i - 1, j, k, itu1) + w(i, j, k, itu1))
1039  nutp = half * (w(i + 1, j, k, itu1) + w(i, j, k, itu1))
1040  nu = rlv(i, j, k) / w(i, j, k, irho)
1041  num = half * (rlv(i - 1, j, k) / w(i - 1, j, k, irho) + nu)
1042  nup = half * (rlv(i + 1, j, k) / w(i + 1, j, k, irho) + nu)
1043  cdm = (num + (one + rsacb2) * nutm) * ttm * cb3inv
1044  cdp = (nup + (one + rsacb2) * nutp) * ttp * cb3inv
1045 
1046  c1m = max(cdm + cam, zero)
1047  c1p = max(cdp + cap, zero)
1048 
1049  bb(i) = -c1m
1050  dd(i) = -c1p
1051 
1052  ! Compute the grid velocity if present.
1053  ! It is taken as the average of i and i-1,
1054 
1055  if (addgridvelocities) &
1056  qs = half * (sfacei(i, j, k) + sfacei(i - 1, j, k)) * voli
1057 
1058  ! Off-diagonal terms due to the advection term in
1059  ! i-direction. First order approximation.
1060 
1061  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
1062  um = zero
1063  up = zero
1064  if (uu < zero) um = uu
1065  if (uu > zero) up = uu
1066 
1067  bb(i) = bb(i) - up
1068  dd(i) = dd(i) + um
1069 
1070  ! Store the central jacobian and rhs in cc and ff.
1071  ! Multiply the off-diagonal terms and rhs by the iblank
1072  ! value so the update determined for iblank = 0 is zero.
1073 
1074  rblank = max(real(iblank(i, j, k), realtype), zero)
1075 
1076  cc(i) = qq(i, j, k)
1077  ff(i) = scratch(i, j, k, idvt) * rblank
1078 
1079  bb(i) = bb(i) * rblank
1080  dd(i) = dd(i) * rblank
1081 
1082  ! Set the off diagonal terms to zero if the wall is flagged.
1083 
1084  if ((i == 2 .and. flagi2(j, k)) .or. &
1085  (i == il .and. flagil(j, k)) .or. &
1086  (j == 2 .and. flagj2(i, k)) .or. &
1087  (j == jl .and. flagjl(i, k)) .or. &
1088  (k == 2 .and. flagk2(i, j)) .or. &
1089  (k == kl .and. flagkl(i, j))) then
1090  bb(i) = zero
1091  dd(i) = zero
1092  end if
1093 
1094  end do
1095 
1096  ! Solve the tri-diagonal system in i-direction.
1097  ! First the backward sweep to eliMinate the upper diagonal dd.
1098 
1099  do i = nx, 2, -1
1100  f = dd(i) / cc(i + 1)
1101  cc(i) = cc(i) - f * bb(i + 1)
1102  ff(i) = ff(i) - f * ff(i + 1)
1103  end do
1104 
1105  ! The matrix is now in lower block bi-diagonal form.
1106  ! Perform a forward sweep to compute the solution.
1107 
1108  ff(2) = ff(2) / cc(2)
1109  do i = 3, il
1110  ff(i) = ff(i) - bb(i) * ff(i - 1)
1111  ff(i) = ff(i) / cc(i)
1112  end do
1113 
1114  ! Determine the new rhs for the next direction.
1115 
1116  do i = 2, il
1117  scratch(i, j, k, idvt) = ff(i) * qq(i, j, k)
1118  end do
1119 
1120  end do
1121  end do
1122  !
1123  ! dd-ADI step in k-direction. As we solve in k-direction, the
1124  ! k-loop is the innermost loop.
1125  !
1126  do j = 2, jl
1127  do i = 2, il
1128  do k = 2, kl
1129 
1130  ! More or less the same code is executed here as above when
1131  ! the residual was built. However, now the off-diagonal
1132  ! terms for the dd-ADI must be built and stored. This could
1133  ! have been done earlier, but then all the coefficients had
1134  ! to be stored. To save memory, they are recomputed.
1135  ! Consequently, see the k-loop to build the residual for
1136  ! the comments.
1137 
1138  voli = one / vol(i, j, k)
1139  volmi = two / (vol(i, j, k) + vol(i, j, k - 1))
1140  volpi = two / (vol(i, j, k) + vol(i, j, k + 1))
1141 
1142  xm = sk(i, j, k - 1, 1) * volmi
1143  ym = sk(i, j, k - 1, 2) * volmi
1144  zm = sk(i, j, k - 1, 3) * volmi
1145  xp = sk(i, j, k, 1) * volpi
1146  yp = sk(i, j, k, 2) * volpi
1147  zp = sk(i, j, k, 3) * volpi
1148 
1149  xa = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1)) * voli
1150  ya = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2)) * voli
1151  za = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3)) * voli
1152  ttm = xm * xa + ym * ya + zm * za
1153  ttp = xp * xa + yp * ya + zp * za
1154 
1155  cnud = -rsacb2 * w(i, j, k, itu1) * cb3inv
1156  cam = ttm * cnud
1157  cap = ttp * cnud
1158 
1159  ! Off-diagonal terms due to the diffusion terms
1160  ! in k-direction.
1161 
1162  nutm = half * (w(i, j, k - 1, itu1) + w(i, j, k, itu1))
1163  nutp = half * (w(i, j, k + 1, itu1) + w(i, j, k, itu1))
1164  nu = rlv(i, j, k) / w(i, j, k, irho)
1165  num = half * (rlv(i, j, k - 1) / w(i, j, k - 1, irho) + nu)
1166  nup = half * (rlv(i, j, k + 1) / w(i, j, k + 1, irho) + nu)
1167  cdm = (num + (one + rsacb2) * nutm) * ttm * cb3inv
1168  cdp = (nup + (one + rsacb2) * nutp) * ttp * cb3inv
1169 
1170  c1m = max(cdm + cam, zero)
1171  c1p = max(cdp + cap, zero)
1172 
1173  bb(k) = -c1m
1174  dd(k) = -c1p
1175 
1176  ! Compute the grid velocity if present.
1177  ! It is taken as the average of k and k-1,
1178 
1179  if (addgridvelocities) &
1180  qs = half * (sfacek(i, j, k) + sfacek(i, j, k - 1)) * voli
1181 
1182  ! Off-diagonal terms due to the advection term in
1183  ! k-direction. First order approximation.
1184 
1185  uu = xa * w(i, j, k, ivx) + ya * w(i, j, k, ivy) + za * w(i, j, k, ivz) - qs
1186  um = zero
1187  up = zero
1188  if (uu < zero) um = uu
1189  if (uu > zero) up = uu
1190 
1191  bb(k) = bb(k) - up
1192  dd(k) = dd(k) + um
1193 
1194  ! Store the central jacobian and rhs in cc and ff.
1195  ! Multiply the off-diagonal terms and rhs by the iblank
1196  ! value so the update determined for iblank = 0 is zero.
1197 
1198  rblank = max(real(iblank(i, j, k), realtype), zero)
1199 
1200  cc(k) = qq(i, j, k)
1201  ff(k) = scratch(i, j, k, idvt) * rblank
1202 
1203  bb(k) = bb(k) * rblank
1204  dd(k) = dd(k) * rblank
1205 
1206  ! Set the off diagonal terms to zero if the wall is flagged.
1207 
1208  if ((i == 2 .and. flagi2(j, k)) .or. &
1209  (i == il .and. flagil(j, k)) .or. &
1210  (j == 2 .and. flagj2(i, k)) .or. &
1211  (j == jl .and. flagjl(i, k)) .or. &
1212  (k == 2 .and. flagk2(i, j)) .or. &
1213  (k == kl .and. flagkl(i, j))) then
1214  bb(k) = zero
1215  dd(k) = zero
1216  end if
1217 
1218  end do
1219 
1220  ! Solve the tri-diagonal system in k-direction.
1221  ! First the backward sweep to eliMinate the upper diagonal dd.
1222 
1223  do k = nz, 2, -1
1224  f = dd(k) / cc(k + 1)
1225  cc(k) = cc(k) - f * bb(k + 1)
1226  ff(k) = ff(k) - f * ff(k + 1)
1227  end do
1228 
1229  ! The matrix is now in lower block bi-diagonal form.
1230  ! Perform a forward sweep to compute the solution.
1231 
1232  ff(2) = ff(2) / cc(2)
1233  do k = 3, kl
1234  ff(k) = ff(k) - bb(k) * ff(k - 1)
1235  ff(k) = ff(k) / cc(k)
1236  end do
1237 
1238  ! Store the update in dvt.
1239 
1240  do k = 2, kl
1241  scratch(i, j, k, idvt) = ff(k)
1242  end do
1243 
1244  end do
1245  end do
1246  !
1247  ! Update the turbulent variables. For explicit relaxation the
1248  ! update must be relaxed; for implicit relaxation this has been
1249  ! done via the time step.
1250  !
1251  factor = one
1252  if (turbrelax == turbrelaxexplicit) factor = alfaturb
1253 
1254  do k = 2, kl
1255  do j = 2, jl
1256  do i = 2, il
1257  w(i, j, k, itu1) = w(i, j, k, itu1) + factor * scratch(i, j, k, idvt)
1258  w(i, j, k, itu1) = max(w(i, j, k, itu1), zero)
1259  end do
1260  end do
1261  end do
1262 
1263  end subroutine sasolve
1264 #endif
1265 end module sa
bcdata
Definition: BCData.F90:1
blockpointers::bmtk2
real(kind=realtype), dimension(:, :, :, :), pointer bmtk2
Definition: blockPointers.F90:162
blockpointers::sfacek
real(kind=realtype), dimension(:, :, :), pointer sfacek
Definition: blockPointers.F90:119
blockpointers::addgridvelocities
logical addgridvelocities
Definition: blockPointers.F90:117
blockpointers::jl
integer(kind=inttype) jl
Definition: blockPointers.F90:38
blockpointers::bmti1
real(kind=realtype), dimension(:, :, :, :), pointer bmti1
Definition: blockPointers.F90:157
blockpointers::nviscbocos
integer(kind=inttype) nviscbocos
Definition: blockPointers.F90:47
blockpointers::bmtj1
real(kind=realtype), dimension(:, :, :, :), pointer bmtj1
Definition: blockPointers.F90:159
blockpointers::nx
integer(kind=inttype) nx
Definition: blockPointers.F90:38
blockpointers::bmti2
real(kind=realtype), dimension(:, :, :, :), pointer bmti2
Definition: blockPointers.F90:158
blockpointers::ny
integer(kind=inttype) ny
Definition: blockPointers.F90:38
blockpointers::w
real(kind=realtype), dimension(:, :, :, :), pointer w
Definition: blockPointers.F90:120
blockpointers::scratch
real(kind=realtype), dimension(:, :, :, :), pointer scratch
Definition: blockPointers.F90:129
blockpointers::sfacei
real(kind=realtype), dimension(:, :, :), pointer sfacei
Definition: blockPointers.F90:119
blockpointers::viscsubface
type(viscsubfacetype), dimension(:), pointer viscsubface
Definition: blockPointers.F90:85
blockpointers::d2wall
real(kind=realtype), dimension(:, :, :), pointer d2wall
Definition: blockPointers.F90:154
blockpointers::iblank
integer(kind=inttype), dimension(:, :, :), pointer iblank
Definition: blockPointers.F90:71
blockpointers::rlv
real(kind=realtype), dimension(:, :, :), pointer rlv
Definition: blockPointers.F90:125
blockpointers::bcfaceid
integer(kind=inttype), dimension(:), pointer bcfaceid
Definition: blockPointers.F90:50
blockpointers::si
real(kind=realtype), dimension(:, :, :, :), pointer si
Definition: blockPointers.F90:96
blockpointers::volref
real(kind=realtype), dimension(:, :, :), pointer volref
Definition: blockPointers.F90:98
blockpointers::sj
real(kind=realtype), dimension(:, :, :, :), pointer sj
Definition: blockPointers.F90:96
blockpointers::sectionid
integer(kind=inttype) sectionid
Definition: blockPointers.F90:30
blockpointers::bmtj2
real(kind=realtype), dimension(:, :, :, :), pointer bmtj2
Definition: blockPointers.F90:160
blockpointers::dw
real(kind=realtype), dimension(:, :, :, :), pointer dw
Definition: blockPointers.F90:128
blockpointers::sk
real(kind=realtype), dimension(:, :, :, :), pointer sk
Definition: blockPointers.F90:96
blockpointers::vol
real(kind=realtype), dimension(:, :, :), pointer vol
Definition: blockPointers.F90:97
blockpointers::sfacej
real(kind=realtype), dimension(:, :, :), pointer sfacej
Definition: blockPointers.F90:119
blockpointers::nz
integer(kind=inttype) nz
Definition: blockPointers.F90:38
blockpointers::bmtk1
real(kind=realtype), dimension(:, :, :, :), pointer bmtk1
Definition: blockPointers.F90:161
blockpointers::kl
integer(kind=inttype) kl
Definition: blockPointers.F90:38
blockpointers::il
integer(kind=inttype) il
Definition: blockPointers.F90:38
constants::strain
integer(kind=inttype), parameter strain
Definition: constants.F90:137
constants::zero
real(kind=realtype), parameter zero
Definition: constants.F90:71
constants::three
real(kind=realtype), parameter three
Definition: constants.F90:74
constants::imax
integer(kind=inttype), parameter imax
Definition: constants.F90:293
constants::idvt
integer, parameter idvt
Definition: constants.F90:51
constants::kmin
integer(kind=inttype), parameter kmin
Definition: constants.F90:296
constants::third
real(kind=realtype), parameter third
Definition: constants.F90:81
constants::jmax
integer(kind=inttype), parameter jmax
Definition: constants.F90:295
constants::six
real(kind=realtype), parameter six
Definition: constants.F90:77
constants::itu1
integer, parameter itu1
Definition: constants.F90:40
constants::irho
integer, parameter irho
Definition: constants.F90:34
constants::vorticity
integer(kind=inttype), parameter vorticity
Definition: constants.F90:137
constants::ivx
integer, parameter ivx
Definition: constants.F90:35
constants::one
real(kind=realtype), parameter one
Definition: constants.F90:72
constants::half
real(kind=realtype), parameter half
Definition: constants.F90:80
constants::turbrelaximplicit
integer(kind=inttype), parameter turbrelaximplicit
Definition: constants.F90:223
constants::turbrelaxexplicit
integer(kind=inttype), parameter turbrelaxexplicit
Definition: constants.F90:223
constants::katolaunder
integer(kind=inttype), parameter katolaunder
Definition: constants.F90:137
constants::imin
integer(kind=inttype), parameter imin
Definition: constants.F90:292
constants::ivz
integer, parameter ivz
Definition: constants.F90:37
constants::two
real(kind=realtype), parameter two
Definition: constants.F90:73
constants::fourth
real(kind=realtype), parameter fourth
Definition: constants.F90:82
constants::kmax
integer(kind=inttype), parameter kmax
Definition: constants.F90:297
constants::ivy
integer, parameter ivy
Definition: constants.F90:36
constants::jmin
integer(kind=inttype), parameter jmin
Definition: constants.F90:294
constants::sixth
real(kind=realtype), parameter sixth
Definition: constants.F90:83
flowvarrefstate::timeref
real(kind=realtype) timeref
Definition: flowVarRefState.F90:40
inputiteration::alfaturb
real(kind=realtype) alfaturb
Definition: inputParam.F90:276
inputiteration::turbrelax
integer(kind=inttype) turbrelax
Definition: inputParam.F90:269
inputphysics::useft2sa
logical useft2sa
Definition: inputParam.F90:587
inputphysics::turbprod
integer(kind=inttype) turbprod
Definition: inputParam.F90:584
inputphysics::wallfunctions
logical wallfunctions
Definition: inputParam.F90:589
inputphysics::userotationsa
logical userotationsa
Definition: inputParam.F90:587
inputtimespectral::ntimeintervalsspectral
integer(kind=inttype) ntimeintervalsspectral
Definition: inputParam.F90:645
iteration::currentlevel
integer(kind=inttype) currentlevel
Definition: iteration.f90:18
paramturb::rsacw1
real(kind=realtype), parameter rsacw1
Definition: paramTurb.F90:18
paramturb::rsak
real(kind=realtype), parameter rsak
Definition: paramTurb.F90:13
paramturb::rsact4
real(kind=realtype), parameter rsact4
Definition: paramTurb.F90:25
paramturb::rsacrot
real(kind=realtype), parameter rsacrot
Definition: paramTurb.F90:26
paramturb::rsacb2
real(kind=realtype), parameter rsacb2
Definition: paramTurb.F90:15
paramturb::rsact3
real(kind=realtype), parameter rsact3
Definition: paramTurb.F90:24
paramturb::rsacw3
real(kind=realtype), parameter rsacw3
Definition: paramTurb.F90:21
paramturb::rsacw2
real(kind=realtype), parameter rsacw2
Definition: paramTurb.F90:20
paramturb::rsacv1
real(kind=realtype), parameter rsacv1
Definition: paramTurb.F90:17
paramturb::rsacb3
real(kind=realtype), parameter rsacb3
Definition: paramTurb.F90:16
paramturb::rsacb1
real(kind=realtype), parameter rsacb1
Definition: paramTurb.F90:14
sa
Definition: sa.F90:5
sa::kar2inv
real(kind=realtype) kar2inv
Definition: sa.F90:8
sa::rrlv
real(kind=realtype), dimension(:, :), pointer rrlv
Definition: sa.F90:11
sa::cb3inv
real(kind=realtype) cb3inv
Definition: sa.F90:8
sa::ddvt
real(kind=realtype), dimension(:, :, :), pointer ddvt
Definition: sa.F90:10
sa::dd2wall
real(kind=realtype), dimension(:, :), pointer dd2wall
Definition: sa.F90:12
sa::sasolve
subroutine sasolve
Definition: sa.F90:715
sa::cw36
real(kind=realtype) cw36
Definition: sa.F90:8
sa::sasource
subroutine sasource
Definition: sa.F90:89
sa::saresscale
subroutine saresscale
Definition: sa.F90:676
sa::cv13
real(kind=realtype) cv13
Definition: sa.F90:8
sa::saviscous
subroutine saviscous
Definition: sa.F90:345
sa::ddw
real(kind=realtype), dimension(:, :, :), pointer ddw
Definition: sa.F90:10
sa::sa_block
subroutine sa_block(resOnly)
Definition: sa.F90:17
sa::qq
real(kind=realtype), dimension(:, :, :), allocatable qq
Definition: sa.F90:9
sa::ww
real(kind=realtype), dimension(:, :, :), pointer ww
Definition: sa.F90:10
section::sections
type(sectiontype), dimension(:), allocatable sections
Definition: section.f90:46
sst
Definition: SST.F90:1
turbbcroutines::bcturbtreatment
subroutine bcturbtreatment
Definition: turbBCRoutines.F90:663
turbbcroutines::applyallturbbcthisblock
subroutine applyallturbbcthisblock(secondHalo)
Definition: turbBCRoutines.F90:50
turbcurvefits::curvetupyp
subroutine curvetupyp(tup, yp, ntu1, ntu2)
Definition: turbCurveFits.F90:104
turbutils::unsteadyturbterm
subroutine unsteadyturbterm(mAdv, nAdv, offset, qq)
Definition: turbUtils.F90:411
turbutils::saeddyviscosity
subroutine saeddyviscosity(iBeg, iEnd, jBeg, jEnd, kBeg, kEnd)
Definition: turbUtils.F90:657
turbutils::turbadvection
subroutine turbadvection(mAdv, nAdv, offset, qq)
Definition: turbUtils.F90:826
Generated by doxygen 1.9.1