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