ADflow  v1.0
ADflow is a finite volume RANS solver tailored for gradient-based aerodynamic design optimization.
cartMesh.F90
Go to the documentation of this file.
1 module cartmesh
2 
3  use oversetdata
4  use communication
5  use utils
6  use haloexchange
8  use su_cgns
9  implicit none
10 
11 contains
12 
13  subroutine createcartmesh(level, sps)
14 
15  use constants
16  use blockpointers
17  use surfacefamilies, only: bcfamgroups
18  use su_cgns
19 #include <petsc/finclude/petsc.h>
20  use petsc
21  implicit none
22 
23  ! Input Params
24  integer(kind=intType), intent(in) :: level, sps
25 
26  ! Working params
27  integer(kind=intType) :: i, j, k, l, ii, jj, kk, iBeg, iEnd, jBeg
28  integer(kind=intType) :: jEnd, kBeg, kEnd, nn, mm, iDim, symOnFace(6)
29  integer(kind=intType) :: count, countLocal, symOnFaceLocal(6), procDims(3)
30  integer(kind=intType) :: cellDims(3), nNodes, nCells, ierr, globalInd
31  integer(kind=intType) :: nSubI, nSubJ, nSeed, iSeed, nChanged, loopIter
32  integer(kind=intType) :: iSize, jSize, kSize, nChangedLocal, stackPointer
33  integer(kind=intType), dimension(:), allocatable :: indices
34  integer(kind=intType), dimension(:, :), allocatable :: lSizes
35  integer(kind=intType), dimension(:, :), allocatable :: stack, floodSeeds
36 
37  real(kind=realtype), dimension(:), pointer :: cartpointer
38  real(kind=realtype), dimension(:, :, :), pointer :: xx
39  real(kind=realtype), dimension(:, :, :), pointer :: arrvals, changed
40  real(kind=realtype), dimension(:), allocatable :: values
41  real(kind=realtype), dimension(3) :: xminlocal, xmaxlocal, xmin, xmax, sss
42  real(kind=realtype), dimension(3) :: pt1, pt2, pt3, pt4, newpt, v1, v2
43  real(kind=realtype) :: arealocal, area, areaavg, err1, err2, h
44  real(kind=realtype) :: cooravg, scalesize, length, u, v
45 
46  dm cartarray
47  ao cartordering
48  vec cartvecglobal, cartveclocal, blankvec, blankveclocal, changedveclocal, changedvecglobal
49  is is1, is2
50  vecscatter blankscatter, blankscatterlocal
51 
52  xminlocal = large
53  xmaxlocal = -large
54  arealocal = zero
55  countlocal = 0
56  ! First we have to determine the bounding box of our surfaces.
57 
58  do nn = 1, ndom
59  call setpointers(nn, level, sps)
60  do mm = 1, nbocos
61  if (bctype(mm) == nswalladiabatic .or. &
62  bctype(mm) == nswallisothermal .or. &
63  bctype(mm) == eulerwall) then
64 
65  jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
66  ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
67 
68  select case (bcfaceid(mm))
69  case (imin)
70  xx => x(1, :, :, :)
71  case (imax)
72  xx => x(il, :, :, :)
73  case (jmin)
74  xx => x(:, 1, :, :)
75  case (jmax)
76  xx => x(:, jl, :, :)
77  case (kmin)
78  xx => x(:, :, 1, :)
79  case (kmax)
80  xx => x(:, :, kl, :)
81  end select
82 
83  do j = jbeg + 1, jend + 1
84  do i = ibeg + 1, iend + 1
85  do idim = 1, 3
86  xminlocal(idim) = min(xminlocal(idim), xx(i, j, idim))
87  xmaxlocal(idim) = max(xmaxlocal(idim), xx(i, j, idim))
88  end do
89  end do
90  end do
91 
92  ! Determine the total area of the patch:
93  do j = jbeg + 1, jend
94  do i = ibeg + 1, iend
95  v1 = xx(i, j, :) - xx(i - 1, j - 1, :)
96  v2 = xx(i - 1, j, :) - xx(i, j - 1, :)
97  ! Cross Product
98  sss(1) = (v1(2) * v2(3) - v1(3) * v2(2))
99  sss(2) = (v1(3) * v2(1) - v1(1) * v2(3))
100  sss(3) = (v1(1) * v2(2) - v1(2) * v2(1))
101 
102  arealocal = arealocal + half * (sss(1)**2 + sss(2)**2 + sss(3)**2)
103  countlocal = countlocal + 1
104  end do
105  end do
106  end if
107  end do
108  end do
109 
110  call mpi_allreduce(xminlocal, xmin, 3, adflow_real, mpi_min, adflow_comm_world, ierr)
111  call echk(ierr, __file__, __line__)
112 
113  call mpi_allreduce(xmaxlocal, xmax, 3, adflow_real, mpi_max, adflow_comm_world, ierr)
114  call echk(ierr, __file__, __line__)
115 
116  call mpi_allreduce(arealocal, area, 1, adflow_real, mpi_sum, adflow_comm_world, ierr)
117  call echk(ierr, __file__, __line__)
118 
119  call mpi_allreduce(countlocal, count, 1, adflow_integer, mpi_sum, adflow_comm_world, ierr)
120  call echk(ierr, __file__, __line__)
121 
122  ! Compute the average area:
123  areaavg = area / count
124 
125  ! Compute the master size 'h':
126  h = sqrt(areaavg)
127  h = h * 6
128  ! Before we setup the cartesian mesh, we will expand the bounds
129  ! slightly to ensure that there isn't a wall cell right at the
130  ! boundary. This is actually vastly tricker than it sounds becuase
131  ! we have to correctly account for symmetry planes. That is we
132  ! CANNOT
133  ! outside-in flood would just "go around the back" of the symmetry
134  ! plane and flood out the inside. Not good. So what we do is now
135  ! that we know the exact bounds of our surface, we loop over the
136  ! symmetry planes and flag values in an array of lenght 6
137  ! corresponding to (iLow, iHigh, jLow, jHigh, kLow, kHigh). We set
138  ! the value to 1 if the symmetry plane I'm looking at corresponds
139  ! to that value. Then we all reduce so that everyone knows which
140  ! of the 6 faces cannot be extended.
141 
142  scalesize = mynorm2(xmax - xmin)
143 
144  symonfacelocal = 0
145  do nn = 1, ndom
146  call setpointers(nn, level, sps)
147  do mm = 1, nbocos
148  if (bctype(mm) == symm) then
149 
150  jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
151  ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
152 
153  select case (bcfaceid(mm))
154  case (imin)
155  xx => x(1, :, :, :)
156  case (imax)
157  xx => x(il, :, :, :)
158  case (jmin)
159  xx => x(:, 1, :, :)
160  case (jmax)
161  xx => x(:, jl, :, :)
162  case (kmin)
163  xx => x(:, :, 1, :)
164  case (kmax)
165  xx => x(:, :, kl, :)
166  end select
167 
168  ! First we have to determine the principle coordinate
169  ! direction of the face. Do this with symNorm which is
170  ! already computed.
171 
172  ! Location, ie coordiante direction of dominate direction
173  idim = maxloc(abs(bcdata(mm)%symNorm), 1)
174 
175  ! Now determine the average value in "iDim" dimension
176  cooravg = sum(xx(ibeg + 1:iend + 1, jbeg + 1:jend + 1, idim)) &
177  / ((iend - ibeg + 1) * (jend - jbeg + 1))
178 
179  ! Check if it is sufficently close to the bounding box:
180  err1 = abs(cooravg - xmin(idim)) / scalesize
181  err2 = abs(cooravg - xmax(idim)) / scalesize
182 
183  if (err1 < 1e-8) then
184  symonfacelocal(2 * idim - 1) = 1
185  end if
186 
187  if (err2 < 1e-8) then
188  symonfacelocal(2 * idim) = 1
189  end if
190  end if
191  end do
192  end do
193 
194  ! Now we all reduce with max. Values that are not zero have
195  ! symmetry planes on them.
196 
197  call mpi_allreduce(symonfacelocal, symonface, 6, adflow_integer, mpi_max, adflow_comm_world, ierr)
198  call echk(ierr, __file__, __line__)
199 
200  ! Now we will expand xmin/xmax by 2h (just to be sure) but only
201  ! on the planes without symmetry conditions
202 
203  do idim = 1, 3
204  if (symonface(2 * idim - 1) == 0) then
205  xmin(idim) = xmin(idim) - 2 * h
206  end if
207 
208  if (symonface(2 * idim) == 0) then
209  xmax(idim) = xmax(idim) + 2 * h
210  end if
211  end do
212 
213  ! Ok. Now we have the final size of our cartesian block. The
214  ! next step is to determine how it is to be partitioned.
215  call triplefactor(nproc, procdims)
216 
217  ! And the dimensions in each direction. Remember the
218  ! distributed array is cell-centered based.
219  celldims = floor((xmax - xmin) / h) + 1
220 
221  ! Now reorder the procDims so that the largest dimension
222  ! corresponds to largest number of procs. Since we know procDims
223  ! are sorted, we can do the 6 cases:
224 
225  if (celldims(1) >= celldims(2) .and. celldims(2) >= celldims(3)) then
226  procdims = (/procdims(3), procdims(2), procdims(1)/)
227  else if (celldims(1) >= celldims(3) .and. celldims(3) >= celldims(2)) then
228  procdims = (/procdims(3), procdims(1), procdims(2)/)
229 
230  else if (celldims(2) >= celldims(1) .and. celldims(1) >= celldims(3)) then
231  procdims = (/procdims(2), procdims(3), procdims(1)/)
232  else if (celldims(2) >= celldims(3) .and. celldims(3) >= celldims(1)) then
233  procdims = (/procdims(1), procdims(3), procdims(2)/)
234 
235  else if (celldims(3) >= celldims(1) .and. celldims(1) >= celldims(2)) then
236  procdims = (/procdims(2), procdims(1), procdims(3)/)
237  else if (celldims(3) >= celldims(2) .and. celldims(2) >= celldims(1)) then
238  procdims = (/procdims(1), procdims(2), procdims(3)/)
239  end if
240 
241  ! Since we can't pass in null array objects to petsc (this is
242  ! still broken) we have to do up all the arguments ourself. In
243  ! particular, the lx, ly and lz one which are quite annonying.
244 
245  allocate (lsizes(maxval(procdims), 3))
246  do idim = 1, 3
247  do j = 1, procdims(idim)
248 
249  ii = celldims(idim) / procdims(idim)
250  if (mod(celldims(idim), procdims(idim)) > j - 1) then
251  ii = ii + 1
252  end if
253  lsizes(j, idim) = ii
254  end do
255  end do
256  if (myid == 0) then
257  print *, 'xmin:', xmin
258  print *, 'xmax:', xmax
259  print *, 'dims:', celldims
260  print *, 'H:', h
261  print *, 'procDims:', procdims
262  print *, 'I sizes:', lsizes(1:procdims(1), 1)
263  print *, 'J sizes:', lsizes(1:procdims(2), 2)
264  print *, 'K sizes:', lsizes(1:procdims(3), 3)
265  end if
266 
267  call dmdacreate3d(adflow_comm_world, dm_boundary_ghosted, dm_boundary_ghosted, &
268  dm_boundary_ghosted, dmda_stencil_star, celldims(1), celldims(2), &
269  celldims(3), procdims(1), procdims(2), procdims(3), 1, 1, &
270  lsizes(1:procdims(1), 1), lsizes(1:procdims(2), 2), lsizes(1:procdims(3), 3), &
271  cartarray, ierr)
272  call echk(ierr, __file__, __line__)
273  deallocate (lsizes)
274 
275  ! Now loop back over the surfaces. For each surface, determine
276  ! it's global index of the point that will be flagged as a cut
277  ! cell. Due to the stupid PETSC ordering crap, we have to store
278  ! all the indices as we go.
279 
280  i = 0
281  do nn = 1, ndom
282  call setpointers(nn, level, sps)
283  call getwallsize(bcfamgroups(ibcgroupwalls)%famList, nnodes, ncells, .false.)
284  i = i + nnodes
285  end do
286  i = i * 10
287  allocate (indices(i))
288  count = 0
289  do nn = 1, ndom
290  call setpointers(nn, level, sps)
291  do mm = 1, nbocos
292  if (bctype(mm) == nswalladiabatic .or. &
293  bctype(mm) == nswallisothermal .or. &
294  bctype(mm) == eulerwall) then
295 
296  jbeg = bcdata(mm)%jnBeg; jend = bcdata(mm)%jnEnd
297  ibeg = bcdata(mm)%inBeg; iend = bcdata(mm)%inEnd
298 
299  select case (bcfaceid(mm))
300  case (imin)
301  xx => x(1, :, :, :)
302  case (imax)
303  xx => x(il, :, :, :)
304  case (jmin)
305  xx => x(:, 1, :, :)
306  case (jmax)
307  xx => x(:, jl, :, :)
308  case (kmin)
309  xx => x(:, :, 1, :)
310  case (kmax)
311  xx => x(:, :, kl, :)
312  end select
313 
314  ! Loop over Raw Nodes
315  do j = jbeg, jend
316  do i = ibeg, iend
317 
318  ! Note that ii, jj, kk are zero based since we are
319  ! working with PETSc indices.
320  ii = int((xx(i + 1, j + 1, 1) - xmin(1)) / h)
321  jj = int((xx(i + 1, j + 1, 2) - xmin(2)) / h)
322  kk = int((xx(i + 1, j + 1, 3) - xmin(3)) / h)
323 
324  globalind = celldims(1) * celldims(2) * kk + celldims(1) * jj + ii
325  count = count + 1
326  indices(count) = globalind
327  end do
328  end do
329 
330  ! Loop over iEdges
331  do j = jbeg, jend
332  do i = ibeg + 1, iend
333 
334  ! Check the Length of this edge
335  pt1 = xx(i, j + 1, :)
336  pt2 = xx(i + 1, j + 1, :)
337 
338  length = mynorm2(pt1 - pt2)
339  nsubi = int(length / h)
340 
341  do k = 1, nsubi
342 
343  newpt = pt1 + dble(k) / (nsubi + 1) * (pt2 - pt1)
344  ii = int((newpt(1) - xmin(1)) / h)
345  jj = int((newpt(2) - xmin(2)) / h)
346  kk = int((newpt(3) - xmin(3)) / h)
347 
348  globalind = celldims(1) * celldims(2) * kk + celldims(1) * jj + ii
349  count = count + 1
350  indices(count) = globalind
351  end do
352  end do
353  end do
354 
355  ! Loop over jEdges
356  do j = jbeg + 1, jend
357  do i = ibeg, iend
358 
359  ! Check the Length of this edge
360  pt1 = xx(i + 1, j, :)
361  pt2 = xx(i + 1, j + 1, :)
362 
363  length = mynorm2(pt1 - pt2)
364  nsubj = int(length / h)
365  do k = 1, nsubj
366 
367  newpt = pt1 + dble(k) / (nsubj + 1) * (pt2 - pt1)
368  ii = int((newpt(1) - xmin(1)) / h)
369  jj = int((newpt(2) - xmin(2)) / h)
370  kk = int((newpt(3) - xmin(3)) / h)
371 
372  globalind = celldims(1) * celldims(2) * kk + celldims(1) * jj + ii
373  count = count + 1
374  indices(count) = globalind
375  end do
376  end do
377  end do
378 
379  ! Loop over faces
380  do j = jbeg + 1, jend
381  do i = ibeg + 1, iend
382 
383  ! Extract the 4 pts. CCW ordering
384  pt1 = xx(i, j, :)
385  pt2 = xx(i + 1, j, :)
386 
387  pt3 = xx(i + 1, j + 1, :)
388  pt4 = xx(i, j + 1, :)
389 
390  ! Sub pts in I
391  length = mynorm2(pt2 - pt1)
392  nsubi = int(length / h)
393 
394  length = mynorm2(pt3 - pt4)
395  nsubi = max(nsubi, int(length / h))
396 
397  ! Sub Pts in J
398  length = mynorm2(pt4 - pt1)
399  nsubj = int(length / h)
400 
401  length = mynorm2(pt3 - pt2)
402  nsubj = max(nsubj, int(length / h))
403 
404  do l = 1, nsubj
405  do k = 1, nsubi
406  u = dble(k) / (nsubi + 1)
407  v = dble(l) / (nsubj + 1)
408 
409  newpt = (one - u) * (one - v) * pt1 + u * (one - v) * pt2 + &
410  u * v * pt3 + (one - u) * v * pt4
411 
412  ii = int((newpt(1) - xmin(1)) / h)
413  jj = int((newpt(2) - xmin(2)) / h)
414  kk = int((newpt(3) - xmin(3)) / h)
415 
416  globalind = celldims(1) * celldims(2) * kk + celldims(1) * jj + ii
417  count = count + 1
418  indices(count) = globalind
419  end do
420  end do
421  end do
422  end do
423  end if
424  end do
425  end do
426 
427  call dmcreateglobalvector(cartarray, cartvecglobal, ierr)
428  call echk(ierr, __file__, __line__)
429 
430  ! Initialize all values to 1 ("Compute")
431  call vecset(cartvecglobal, one, ierr)
432  call echk(ierr, __file__, __line__)
433 
434  call dmdagetao(cartarray, cartordering, ierr)
435  call echk(ierr, __file__, __line__)
436 
437  ! Convert the indices from application ordering to petsc
438  call aoapplicationtopetsc(cartordering, size(indices), indices, ierr)
439  call echk(ierr, __file__, __line__)
440 
441  ! Now we set all the values, with a simple single vecSet call
442  allocate (values(count))
443  values = -three ! -3 is flood seed. Use the same notation here.
444 
445  call vecsetvalues(cartvecglobal, count, indices, values, insert_values, ierr)
446  call echk(ierr, __file__, __line__)
447  deallocate (values, indices)
448 
449  ! Don't forget to assemble
450  call vecassemblybegin(cartvecglobal, ierr)
451  call echk(ierr, __file__, __line__)
452 
453  call vecassemblyend(cartvecglobal, ierr)
454  call echk(ierr, __file__, __line__)
455 
456  ! These are "get" vecs. We have to restore them later.
457  call dmgetlocalvector(cartarray, cartveclocal, ierr)
458  call echk(ierr, __file__, __line__)
459 
460  call dmgetlocalvector(cartarray, changedveclocal, ierr)
461  call echk(ierr, __file__, __line__)
462 
463  call vecset(changedveclocal, zero, ierr)
464  call echk(ierr, __file__, __line__)
465 
466  call dmgetglobalvector(cartarray, changedvecglobal, ierr)
467  call echk(ierr, __file__, __line__)
468 
469  ! Now the next step is to perform the flooding from the outside in.
470  call dmglobaltolocalbegin(cartarray, cartvecglobal, insert_values, cartveclocal, ierr)
471  call echk(ierr, __file__, __line__)
472 
473  call dmglobaltolocalend(cartarray, cartvecglobal, insert_values, cartveclocal, ierr)
474  call echk(ierr, __file__, __line__)
475 
476  ! Determine the bounds of the arrays we will be getting back.
477  call dmdavecgetarrayf90(cartarray, cartveclocal, arrvals, ierr)
478  call echk(ierr, __file__, __line__)
479 
480  ibeg = lbound(arrvals, 1)
481  jbeg = lbound(arrvals, 2)
482  kbeg = lbound(arrvals, 3)
483 
484  iend = ubound(arrvals, 1)
485  jend = ubound(arrvals, 2)
486  kend = ubound(arrvals, 3)
487 
488  isize = iend - ibeg + 1
489  jsize = jend - jbeg + 1
490  ksize = kend - kbeg + 1
491 
492  call dmdavecrestorearrayf90(cartarray, cartveclocal, arrvals, ierr)
493  call echk(ierr, __file__, __line__)
494 
495  loopiter = 1
496  allocate (stack(3, 6 * isize * jsize * ksize + 1))
497  allocate (floodseeds(3, 2 * isize * jsize + 2 * isize * ksize + 2 * jsize * ksize))
498 
499  parallelsyncloop: do
500 
501  call dmdavecgetarrayf90(cartarray, cartveclocal, arrvals, ierr)
502  call echk(ierr, __file__, __line__)
503 
504  call dmdavecgetarrayf90(cartarray, changedveclocal, changed, ierr)
505  call echk(ierr, __file__, __line__)
506 
507  ! Keep track of the total number of fringes we've modified
508  nchangedlocal = 0
509 
510  ! Allocate space for our queue (stack). It needs to be 6*nx*ny*nz + 1:
511  ! 6 for each of the 6 coordinate directions plus our extra
512  ! seed. It should never come close to this unless the entire
513  ! block will be blanked.
514 
515  nseed = 0
516  if (loopiter == 1) then
517 
518  if (myid == 0) then
519 
520  ! Set the single seed on the bottom corner of the root proc
521 
522  call addseed(0, 0, 0)
523 
524  end if
525 
526  else
527 
528  ! On the second and subsequent passes, check each 1st
529  ! non-corner halos in the 6 faces to see if we received
530  ! "changed" info from neighbour proc. This will allow us to
531  ! continue the flooding on this processor/block. Note that
532  ! even in a single processor case, the halo exchange in
533  ! necessary to communicate between two local blocks
534 
535  ! iMin/iMax
536  do k = kbeg + 1, kend - 1
537  do j = jbeg + 1, jend - 1
538  if (int(changed(ibeg, j, k)) == 1) then
539  call addseed(ibeg + 1, j, k)
540  end if
541  if (int(changed(iend, j, k)) == 1) then
542  call addseed(iend - 1, j, k)
543  end if
544  end do
545  end do
546 
547  ! jMin/jMax
548  do k = kbeg + 1, kend - 1
549  do i = ibeg + 1, iend - 1
550  if (int(changed(i, jbeg, k)) == 1) then
551  call addseed(i, jbeg + 1, k)
552  end if
553  if (int(changed(i, jend, k)) == 1) then
554  call addseed(i, jend - 1, k)
555  end if
556  end do
557  end do
558 
559  ! kMin:
560  do j = jbeg + 1, jend - 1
561  do i = ibeg + 1, iend - 1
562  if (int(changed(i, j, kbeg)) == 1) then
563  call addseed(i, j, kbeg + 1)
564  end if
565  if (int(changed(i, j, kend)) == 1) then
566  call addseed(i, j, kend - 1)
567  end if
568  end do
569  end do
570  end if
571 
572  ! Loop over our seeds we currently have
573  do iseed = 1, nseed
574 
575  ! Put the particular seed in the first slot of the stack
576  stack(:, 1) = floodseeds(:, iseed)
577 
578  ! Reset the stack pointer length back to just the one seed
579  ! we have
580  stackpointer = 1
581 
582  ! Start the flooding (stacked based, not recursive)
583  do while (stackpointer > 0)
584 
585  ! 'Pop' the current point off the stack
586  i = stack(1, stackpointer)
587  j = stack(2, stackpointer)
588  k = stack(3, stackpointer)
589  stackpointer = stackpointer - 1
590 
591  if (int(arrvals(i, j, k)) == 1) then
592 
593  ! Flag the cell (using changed) as being changed
594  changed(i, j, k) = one
595 
596  ! Keep track of the total number we've changed. For
597  ! reporting purposes...only count the ones that are
598  ! on actual compute cells:
599  if (onblock(i, j, k)) then
600  nchangedlocal = nchangedlocal + 1
601  end if
602 
603  ! Set the value as flooded. Use 0 as per usual.
604  arrvals(i, j, k) = zero
605 
606  ! Now add the six nearest neighbours to the stack
607  ! provided they are in the owned cell range:
608 
609  if (i - 1 >= ibeg) then
610  stackpointer = stackpointer + 1
611  stack(:, stackpointer) = (/i - 1, j, k/)
612  end if
613 
614  if (i + 1 <= iend) then
615  stackpointer = stackpointer + 1
616  stack(:, stackpointer) = (/i + 1, j, k/)
617  end if
618 
619  if (j - 1 >= jbeg) then
620  stackpointer = stackpointer + 1
621  stack(:, stackpointer) = (/i, j - 1, k/)
622  end if
623 
624  if (j + 1 <= jend) then
625  stackpointer = stackpointer + 1
626  stack(:, stackpointer) = (/i, j + 1, k/)
627  end if
628 
629  if (k - 1 >= kbeg) then
630  stackpointer = stackpointer + 1
631  stack(:, stackpointer) = (/i, j, k - 1/)
632  end if
633 
634  if (k + 1 <= kend) then
635  stackpointer = stackpointer + 1
636  stack(:, stackpointer) = (/i, j, k + 1/)
637  end if
638  end if
639  end do
640  end do
641 
642  call dmdavecrestorearrayf90(cartarray, cartveclocal, arrvals, ierr)
643  call echk(ierr, __file__, __line__)
644 
645  call dmdavecrestorearrayf90(cartarray, changedveclocal, changed, ierr)
646  call echk(ierr, __file__, __line__)
647 
648  ! Exchange "changed"
649  call dmlocaltoglobalbegin(cartarray, changedveclocal, insert_values, changedvecglobal, ierr)
650  call echk(ierr, __file__, __line__)
651 
652  call dmlocaltoglobalend(cartarray, changedveclocal, insert_values, changedvecglobal, ierr)
653  call echk(ierr, __file__, __line__)
654 
655  call dmglobaltolocalbegin(cartarray, changedvecglobal, insert_values, changedveclocal, ierr)
656  call echk(ierr, __file__, __line__)
657 
658  call dmglobaltolocalend(cartarray, changedvecglobal, insert_values, changedveclocal, ierr)
659  call echk(ierr, __file__, __line__)
660 
661  ! Determine if cells got changd. If so do another loop.
662  call mpi_allreduce(nchangedlocal, nchanged, 1, adflow_integer, mpi_sum, &
663  adflow_comm_world, ierr)
664  call echk(ierr, __file__, __line__)
665  if (myid == 0) then
666  print *, 'Cart Flood Iteration:', loopiter, 'Blanked ', nchanged, 'Interior Cells.'
667  end if
668 
669  if (nchanged == 0) then
670  exit parallelsyncloop
671  end if
672 
673  loopiter = loopiter + 1
674  end do parallelsyncloop
675 
676  deallocate (stack, floodseeds)
677 
678  ! Now that we have flooded everything, any cells left over must
679  ! be *inside. Do one last pass through and flip those.
680 
681  call dmdavecgetarrayf90(cartarray, cartveclocal, arrvals, ierr)
682  call echk(ierr, __file__, __line__)
683 
684  do k = kbeg + 1, kend - 1
685  do j = jbeg + 1, jend - 1
686  do i = ibeg + 1, iend - 1
687  if (int(arrvals(i, j, k)) == 1) then
688  arrvals(i, j, k) = -three
689  end if
690  end do
691  end do
692  end do
693 
694  call dmdavecrestorearrayf90(cartarray, cartveclocal, arrvals, ierr)
695  call echk(ierr, __file__, __line__)
696 
697  ! Restore the vectors obtained with "get"
698  call dmrestorelocalvector(cartarray, cartveclocal, ierr)
699  call echk(ierr, __file__, __line__)
700 
701  call dmrestorelocalvector(cartarray, changedveclocal, ierr)
702  call echk(ierr, __file__, __line__)
703 
704  call dmrestoreglobalvector(cartarray, changedvecglobal, ierr)
705  call echk(ierr, __file__, __line__)
706 
707  ! Now update the final global variables in cartArray.
708  call dmlocaltoglobalbegin(cartarray, cartveclocal, insert_values, cartvecglobal, ierr)
709  call echk(ierr, __file__, __line__)
710 
711  call dmlocaltoglobalend(cartarray, cartveclocal, insert_values, cartvecglobal, ierr)
712  call echk(ierr, __file__, __line__)
713 
714  ! Now that we are done with the flooding and any local
715  ! communication, we can create a global vector that is in the
716  ! ordering that we actually want.
717 
718  i = celldims(1) * celldims(2) * celldims(3)
719  call veccreatempi(adflow_comm_world, petsc_decide, i, blankvec, ierr)
720  call echk(ierr, __file__, __line__)
721 
722  ! Create the index for the real global Vector
723  call vecgetownershiprange(blankvec, i, j, ierr)
724  call echk(ierr, __file__, __line__)
725  call iscreatestride(adflow_comm_world, j - i, i, 1, is2, ierr)
726  call echk(ierr, __file__, __line__)
727 
728  call isduplicate(is2, is1, ierr)
729  call aoapplicationtopetscis(cartordering, is1, ierr)
730  call echk(ierr, __file__, __line__)
731 
732  call vecscattercreate(cartvecglobal, is1, blankvec, is2, blankscatter, ierr)
733  call echk(ierr, __file__, __line__)
734 
735  call vecscatterbegin(blankscatter, cartvecglobal, blankvec, &
736  insert_values, scatter_forward, ierr)
737  call echk(ierr, __file__, __line__)
738 
739  call vecscatterend(blankscatter, cartvecglobal, blankvec, &
740  insert_values, scatter_forward, ierr)
741  call echk(ierr, __file__, __line__)
742 
743  ! We are now done with everything petsc related except for
744  ! blankVec. That's all we need now.
745  call vecscatterdestroy(blankscatter, ierr)
746  call echk(ierr, __file__, __line__)
747 
748  call isdestroy(is1, ierr)
749  call echk(ierr, __file__, __line__)
750 
751  call isdestroy(is2, ierr)
752  call echk(ierr, __file__, __line__)
753 
754  call vecdestroy(cartvecglobal, ierr)
755  call echk(ierr, __file__, __line__)
756 
757  call dmdestroy(cartarray, ierr)
758  call echk(ierr, __file__, __line__)
759 
760  ! Write our mesh
761  call writecartmesh(blankvec, celldims, xmin, h)
762 
763  ! Now what we have to do is to distribute parts of the cart mesh
764  ! to the processors that need it. For now, just do a create to all.
765 
766  call vecscattercreatetoall(blankvec, blankscatterlocal, blankveclocal, ierr)
767  call echk(ierr, __file__, __line__)
768 
769  call vecscatterbegin(blankscatterlocal, blankvec, blankveclocal, insert_values, scatter_forward, ierr)
770  call echk(ierr, __file__, __line__)
771 
772  call vecscatterend(blankscatterlocal, blankvec, blankveclocal, insert_values, scatter_forward, ierr)
773  call echk(ierr, __file__, __line__)
774 
775  call vecscatterdestroy(blankscatterlocal, ierr)
776  call echk(ierr, __file__, __line__)
777 
778  call vecgetarrayf90(blankveclocal, cartpointer, ierr)
779  call echk(ierr, __file__, __line__)
780 
781  do nn = 1, ndom
782  call setpointers(nn, level, sps)
783  do k = 2, kl
784  do j = 2, jl
785  do i = 2, il
786  ii = ii + 1
787  do idim = 1, 3
788  pt1(idim) = eighth * ( &
789  x(i - 1, j - 1, k - 1, idim) + &
790  x(i, j - 1, k - 1, idim) + &
791  x(i - 1, j, k - 1, idim) + &
792  x(i, j, k - 1, idim) + &
793  x(i - 1, j - 1, k, idim) + &
794  x(i, j - 1, k, idim) + &
795  x(i - 1, j, k, idim) + &
796  x(i, j, k, idim))
797  end do
798 
799  ! Now Simply check if the cell this point is in has a value of -3
800  ii = int((pt1(1) - xmin(1)) / h)
801  jj = int((pt1(2) - xmin(2)) / h)
802  kk = int((pt1(3) - xmin(3)) / h)
803 
804  ! Clip the bounds to the actual ranges:
805  ii = min(max(0, ii), celldims(1) - 1)
806  jj = min(max(0, jj), celldims(2) - 1)
807  kk = min(max(0, kk), celldims(3) - 1)
808 
809  ! GlobalInd is 0 based
810  globalind = celldims(1) * celldims(2) * kk + celldims(1) * jj + ii
811 
812  ! Only blank if we are more than sqrt(3)*h from our own wall:
813  if (xseed(i, j, k, 1) < large) then
814  ! We have a wall: See how far we are away form it
815  length = mynorm2(pt1 - xseed(i, j, k, :))
816  if (length > 1.73205080 * h) then
817  if (int(cartpointer(globalind + 1)) == -3) then
818  iblank(i, j, k) = -3
819  end if
820  end if
821  else
822  ! No wall...no wall check.
823  if (int(cartpointer(globalind + 1)) == -3) then
824  iblank(i, j, k) = -3
825  end if
826  end if
827  end do
828  end do
829  end do
830  end do
831 
832  call vecrestorearrayf90(blankveclocal, cartpointer, ierr)
833  call echk(ierr, __file__, __line__)
834 
835  ! Clean up the remaining PETScmemory
836  call vecdestroy(blankvec, ierr)
837  call echk(ierr, __file__, __line__)
838 
839  call vecdestroy(blankveclocal, ierr)
840  call echk(ierr, __file__, __line__)
841 
842  ! Update the iblank info.
843  domainloop: do nn = 1, ndom
844  flowdoms(nn, level, sps)%intCommVars(1)%var => &
845  flowdoms(nn, level, sps)%iblank(:, :, :)
846  end do domainloop
847 
848  ! Run the generic integer exchange
850 
851  contains
852  ! Simple routine to make code easier to read above
853  subroutine addseed(i, j, k)
854  use constants
855  implicit none
856  integer(kind=intType), intent(in) :: i, j, k
857  nseed = nseed + 1
858  floodseeds(:, nseed) = (/i, j, k/)
859  end subroutine addseed
860 
861  function onblock(i, j, k)
862 
863  use constants
864  implicit none
865 
866  integer(kind=intType), intent(in) :: i, j, k
867  logical :: onblock
868 
869  if (i >= 0 .and. i <= iend - 1 .and. j >= 0 .and. j <= jend - 1 .and. k >= 0 .and. k <= kend - 1) then
870  onblock = .true.
871  else
872  onblock = .false.
873  end if
874  end function onblock
875  end subroutine createcartmesh
876 
877  subroutine writecartmesh(blankVec, cellDims, xMin, h)
878 
879 #include <petsc/finclude/petsc.h>
880  use petsc
881  implicit none
882 
883  ! Input
884  integer(kind=intType), intent(in), dimension(3) :: cellDims
885  real(kind=realtype), intent(in), dimension(3) :: xmin
886  real(kind=realtype), intent(in) :: h
887  vec blankvec
888  ! Working
889  integer(kind=intType) :: ierr
890 
891  ! CGNS
892  character(len=32) :: coorNames(3)
893  integer(kind=intType) :: base, zoneID, coordID, cg, zone, iField, iSol, i, j, k, iDim
894  real(kind=realtype), dimension(:, :, :, :), allocatable :: xtmp
895  real(kind=realtype), dimension(:), pointer :: cartpointer
896 
897  vec blankveclocal
898  is is1, is2
899  vecscatter blankscatterlocal
900 
901  coornames(1) = "CoordinateX"
902  coornames(2) = "CoordinateY"
903  coornames(3) = "CoordinateZ"
904 
905  call vecscattercreatetozero(blankvec, blankscatterlocal, blankveclocal, ierr)
906  call echk(ierr, __file__, __line__)
907 
908  call vecscatterbegin(blankscatterlocal, blankvec, blankveclocal, &
909  insert_values, scatter_forward, ierr)
910  call echk(ierr, __file__, __line__)
911 
912  call vecscatterend(blankscatterlocal, blankvec, blankveclocal, &
913  insert_values, scatter_forward, ierr)
914  call echk(ierr, __file__, __line__)
915 
916  if (myid == 0) then
917  ! Open the CGNS File
918  call cg_open_f("cartblock.cgns", mode_write, cg, ierr)
919  base = 1
920  call cg_base_write_f(cg, "Base#1", 3, 3, base, ierr)
921 
922  call cg_zone_write_f(cg, base, "cartblock", int((/celldims(1) + 1, celldims(2) + 1, celldims(3) + 1, &
923  celldims(1), celldims(2), celldims(3), 0, 0, 0/), &
924  cgsize_t), structured, zoneid, ierr)
925 
926  allocate (xtmp(celldims(1) + 1, celldims(2) + 1, celldims(3) + 1, 3))
927 
928  do k = 1, celldims(3) + 1
929  do j = 1, celldims(2) + 1
930  do i = 1, celldims(1) + 1
931  xtmp(i, j, k, 1) = xmin(1) + (i - 1) * h
932  end do
933  end do
934  end do
935 
936  do k = 1, celldims(3) + 1
937  do j = 1, celldims(2) + 1
938  do i = 1, celldims(1) + 1
939  xtmp(i, j, k, 2) = xmin(2) + (j - 1) * h
940  end do
941  end do
942  end do
943 
944  do k = 1, celldims(3) + 1
945  do j = 1, celldims(2) + 1
946  do i = 1, celldims(1) + 1
947  xtmp(i, j, k, 3) = xmin(3) + (k - 1) * h
948  end do
949  end do
950  end do
951 
952  do idim = 1, 3
953  call cg_coord_write_f(cg, base, zoneid, realdouble, coornames(idim), &
954  xtmp(:, :, :, idim), coordid, ierr)
955  end do
956  deallocate (xtmp)
957 
958  call cg_sol_write_f(cg, base, zoneid, "flowSolution", cellcenter, isol, ierr)
959  call vecgetarrayf90(blankveclocal, cartpointer, ierr)
960  call echk(ierr, __file__, __line__)
961 
962  call cg_field_write_f(cg, base, zoneid, isol, realdouble, "iBlank", &
963  cartpointer, ifield, ierr)
964 
965  call vecrestorearrayf90(blankveclocal, cartpointer, ierr)
966  call echk(ierr, __file__, __line__)
967  call cg_close_f(cg, ierr)
968  end if
969 
970  call vecdestroy(blankveclocal, ierr)
971  call echk(ierr, __file__, __line__)
972 
973  call vecscatterdestroy(blankscatterlocal, ierr)
974  call echk(ierr, __file__, __line__)
975 
976  end subroutine writecartmesh
977 
978  subroutine triplefactor(N, s)
979  use constants
980  use sorting, only: qsortintegers
981  implicit none
982  ! Input/Output
983  integer(kind=intType), intent(in) :: N
984  integer(kind=intType), intent(out) :: s(3)
985 
986  ! Working
987  integer(kind=intType) :: a, b, a1, b1, a2, b2, s1(3), s2(3)
988 
989  ! Determine a set of triple factors for integer N
990 
991  call largefactor(n, a, b)
992  call largefactor(b, a1, a2)
993  call largefactor(a, b1, b2)
994 
995  ! Our options are a, a1 and a2 OR b, b1, and b2
996  s1 = (/a, a1, a2/)
997  s2 = (/b, b1, b2/)
998 
999  ! Sort them
1000  call qsortintegers(s1, 3)
1001  call qsortintegers(s2, 3)
1002 
1003  ! And take the set that has the largest, smallest value.
1004  if (s1(1) > s2(1)) then
1005  s = s1
1006  else
1007  s = s2
1008  end if
1009  end subroutine triplefactor
1010 
1011  subroutine largefactor(N, f1, f2)
1012 
1013  implicit none
1014  integer(kind=intType) :: N, f1, f2, i, j, s
1015  ! Return the two factors that are closest to the sqrt(N)
1016 
1017  s = int(sqrt(dble(n)))
1018  do j = s, 1, -1
1019  if (mod(n, j) == 0) then
1020  f1 = n / j
1021  f2 = j
1022  exit
1023  end if
1024  end do
1025  end subroutine largefactor
1026 
1027 end module cartmesh
subroutine addseed(i, j, k)
Definition: cartMesh.F90:854
logical function onblock(i, j, k)
Definition: fortranPC.F90:450
Definition: BCData.F90:1
integer(kind=inttype) jl
integer(kind=inttype), dimension(:, :, :), pointer iblank
integer(kind=inttype), dimension(:), pointer bcfaceid
integer(kind=inttype) nbocos
integer(kind=inttype), dimension(:), pointer bctype
real(kind=realtype), dimension(:, :, :, :), pointer x
real(kind=realtype), dimension(:, :, :, :), pointer xseed
integer(kind=inttype) kl
integer(kind=inttype) il
subroutine triplefactor(N, s)
Definition: cartMesh.F90:979
subroutine largefactor(N, f1, f2)
Definition: cartMesh.F90:1012
subroutine writecartmesh(blankVec, cellDims, xMin, h)
Definition: cartMesh.F90:878
subroutine createcartmesh(level, sps)
Definition: cartMesh.F90:14
type(internalcommtype), dimension(:), allocatable internalcell_2nd
type(commtype), dimension(:), allocatable commpatterncell_2nd
integer adflow_comm_world
integer(kind=inttype), parameter eulerwall
Definition: constants.F90:262
real(kind=realtype), parameter zero
Definition: constants.F90:71
real(kind=realtype), parameter three
Definition: constants.F90:74
integer(kind=inttype), parameter imax
Definition: constants.F90:293
integer(kind=inttype), parameter kmin
Definition: constants.F90:296
integer(kind=inttype), parameter jmax
Definition: constants.F90:295
real(kind=realtype), parameter eighth
Definition: constants.F90:84
integer(kind=inttype), parameter nswalladiabatic
Definition: constants.F90:260
integer(kind=inttype), parameter symm
Definition: constants.F90:258
real(kind=realtype), parameter one
Definition: constants.F90:72
real(kind=realtype), parameter half
Definition: constants.F90:80
integer(kind=inttype), parameter imin
Definition: constants.F90:292
integer(kind=inttype), parameter nswallisothermal
Definition: constants.F90:261
integer(kind=inttype), parameter ibcgroupwalls
Definition: constants.F90:309
real(kind=realtype), parameter large
Definition: constants.F90:24
integer(kind=inttype), parameter kmax
Definition: constants.F90:297
integer(kind=inttype), parameter jmin
Definition: constants.F90:294
subroutine whalo1to1intgeneric(nVar, level, sps, commPattern, internal)
subroutine getwallsize(famList, nNodes, nCells, dualMesh)
subroutine qsortintegers(arr, nn)
Definition: sorting.F90:101
type(bcgrouptype), dimension(nfamexchange) bcfamgroups
Definition: utils.F90:1
real(kind=realtype) function mynorm2(x)
Definition: utils.F90:1697
subroutine echk(errorcode, file, line)
Definition: utils.F90:5722
subroutine setpointers(nn, mm, ll)
Definition: utils.F90:3237