ADflow  v1.0
ADflow is a finite volume RANS solver tailored for gradient-based aerodynamic design optimization.
adjointExtra.F90
Go to the documentation of this file.
1 module adjointextra
2 
3 contains
4 
5  subroutine volume_block
6 
7  ! This is COPY of metric.f90. It was necessary to copy this file
8  ! since there is debugging stuff in the original that is not
9  ! necessary for AD.
10  use constants
11  use blockpointers
12  use cgnsgrid
13  use communication
14  use inputtimespectral
15 
16  implicit none
17  !
18  ! Local parameter.
19  !
20  real(kind=realtype), parameter :: thresvolume = 1.e-2_realtype
21  real(kind=realtype), parameter :: halocellratio = 1e-10_realtype
22  !
23  ! Local variables.
24  !
25  integer(kind=intType) :: i, j, k, n, m, l, ii
26  integer(kind=intType) :: mm
27  real(kind=realtype) :: fact, mult
28  real(kind=realtype) :: xp, yp, zp, vp1, vp2, vp3, vp4, vp5, vp6
29  real(kind=realtype) :: xxp, yyp, zzp
30  real(kind=realtype), dimension(3) :: v1, v2
31 
32  ! Compute the volumes. The hexahedron is split into 6 pyramids
33  ! whose volumes are computed. The volume is positive for a
34  ! right handed block.
35  ! Initialize the volumes to zero. The reasons is that the second
36  ! level halo's must be initialized to zero and for convenience
37  ! all the volumes are set to zero.
38 
39  vol = zero
40 
41  do k = 1, ke
42  n = k - 1
43  do j = 1, je
44  m = j - 1
45  do i = 1, ie
46 
47  l = i - 1
48 
49  ! Compute the coordinates of the center of gravity.
50 
51  xp = eighth * (x(i, j, k, 1) + x(i, m, k, 1) &
52  + x(i, m, n, 1) + x(i, j, n, 1) &
53  + x(l, j, k, 1) + x(l, m, k, 1) &
54  + x(l, m, n, 1) + x(l, j, n, 1))
55  yp = eighth * (x(i, j, k, 2) + x(i, m, k, 2) &
56  + x(i, m, n, 2) + x(i, j, n, 2) &
57  + x(l, j, k, 2) + x(l, m, k, 2) &
58  + x(l, m, n, 2) + x(l, j, n, 2))
59  zp = eighth * (x(i, j, k, 3) + x(i, m, k, 3) &
60  + x(i, m, n, 3) + x(i, j, n, 3) &
61  + x(l, j, k, 3) + x(l, m, k, 3) &
62  + x(l, m, n, 3) + x(l, j, n, 3))
63 
64  ! Compute the volumes of the 6 sub pyramids. The
65  ! arguments of volpym must be such that for a (regular)
66  ! right handed hexahedron all volumes are positive.
67 
68  call volpym(x(i, j, k, 1), x(i, j, k, 2), x(i, j, k, 3), &
69  x(i, j, n, 1), x(i, j, n, 2), x(i, j, n, 3), &
70  x(i, m, n, 1), x(i, m, n, 2), x(i, m, n, 3), &
71  x(i, m, k, 1), x(i, m, k, 2), x(i, m, k, 3), vp1)
72 
73  call volpym(x(l, j, k, 1), x(l, j, k, 2), x(l, j, k, 3), &
74  x(l, m, k, 1), x(l, m, k, 2), x(l, m, k, 3), &
75  x(l, m, n, 1), x(l, m, n, 2), x(l, m, n, 3), &
76  x(l, j, n, 1), x(l, j, n, 2), x(l, j, n, 3), vp2)
77 
78  call volpym(x(i, j, k, 1), x(i, j, k, 2), x(i, j, k, 3), &
79  x(l, j, k, 1), x(l, j, k, 2), x(l, j, k, 3), &
80  x(l, j, n, 1), x(l, j, n, 2), x(l, j, n, 3), &
81  x(i, j, n, 1), x(i, j, n, 2), x(i, j, n, 3), vp3)
82 
83  call volpym(x(i, m, k, 1), x(i, m, k, 2), x(i, m, k, 3), &
84  x(i, m, n, 1), x(i, m, n, 2), x(i, m, n, 3), &
85  x(l, m, n, 1), x(l, m, n, 2), x(l, m, n, 3), &
86  x(l, m, k, 1), x(l, m, k, 2), x(l, m, k, 3), vp4)
87 
88  call volpym(x(i, j, k, 1), x(i, j, k, 2), x(i, j, k, 3), &
89  x(i, m, k, 1), x(i, m, k, 2), x(i, m, k, 3), &
90  x(l, m, k, 1), x(l, m, k, 2), x(l, m, k, 3), &
91  x(l, j, k, 1), x(l, j, k, 2), x(l, j, k, 3), vp5)
92 
93  call volpym(x(i, j, n, 1), x(i, j, n, 2), x(i, j, n, 3), &
94  x(l, j, n, 1), x(l, j, n, 2), x(l, j, n, 3), &
95  x(l, m, n, 1), x(l, m, n, 2), x(l, m, n, 3), &
96  x(i, m, n, 1), x(i, m, n, 2), x(i, m, n, 3), vp6)
97 
98  ! Set the volume to 1/6 of the sum of the volumes of the
99  ! pyramid. Remember that volpym computes 6 times the
100  ! volume.
101 
102  vol(i, j, k) = sixth * (vp1 + vp2 + vp3 + vp4 + vp5 + vp6)
103 
104  ! Set the volume to the absolute value.
105  vol(i, j, k) = abs(vol(i, j, k))
106  end do
107  end do
108  end do
109 
110  ! Some additional safety stuff for halo volumes.
111 
112  do k = 2, kl
113  do j = 2, jl
114  if (vol(1, j, k) / vol(2, j, k) < halocellratio) then
115  vol(1, j, k) = vol(2, j, k)
116  end if
117  if (vol(ie, j, k) / vol(il, j, k) < halocellratio) then
118  vol(ie, j, k) = vol(il, j, k)
119  end if
120  end do
121  end do
122 
123  do k = 2, kl
124  do i = 1, ie
125  if (vol(i, 1, k) / vol(i, 2, k) < halocellratio) then
126  vol(i, 1, k) = vol(i, 2, k)
127  end if
128  if (vol(i, je, k) / vol(i, jl, k) < halocellratio) then
129  vol(i, je, k) = vol(i, jl, k)
130  end if
131  end do
132  end do
133 
134  do j = 1, je
135  do i = 1, ie
136  if (vol(i, j, 1) / vol(i, j, 2) < halocellratio) then
137  vol(i, j, 1) = vol(i, j, 2)
138  end if
139  if (vol(i, j, ke) / vol(i, j, kl) < halocellratio) then
140  vol(i, j, ke) = vol(i, j, kl)
141  end if
142  end do
143  end do
144 
145  contains
146 
147  subroutine volpym(xa, ya, za, xb, yb, zb, xc, yc, zc, xd, yd, zd, volume)
148  !
149  ! volpym computes 6 times the volume of a pyramid. Node p,
150  ! whose coordinates are set in the subroutine metric itself,
151  ! is the top node and a-b-c-d is the quadrilateral surface.
152  ! It is assumed that the cross product vCa * vDb points in
153  ! the direction of the top node. Here vCa is the diagonal
154  ! running from node c to node a and vDb the diagonal from
155  ! node d to node b.
156  !
157  use precision
158  implicit none
159  !
160  ! Function type.
161  !
162  real(kind=realtype) :: volume
163  !
164  ! Function arguments.
165  !
166  real(kind=realtype), intent(in) :: xa, ya, za, xb, yb, zb
167  real(kind=realtype), intent(in) :: xc, yc, zc, xd, yd, zd
168 
169  volume = (xp - fourth * (xa + xb + xc + xd)) &
170  * ((ya - yc) * (zb - zd) - (za - zc) * (yb - yd)) + &
171  (yp - fourth * (ya + yb + yc + yd)) &
172  * ((za - zc) * (xb - xd) - (xa - xc) * (zb - zd)) + &
173  (zp - fourth * (za + zb + zc + zd)) &
174  * ((xa - xc) * (yb - yd) - (ya - yc) * (xb - xd))
175 
176  end subroutine volpym
177  end subroutine volume_block
178 
179  subroutine metric_block
180  use constants
181  use blockpointers
182  implicit none
183 
184  ! Local variables.
185  integer(kind=intType) :: i, j, k, n, m, l, ii
186  real(kind=realtype) :: fact
187  real(kind=realtype) :: xxp, yyp, zzp
188  real(kind=realtype), dimension(3) :: v1, v2
189 
190  ! Set the factor in the surface normals computation. For a
191  ! left handed block this factor is negative, such that the
192  ! normals still point in the direction of increasing index.
193  ! The formulae used later on assume a right handed block
194  ! and fact is used to correct this for a left handed block,
195  ! as well as the scaling factor of 0.5
196 
197  if (righthanded) then
198  fact = half
199  else
200  fact = -half
201  end if
202 
203  !
204  ! Computation of the face normals in i-, j- and k-direction.
205  ! Formula's are valid for a right handed block; for a left
206  ! handed block the correct orientation is obtained via fact.
207  ! The normals point in the direction of increasing index.
208  ! The absolute value of fact is 0.5, because the cross
209  ! product of the two diagonals is twice the normal vector.
210  ! Note that also the normals of the first level halo cells
211  ! are computed. These are needed for the viscous fluxes.
212  !
213  ! Projected areas of cell faces in the i direction.
214  !$AD II-LOOP
215  do ii = 0, ke * je * (ie + 1) - 1
216  i = mod(ii, ie + 1) + 0 ! 0:ie
217  j = mod(ii / (ie + 1), je) + 1 !1:je
218  k = ii / ((ie + 1) * je) + 1 !1:ke
219 
220  n = k - 1
221  m = j - 1
222 
223  ! Determine the two diagonal vectors of the face.
224 
225  v1(1) = x(i, j, n, 1) - x(i, m, k, 1)
226  v1(2) = x(i, j, n, 2) - x(i, m, k, 2)
227  v1(3) = x(i, j, n, 3) - x(i, m, k, 3)
228 
229  v2(1) = x(i, j, k, 1) - x(i, m, n, 1)
230  v2(2) = x(i, j, k, 2) - x(i, m, n, 2)
231  v2(3) = x(i, j, k, 3) - x(i, m, n, 3)
232 
233  ! The face normal, which is the cross product of the two
234  ! diagonal vectors times fact; remember that fact is
235  ! either -0.5 or 0.5.
236 
237  si(i, j, k, 1) = fact * (v1(2) * v2(3) - v1(3) * v2(2))
238  si(i, j, k, 2) = fact * (v1(3) * v2(1) - v1(1) * v2(3))
239  si(i, j, k, 3) = fact * (v1(1) * v2(2) - v1(2) * v2(1))
240  end do
241 
242  ! Projected areas of cell faces in the j direction
243  !$AD II-LOOP
244  do ii = 0, ke * (je + 1) * ie - 1
245  i = mod(ii, ie) + 1 ! 1:ie
246  j = mod(ii / ie, je + 1) + 0 !0:je
247  k = ii / (ie * (je + 1)) + 1 !1:ke
248  n = k - 1
249  l = i - 1
250 
251  ! Determine the two diagonal vectors of the face.
252 
253  v1(1) = x(i, j, n, 1) - x(l, j, k, 1)
254  v1(2) = x(i, j, n, 2) - x(l, j, k, 2)
255  v1(3) = x(i, j, n, 3) - x(l, j, k, 3)
256 
257  v2(1) = x(l, j, n, 1) - x(i, j, k, 1)
258  v2(2) = x(l, j, n, 2) - x(i, j, k, 2)
259  v2(3) = x(l, j, n, 3) - x(i, j, k, 3)
260 
261  ! The face normal, which is the cross product of the two
262  ! diagonal vectors times fact; remember that fact is
263  ! either -0.5 or 0.5.
264 
265  sj(i, j, k, 1) = fact * (v1(2) * v2(3) - v1(3) * v2(2))
266  sj(i, j, k, 2) = fact * (v1(3) * v2(1) - v1(1) * v2(3))
267  sj(i, j, k, 3) = fact * (v1(1) * v2(2) - v1(2) * v2(1))
268 
269  end do
270 
271  ! Projected areas of cell faces in the k direction.
272  !$AD II-LOOP
273  do ii = 0, (ke + 1) * je * ie - 1
274  i = mod(ii, ie) + 1 ! 1:ie
275  j = mod(ii / ie, je) + 1 !1:je
276  k = ii / (ie * je) + 0 !0:ke
277  m = j - 1
278  l = i - 1
279 
280  ! Determine the two diagonal vectors of the face.
281 
282  v1(1) = x(i, j, k, 1) - x(l, m, k, 1)
283  v1(2) = x(i, j, k, 2) - x(l, m, k, 2)
284  v1(3) = x(i, j, k, 3) - x(l, m, k, 3)
285 
286  v2(1) = x(l, j, k, 1) - x(i, m, k, 1)
287  v2(2) = x(l, j, k, 2) - x(i, m, k, 2)
288  v2(3) = x(l, j, k, 3) - x(i, m, k, 3)
289 
290  ! The face normal, which is the cross product of the two
291  ! diagonal vectors times fact; remember that fact is
292  ! either -0.5 or 0.5.
293 
294  sk(i, j, k, 1) = fact * (v1(2) * v2(3) - v1(3) * v2(2))
295  sk(i, j, k, 2) = fact * (v1(3) * v2(1) - v1(1) * v2(3))
296  sk(i, j, k, 3) = fact * (v1(1) * v2(2) - v1(2) * v2(1))
297  end do
298  end subroutine metric_block
299 
300  subroutine boundarynormals
301 
302  ! The unit normals on the boundary faces. These always point
303  ! out of the domain, so a multiplication by -1 is needed for
304  ! the iMin, jMin and kMin boundaries.
305  !
306  use constants
307  use blockpointers
308  use cgnsgrid
309  use communication
310  use inputtimespectral
311  implicit none
312 
313  ! Local variables.
314  integer(kind=intType) :: i, j, ii
315  integer(kind=intType) :: mm
316  real(kind=realtype) :: fact, mult
317  real(kind=realtype) :: xxp, yyp, zzp
318 
319  !Loop over the boundary subfaces of this block.
320  !$AD II-LOOP
321  bocoloop: do mm = 1, nbocos
322 
323  ! Loop over the boundary faces of the subface.
324  !$AD II-LOOP
325  do ii = 0, (bcdata(mm)%jcEnd - bcdata(mm)%jcBeg + 1) * (bcdata(mm)%icEnd - bcdata(mm)%icBeg + 1) - 1
326  i = mod(ii, (bcdata(mm)%icEnd - bcdata(mm)%icBeg + 1)) + bcdata(mm)%icBeg
327  j = ii / (bcdata(mm)%icEnd - bcdata(mm)%icBeg + 1) + bcdata(mm)%jcBeg
328 
329  select case (bcfaceid(mm))
330  case (imin)
331  mult = -one
332  xxp = si(1, i, j, 1); yyp = si(1, i, j, 2); zzp = si(1, i, j, 3)
333  case (imax)
334  mult = one
335  xxp = si(il, i, j, 1); yyp = si(il, i, j, 2); zzp = si(il, i, j, 3)
336  case (jmin)
337  mult = -one
338  xxp = sj(i, 1, j, 1); yyp = sj(i, 1, j, 2); zzp = sj(i, 1, j, 3)
339  case (jmax)
340  mult = one
341  xxp = sj(i, jl, j, 1); yyp = sj(i, jl, j, 2); zzp = sj(i, jl, j, 3)
342  case (kmin)
343  mult = -one
344  xxp = sk(i, j, 1, 1); yyp = sk(i, j, 1, 2); zzp = sk(i, j, 1, 3)
345  case (kmax)
346  mult = one
347  xxp = sk(i, j, kl, 1); yyp = sk(i, j, kl, 2); zzp = sk(i, j, kl, 3)
348  end select
349 
350  ! Compute the inverse of the length of the normal vector
351  ! and possibly correct for inward pointing.
352 
353  fact = sqrt(xxp * xxp + yyp * yyp + zzp * zzp)
354  if (fact > zero) fact = mult / fact
355 
356  ! Compute the unit normal.
357 
358  bcdata(mm)%norm(i, j, 1) = fact * xxp
359  bcdata(mm)%norm(i, j, 2) = fact * yyp
360  bcdata(mm)%norm(i, j, 3) = fact * zzp
361  end do
362  end do bocoloop
363  end subroutine boundarynormals
364 
365  subroutine xhalo_block
366  !
367  ! xhalo determines the coordinates of the nodal halo's.
368  ! First it sets all halo coordinates by simple extrapolation,
369  ! then the symmetry planes are treated (also the unit normal of
370  ! symmetry planes are determined) and finally an exchange is
371  ! made for the internal halo's.
372  !
373  use constants
374  use blockpointers
375  use communication
376  use inputtimespectral
377  implicit none
378  !
379  ! Local variables.
380  !
381  integer(kind=intType) :: mm, i, j, k
382  integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, iiMax, jjMax
383  logical err
384  real(kind=realtype) :: length, dot
385  real(kind=realtype), dimension(3) :: v1, v2, norm
386 
387  ! Extrapolation in i-direction.
388 
389  do k = 1, kl
390  do j = 1, jl
391  x(0, j, k, 1) = two * x(1, j, k, 1) - x(2, j, k, 1)
392  x(0, j, k, 2) = two * x(1, j, k, 2) - x(2, j, k, 2)
393  x(0, j, k, 3) = two * x(1, j, k, 3) - x(2, j, k, 3)
394 
395  x(ie, j, k, 1) = two * x(il, j, k, 1) - x(nx, j, k, 1)
396  x(ie, j, k, 2) = two * x(il, j, k, 2) - x(nx, j, k, 2)
397  x(ie, j, k, 3) = two * x(il, j, k, 3) - x(nx, j, k, 3)
398  end do
399  end do
400 
401  ! Extrapolation in j-direction.
402 
403  do k = 1, kl
404  do i = 0, ie
405  x(i, 0, k, 1) = two * x(i, 1, k, 1) - x(i, 2, k, 1)
406  x(i, 0, k, 2) = two * x(i, 1, k, 2) - x(i, 2, k, 2)
407  x(i, 0, k, 3) = two * x(i, 1, k, 3) - x(i, 2, k, 3)
408 
409  x(i, je, k, 1) = two * x(i, jl, k, 1) - x(i, ny, k, 1)
410  x(i, je, k, 2) = two * x(i, jl, k, 2) - x(i, ny, k, 2)
411  x(i, je, k, 3) = two * x(i, jl, k, 3) - x(i, ny, k, 3)
412  end do
413  end do
414 
415  ! Extrapolation in k-direction.
416 
417  do j = 0, je
418  do i = 0, ie
419  x(i, j, 0, 1) = two * x(i, j, 1, 1) - x(i, j, 2, 1)
420  x(i, j, 0, 2) = two * x(i, j, 1, 2) - x(i, j, 2, 2)
421  x(i, j, 0, 3) = two * x(i, j, 1, 3) - x(i, j, 2, 3)
422 
423  x(i, j, ke, 1) = two * x(i, j, kl, 1) - x(i, j, nz, 1)
424  x(i, j, ke, 2) = two * x(i, j, kl, 2) - x(i, j, nz, 2)
425  x(i, j, ke, 3) = two * x(i, j, kl, 3) - x(i, j, nz, 3)
426  end do
427  end do
428  !
429  ! Mirror the halo coordinates adjacent to the symmetry
430  ! planes
431  !
432  ! Loop over boundary subfaces.
433 
434  loopbocos: do mm = 1, nbocos
435 
436  ! The actual correction of the coordinates only takes
437  ! place for symmetry planes.
438 
439  testsymmetry: if (bctype(mm) == symm) then
440 
441  ! Set some variables, depending on the block face on
442  ! which the subface is located.
443  norm(1) = bcdata(mm)%symNorm(1)
444  norm(2) = bcdata(mm)%symNorm(2)
445  norm(3) = bcdata(mm)%symNorm(3)
446 
447  length = sqrt(norm(1)**2 + norm(2)**2 + norm(3)**2)
448 
449  ! Compute the unit normal of the subface.
450 
451  norm(1) = norm(1) / length
452  norm(2) = norm(2) / length
453  norm(3) = norm(3) / length
454  ! See xhalo_block for comments for below:
455  testsingular: if (length > eps) then
456 
457  select case (bcfaceid(mm))
458  case (imin)
459  ibeg = jnbeg(mm); iend = jnend(mm); iimax = jl
460  jbeg = knbeg(mm); jend = knend(mm); jjmax = kl
461 
462  if (ibeg == 1) ibeg = 0
463  if (iend == iimax) iend = iimax + 1
464 
465  if (jbeg == 1) jbeg = 0
466  if (jend == jjmax) jend = jjmax + 1
467 
468  do j = jbeg, jend
469  do i = ibeg, iend
470  v1(1) = x(1, i, j, 1) - x(2, i, j, 1)
471  v1(2) = x(1, i, j, 2) - x(2, i, j, 2)
472  v1(3) = x(1, i, j, 3) - x(2, i, j, 3)
473  dot = two * (v1(1) * norm(1) + v1(2) * norm(2) &
474  + v1(3) * norm(3))
475  x(0, i, j, 1) = x(2, i, j, 1) + dot * norm(1)
476  x(0, i, j, 2) = x(2, i, j, 2) + dot * norm(2)
477  x(0, i, j, 3) = x(2, i, j, 3) + dot * norm(3)
478  end do
479  end do
480 
481  case (imax)
482  ibeg = jnbeg(mm); iend = jnend(mm); iimax = jl
483  jbeg = knbeg(mm); jend = knend(mm); jjmax = kl
484 
485  if (ibeg == 1) ibeg = 0
486  if (iend == iimax) iend = iimax + 1
487 
488  if (jbeg == 1) jbeg = 0
489  if (jend == jjmax) jend = jjmax + 1
490 
491  do j = jbeg, jend
492  do i = ibeg, iend
493  v1(1) = x(il, i, j, 1) - x(nx, i, j, 1)
494  v1(2) = x(il, i, j, 2) - x(nx, i, j, 2)
495  v1(3) = x(il, i, j, 3) - x(nx, i, j, 3)
496  dot = two * (v1(1) * norm(1) + v1(2) * norm(2) &
497  + v1(3) * norm(3))
498  x(ie, i, j, 1) = x(nx, i, j, 1) + dot * norm(1)
499  x(ie, i, j, 2) = x(nx, i, j, 2) + dot * norm(2)
500  x(ie, i, j, 3) = x(nx, i, j, 3) + dot * norm(3)
501  end do
502  end do
503 
504  case (jmin)
505  ibeg = inbeg(mm); iend = inend(mm); iimax = il
506  jbeg = knbeg(mm); jend = knend(mm); jjmax = kl
507 
508  if (ibeg == 1) ibeg = 0
509  if (iend == iimax) iend = iimax + 1
510 
511  if (jbeg == 1) jbeg = 0
512  if (jend == jjmax) jend = jjmax + 1
513 
514  do j = jbeg, jend
515  do i = ibeg, iend
516  v1(1) = x(i, 1, j, 1) - x(i, 2, j, 1)
517  v1(2) = x(i, 1, j, 2) - x(i, 2, j, 2)
518  v1(3) = x(i, 1, j, 3) - x(i, 2, j, 3)
519  dot = two * (v1(1) * norm(1) + v1(2) * norm(2) &
520  + v1(3) * norm(3))
521  x(i, 0, j, 1) = x(i, 2, j, 1) + dot * norm(1)
522  x(i, 0, j, 2) = x(i, 2, j, 2) + dot * norm(2)
523  x(i, 0, j, 3) = x(i, 2, j, 3) + dot * norm(3)
524  end do
525  end do
526 
527  case (jmax)
528  ibeg = inbeg(mm); iend = inend(mm); iimax = il
529  jbeg = knbeg(mm); jend = knend(mm); jjmax = kl
530 
531  if (ibeg == 1) ibeg = 0
532  if (iend == iimax) iend = iimax + 1
533 
534  if (jbeg == 1) jbeg = 0
535  if (jend == jjmax) jend = jjmax + 1
536 
537  do j = jbeg, jend
538  do i = ibeg, iend
539  v1(1) = x(i, jl, j, 1) - x(i, ny, j, 1)
540  v1(2) = x(i, jl, j, 2) - x(i, ny, j, 2)
541  v1(3) = x(i, jl, j, 3) - x(i, ny, j, 3)
542  dot = two * (v1(1) * norm(1) + v1(2) * norm(2) &
543  + v1(3) * norm(3))
544  x(i, je, j, 1) = x(i, ny, j, 1) + dot * norm(1)
545  x(i, je, j, 2) = x(i, ny, j, 2) + dot * norm(2)
546  x(i, je, j, 3) = x(i, ny, j, 3) + dot * norm(3)
547  end do
548  end do
549 
550  case (kmin)
551  ibeg = inbeg(mm); iend = inend(mm); iimax = il
552  jbeg = jnbeg(mm); jend = jnend(mm); jjmax = jl
553 
554  if (ibeg == 1) ibeg = 0
555  if (iend == iimax) iend = iimax + 1
556 
557  if (jbeg == 1) jbeg = 0
558  if (jend == jjmax) jend = jjmax + 1
559 
560  do j = jbeg, jend
561  do i = ibeg, iend
562  v1(1) = x(i, j, 1, 1) - x(i, j, 2, 1)
563  v1(2) = x(i, j, 1, 2) - x(i, j, 2, 2)
564  v1(3) = x(i, j, 1, 3) - x(i, j, 2, 3)
565  dot = two * (v1(1) * norm(1) + v1(2) * norm(2) &
566  + v1(3) * norm(3))
567  x(i, j, 0, 1) = x(i, j, 2, 1) + dot * norm(1)
568  x(i, j, 0, 2) = x(i, j, 2, 2) + dot * norm(2)
569  x(i, j, 0, 3) = x(i, j, 2, 3) + dot * norm(3)
570  end do
571  end do
572 
573  case (kmax)
574  ibeg = inbeg(mm); iend = inend(mm); iimax = il
575  jbeg = jnbeg(mm); jend = jnend(mm); jjmax = jl
576 
577  if (ibeg == 1) ibeg = 0
578  if (iend == iimax) iend = iimax + 1
579 
580  if (jbeg == 1) jbeg = 0
581  if (jend == jjmax) jend = jjmax + 1
582 
583  do j = jbeg, jend
584  do i = ibeg, iend
585  v1(1) = x(i, j, kl, 1) - x(i, j, nz, 1)
586  v1(2) = x(i, j, kl, 2) - x(i, j, nz, 2)
587  v1(3) = x(i, j, kl, 3) - x(i, j, nz, 3)
588  dot = two * (v1(1) * norm(1) + v1(2) * norm(2) &
589  + v1(3) * norm(3))
590  x(i, j, ke, 1) = x(i, j, nz, 1) + dot * norm(1)
591  x(i, j, ke, 2) = x(i, j, nz, 2) + dot * norm(2)
592  x(i, j, ke, 3) = x(i, j, nz, 3) + dot * norm(3)
593  end do
594  end do
595  end select
596  end if testsingular
597  end if testsymmetry
598  end do loopbocos
599  end subroutine xhalo_block
600 
601  subroutine resscale
602 
603  use constants
604  use blockpointers, only: il, jl, kl, nx, ny, nz, volref, dw
605  use flowvarrefstate, only: nwf, nt1, nt2
606  use inputiteration, only: turbresscale
607  implicit none
608 
609  ! Local Variables
610  integer(kind=intType) :: i, j, k, ii, nTurb
611  real(kind=realtype) :: ovol
612 
613  ! Divide through by the reference volume
614  nturb = nt2 - nt1 + 1
615 #ifdef TAPENADE_REVERSE
616  !$AD II-LOOP
617  do ii = 0, nx * ny * nz - 1
618  i = mod(ii, nx) + 2
619  j = mod(ii / nx, ny) + 2
620  k = ii / (nx * ny) + 2
621 #else
622  do k = 2, kl
623  do j = 2, jl
624  do i = 2, il
625 #endif
626  ovol = one / volref(i, j, k)
627  dw(i, j, k, 1:nwf) = dw(i, j, k, 1:nwf) * ovol
628  dw(i, j, k, nt1:nt2) = dw(i, j, k, nt1:nt2) * ovol * turbresscale(1:nturb)
629 #ifdef TAPENADE_REVERSE
630  end do
631 #else
632  end do
633  end do
634  end do
635 #endif
636  end subroutine resscale
637 
638  subroutine sumdwandfw
639 
640  use constants
641  use blockpointers, only: il, jl, kl, dw, fw, iblank
642  use flowvarrefstate, only: nwf
643 
644  implicit none
645 
646  ! Local Variables
647  integer(kind=intType) :: i, j, k, l
648 
649  do l = 1, nwf
650  do k = 2, kl
651  do j = 2, jl
652  do i = 2, il
653  dw(i, j, k, l) = (dw(i, j, k, l) + fw(i, j, k, l)) &
654  * max(real(iblank(i, j, k), realtype), zero)
655  end do
656  end do
657  end do
658  end do
659  end subroutine sumdwandfw
660 
661 end module adjointextra
volpym
subroutine volpym(xa, ya, za, xb, yb, zb, xc, yc, zc, xd, yd, zd, volume)
Definition: adjointExtra.F90:148
adjointextra::resscale
subroutine resscale
Definition: adjointExtra.F90:602
adjointextra::volume_block
subroutine volume_block
Definition: adjointExtra.F90:6
adjointextra::xhalo_block
subroutine xhalo_block
Definition: adjointExtra.F90:366
adjointextra::sumdwandfw
subroutine sumdwandfw
Definition: adjointExtra.F90:639
adjointextra::boundarynormals
subroutine boundarynormals
Definition: adjointExtra.F90:301
adjointextra::metric_block
subroutine metric_block
Definition: adjointExtra.F90:180
bcdata
Definition: BCData.F90:1
blockpointers::jl
integer(kind=inttype) jl
Definition: blockPointers.F90:38
blockpointers::knend
integer(kind=inttype), dimension(:), pointer knend
Definition: blockPointers.F90:56
blockpointers::inend
integer(kind=inttype), dimension(:), pointer inend
Definition: blockPointers.F90:54
blockpointers::knbeg
integer(kind=inttype), dimension(:), pointer knbeg
Definition: blockPointers.F90:56
blockpointers::righthanded
logical righthanded
Definition: blockPointers.F90:42
blockpointers::jnbeg
integer(kind=inttype), dimension(:), pointer jnbeg
Definition: blockPointers.F90:55
blockpointers::nx
integer(kind=inttype) nx
Definition: blockPointers.F90:38
blockpointers::ny
integer(kind=inttype) ny
Definition: blockPointers.F90:38
blockpointers::ie
integer(kind=inttype) ie
Definition: blockPointers.F90:39
blockpointers::iblank
integer(kind=inttype), dimension(:, :, :), pointer iblank
Definition: blockPointers.F90:71
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::nbocos
integer(kind=inttype) nbocos
Definition: blockPointers.F90:47
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::jnend
integer(kind=inttype), dimension(:), pointer jnend
Definition: blockPointers.F90:55
blockpointers::bctype
integer(kind=inttype), dimension(:), pointer bctype
Definition: blockPointers.F90:49
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::xd
real(kind=realtype), dimension(:, :, :, :), pointer xd
Definition: blockPointers.F90:185
blockpointers::fw
real(kind=realtype), dimension(:, :, :, :), pointer fw
Definition: blockPointers.F90:128
blockpointers::nz
integer(kind=inttype) nz
Definition: blockPointers.F90:38
blockpointers::je
integer(kind=inttype) je
Definition: blockPointers.F90:39
blockpointers::ke
integer(kind=inttype) ke
Definition: blockPointers.F90:39
blockpointers::x
real(kind=realtype), dimension(:, :, :, :), pointer x
Definition: blockPointers.F90:94
blockpointers::inbeg
integer(kind=inttype), dimension(:), pointer inbeg
Definition: blockPointers.F90:54
blockpointers::kl
integer(kind=inttype) kl
Definition: blockPointers.F90:38
blockpointers::il
integer(kind=inttype) il
Definition: blockPointers.F90:38
constants::zero
real(kind=realtype), parameter zero
Definition: constants.F90:71
constants::imax
integer(kind=inttype), parameter imax
Definition: constants.F90:293
constants::kmin
integer(kind=inttype), parameter kmin
Definition: constants.F90:296
constants::jmax
integer(kind=inttype), parameter jmax
Definition: constants.F90:295
constants::eps
real(kind=realtype), parameter eps
Definition: constants.F90:23
constants::eighth
real(kind=realtype), parameter eighth
Definition: constants.F90:84
constants::symm
integer(kind=inttype), parameter symm
Definition: constants.F90:258
constants::one
real(kind=realtype), parameter one
Definition: constants.F90:72
constants::half
real(kind=realtype), parameter half
Definition: constants.F90:80
constants::imin
integer(kind=inttype), parameter imin
Definition: constants.F90:292
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::jmin
integer(kind=inttype), parameter jmin
Definition: constants.F90:294
constants::sixth
real(kind=realtype), parameter sixth
Definition: constants.F90:83
flowvarrefstate::nt1
integer(kind=inttype) nt1
Definition: flowVarRefState.F90:17
flowvarrefstate::nwf
integer(kind=inttype) nwf
Definition: flowVarRefState.F90:17
flowvarrefstate::nt2
integer(kind=inttype) nt2
Definition: flowVarRefState.F90:17
inputiteration::turbresscale
real(kind=realtype), dimension(4) turbresscale
Definition: inputParam.F90:293
precision::realtype
integer, parameter realtype
Definition: precision.F90:112
Generated by doxygen 1.9.1