ADflow  v1.0
ADflow is a finite volume RANS solver tailored for gradient-based aerodynamic design optimization.
loadBalance.F90
Go to the documentation of this file.
1 module loadbalance
2 
3 contains
4 
5  subroutine loadbalancegrid
6  !
7  ! loadBalance determines the mapping of the blocks onto the
8  ! processors. If the user allows so blocks my be split to obtain
9  ! a better load balance.
10  !
11  use constants
12  use block, only: flowdoms, ndom
13  use cgnsgrid, only: cgnsdoms, cgnsndom
14  use communication, only: adflow_comm_world, nproc, myid
15  use inputmotion, only: gridmotionspecified
16  use inputparallel, only: partitionlikenproc
17  use inputphysics, only: equationmode
18  use inputtimespectral, only: ntimeintervalsspectral
19  use iteration, only: deforming_grid
20  use partitionmod, only: subblocksofcgnstype, blocks, part, nblocks, &
21  sortrangessplitinfo, qsortsubblocksofcgnstype
22  use utils, only: terminate
23  implicit none
24  !
25  ! Local variables.
26  !
27  integer :: ierr
28 
29  integer(kind=intType) :: i, j, k, nn, mm, ii, jj, kk
30  integer(kind=intType) :: nViscBocos
31 
32  integer(kind=intType), dimension(0:nProc - 1) :: nBlockPerProc
33 
34  integer(kind=intType), dimension(:), allocatable :: oldSubfaceID
35 
36  type(subblocksofcgnstype), dimension(:), allocatable :: &
37  subblocksOfCGNS
38 
39  ! Determine the block distribution over the processors.
40 
41  if (partitionlikenproc > nproc) then
42  nproc = partitionlikenproc
43  end if
44 
45  call blockdistribution
46 
47  ! Restore the size of the comm
48  call mpi_comm_size(adflow_comm_world, nproc, ierr)
49 
50  ! We ned to modify what comes back if we are using the
51  ! partitionLikenProc option:
52  if (partitionlikenproc > nproc) then
53  do i = 1, nblocks
54  part(i) = mod(part(i), nproc)
55  end do
56  end if
57 
58  !
59  ! Determine the local block info.
60  !
61  ! Initialize nBlockPerProc to 0.
62 
63  nblockperproc = 0
64 
65  ! Determine the number of blocks the current processor will store
66  ! and the local block number for every block.
67 
68  ndom = 0
69  do i = 1, nblocks
70  if (part(i) == myid) ndom = ndom + 1
71 
72  nblockperproc(part(i)) = nblockperproc(part(i)) + 1
73  blocks(i)%blockID = nblockperproc(part(i))
74  end do
75 
76  ! Allocate the memory for flowDoms and initialize its pointers
77  ! to null pointers.
78 
79  call initflowdoms
80 
81  ! Repeat the loop, but now store the info of the blocks
82  ! in flowDoms. Store the number of time intervals for the spectral
83  ! method a bit easier in mm. Note that this number is 1 for the
84  ! steady and unsteady modes.
85 
86  nn = 0
87  mm = ntimeintervalsspectral
88  domains: do i = 1, nblocks
89  myblock: if (part(i) == myid) then
90 
91  ! Update the counter nn.
92 
93  nn = nn + 1
94 
95  ! Copy the dimensions of the block.
96 
97  flowdoms(nn, 1, 1:mm)%nx = blocks(i)%nx
98  flowdoms(nn, 1, 1:mm)%ny = blocks(i)%ny
99  flowdoms(nn, 1, 1:mm)%nz = blocks(i)%nz
100 
101  flowdoms(nn, 1, 1:mm)%il = blocks(i)%il
102  flowdoms(nn, 1, 1:mm)%jl = blocks(i)%jl
103  flowdoms(nn, 1, 1:mm)%kl = blocks(i)%kl
104 
105  ! The number of single halo quantities.
106 
107  flowdoms(nn, 1, 1:mm)%ie = blocks(i)%il + 1
108  flowdoms(nn, 1, 1:mm)%je = blocks(i)%jl + 1
109  flowdoms(nn, 1, 1:mm)%ke = blocks(i)%kl + 1
110 
111  ! The number of double halo quantities.
112 
113  flowdoms(nn, 1, 1:mm)%ib = blocks(i)%il + 2
114  flowdoms(nn, 1, 1:mm)%jb = blocks(i)%jl + 2
115  flowdoms(nn, 1, 1:mm)%kb = blocks(i)%kl + 2
116 
117  ! Relation to the original cgns grid.
118 
119  flowdoms(nn, 1, 1:mm)%cgnsBlockID = blocks(i)%cgnsBlockID
120 
121  flowdoms(nn, 1, 1:mm)%iBegOr = blocks(i)%iBegOr
122  flowdoms(nn, 1, 1:mm)%jBegOr = blocks(i)%jBegOr
123  flowdoms(nn, 1, 1:mm)%kBegOr = blocks(i)%kBegOr
124 
125  flowdoms(nn, 1, 1:mm)%iEndOr = blocks(i)%iEndOr
126  flowdoms(nn, 1, 1:mm)%jEndOr = blocks(i)%jEndOr
127  flowdoms(nn, 1, 1:mm)%kEndOr = blocks(i)%kEndOr
128 
129  ! Determine whether or not the block is moving.
130  ! First initialize it to gridMotionSpecified. This is
131  ! .true. if a rigid body motion was specified for the
132  ! entire grid; otherwise it is .false.
133 
134  flowdoms(nn, 1, 1:mm)%blockIsMoving = gridmotionspecified
135 
136  ! Check whether the corresponding cgns block is moving.
137  ! Although it is possible that boundaries of a block rotate
138  ! differently than the block itself, this should not be
139  ! taken into account here; that's a matter of BC's.
140  ! Here only the internal block structure is looked at.
141 
142  k = flowdoms(nn, 1, 1)%cgnsBlockID
143  if (cgnsdoms(k)%rotatingFrameSpecified) &
144  flowdoms(nn, 1, 1:mm)%blockIsMoving = .true.
145 
146  ! For an unsteady computation on a deforming mesh
147  ! blockIsMoving is always .true. Note that the time spectral
148  ! method is also an unsteady computation.
149 
150  if (deforming_grid .and. &
151  (equationmode == unsteady .or. &
152  equationmode == timespectral)) &
153  flowdoms(nn, 1, 1:mm)%blockIsMoving = .true.
154 
155  ! Set addGridVelocities to blockIsMoving. This could be
156  ! overwritten later when the code is running in python mode.
157 
158  flowdoms(nn, 1, 1:mm)%addGridVelocities = &
159  flowdoms(nn, 1, 1:mm)%blockIsMoving
160 
161  ! Set the number of subfaces and allocate the memory for the
162  ! subface info. Note that this memory is only allocated for
163  ! the first spectral time value; the other ones are identical.
164 
165  flowdoms(nn, 1, 1:mm)%nBocos = blocks(i)%nBocos
166  flowdoms(nn, 1, 1:mm)%n1to1 = blocks(i)%n1to1
167  flowdoms(nn, 1, 1:mm)%nSubface = blocks(i)%nSubface
168  j = blocks(i)%nSubface
169 
170  allocate (flowdoms(nn, 1, 1)%BCType(j), &
171  flowdoms(nn, 1, 1)%BCFaceID(j), &
172  flowdoms(nn, 1, 1)%cgnsSubface(j), &
173  flowdoms(nn, 1, 1)%inBeg(j), &
174  flowdoms(nn, 1, 1)%jnBeg(j), &
175  flowdoms(nn, 1, 1)%knBeg(j), &
176  flowdoms(nn, 1, 1)%inEnd(j), &
177  flowdoms(nn, 1, 1)%jnEnd(j), &
178  flowdoms(nn, 1, 1)%knEnd(j), &
179  flowdoms(nn, 1, 1)%dinBeg(j), &
180  flowdoms(nn, 1, 1)%djnBeg(j), &
181  flowdoms(nn, 1, 1)%dknBeg(j), &
182  flowdoms(nn, 1, 1)%dinEnd(j), &
183  flowdoms(nn, 1, 1)%djnEnd(j), &
184  flowdoms(nn, 1, 1)%dknEnd(j), &
185  flowdoms(nn, 1, 1)%neighProc(j), &
186  flowdoms(nn, 1, 1)%neighBlock(j), &
187  flowdoms(nn, 1, 1)%l1(j), &
188  flowdoms(nn, 1, 1)%l2(j), &
189  flowdoms(nn, 1, 1)%l3(j), &
190  flowdoms(nn, 1, 1)%groupNum(j), &
191  stat=ierr)
192 
193  if (ierr /= 0) &
194  call terminate("loadBalance", &
195  "Memory allocation failure for subface info")
196 
197  ! Determine the new numbering of the boundary subfaces, such
198  ! that the viscous subfaces are numbered first, followed by
199  ! the inViscid subfaces, etc.
200 
201  allocate (oldsubfaceid(blocks(i)%nBocos), stat=ierr)
202  if (ierr /= 0) &
203  call terminate("loadBalance", &
204  "Memory allocation failure for oldSubfaceID")
205 
206  call sortsubfaces(oldsubfaceid, blocks(i))
207 
208  ! Initialize the number of viscous boundary subfaces to 0.
209 
210  nviscbocos = 0
211 
212  ! Copy the info. Set the neighboring proc and block id to -1
213  ! and 0 repectively in this loop. This is okay for boundary
214  ! faces, but must be corrected for the internal block
215  ! boundaries.
216 
217  do j = 1, blocks(i)%nSubface
218 
219  ! Store the old subface id in k. For boundary faces the
220  ! sorting is taken into account; for 1 to 1 subfaces the
221  ! number is identical to the subface id in block.
222 
223  k = j
224  if (j <= blocks(i)%nBocos) k = oldsubfaceid(j)
225 
226  ! Copy the info.
227 
228  flowdoms(nn, 1, 1)%BCType(j) = blocks(i)%BCType(k)
229  flowdoms(nn, 1, 1)%BCFaceID(j) = blocks(i)%BCFaceID(k)
230  flowdoms(nn, 1, 1)%cgnsSubface(j) = blocks(i)%cgnsSubface(k)
231 
232  flowdoms(nn, 1, 1)%inBeg(j) = blocks(i)%inBeg(k)
233  flowdoms(nn, 1, 1)%jnBeg(j) = blocks(i)%jnBeg(k)
234  flowdoms(nn, 1, 1)%knBeg(j) = blocks(i)%knBeg(k)
235  flowdoms(nn, 1, 1)%inEnd(j) = blocks(i)%inEnd(k)
236  flowdoms(nn, 1, 1)%jnEnd(j) = blocks(i)%jnEnd(k)
237  flowdoms(nn, 1, 1)%knEnd(j) = blocks(i)%knEnd(k)
238 
239  flowdoms(nn, 1, 1)%dinBeg(j) = blocks(i)%dinBeg(k)
240  flowdoms(nn, 1, 1)%djnBeg(j) = blocks(i)%djnBeg(k)
241  flowdoms(nn, 1, 1)%dknBeg(j) = blocks(i)%dknBeg(k)
242  flowdoms(nn, 1, 1)%dinEnd(j) = blocks(i)%dinEnd(k)
243  flowdoms(nn, 1, 1)%djnEnd(j) = blocks(i)%djnEnd(k)
244  flowdoms(nn, 1, 1)%dknEnd(j) = blocks(i)%dknEnd(k)
245 
246  flowdoms(nn, 1, 1)%neighProc(j) = -1
247  flowdoms(nn, 1, 1)%neighBlock(j) = 0
248 
249  flowdoms(nn, 1, 1)%l1(j) = blocks(i)%l1(k)
250  flowdoms(nn, 1, 1)%l2(j) = blocks(i)%l2(k)
251  flowdoms(nn, 1, 1)%l3(j) = blocks(i)%l3(k)
252 
253  flowdoms(nn, 1, 1)%groupNum(j) = blocks(i)%groupNum(k)
254 
255  ! Update the number of viscous boundaries if this
256  ! is a viscous subface.
257 
258  if (flowdoms(nn, 1, 1)%BCType(j) == nswalladiabatic .or. &
259  flowdoms(nn, 1, 1)%BCType(j) == nswallisothermal) &
260  nviscbocos = nviscbocos + 1
261  end do
262 
263  flowdoms(nn, 1, 1:mm)%nViscBocos = nviscbocos
264 
265  ! Correct the neighboring block and proc ID for internal
266  ! block boundaries.
267 
268  do k = 1, blocks(i)%n1to1
269  j = blocks(i)%nBocos + k
270 
271  flowdoms(nn, 1, 1)%neighProc(j) = part(blocks(i)%neighBlock(j))
272  flowdoms(nn, 1, 1)%neighBlock(j) = &
273  blocks(blocks(i)%neighBlock(j))%blockID
274  end do
275 
276  ! Release the memory of oldSubfaceID.
277 
278  deallocate (oldsubfaceid, stat=ierr)
279  if (ierr /= 0) &
280  call terminate("loadBalance", &
281  "Deallocation error for oldSubfaceID")
282 
283  end if myblock
284 
285  ! Release the memory of the subface on this block.
286 
287  deallocate (blocks(i)%bcType, blocks(i)%bcFaceid, &
288  blocks(i)%cgnsSubface, blocks(i)%inBeg, &
289  blocks(i)%jnBeg, blocks(i)%knBeg, &
290  blocks(i)%inEnd, blocks(i)%jnEnd, &
291  blocks(i)%knEnd, blocks(i)%dinBeg, &
292  blocks(i)%djnBeg, blocks(i)%dknBeg, &
293  blocks(i)%dinEnd, blocks(i)%djnEnd, &
294  blocks(i)%dknEnd, blocks(i)%neighBlock, &
295  blocks(i)%l1, blocks(i)%l2, &
296  blocks(i)%l3, blocks(i)%groupNum, &
297  stat=ierr)
298  if (ierr /= 0) &
299  call terminate("loadBalance", &
300  "Deallocation error for boundary info")
301  end do domains
302 
303  ! Determine the number of processors, the processor ID's on
304  ! which the original cgns blocks are stored, the local
305  ! block ID's and the nodal ranges of the subblocks. As blocks
306  ! can be split during run-time, multiple processors can store a
307  ! part of original block.
308  !
309  ! Allocate the memory for subblocksOfCGNS.
310 
311  allocate (subblocksofcgns(nblocks), stat=ierr)
312  if (ierr /= 0) &
313  call terminate("loadBalance", &
314  "Memory allocation failure for subblocksOfCGNS")
315 
316  ! Copy the data into subblocksOfCGNS.
317 
318  do nn = 1, nblocks
319  subblocksofcgns(nn)%cgnsBlockID = blocks(nn)%cgnsBlockID
320  subblocksofcgns(nn)%procID = part(nn)
321  subblocksofcgns(nn)%blockID = blocks(nn)%blockID
322 
323  subblocksofcgns(nn)%iBegOr = blocks(nn)%iBegOr
324  subblocksofcgns(nn)%iEndOr = blocks(nn)%iEndOr
325  subblocksofcgns(nn)%jBegOr = blocks(nn)%jBegOr
326  subblocksofcgns(nn)%jEndOr = blocks(nn)%jEndOr
327  subblocksofcgns(nn)%kBegOr = blocks(nn)%kBegOr
328  subblocksofcgns(nn)%kEndOr = blocks(nn)%kEndOr
329  end do
330 
331  ! Sort subblocksOfCGNS in increasing order.
332 
333  call qsortsubblocksofcgnstype(subblocksofcgns, nblocks)
334 
335  ! Loop over the number of cgns blocks and find out the number of
336  ! subblocks it contains.
337 
338  ii = 1
339  subblockloop: do nn = 1, cgnsndom
340 
341  ! Determine the ending index jj in subblocksOfCGNS for this
342  ! CGNS block. The starting index is ii.
343 
344  if (nn == cgnsndom) then
345  jj = nblocks
346  else
347  jj = ii
348  do
349  if (subblocksofcgns(jj + 1)%cgnsBlockID > nn) exit
350  jj = jj + 1
351  end do
352  end if
353 
354  ! Set nSubBlocks and allocate the memory for procStored,
355  ! localBlockID, iBegOr, iEndOr, etc.
356 
357  cgnsdoms(nn)%nSubBlocks = jj - ii + 1
358  k = cgnsdoms(nn)%nSubBlocks
359  allocate (cgnsdoms(nn)%iBegOr(k), cgnsdoms(nn)%iEndOr(k), &
360  cgnsdoms(nn)%jBegOr(k), cgnsdoms(nn)%jEndOr(k), &
361  cgnsdoms(nn)%kBegOr(k), cgnsdoms(nn)%kEndOr(k), &
362  cgnsdoms(nn)%procStored(k), &
363  cgnsdoms(nn)%localBlockID(k), stat=ierr)
364  if (ierr /= 0) &
365  call terminate("loadBalance", &
366  "Memory allocation failure for procStored, &
367  &localBlockID, iBegOr, iEndOr, etc.")
368 
369  ! Copy the processor ID's, the local block ID's
370  ! and the subranges.
371 
372  do i = 1, cgnsdoms(nn)%nSubBlocks
373  j = i + ii - 1
374  cgnsdoms(nn)%procStored(i) = subblocksofcgns(j)%procID
375  cgnsdoms(nn)%localBlockID(i) = subblocksofcgns(j)%blockID
376 
377  cgnsdoms(nn)%iBegOr(i) = subblocksofcgns(j)%iBegOr
378  cgnsdoms(nn)%iEndOr(i) = subblocksofcgns(j)%iEndOr
379  cgnsdoms(nn)%jBegOr(i) = subblocksofcgns(j)%jBegOr
380  cgnsdoms(nn)%jEndOr(i) = subblocksofcgns(j)%jEndOr
381  cgnsdoms(nn)%kBegOr(i) = subblocksofcgns(j)%kBegOr
382  cgnsdoms(nn)%kEndOr(i) = subblocksofcgns(j)%kEndOr
383  end do
384 
385  ! Set ii for the next CGNS block.
386 
387  ii = jj + 1
388 
389  end do subblockloop
390 
391  ! Release the memory of blocks, part and subblocksOfCGNS.
392 
393  deallocate (blocks, part, subblocksofcgns, stat=ierr)
394  if (ierr /= 0) &
395  call terminate("loadBalance", &
396  "Deallocation error for blocks, part and &
397  &subblocksOfCGNS")
398 
399  !j = 20+myID
400  !do nn=1,ndom
401  ! write(j,"(8I4)") nn, flowDoms(nn,1,1)%cgnsBlockID, &
402  ! flowDoms(nn,1,1)%iBegOr, flowDoms(nn,1,1)%iEndOr, &
403  ! flowDoms(nn,1,1)%jBegOr, flowDoms(nn,1,1)%jEndOr, &
404  ! flowDoms(nn,1,1)%kBegOr, flowDoms(nn,1,1)%kEndOr
405  !enddo
406 
407  end subroutine loadbalancegrid
408 
409  subroutine blockdistribution
410  !
411  ! blockDistribution determines the distribution of the blocks
412  ! over the processors. If blocks must be split to obtain a good
413  ! load balance an iterative algorithm is used to determine the
414  ! best way to split them.
415  !
416  use constants
417  use cgnsgrid, only: cgnsdoms, cgnsndom
418  use communication, only: myid, adflow_comm_world, nproc
419  use inputparallel, only: loadbalanceiter, loadimbalance, splitblocks
420  use partitionmod, onlY: splitcgnstype, distributionblocktype, ubvec, &
421  blocks, part, nblocks, sortrangessplitinfo
422  use utils, only: terminate
423  implicit none
424  !
425  ! Local variables.
426  !
427  integer :: ierr
428 
429  integer(kind=intType) :: nn, nx, ny, nz
430  integer(kind=intType) :: nCellsTot, nCellsEven, nCellsUpper
431  integer(kind=intType) :: nCellsPerProcMax
432  integer(kind=intType) :: iter, iterMax
433 
434  logical :: cellsBalanced, facesBalanced
435  logical :: emptyPartitions, commNeglected
436 
437  type(splitcgnstype), dimension(cgnsNDom) :: splitInfo
438 
439  character(len=maxStringLen) :: loadFormat
440 
441  ! If it is not allowed to split the blocks, check that the
442  ! number of blocks is equal or larger than the number of
443  ! processors. If not, print an error message and exit.
444 
445  if (.not. splitblocks .and. nproc > cgnsndom) then
446 
447  ! Processor 0 prints the error message, while the others wait
448  ! to get killed.
449 
450  if (myid == 0) &
451  call terminate("blockDistribution", &
452  "Number of processors is larger than number &
453  &of blocks, but it is not allowed to split &
454  &blocks")
455  call mpi_barrier(adflow_comm_world, ierr)
456  end if
457  !
458  ! Determine how the blocks must be split (if allowed) for
459  ! load balancing reasons.
460  !
461  ! Initialize splitInfo to the original cgns blocks and determine
462  ! the total number of cells.
463 
464  ncellstot = 0
465  nblocks = cgnsndom
466 
467  do nn = 1, cgnsndom
468  splitinfo(nn)%nSubBlocks = 1
469 
470  allocate (splitinfo(nn)%ranges(1, 3, 2), stat=ierr)
471  if (ierr /= 0) &
472  call terminate("blockDistribution", &
473  "Memory allocation failure for ranges")
474 
475  splitinfo(nn)%ranges(1, 1, 1) = 1
476  splitinfo(nn)%ranges(1, 2, 1) = 1
477  splitinfo(nn)%ranges(1, 3, 1) = 1
478 
479  splitinfo(nn)%ranges(1, 1, 2) = cgnsdoms(nn)%il
480  splitinfo(nn)%ranges(1, 2, 2) = cgnsdoms(nn)%jl
481  splitinfo(nn)%ranges(1, 3, 2) = cgnsdoms(nn)%kl
482 
483  nx = cgnsdoms(nn)%nx
484  ny = cgnsdoms(nn)%ny
485  nz = cgnsdoms(nn)%nz
486 
487  ncellstot = ncellstot + nx * ny * nz
488  end do
489 
490  ! Determine the desirable number of cells per processor and
491  ! determine the upper bound based on the imbalance tolerance.
492  ! Initialize nCellsPerProcMax, which is used to decide whether
493  ! or not a block should be split, to this upper limit.
494 
495  ncellseven = ncellstot / nproc
496  ncellsupper = ncellseven + loadimbalance * ncellseven
497  ncellsperprocmax = ncellsupper
498 
499  ! Start the loop to determine the block topology to be used
500  ! in the computation.
501 
502  initsplitloop: do iter = 1, 2
503 
504  ! Loop over the number of original cgns blocks and determine
505  ! whether or not they must be split further.
506 
507  do nn = 1, cgnsndom
508  if (.not. splittingisokay(nn)) &
509  call splitblockinitialization(nn)
510  end do
511 
512  ! Exit the loop if the number of processors is smaller than or
513  ! equal to the number of blocks.
514 
515  if (nblocks >= nproc) exit
516 
517  ! The number of blocks is smaller than the number of processors.
518  ! Set the tolerance for splitting to the desired number of
519  ! cells on a processor.
520 
521  ncellsperprocmax = ncellseven
522 
523  end do initsplitloop
524 
525  ! Set the number of iterations to determine the distribution over
526  ! the processors. If blocks cannot be split this is set to 1;
527  ! otherwise it is set to 2.
528 
529  itermax = 1
530  if (splitblocks) itermax = loadbalanceiter
531 
532  ! Loop to determine a good load balance.
533 
534  distributionloop: do iter = 1, itermax
535 
536  ! Determine the computational blocks from the splitting info of
537  ! the original blocks.
538 
539  call determinecomputeblocks(splitinfo)
540  ! Apply the graph partitioning to the computational blocks.
541 
542  call graphpartitioning(emptypartitions, commneglected)
543 
544  ! Determine whether the load balance is okay. If empty
545  ! partitions are present the load balance is per definition
546  ! not okay and there is no need to call checkLoadBalance.
547 
548  if (emptypartitions) then
549  cellsbalanced = .false.
550  facesbalanced = .false.
551  else
552  call checkloadbalance(cellsbalanced, facesbalanced)
553  end if
554 
555  ! Exit the loop if the cells or the faces are load balanced
556  ! or if the maximum number of iterations have been reached.
557 
558  if (cellsbalanced .or. facesbalanced .or. &
559  iter == itermax) exit
560 
561  ! exit
562 
563  ! Split some blocks on the processors with too many cells/faces.
564  call splitblocksloadbalance
565 
566  end do distributionloop
567 
568  ! Deallocate the memory for ranges.
569 
570  do nn = 1, cgnsndom
571  deallocate (splitinfo(nn)%ranges, stat=ierr)
572  if (ierr /= 0) &
573  call terminate("blockDistribution", &
574  "Deallocation failure for ranges")
575  end do
576 
577  ! If empty partitions are present print an error message and
578  ! exit. Only processor 0 prints the message to avoid a mess.
579 
580  if (emptypartitions) then
581  if (myid == 0) &
582  call terminate("blockDistribution", &
583  "Empty partitions present")
584  call mpi_barrier(adflow_comm_world, ierr)
585  end if
586 
587  ! Processor 0 prints a warning if the communication was
588  ! neglected to obtain a valid partitioning.
589 
590  if (commneglected .and. myid == 0) then
591  print "(a)", "#"
592  print "(a)", "# Warning"
593  print "(a)", "# Communication costs neglected to obtain a valid partitioning."
594  end if
595 
596  ! If the load imbalance tolerance was not met, print a warning
597  ! message if I am processor 0.
598 
599  loadformat = "(*(A, F6.3))"
600 
601  if (myid == 0) then
602  if (.not. (cellsbalanced .and. facesbalanced)) then
603  print "(a)", "#"
604  print "(a)", "# Warning"
605  print loadformat, "# Specified load imbalance tolerance", real(loadimbalance), " not achieved."
606  print loadformat, "# Continuing with", real(ubvec(1)), &
607  " load imbalance for the cells and", real(ubvec(2)), " for the faces"
608  print "(a)", "#"
609  else
610  print "(a)", "#"
611  print loadformat, "# Specified load imbalance tolerance", real(loadimbalance), " achieved."
612  print loadformat, "# Continuing with", real(ubvec(1)), &
613  " load imbalance for the cells and", real(ubvec(2)), " for the faces"
614  print "(a)", "#"
615  end if
616 
617  end if
618 
619  !=================================================================
620 
621  contains
622 
623  !===============================================================
624 
625  logical function splittingisokay(cgnsID)
626  !
627  ! splittingIsOkay determines whether or not the splitting of
628  ! the given cgns block is okay in the sense that all subblocks
629  ! are smaller than the allowed number of cells and faces.
630  !
631  implicit none
632  !
633  ! Function argument
634  !
635  integer(kind=intType), intent(in) :: cgnsid
636  !
637  ! Local variables.
638  !
639  integer(kind=intType) :: i, nx, ny, nz, ncells
640 
641  ! Initialize splittingIsOkay to .true.
642 
643  splittingisokay = .true.
644 
645  ! Loop over the subblocks.
646 
647  do i = 1, splitinfo(cgnsid)%nSubBlocks
648 
649  ! Determine the number of cells in the three directions.
650 
651  nx = splitinfo(cgnsid)%ranges(i, 1, 2) &
652  - splitinfo(cgnsid)%ranges(i, 1, 1)
653  ny = splitinfo(cgnsid)%ranges(i, 2, 2) &
654  - splitinfo(cgnsid)%ranges(i, 2, 1)
655  nz = splitinfo(cgnsid)%ranges(i, 3, 2) &
656  - splitinfo(cgnsid)%ranges(i, 3, 1)
657 
658  ! Determine the number of cells for this subblock.
659 
660  ncells = nx * ny * nz
661 
662  ! Check whether this number is smaller or equal to the
663  ! maximum allowed number. If not, set splittingIsOkay to
664  ! .false. and exit the loop.
665 
666  if (ncells > ncellsperprocmax) then
667  splittingisokay = .false.
668  exit
669  end if
670 
671  end do
672 
673  end function splittingisokay
674 
675  !===============================================================
676 
677  subroutine splitblockinitialization(cgnsID)
678  !
679  ! splitBlockInitialization splits the given cgns block ID
680  ! into a number of subbocks during the initialization phase.
681  !
682  implicit none
683  !
684  ! Subroutine argument.
685  !
686  integer(kind=intType), intent(in) :: cgnsID
687  !
688  ! Local variables.
689  !
690  integer :: ierr
691 
692  integer(kind=intType) :: nn, mm
693  integer(kind=intType) :: nx, ny, nz, nCells, nSub
694 
695  integer(kind=intType), dimension(:, :, :), allocatable :: tmpRange
696 
697  type(distributionblocktype) :: tmpBlock
698 
699  ! Check whether it is allowed to split blocks.
700 
701  if (.not. splitblocks) &
702  call terminate("splitBlockInitialization", &
703  "Block must be split for load balance, &
704  &but I am not allowed to do so")
705 
706  ! Store the number of cells of this block a bit easier.
707 
708  nx = cgnsdoms(cgnsid)%nx
709  ny = cgnsdoms(cgnsid)%ny
710  nz = cgnsdoms(cgnsid)%nz
711  ncells = nx * ny * nz
712 
713  ! Copy the information from the given cgns block into
714  ! tmpBlock, such that the routine splitBlock can be used.
715 
716  ! First the scalar info.
717 
718  tmpblock%nx = nx; tmpblock%il = nx + 1
719  tmpblock%ny = ny; tmpblock%jl = ny + 1
720  tmpblock%nz = nz; tmpblock%kl = nz + 1
721 
722  tmpblock%nCell = ncells
723  tmpblock%nface = (nx + 1) * ny * nz + (ny + 1) * nx * nz &
724  + (nz + 1) * nx * ny
725 
726  tmpblock%cgnsBlockID = cgnsid
727 
728  tmpblock%iBegor = 1; tmpblock%iEndor = tmpblock%il
729  tmpblock%jBegor = 1; tmpblock%jEndor = tmpblock%jl
730  tmpblock%kBegor = 1; tmpblock%kEndor = tmpblock%kl
731 
732  tmpblock%nBocos = cgnsdoms(cgnsid)%nBocos
733  tmpblock%n1to1 = cgnsdoms(cgnsid)%n1to1
734  tmpblock%nSubface = tmpblock%nBocos + tmpblock%n1to1
735 
736  ! Allocate the memory for the nodal ranges of the subfaces
737  ! and nullify the other pointers.
738 
739  nsub = tmpblock%nSubface
740 
741  allocate (tmpblock%BCType(nsub), tmpblock%BCFaceID(nsub), &
742  tmpblock%inBeg(nsub), tmpblock%inEnd(nsub), &
743  tmpblock%jnBeg(nsub), tmpblock%jnEnd(nsub), &
744  tmpblock%knBeg(nsub), tmpblock%knEnd(nsub), &
745  stat=ierr)
746  if (ierr /= 0) &
747  call terminate("splitBlockInitialization", &
748  "Deallocation failure for the subface &
749  &info in tmpBlock")
750 
751  nullify (tmpblock%cgnsSubface, tmpblock%neighBlock, &
752  tmpblock%dinBeg, tmpblock%dinEnd, &
753  tmpblock%djnBeg, tmpblock%djnEnd, &
754  tmpblock%dknBeg, tmpblock%dknEnd, &
755  tmpblock%l1, tmpblock%l2, &
756  tmpblock%l3, tmpblock%groupNum)
757 
758  ! Copy the range of the subfaces into tmpBlock and set the
759  ! corresponding boundary condition.
760  ! First the boundary faces.
761 
762  do nn = 1, cgnsdoms(cgnsid)%nBocos
763  tmpblock%inBeg(nn) = cgnsdoms(cgnsid)%bocoInfo(nn)%iBeg
764  tmpblock%jnBeg(nn) = cgnsdoms(cgnsid)%bocoInfo(nn)%jBeg
765  tmpblock%knBeg(nn) = cgnsdoms(cgnsid)%bocoInfo(nn)%kBeg
766 
767  tmpblock%inEnd(nn) = cgnsdoms(cgnsid)%bocoInfo(nn)%iEnd
768  tmpblock%jnEnd(nn) = cgnsdoms(cgnsid)%bocoInfo(nn)%jEnd
769  tmpblock%knEnd(nn) = cgnsdoms(cgnsid)%bocoInfo(nn)%kEnd
770 
771  tmpblock%BCType(nn) = cgnsdoms(cgnsid)%bocoInfo(nn)%BCType
772  end do
773 
774  ! And the internal block faces; set the boundary condition to
775  ! B2BMatch.
776 
777  do mm = 1, cgnsdoms(cgnsid)%n1to1
778  nn = mm + cgnsdoms(cgnsid)%nBocos
779 
780  tmpblock%inBeg(nn) = cgnsdoms(cgnsid)%conn1to1(mm)%iBeg
781  tmpblock%jnBeg(nn) = cgnsdoms(cgnsid)%conn1to1(mm)%jBeg
782  tmpblock%knBeg(nn) = cgnsdoms(cgnsid)%conn1to1(mm)%kBeg
783 
784  tmpblock%inEnd(nn) = cgnsdoms(cgnsid)%conn1to1(mm)%iEnd
785  tmpblock%jnEnd(nn) = cgnsdoms(cgnsid)%conn1to1(mm)%jEnd
786  tmpblock%knEnd(nn) = cgnsdoms(cgnsid)%conn1to1(mm)%kEnd
787 
788  tmpblock%BCType(nn) = b2bmatch
789  end do
790 
791  ! Determine for the block face on which the subface is located.
792  ! Assume that the subface connectivity is correct. This will be
793  ! tested later on in determineComputeBlocks.
794 
795  do nn = 1, tmpblock%nSubface
796  if (tmpblock%inBeg(nn) == tmpblock%inEnd(nn)) then
797  tmpblock%BCFaceID(nn) = imax
798  if (tmpblock%inBeg(nn) == 1) tmpblock%BCFaceID(nn) = imin
799  else if (tmpblock%jnBeg(nn) == tmpblock%jnEnd(nn)) then
800  tmpblock%BCFaceID(nn) = jmax
801  if (tmpblock%jnBeg(nn) == 1) tmpblock%BCFaceID(nn) = jmin
802  else
803  tmpblock%BCFaceID(nn) = kmax
804  if (tmpblock%knBeg(nn) == 1) tmpblock%BCFaceID(nn) = kmin
805  end if
806  end do
807 
808  ! Determine the number of subblocks into which this block is
809  ! to be split.
810 
811  nsub = nint(real(ncells, realtype) / real(ncellseven, realtype))
812  nsub = max(nsub, 2_inttype)
813 
814  ! Deallocate the memory of the splitting info and allocate
815  ! the memory of tmpRange.
816 
817  deallocate (splitinfo(cgnsid)%ranges, stat=ierr)
818  if (ierr /= 0) &
819  call terminate("splitBlockInitialization", &
820  "Deallocation failure for ranges")
821 
822  allocate (tmprange(nsub, 3, 2), stat=ierr)
823  if (ierr /= 0) &
824  call terminate("splitBlockInitialization", &
825  "Memory allocation failure for tmpRange")
826 
827  ! Split the block into the subblocks. It is possible that
828  ! the block is split into less subblocks than the desired
829  ! number. The actual number of subblocks is returned in nSub.
830 
831  call splitblock(tmpblock, nsub, ncellseven, tmprange)
832 
833  ! Allocate the memory for ranges.
834 
835  allocate (splitinfo(cgnsid)%ranges(nsub, 3, 2), stat=ierr)
836  if (ierr /= 0) &
837  call terminate("splitBlockInitialization", &
838  "Memory allocation failure for ranges")
839 
840  ! Determine the new number of computational blocks. This is
841  ! the old number plus the number or extra blocks created
842  ! by the splitting.
843 
844  nblocks = nblocks + nsub - splitinfo(cgnsid)%nSubBlocks
845 
846  ! Copy tmpRange into ranges of splitInfo.
847 
848  splitinfo(cgnsid)%nSubBlocks = nsub
849 
850  do nn = 1, nsub
851  splitinfo(cgnsid)%ranges(nn, 1, 1) = tmprange(nn, 1, 1)
852  splitinfo(cgnsid)%ranges(nn, 2, 1) = tmprange(nn, 2, 1)
853  splitinfo(cgnsid)%ranges(nn, 3, 1) = tmprange(nn, 3, 1)
854 
855  splitinfo(cgnsid)%ranges(nn, 1, 2) = tmprange(nn, 1, 2)
856  splitinfo(cgnsid)%ranges(nn, 2, 2) = tmprange(nn, 2, 2)
857  splitinfo(cgnsid)%ranges(nn, 3, 2) = tmprange(nn, 3, 2)
858  end do
859 
860  ! Deallocate the memory allocated for tmpBlock and tmpRange.
861 
862  deallocate (tmpblock%BCType, tmpblock%BCFaceID, &
863  tmpblock%inBeg, tmpblock%inEnd, &
864  tmpblock%jnBeg, tmpblock%jnEnd, &
865  tmpblock%knBeg, tmpblock%knEnd, &
866  tmprange, stat=ierr)
867  if (ierr /= 0) &
868  call terminate("splitBlockInitialization", &
869  "Deallocation failure for tmpBlock &
870  &and tmpRange")
871 
872  ! Sort the subranges of this block in increasing order.
873 
874  call sortrangessplitinfo(splitinfo(cgnsid))
875 
876  end subroutine splitblockinitialization
877 
878  !===============================================================
879 
880  subroutine splitblocksloadbalance
881  !
882  ! splitBlocksLoadBalance splits some (sub)blocks even
883  ! further to obtain a better load balance.
884  !
885  use sorting, only: bsearchintegers, qsortintegers
886  implicit none
887  !
888  ! Local variables.
889  !
890  integer :: ierr
891 
892  integer(kind=intType) :: i, j, k, kk, jj
893  integer(kind=intType) :: nCellDiff, proc, nSub, cgnsID
894  integer(kind=intType) :: ncgnsSort
895 
896  integer(kind=intType), dimension(0:cgnsNDom) :: nSubBlocks
897  integer(kind=intType), dimension(cgnsNDom) :: cgnsSort
898 
899  integer(kind=intType), dimension(nProc) :: nCell, nCellOr
900  integer(kind=intType), dimension(nProc) :: tmp, procNCell
901  integer(kind=intType), dimension(0:nProc) :: nBlockProc
902  integer(kind=intType), dimension(0:nProc) :: multNCell
903  integer(kind=intType), dimension(nBlocks) :: blockProc
904 
905  integer(kind=intType), dimension(:, :, :), allocatable :: tmpRange
906  integer(kind=intType), dimension(:, :, :), pointer :: oldRanges
907 
908  logical, dimension(cgnsNDom) :: cgnsBlockFlagged
909 
910  ! Determine the number of blocks per cgns block in cumulative
911  ! storage format. These values will serve as an offset to
912  ! determine the local subblock ID.
913  ! Initialize in the same loop cgnsBlockFlagged to .false.
914 
915  nsubblocks(0) = 0
916  do i = 1, cgnsndom
917  nsubblocks(i) = nsubblocks(i - 1) + splitinfo(i)%nSubBlocks
918  cgnsblockflagged(i) = .false.
919  end do
920 
921  ! Determine the number of cells and blocks per processor.
922  ! Note that part(i) contains the processor ID, which start at 0.
923 
924  ncell = 0
925  nblockproc = 0
926 
927  do i = 1, nblocks
928  j = part(i) + 1
929  ncell(j) = ncell(j) + blocks(i)%nCell
930  nblockproc(j) = nblockproc(j) + 1
931  end do
932 
933  ! Put nBlockProc in cumulative storage format and store the
934  ! starting entries in tmp, which will be used as a counter to
935  ! put blockProc in the correct place.
936 
937  do i = 1, nproc
938  tmp(i) = nblockproc(i - 1)
939  nblockproc(i) = nblockproc(i) + tmp(i)
940  end do
941 
942  ! Determine the block ID's for every processor. Again note
943  ! that part(i) starts at 0; tmp is used as a counter variable.
944 
945  do i = 1, nblocks
946  j = part(i) + 1
947  tmp(j) = tmp(j) + 1
948  blockproc(tmp(j)) = i
949  end do
950 
951  ! Sort nCell in increasing order. Store the original values.
952 
953  ncellor = ncell
954  ncelldiff = nproc
955  call qsortintegers(ncell, ncelldiff)
956 
957  ! Determine the different number of cells per proc and store
958  ! the multiplicity.
959 
960  multncell(0) = 0
961  multncell(1) = 1
962  ncelldiff = 1
963 
964  do i = 2, nproc
965  if (ncell(i) == ncell(ncelldiff)) then
966  multncell(ncelldiff) = multncell(ncelldiff) + 1
967  else
968  ncelldiff = ncelldiff + 1
969  multncell(ncelldiff) = multncell(ncelldiff - 1) + 1
970  ncell(ncelldiff) = ncell(i)
971  end if
972  end do
973 
974  ! Set tmp to multNCell; it will serve as a counter to determine
975  ! the corresponding processor ID's of the sorted cell numbers
976 
977  do i = 1, ncelldiff
978  tmp(i) = multncell(i - 1)
979  end do
980 
981  ! Determine the processor ID's corresponding to the sorted
982  ! number of cells. Note that the processor ID's start at 0.
983 
984  do i = 1, nproc
985  j = bsearchintegers(ncellor(i), ncell(1:ncelldiff))
986  tmp(j) = tmp(j) + 1
987  procncell(tmp(j)) = i - 1
988  end do
989 
990  ! Loop to subdivide blocks. As nCell is sorted in increasing
991  ! order, the largest number of cells are at the end of this
992  ! array. Therefore this loop starts at the back.
993 
994  i = ncelldiff
995  ncgnssort = 0
996  divisionloop: do
997 
998  ! Condition to exit the loop. The number of cells per
999  ! processor is allowed in the load balance.
1000 
1001  if (ncell(i) <= ncellsupper) exit
1002 
1003  ! Loop over the processors with this number of cells.
1004 
1005  processorloop: do j = (multncell(i - 1) + 1), multncell(i)
1006 
1007  ! Store the processor ID a bit easier.
1008 
1009  proc = procncell(j)
1010 
1011  ! Determine the largest block of this processor, which will
1012  ! be stored in index jj. As only the largest is needed a
1013  ! linear search algorithm is okay.
1014 
1015  k = nblockproc(proc) + 1
1016  jj = blockproc(k)
1017 
1018  do k = (nblockproc(proc) + 2), nblockproc(proc + 1)
1019  kk = blockproc(k)
1020  if (blocks(kk)%nCell > blocks(jj)%nCell) jj = kk
1021  end do
1022 
1023  ! Determine the local subblock number of computational block
1024  ! jj in the original cgns block. This will be stored in kk.
1025 
1026  cgnsid = blocks(jj)%cgnsBlockID
1027  kk = jj - nsubblocks(cgnsid - 1)
1028 
1029  ! Determine the number of cells, which should be split from
1030  ! block jj, such that the load balance of processor proc
1031  ! would be optimal (in terms of cells).
1032 
1033  ncelldiff = ncell(i) - ncellseven
1034 
1035  ! A very theoretical case would be that nCellDiff is larger
1036  ! than the number of cells of block jj. This could occur,
1037  ! because the partitioning is a multi-constraint
1038  ! partitioning; however, it is not likely. In this case
1039  ! nothing is done and the next processor is treated.
1040 
1041  if (ncelldiff >= blocks(jj)%nCell) cycle
1042 
1043  ! Determine the situation we are having here. If nCellDiff
1044  ! is less or equal to the number of cells to obtain a good
1045  ! load balance for the least loaded processor the block will
1046  ! be split into 2. Test this.
1047 
1048  if (ncelldiff <= (ncellsupper - ncell(1))) then
1049 
1050  nsub = 2
1051 
1052  else
1053 
1054  ! nCellDiff is larger than the number of cells to create
1055  ! a good load balance for the least loaded processor.
1056  ! This probably means that the block should be split in
1057  ! more than 2 subblocks.
1058  ! Store the number of cells which should be kept in block
1059  ! jj for an optimal load balance in nCellDiff.
1060 
1061  ncelldiff = blocks(jj)%nCell - ncelldiff
1062 
1063  ! Determine the number of subblocks.
1064 
1065  nsub = nint(real(blocks(jj)%nCell, realtype) &
1066  / real(ncelldiff, realtype))
1067  nsub = max(nsub, 2_inttype)
1068 
1069  end if
1070 
1071  ! Allocate the memory for tmpRange, which will store the
1072  ! block splitting of block jj.
1073 
1074  allocate (tmprange(nsub, 3, 2), stat=ierr)
1075  if (ierr /= 0) &
1076  call terminate("splitBlocksLoadBalance", &
1077  "Memory allocation failure for tmpRange")
1078 
1079  ! Split computational block into the desired number of
1080  ! subblocks, if possible. The actual number of splittings
1081  ! is returned in nSub.
1082 
1083  call splitblock(blocks(jj), nsub, ncelldiff, tmprange)
1084 
1085  ! Store the old ranges, store the old number of subblocks
1086  ! in jj, determine the new number of subblocks and allocate
1087  ! the memory for ranges.
1088 
1089  jj = splitinfo(cgnsid)%nSubBlocks
1090  oldranges => splitinfo(cgnsid)%ranges
1091 
1092  splitinfo(cgnsid)%nSubBlocks = jj + nsub - 1
1093 
1094  k = splitinfo(cgnsid)%nSubBlocks
1095  allocate (splitinfo(cgnsid)%ranges(k, 3, 2), stat=ierr)
1096  if (ierr /= 0) &
1097  call terminate("splitBlocksLoadBalance", &
1098  "Memory allocation failure for ranges")
1099 
1100  ! Determine the new number of computational blocks.
1101  ! This is the old number plus the number or extra blocks
1102  ! created by the splitting.
1103 
1104  nblocks = nblocks + splitinfo(cgnsid)%nSubBlocks - jj
1105 
1106  ! Copy the original info back into ranges.
1107 
1108  do k = 1, jj
1109  splitinfo(cgnsid)%ranges(k, 1, 1) = oldranges(k, 1, 1)
1110  splitinfo(cgnsid)%ranges(k, 2, 1) = oldranges(k, 2, 1)
1111  splitinfo(cgnsid)%ranges(k, 3, 1) = oldranges(k, 3, 1)
1112 
1113  splitinfo(cgnsid)%ranges(k, 1, 2) = oldranges(k, 1, 2)
1114  splitinfo(cgnsid)%ranges(k, 2, 2) = oldranges(k, 2, 2)
1115  splitinfo(cgnsid)%ranges(k, 3, 2) = oldranges(k, 3, 2)
1116  end do
1117 
1118  ! The splitting kk has disappeared. Replace it by the
1119  ! first splitting in tmpRange.
1120 
1121  splitinfo(cgnsid)%ranges(kk, 1, 1) = tmprange(1, 1, 1)
1122  splitinfo(cgnsid)%ranges(kk, 2, 1) = tmprange(1, 2, 1)
1123  splitinfo(cgnsid)%ranges(kk, 3, 1) = tmprange(1, 3, 1)
1124 
1125  splitinfo(cgnsid)%ranges(kk, 1, 2) = tmprange(1, 1, 2)
1126  splitinfo(cgnsid)%ranges(kk, 2, 2) = tmprange(1, 2, 2)
1127  splitinfo(cgnsid)%ranges(kk, 3, 2) = tmprange(1, 3, 2)
1128 
1129  ! Add the other splittings to the end of ranges.
1130 
1131  do k = 2, nsub
1132  jj = jj + 1
1133 
1134  splitinfo(cgnsid)%ranges(jj, 1, 1) = tmprange(k, 1, 1)
1135  splitinfo(cgnsid)%ranges(jj, 2, 1) = tmprange(k, 2, 1)
1136  splitinfo(cgnsid)%ranges(jj, 3, 1) = tmprange(k, 3, 1)
1137 
1138  splitinfo(cgnsid)%ranges(jj, 1, 2) = tmprange(k, 1, 2)
1139  splitinfo(cgnsid)%ranges(jj, 2, 2) = tmprange(k, 2, 2)
1140  splitinfo(cgnsid)%ranges(jj, 3, 2) = tmprange(k, 3, 2)
1141  end do
1142 
1143  ! Deallocate the memory of oldRanges and tmpRange.
1144 
1145  deallocate (oldranges, tmprange, stat=ierr)
1146  if (ierr /= 0) &
1147  call terminate("splitBlocksLoadBalance", &
1148  "Deallocation failure for oldRanges and &
1149  &tmpRange")
1150 
1151  ! Store the cgns ID, if not already stored. In this way
1152  ! it is known for which cgns blocks the subranges must be
1153  ! sorted.
1154 
1155  if (.not. cgnsblockflagged(cgnsid)) then
1156  cgnsblockflagged(cgnsid) = .true.
1157  ncgnssort = ncgnssort + 1
1158  cgnssort(ncgnssort) = cgnsid
1159  end if
1160 
1161  end do processorloop
1162 
1163  ! Decrease the counter i for the next number of cells.
1164 
1165  i = i - 1
1166 
1167  end do divisionloop
1168 
1169  ! Sort the subranges of the CGNS blocks that were split in this
1170  ! routine in increasing order.
1171 
1172  do i = 1, ncgnssort
1173  call sortrangessplitinfo(splitinfo(cgnssort(i)))
1174  end do
1175 
1176  end subroutine splitblocksloadbalance
1177 
1178  end subroutine blockdistribution
1179 
1180  subroutine initflowdoms
1181  !
1182  ! initFlowDoms allocates the memory for flowDoms and initializes
1183  ! its pointers to null pointers, such that they do not have
1184  ! random targets.
1185  !
1186  use constants
1187  use block, only: flowdoms, ndom
1188  use inputiteration, only: nmglevels, mgstartlevel
1189  use inputtimespectral, only: ntimeintervalsspectral
1190  use utils, only: terminate, nullifyflowdompointers
1191  implicit none
1192  !
1193  ! Local variables.
1194  !
1195  integer :: ierr
1196 
1197  integer(kind=intType) :: i, j, k, nn
1198 
1199  ! Allocate the memory for flowDoms. Set nn to the maximum of the
1200  ! number of mg levels needed in the cycle and mg start level.
1201  ! This is namely the amount of grid levels the solver needs.
1202 
1203  nn = max(nmglevels, mgstartlevel)
1204  allocate (flowdoms(ndom, nn, ntimeintervalsspectral), stat=ierr)
1205  if (ierr /= 0) &
1206  call terminate("initFlowDoms", &
1207  "Memory allocation failure for flowDoms")
1208 
1209  ! Loop over all the blocks and initialize its pointers to the
1210  ! null-pointer.
1211 
1212  do k = 1, ntimeintervalsspectral
1213  do j = 1, nn
1214  do i = 1, ndom
1215  call nullifyflowdompointers(i, j, k)
1216  end do
1217  end do
1218  end do
1219 
1220  end subroutine initflowdoms
1221  subroutine sortsubfaces(oldSubfaceID, blockID)
1222  !
1223  ! sortSubfaces sorts the boundary subfaces of the given block
1224  ! such that viscous subfaces are numbered first, followed by
1225  ! inviscid, etc.
1226  !
1227  use constants
1228  use partitionmod, only: distributionblocktype
1229  use sorting, only: qsortintegers, bsearchintegers
1230  implicit none
1231  !
1232  ! Subroutine arguments
1233  !
1234  integer(kind=intType), dimension(*), intent(out) :: oldSubfaceID
1235  type(distributionblocktype), intent(in) :: blockID
1236  !
1237  ! Local variables.
1238  !
1239  integer(kind=intType) :: i, ii, nDiff
1240 
1241  integer(kind=intType), dimension(blockID%nBocos) :: bcPrior
1242  integer(kind=intType), dimension(blockID%nBocos) :: bcPriorSort
1243 
1244  integer(kind=intType), dimension(0:blockID%nBocos) :: mult
1245 
1246  ! Loop over the boundary subfaces and determine the priorities.
1247  ! Store the priorities in both bcPrior and bcPriorSort.
1248 
1249  do i = 1, blockid%nBocos
1250 
1251  select case (blockid%BCType(i))
1252 
1253  case (nswalladiabatic)
1254  bcprior(i) = 1
1255 
1256  case (nswallisothermal)
1257  bcprior(i) = 2
1258 
1259  case (eulerwall)
1260  bcprior(i) = 3
1261 
1262  case (symm)
1263  bcprior(i) = 4
1264 
1265  case (symmpolar)
1266  bcprior(i) = 5
1267 
1268  case (farfield)
1269  bcprior(i) = 6
1270 
1271  case (supersonicinflow)
1272  bcprior(i) = 7
1273 
1274  case (subsonicinflow)
1275  bcprior(i) = 8
1276 
1277  case (supersonicoutflow)
1278  bcprior(i) = 9
1279 
1280  case (subsonicoutflow)
1281  bcprior(i) = 10
1282 
1283  case (massbleedinflow)
1284  bcprior(i) = 11
1285 
1286  case (massbleedoutflow)
1287  bcprior(i) = 12
1288 
1289  case (mdot)
1290  bcprior(i) = 13
1291 
1292  case (bcthrust)
1293  bcprior(i) = 14
1294 
1295  case (extrap)
1296  bcprior(i) = 15
1297 
1298  case (slidinginterface)
1299  bcprior(i) = 19
1300 
1301  case (oversetouterbound)
1302  bcprior(i) = 20
1303 
1304  case (domaininterfaceall)
1305  bcprior(i) = 21
1306 
1307  case (domaininterfacerhouvw)
1308  bcprior(i) = 22
1309 
1310  case (domaininterfacep)
1311  bcprior(i) = 23
1312 
1313  case (domaininterfacerho)
1314  bcprior(i) = 24
1315 
1316  case (domaininterfacetotal)
1317  bcprior(i) = 25
1318 
1319  end select
1320 
1321  bcpriorsort(i) = bcprior(i)
1322 
1323  end do
1324 
1325  ! Sort bcPriorSort in increasing order.
1326 
1327  call qsortintegers(bcpriorsort, blockid%nBocos)
1328 
1329  ! Get rid of the multiple entries and store the multiplicity in
1330  ! cumulative storage format. nDiff contains the number of
1331  ! different boundary conditions for this block. The initialization
1332  ! with the min function is necessary to be able to treat blocks
1333  ! without a boundary condition correctly.
1334 
1335  ndiff = min(1_inttype, blockid%nBocos)
1336  mult(0) = 0
1337  mult(ndiff) = 1
1338 
1339  do i = 2, blockid%nBocos
1340  if (bcpriorsort(i) == bcpriorsort(ndiff)) then
1341  mult(ndiff) = mult(ndiff) + 1
1342  else
1343  ndiff = ndiff + 1
1344  mult(ndiff) = mult(ndiff - 1) + 1
1345  bcpriorsort(ndiff) = bcpriorsort(i)
1346  end if
1347  end do
1348 
1349  ! Determine the old subface ID by searching in the sorted
1350  ! priorities.
1351 
1352  do i = 1, blockid%nBocos
1353 
1354  ii = bsearchintegers(bcprior(i), bcpriorsort(1:ndiff))
1355 
1356  ! Update the lower boundary in the multiplicity for this
1357  ! boundary condition and store the entry a bit easier.
1358 
1359  mult(ii - 1) = mult(ii - 1) + 1
1360  ii = mult(ii - 1)
1361 
1362  ! Store the mapping of the new to old subface id.
1363 
1364  oldsubfaceid(ii) = i
1365 
1366  end do
1367 
1368  end subroutine sortsubfaces
1369 
1370  subroutine determinecomputeblocks(splitInfo)
1371  !
1372  ! determineComputeBlocks determines the computational blocks
1373  ! from the original grid and the given information how to split
1374  ! these blocks.
1375  !
1376  use constants
1377  use cgnsgrid, only: cgnsndom, cgnsdoms
1378  use partitionmod, only: blocks, nblocks, splitcgnstype
1379  use utils, only: terminate
1380  implicit none
1381  !
1382  ! Subroutine arguments.
1383  !
1384  type(splitcgnstype), dimension(cgnsNDom), intent(in) :: splitInfo
1385  !
1386  ! Local variables.
1387  !
1388  integer :: ierr
1389 
1390  integer(kind=intType) :: i, ii, jj, mm
1391  integer(kind=intType) :: nx, ny, nz, nAlloc
1392 
1393  integer(kind=intType), dimension(0:cgnsNDom) :: nSubPerCGNS
1394 
1395  ! Determine the number of subblocks per cgns block in
1396  ! cumulative storage format.
1397 
1398  nsubpercgns(0) = 0
1399  do i = 1, cgnsndom
1400  nsubpercgns(i) = nsubpercgns(i - 1) + splitinfo(i)%nSubblocks
1401  end do
1402 
1403  ! Check whether blocks are already allocated. If so, release
1404  ! the memory.
1405  if (allocated(blocks)) then
1406 
1407  ! Loop over the number of old blocks and release the memory.
1408 
1409  mm = ubound(blocks, 1)
1410  do i = 1, mm
1411  deallocate (blocks(i)%BCType, blocks(i)%BCFaceID, &
1412  blocks(i)%cgnsSubface, blocks(i)%inBeg, &
1413  blocks(i)%jnBeg, blocks(i)%knBeg, &
1414  blocks(i)%inEnd, blocks(i)%jnEnd, &
1415  blocks(i)%knEnd, blocks(i)%dinBeg, &
1416  blocks(i)%djnBeg, blocks(i)%dknBeg, &
1417  blocks(i)%dinEnd, blocks(i)%djnEnd, &
1418  blocks(i)%dknEnd, blocks(i)%neighBlock, &
1419  blocks(i)%l1, blocks(i)%l2, &
1420  blocks(i)%l3, blocks(i)%groupNum, &
1421  stat=ierr)
1422  if (ierr /= 0) &
1423  call terminate("determineComputeBlocks", &
1424  "Deallocation error for subface info")
1425  end do
1426 
1427  ! Release the memory of blocks itself.
1428 
1429  deallocate (blocks, stat=ierr)
1430  if (ierr /= 0) &
1431  call terminate("determineComputeBlocks", &
1432  "Deallocation error for blocks")
1433  end if
1434 
1435  ! Allocate the memory for blocks.
1436 
1437  allocate (blocks(nblocks), stat=ierr)
1438  if (ierr /= 0) &
1439  call terminate("determineComputeBlocks", &
1440  "Memory allocation failure for blocks")
1441 
1442  ! Set the counter ii for the global computational block number
1443  ! and loop over the cgns blocks.
1444 
1445  ii = 0
1446  cgnsloop: do i = 1, cgnsndom
1447 
1448  ! Loop over the number of subblocks of this cgns block.
1449 
1450  subblockloop: do mm = 1, splitinfo(i)%nSubblocks
1451 
1452  ! Update the counter ii and store the number of cells in
1453  ! the three directions in nx, ny and nz.
1454 
1455  ii = ii + 1
1456  nx = splitinfo(i)%ranges(mm, 1, 2) &
1457  - splitinfo(i)%ranges(mm, 1, 1)
1458  ny = splitinfo(i)%ranges(mm, 2, 2) &
1459  - splitinfo(i)%ranges(mm, 2, 1)
1460  nz = splitinfo(i)%ranges(mm, 3, 2) &
1461  - splitinfo(i)%ranges(mm, 3, 1)
1462 
1463  ! Initialize the scalar variables of blocks(ii).
1464 
1465  blocks(ii)%nx = nx; blocks(ii)%il = nx + 1
1466  blocks(ii)%ny = ny; blocks(ii)%jl = ny + 1
1467  blocks(ii)%nz = nz; blocks(ii)%kl = nz + 1
1468 
1469  blocks(ii)%ncell = nx * ny * nz
1470  blocks(ii)%nface = (nx + 1) * ny * nz + (ny + 1) * nx * nz &
1471  + (nz + 1) * nx * ny
1472 
1473  blocks(ii)%cgnsBlockID = i
1474 
1475  blocks(ii)%iBegor = splitinfo(i)%ranges(mm, 1, 1)
1476  blocks(ii)%jBegor = splitinfo(i)%ranges(mm, 2, 1)
1477  blocks(ii)%kBegor = splitinfo(i)%ranges(mm, 3, 1)
1478 
1479  blocks(ii)%iEndor = splitinfo(i)%ranges(mm, 1, 2)
1480  blocks(ii)%jEndor = splitinfo(i)%ranges(mm, 2, 2)
1481  blocks(ii)%kEndor = splitinfo(i)%ranges(mm, 3, 2)
1482 
1483  blocks(ii)%nBocos = 0
1484  blocks(ii)%nSubface = 0
1485  blocks(ii)%n1to1 = 0
1486 
1487  ! Do an allocation for the subface info. NAlloc is such
1488  ! that no reallocation is needed for the boundary info.
1489 
1490  nalloc = cgnsdoms(i)%nBocos + cgnsdoms(i)%n1to1
1491 
1492  allocate (blocks(ii)%BCType(nalloc), &
1493  blocks(ii)%BCFaceID(nalloc), &
1494  blocks(ii)%cgnsSubface(nalloc), &
1495  blocks(ii)%inBeg(nalloc), &
1496  blocks(ii)%jnBeg(nalloc), &
1497  blocks(ii)%knBeg(nalloc), &
1498  blocks(ii)%inEnd(nalloc), &
1499  blocks(ii)%jnEnd(nalloc), &
1500  blocks(ii)%knEnd(nalloc), &
1501  blocks(ii)%dinBeg(nalloc), &
1502  blocks(ii)%djnBeg(nalloc), &
1503  blocks(ii)%dknBeg(nalloc), &
1504  blocks(ii)%dinEnd(nalloc), &
1505  blocks(ii)%djnEnd(nalloc), &
1506  blocks(ii)%dknEnd(nalloc), &
1507  blocks(ii)%neighBlock(nalloc), &
1508  blocks(ii)%l1(nalloc), &
1509  blocks(ii)%l2(nalloc), &
1510  blocks(ii)%l3(nalloc), &
1511  blocks(ii)%groupNum(nalloc), &
1512  stat=ierr)
1513  if (ierr /= 0) &
1514  call terminate("determineComputeBlocks", &
1515  "Memory allocation failure for &
1516  &subface info")
1517 
1518  ! Determine the boundary condition subfaces and the subfaces
1519  ! of the subblock, which are located on the block boundaries
1520  ! of the original cgns block.
1521 
1522  jj = 0
1523  call bcfacessubblock(i, ii, jj)
1524  blocks(ii)%nBocos = jj
1525 
1526  call externalfacessubblock(i, ii, jj, nsubpercgns, &
1527  nalloc, splitinfo)
1528 
1529  ! Determine the subfaces of the subblock created by the
1530  ! splitting of the original block.
1531 
1532  call internalfacessubblock(i, ii, jj, nsubpercgns, &
1533  nalloc, splitinfo(i))
1534  blocks(ii)%nSubface = jj
1535  blocks(ii)%n1to1 = jj - blocks(ii)%nBocos
1536 
1537  end do subblockloop
1538  end do cgnsloop
1539 
1540  end subroutine determinecomputeblocks
1541 
1542  !========================================================================
1543 
1544  subroutine bcfacessubblock(cgnsID, ii, jj)
1545  !
1546  ! BCFacesSubblock determines the boundary subfaces of compute
1547  ! block ii, which is a subblock of the given cgns block.
1548  ! Jj is the counter for the number of subfaces.
1549  !
1550  use cgnsgrid
1551  use communication
1552  use inputphysics
1553  use partitionmod
1554  use utils, only: terminate
1555  use commonformats, only: strings
1556  implicit none
1557  !
1558  ! Subroutine arguments.
1559  !
1560  integer(kind=intType), intent(in) :: cgnsID, ii
1561  integer(kind=intType), intent(inout) :: jj
1562  !
1563  ! Local variables.
1564  !
1565  integer :: ierr
1566 
1567  integer(kind=intType) :: j, mm
1568  integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
1569 
1570  character(len=maxCGNSNameLen) :: zoneName, subName
1571  character(len=2*maxStringLen) :: errorMessage
1572 
1573  ! Loop over the physical boundaries of the original block.
1574 
1575  bocoloop: do j = 1, cgnsdoms(cgnsid)%nBocos
1576 
1577  ! Continue with the next boundary subface if this is a
1578  ! degenerated subface.
1579 
1580  if (.not. cgnsdoms(cgnsid)%bocoInfo(j)%actualFace) cycle
1581 
1582  ! Store the subface range a bit easier. Make sure that the
1583  ! indices run from low to high.
1584 
1585  ibeg = min(cgnsdoms(cgnsid)%bocoInfo(j)%iBeg, &
1586  cgnsdoms(cgnsid)%bocoInfo(j)%iEnd)
1587  iend = max(cgnsdoms(cgnsid)%bocoInfo(j)%iBeg, &
1588  cgnsdoms(cgnsid)%bocoInfo(j)%iEnd)
1589 
1590  jbeg = min(cgnsdoms(cgnsid)%bocoInfo(j)%jBeg, &
1591  cgnsdoms(cgnsid)%bocoInfo(j)%jEnd)
1592  jend = max(cgnsdoms(cgnsid)%bocoInfo(j)%jBeg, &
1593  cgnsdoms(cgnsid)%bocoInfo(j)%jEnd)
1594 
1595  kbeg = min(cgnsdoms(cgnsid)%bocoInfo(j)%kBeg, &
1596  cgnsdoms(cgnsid)%bocoInfo(j)%kEnd)
1597  kend = max(cgnsdoms(cgnsid)%bocoInfo(j)%kBeg, &
1598  cgnsdoms(cgnsid)%bocoInfo(j)%kEnd)
1599 
1600  ! Check for a possible overlap between the current boundary
1601  ! subface and subblock ii.
1602 
1603  overlap: if (ibeg <= blocks(ii)%iEndor .and. &
1604  iend >= blocks(ii)%iBegor .and. &
1605  jbeg <= blocks(ii)%jEndor .and. &
1606  jend >= blocks(ii)%jBegor .and. &
1607  kbeg <= blocks(ii)%kEndor .and. &
1608  kend >= blocks(ii)%kBegor) then
1609 
1610  ! Determine the overlap region between the current boundary
1611  ! face and subblock ii.
1612 
1613  ibeg = max(blocks(ii)%iBegor, ibeg)
1614  iend = min(blocks(ii)%iEndor, iend)
1615 
1616  jbeg = max(blocks(ii)%jBegor, jbeg)
1617  jend = min(blocks(ii)%jEndor, jend)
1618 
1619  kbeg = max(blocks(ii)%kBegor, kbeg)
1620  kend = min(blocks(ii)%kEndor, kend)
1621 
1622  ! Check the number of equal indices, which is stored in mm.
1623 
1624  mm = 0
1625  if (ibeg == iend) mm = mm + 1
1626  if (jbeg == jend) mm = mm + 1
1627  if (kbeg == kend) mm = mm + 1
1628 
1629  ! If no constant index is found something is wrong with the
1630  ! grid. Processor 0 prints an error message, while the
1631  ! others wait until they are killed.
1632 
1633  if (mm == 0) then
1634  if (myid == 0) then
1635  zonename = cgnsdoms(cgnsid)%zoneName
1636  subname = cgnsdoms(cgnsid)%bocoInfo(j)%bocoName
1637  write (errormessage, strings) "Zone ", trim(zonename), ", boundary subface ", trim(subname), &
1638  ": No constant index found for subface"
1639  call terminate("BCFacesSubblock", errormessage)
1640  end if
1641 
1642  call mpi_barrier(adflow_comm_world, ierr)
1643  end if
1644 
1645  ! Continue with the next subface if there is more than
1646  ! one constant index. This means that there is no overlap,
1647  ! but just an adjacency.
1648 
1649  if (mm > 1) cycle
1650 
1651  ! Update the counter jj and determine the range of the
1652  ! subface in the subblock.
1653 
1654  jj = jj + 1
1655 
1656  blocks(ii)%inBeg(jj) = ibeg - blocks(ii)%iBegor + 1
1657  blocks(ii)%inEnd(jj) = iend - blocks(ii)%iBegor + 1
1658 
1659  blocks(ii)%jnBeg(jj) = jbeg - blocks(ii)%jBegor + 1
1660  blocks(ii)%jnEnd(jj) = jend - blocks(ii)%jBegor + 1
1661 
1662  blocks(ii)%knBeg(jj) = kbeg - blocks(ii)%kBegor + 1
1663  blocks(ii)%knEnd(jj) = kend - blocks(ii)%kBegor + 1
1664 
1665  ! Determine the block face id on which this subface
1666  ! is located.
1667 
1668  if (ibeg == iend) then
1669 
1670  blocks(ii)%BCFaceID(jj) = imax
1671  if (ibeg == blocks(ii)%iBegor) blocks(ii)%BCFaceID(jj) = imin
1672 
1673  else if (jbeg == jend) then
1674 
1675  blocks(ii)%BCFaceID(jj) = jmax
1676  if (jbeg == blocks(ii)%jBegor) blocks(ii)%BCFaceID(jj) = jmin
1677 
1678  else
1679 
1680  blocks(ii)%BCFaceID(jj) = kmax
1681  if (kbeg == blocks(ii)%kBegor) blocks(ii)%BCFaceID(jj) = kmin
1682 
1683  end if
1684 
1685  ! Set some variables to 0, which are not relevant
1686  ! for boundary subfaces.
1687 
1688  blocks(ii)%dinBeg(jj) = 0; blocks(ii)%dinEnd(jj) = 0
1689  blocks(ii)%djnBeg(jj) = 0; blocks(ii)%djnEnd(jj) = 0
1690  blocks(ii)%dknBeg(jj) = 0; blocks(ii)%dknEnd(jj) = 0
1691 
1692  blocks(ii)%neighBlock(jj) = 0
1693 
1694  blocks(ii)%l1(jj) = 0
1695  blocks(ii)%l2(jj) = 0
1696  blocks(ii)%l3(jj) = 0
1697 
1698  ! Set the boundary condition and store to which original cgns
1699  ! subface this subface belongs.
1700 
1701  blocks(ii)%BCType(jj) = cgnsdoms(cgnsid)%bocoInfo(j)%BCType
1702 
1703  blocks(ii)%cgnsSubface(jj) = j
1704 
1705  ! Check whether this is a valid boundary condition for
1706  ! the current simulation.
1707 
1708  if (blocks(ii)%BCType(jj) == bcnotvalid) then
1709 
1710  ! To avoid a messy output only processor 0 calls
1711  ! terminate. The other processors will wait until
1712  ! they are killed.
1713 
1714  if (myid == 0) then
1715  zonename = cgnsdoms(cgnsid)%zoneName
1716  subname = cgnsdoms(cgnsid)%bocoInfo(j)%bocoName
1717 
1718  ! Check whether this is an internal or an external
1719  ! flow problem and create the error message
1720  ! accordingly.
1721 
1722  if (flowtype == internalflow) then
1723  write (errormessage, strings) "Zone ", trim(zonename), &
1724  ", boundary subface ", trim(subname), &
1725  ": Not a valid boundary condition for internal flow"
1726  else
1727  write (errormessage, strings) "Zone ", trim(zonename), &
1728  ", boundary subface ", trim(subname), &
1729  ": Not a valid boundary condition for external flow"
1730  end if
1731 
1732  call terminate("BCFacesSubblock", errormessage)
1733  end if
1734 
1735  call mpi_barrier(adflow_comm_world, ierr)
1736  end if
1737 
1738  ! Store the corresponding family a bit easier.
1739 
1740  mm = cgnsdoms(cgnsid)%bocoInfo(j)%familyID
1741 
1742  ! Check if this is either a sliding mesh interface or a
1743  ! bleed flow region. If so the group nummer is set to the
1744  ! sliding interface ID or bleed flow region ID respectivily.
1745  ! Otherwise the group nummer is the family nummer, which
1746  ! is 0 if the subface does not belong to a family.
1747 
1748  select case (blocks(ii)%BCType(jj))
1749  case (slidinginterface)
1750  blocks(ii)%groupNum(jj) = &
1751  cgnsdoms(cgnsid)%bocoInfo(j)%slidingID
1752 
1753  case (massbleedinflow, massbleedoutflow)
1754  blocks(ii)%groupNum(jj) = cgnsfamilies(mm)%bleedRegionID
1755 
1756  case default
1757  blocks(ii)%groupNum(jj) = mm
1758  end select
1759 
1760  end if overlap
1761 
1762  end do bocoloop
1763 
1764  end subroutine bcfacessubblock
1765 
1766  !========================================================================
1767 
1768  subroutine externalfacessubblock(cgnsID, ii, jj, nSubPerCGNS, &
1769  nAlloc, splitInfo)
1770  !
1771  ! externalFacesSubblock determines the block boundaries of
1772  ! the compute block ii which are located on the boundaries of
1773  ! the given original cgns block. As it is possible that due to
1774  ! a splitting of a neighboring block the number of block
1775  ! boundaries is larger than the original number, it must be
1776  ! checked whether enough memory has been allocated.
1777  ! jj is the counter for the number of subfaces.
1778  !
1779  use cgnsgrid
1780  use communication
1781  use partitionmod
1782  use utils, only: delta, terminate
1783  use commonformats, only: strings
1784  implicit none
1785  !
1786  ! Subroutine arguments
1787  !
1788  integer(kind=intType), intent(in) :: cgnsID, ii
1789  integer(kind=intType), intent(inout) :: jj, nAlloc
1790 
1791  integer(kind=intType), dimension(0:cgnsNDom), intent(in) :: &
1792  nSubPerCGNS
1793  type(splitcgnstype), dimension(cgnsNDom), intent(in) :: splitInfo
1794  !
1795  ! Local variables.
1796  !
1797  integer :: ierr
1798 
1799  integer(kind=intType) :: j, k, kk, mm
1800  integer(kind=intType) :: l1, L2, l3
1801  integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd, kBeg, kEnd
1802  integer(kind=intType) :: diBeg, diEnd, djBeg, djEnd
1803  integer(kind=intType) :: dkBeg, dkEnd
1804 
1805  integer(kind=intType), dimension(3, 3) :: trMat
1806 
1807  integer(kind=intType), dimension(:, :, :), pointer :: ranges
1808 
1809  character(len=maxCGNSNameLen) :: zoneName, subName
1810  character(len=2*maxStringLen) :: errorMessage
1811 
1812  logical :: diSwap, djSwap, dkSwap
1813  !
1814  ! Loop over the 1 to 1 block boundaries of the original block.
1815 
1816  n1to1loop: do j = 1, cgnsdoms(cgnsid)%n1to1
1817 
1818  ! Store the subface range a bit easier. Make sure that the
1819  ! indices run from low to high.
1820 
1821  ibeg = min(cgnsdoms(cgnsid)%conn1to1(j)%iBeg, &
1822  cgnsdoms(cgnsid)%conn1to1(j)%iEnd)
1823  iend = max(cgnsdoms(cgnsid)%conn1to1(j)%iBeg, &
1824  cgnsdoms(cgnsid)%conn1to1(j)%iEnd)
1825 
1826  jbeg = min(cgnsdoms(cgnsid)%conn1to1(j)%jBeg, &
1827  cgnsdoms(cgnsid)%conn1to1(j)%jEnd)
1828  jend = max(cgnsdoms(cgnsid)%conn1to1(j)%jBeg, &
1829  cgnsdoms(cgnsid)%conn1to1(j)%jEnd)
1830 
1831  kbeg = min(cgnsdoms(cgnsid)%conn1to1(j)%kBeg, &
1832  cgnsdoms(cgnsid)%conn1to1(j)%kEnd)
1833  kend = max(cgnsdoms(cgnsid)%conn1to1(j)%kBeg, &
1834  cgnsdoms(cgnsid)%conn1to1(j)%kEnd)
1835 
1836  ! Check for a possible overlap between the current boundary
1837  ! subface and subblock ii.
1838 
1839  overlap: if (ibeg <= blocks(ii)%iEndor .and. &
1840  iend >= blocks(ii)%iBegor .and. &
1841  jbeg <= blocks(ii)%jEndor .and. &
1842  jend >= blocks(ii)%jBegor .and. &
1843  kbeg <= blocks(ii)%kEndor .and. &
1844  kend >= blocks(ii)%kBegor) then
1845 
1846  ! Determine the overlap region between the current boundary
1847  ! face and subblock ii.
1848 
1849  ibeg = max(blocks(ii)%iBegor, ibeg)
1850  iend = min(blocks(ii)%iEndor, iend)
1851 
1852  jbeg = max(blocks(ii)%jBegor, jbeg)
1853  jend = min(blocks(ii)%jEndor, jend)
1854 
1855  kbeg = max(blocks(ii)%kBegor, kbeg)
1856  kend = min(blocks(ii)%kEndor, kend)
1857 
1858  ! Check the number of equal indices, which is stored in kk.
1859 
1860  kk = 0
1861  if (ibeg == iend) kk = kk + 1
1862  if (jbeg == jend) kk = kk + 1
1863  if (kbeg == kend) kk = kk + 1
1864 
1865  ! If no constant index is found something is wrong with the
1866  ! grid. Processor 0 prints an error message, while the
1867  ! others wait until they are killed.
1868 
1869  if (kk == 0) then
1870  if (myid == 0) then
1871  zonename = cgnsdoms(cgnsid)%zoneName
1872  subname = cgnsdoms(cgnsid)%bocoInfo(j)%bocoName
1873  write (errormessage, strings) "Zone ", trim(zonename), &
1874  ", 1 to 1 block connectivity ", trim(subname), &
1875  ": No constant index found for subface"
1876  call terminate("externalFacesSubblock", errormessage)
1877  end if
1878 
1879  call mpi_barrier(adflow_comm_world, ierr)
1880  end if
1881 
1882  ! Continue with the next subface if there is more than
1883  ! one constant index. This means that there is no overlap,
1884  ! but just an adjacency.
1885 
1886  if (kk > 1) cycle
1887 
1888  ! Preserve negative running indices of the subface.
1889 
1890  if (cgnsdoms(cgnsid)%conn1to1(j)%iEnd < &
1891  cgnsdoms(cgnsid)%conn1to1(j)%iBeg) then
1892  mm = ibeg
1893  ibeg = iend
1894  iend = mm
1895  end if
1896 
1897  if (cgnsdoms(cgnsid)%conn1to1(j)%jEnd < &
1898  cgnsdoms(cgnsid)%conn1to1(j)%jBeg) then
1899  mm = jbeg
1900  jbeg = jend
1901  jend = mm
1902  end if
1903 
1904  if (cgnsdoms(cgnsid)%conn1to1(j)%kEnd < &
1905  cgnsdoms(cgnsid)%conn1to1(j)%kBeg) then
1906  mm = kbeg
1907  kbeg = kend
1908  kend = mm
1909  end if
1910 
1911  ! Determine the transformation matrix between the
1912  ! current subface and the donor subface.
1913 
1914  l1 = cgnsdoms(cgnsid)%conn1to1(j)%l1
1915  l2 = cgnsdoms(cgnsid)%conn1to1(j)%l2
1916  l3 = cgnsdoms(cgnsid)%conn1to1(j)%l3
1917 
1918  trmat(1, 1) = sign(1_inttype, l1) * delta(l1, 1_inttype)
1919  trmat(2, 1) = sign(1_inttype, l1) * delta(l1, 2_inttype)
1920  trmat(3, 1) = sign(1_inttype, l1) * delta(l1, 3_inttype)
1921 
1922  trmat(1, 2) = sign(1_inttype, l2) * delta(l2, 1_inttype)
1923  trmat(2, 2) = sign(1_inttype, l2) * delta(l2, 2_inttype)
1924  trmat(3, 2) = sign(1_inttype, l2) * delta(l2, 3_inttype)
1925 
1926  trmat(1, 3) = sign(1_inttype, l3) * delta(l3, 1_inttype)
1927  trmat(2, 3) = sign(1_inttype, l3) * delta(l3, 2_inttype)
1928  trmat(3, 3) = sign(1_inttype, l3) * delta(l3, 3_inttype)
1929 
1930  ! Determine the corresponding donor range of the subface
1931  ! iBeg, iEnd; jBeg, jEnd; kBeg, kEnd.
1932 
1933  l1 = ibeg - cgnsdoms(cgnsid)%conn1to1(j)%iBeg
1934  l2 = jbeg - cgnsdoms(cgnsid)%conn1to1(j)%jBeg
1935  l3 = kbeg - cgnsdoms(cgnsid)%conn1to1(j)%kBeg
1936 
1937  dibeg = cgnsdoms(cgnsid)%conn1to1(j)%diBeg &
1938  + trmat(1, 1) * l1 + trmat(1, 2) * l2 + trmat(1, 3) * l3
1939  djbeg = cgnsdoms(cgnsid)%conn1to1(j)%djBeg &
1940  + trmat(2, 1) * l1 + trmat(2, 2) * l2 + trmat(2, 3) * l3
1941  dkbeg = cgnsdoms(cgnsid)%conn1to1(j)%dkBeg &
1942  + trmat(3, 1) * l1 + trmat(3, 2) * l2 + trmat(3, 3) * l3
1943 
1944  l1 = iend - cgnsdoms(cgnsid)%conn1to1(j)%iBeg
1945  l2 = jend - cgnsdoms(cgnsid)%conn1to1(j)%jBeg
1946  l3 = kend - cgnsdoms(cgnsid)%conn1to1(j)%kBeg
1947 
1948  diend = cgnsdoms(cgnsid)%conn1to1(j)%diBeg &
1949  + trmat(1, 1) * l1 + trmat(1, 2) * l2 + trmat(1, 3) * l3
1950  djend = cgnsdoms(cgnsid)%conn1to1(j)%djBeg &
1951  + trmat(2, 1) * l1 + trmat(2, 2) * l2 + trmat(2, 3) * l3
1952  dkend = cgnsdoms(cgnsid)%conn1to1(j)%dkBeg &
1953  + trmat(3, 1) * l1 + trmat(3, 2) * l2 + trmat(3, 3) * l3
1954 
1955  ! Make sure that the donor indices are positive running
1956  ! indices. If they must be swapped, the corresponding
1957  ! logical is set to .true.
1958 
1959  diswap = .false.
1960  if (dibeg > diend) then
1961  mm = dibeg; dibeg = diend; diend = mm; diswap = .true.
1962  end if
1963 
1964  djswap = .false.
1965  if (djbeg > djend) then
1966  mm = djbeg; djbeg = djend; djend = mm; djswap = .true.
1967  end if
1968 
1969  dkswap = .false.
1970  if (dkbeg > dkend) then
1971  mm = dkbeg; dkbeg = dkend; dkend = mm; dkswap = .true.
1972  end if
1973 
1974  ! Store the index of the donor block a bit easier and loop
1975  ! over its subblocks to find the donor range.
1976 
1977  mm = cgnsdoms(cgnsid)%conn1to1(j)%donorBlock
1978  ranges => splitinfo(mm)%ranges
1979 
1980  donorloop: do k = 1, splitinfo(mm)%nSubblocks
1981 
1982  ! Check whether this subblock and the given donor range
1983  ! overlap.
1984 
1985  donoroverlap: if (dibeg <= ranges(k, 1, 2) .and. &
1986  diend >= ranges(k, 1, 1) .and. &
1987  djbeg <= ranges(k, 2, 2) .and. &
1988  djend >= ranges(k, 2, 1) .and. &
1989  dkbeg <= ranges(k, 3, 2) .and. &
1990  dkend >= ranges(k, 3, 1)) then
1991 
1992  ! Determine the range of the donor face, which is
1993  ! stored in iBeg, iEnd, etc.
1994 
1995  ibeg = max(dibeg, ranges(k, 1, 1))
1996  iend = min(diend, ranges(k, 1, 2))
1997 
1998  jbeg = max(djbeg, ranges(k, 2, 1))
1999  jend = min(djend, ranges(k, 2, 2))
2000 
2001  kbeg = max(dkbeg, ranges(k, 3, 1))
2002  kend = min(dkend, ranges(k, 3, 2))
2003 
2004  ! Check whether the subfaces are truely overlapping
2005  ! or just adjacent.
2006 
2007  kk = 0
2008  if (ibeg == iend) kk = kk + 1
2009  if (jbeg == jend) kk = kk + 1
2010  if (kbeg == kend) kk = kk + 1
2011 
2012  if (kk > 1) cycle
2013 
2014  ! Update the counter jj and check whether enough memory
2015  ! has been allocated. If not, reallocate.
2016 
2017  jj = jj + 1
2018  if (jj > nalloc) call reallocsubfacememory(ii, nalloc)
2019 
2020  ! Set some info for this subface, which can be
2021  ! determined relatively easily.
2022 
2023  blocks(ii)%BCType(jj) = b2bmatch
2024  blocks(ii)%cgnsSubface(jj) = j
2025  blocks(ii)%groupNum(jj) = 0
2026 
2027  blocks(ii)%l1(jj) = cgnsdoms(cgnsid)%conn1to1(j)%l1
2028  blocks(ii)%l2(jj) = cgnsdoms(cgnsid)%conn1to1(j)%l2
2029  blocks(ii)%l3(jj) = cgnsdoms(cgnsid)%conn1to1(j)%l3
2030 
2031  ! Determine the neighboring block id.
2032 
2033  blocks(ii)%neighBlock(jj) = nsubpercgns(mm - 1) + k
2034 
2035  ! Determine the range of the donor. First switch the
2036  ! indices if the original indices were swapped.
2037 
2038  if (diswap) then
2039  kk = ibeg; ibeg = iend; iend = kk
2040  end if
2041 
2042  if (djswap) then
2043  kk = jbeg; jbeg = jend; jend = kk
2044  end if
2045 
2046  if (dkswap) then
2047  kk = kbeg; kbeg = kend; kend = kk
2048  end if
2049 
2050  ! Determine the local range of the donor, i.e. the
2051  ! offset in the original block must be substracted.
2052 
2053  blocks(ii)%dinBeg(jj) = ibeg - ranges(k, 1, 1) + 1
2054  blocks(ii)%djnBeg(jj) = jbeg - ranges(k, 2, 1) + 1
2055  blocks(ii)%dknBeg(jj) = kbeg - ranges(k, 3, 1) + 1
2056 
2057  blocks(ii)%dinEnd(jj) = iend - ranges(k, 1, 1) + 1
2058  blocks(ii)%djnEnd(jj) = jend - ranges(k, 2, 1) + 1
2059  blocks(ii)%dknEnd(jj) = kend - ranges(k, 3, 1) + 1
2060 
2061  ! Transform the donor range in the original donor block
2062  ! back the a subface range in the original cgns block.
2063  ! The inverse of the transformation matrix trMat is
2064  ! the transpose.
2065 
2066  l1 = ibeg - cgnsdoms(cgnsid)%conn1to1(j)%diBeg
2067  l2 = jbeg - cgnsdoms(cgnsid)%conn1to1(j)%djBeg
2068  l3 = kbeg - cgnsdoms(cgnsid)%conn1to1(j)%dkBeg
2069 
2070  ibeg = cgnsdoms(cgnsid)%conn1to1(j)%iBeg &
2071  + trmat(1, 1) * l1 + trmat(2, 1) * l2 + trmat(3, 1) * l3
2072  jbeg = cgnsdoms(cgnsid)%conn1to1(j)%jBeg &
2073  + trmat(1, 2) * l1 + trmat(2, 2) * l2 + trmat(3, 2) * l3
2074  kbeg = cgnsdoms(cgnsid)%conn1to1(j)%kBeg &
2075  + trmat(1, 3) * l1 + trmat(2, 3) * l2 + trmat(3, 3) * l3
2076 
2077  l1 = iend - cgnsdoms(cgnsid)%conn1to1(j)%diBeg
2078  l2 = jend - cgnsdoms(cgnsid)%conn1to1(j)%djBeg
2079  l3 = kend - cgnsdoms(cgnsid)%conn1to1(j)%dkBeg
2080 
2081  iend = cgnsdoms(cgnsid)%conn1to1(j)%iBeg &
2082  + trmat(1, 1) * l1 + trmat(2, 1) * l2 + trmat(3, 1) * l3
2083  jend = cgnsdoms(cgnsid)%conn1to1(j)%jBeg &
2084  + trmat(1, 2) * l1 + trmat(2, 2) * l2 + trmat(3, 2) * l3
2085  kend = cgnsdoms(cgnsid)%conn1to1(j)%kBeg &
2086  + trmat(1, 3) * l1 + trmat(2, 3) * l2 + trmat(3, 3) * l3
2087 
2088  ! Store the subface range of the new block, i.e.
2089  ! An offset must be subtracted.
2090 
2091  blocks(ii)%inBeg(jj) = ibeg - blocks(ii)%iBegor + 1
2092  blocks(ii)%jnBeg(jj) = jbeg - blocks(ii)%jBegor + 1
2093  blocks(ii)%knBeg(jj) = kbeg - blocks(ii)%kBegor + 1
2094 
2095  blocks(ii)%inEnd(jj) = iend - blocks(ii)%iBegor + 1
2096  blocks(ii)%jnEnd(jj) = jend - blocks(ii)%jBegor + 1
2097  blocks(ii)%knEnd(jj) = kend - blocks(ii)%kBegor + 1
2098 
2099  ! Determine the block face id on which this subface
2100  ! is located.
2101 
2102  if (ibeg == iend) then
2103 
2104  blocks(ii)%BCFaceID(jj) = imax
2105  if (ibeg == blocks(ii)%iBegor) &
2106  blocks(ii)%BCFaceID(jj) = imin
2107 
2108  else if (jbeg == jend) then
2109 
2110  blocks(ii)%BCFaceID(jj) = jmax
2111  if (jbeg == blocks(ii)%jBegor) &
2112  blocks(ii)%BCFaceID(jj) = jmin
2113 
2114  else
2115 
2116  blocks(ii)%BCFaceID(jj) = kmax
2117  if (kbeg == blocks(ii)%kBegor) &
2118  blocks(ii)%BCFaceID(jj) = kmin
2119 
2120  end if
2121 
2122  end if donoroverlap
2123  end do donorloop
2124  end if overlap
2125  end do n1to1loop
2126 
2127  end subroutine externalfacessubblock
2128 
2129  !========================================================================
2130 
2131  subroutine internalfacessubblock(cgnsID, ii, jj, nSubPerCGNS, &
2132  nAlloc, splitInfo)
2133  !
2134  ! internalFacesSubblock determines the block boundaries of
2135  ! the compute block ii which are created due to the splitting of
2136  ! the original block into subblock. As the number of these
2137  ! internal boundaries is not known, it must be checked whether
2138  ! enough memory has been allocated. jj is the counter for the
2139  ! number of subfaces.
2140  !
2141  use cgnsgrid
2142  use communication
2143  use partitionmod
2144  implicit none
2145  !
2146  ! Subroutine arguments
2147  !
2148  integer(kind=intType), intent(in) :: cgnsID, ii
2149  integer(kind=intType), intent(inout) :: jj, nAlloc
2150 
2151  integer(kind=intType), dimension(0:cgnsNDom), intent(in) :: &
2152  nSubPerCGNS
2153  type(splitcgnstype), intent(in) :: splitInfo
2154  !
2155  ! Local variables.
2156  !
2157  integer(kind=intType) :: indFace, jBeg, jEnd, kBeg, kEnd
2158  integer(kind=intType) :: i, i2, j, k, faceID
2159 
2160  ! iMin face.
2161 
2162  if (blocks(ii)%iBegor > 1) then
2163 
2164  ! Imin face is created through splitting. Set some variables
2165  ! for the general treatment.
2166 
2167  indface = blocks(ii)%iBegor
2168  jbeg = blocks(ii)%jBegor
2169  jend = blocks(ii)%jEndor
2170  kbeg = blocks(ii)%kBegor
2171  kend = blocks(ii)%kEndor
2172 
2173  i = 1; j = 2; k = 3
2174  i2 = 2; faceid = imin
2175 
2176  ! Search for neighbors in the subblocks of the given
2177  ! cgns blocks.
2178 
2179  call searchinternalneighbors
2180 
2181  end if
2182 
2183  ! iMax face.
2184 
2185  if (blocks(ii)%iEndor < cgnsdoms(cgnsid)%il) then
2186 
2187  ! Imax face is created through splitting. Set some variables
2188  ! for the general treatment.
2189 
2190  indface = blocks(ii)%iEndor
2191  jbeg = blocks(ii)%jBegor
2192  jend = blocks(ii)%jEndor
2193  kbeg = blocks(ii)%kBegor
2194  kend = blocks(ii)%kEndor
2195 
2196  i = 1; j = 2; k = 3
2197  i2 = 1; faceid = imax
2198 
2199  ! Search for neighbors in the subblocks of the given
2200  ! cgns blocks.
2201 
2202  call searchinternalneighbors
2203 
2204  end if
2205 
2206  ! jMin face.
2207 
2208  if (blocks(ii)%jBegor > 1) then
2209 
2210  ! Jmin face is created through splitting. Set some variables
2211  ! for the general treatment.
2212 
2213  indface = blocks(ii)%jBegor
2214  jbeg = blocks(ii)%iBegor
2215  jend = blocks(ii)%iEndor
2216  kbeg = blocks(ii)%kBegor
2217  kend = blocks(ii)%kEndor
2218 
2219  i = 2; j = 1; k = 3
2220  i2 = 2; faceid = jmin
2221 
2222  ! Search for neighbors in the subblocks of the given
2223  ! cgns blocks.
2224 
2225  call searchinternalneighbors
2226 
2227  end if
2228 
2229  ! jMax face.
2230 
2231  if (blocks(ii)%jEndor < cgnsdoms(cgnsid)%jl) then
2232 
2233  ! Jmax face is created through splitting. Set some variables
2234  ! for the general treatment.
2235 
2236  indface = blocks(ii)%jEndor
2237  jbeg = blocks(ii)%iBegor
2238  jend = blocks(ii)%iEndor
2239  kbeg = blocks(ii)%kBegor
2240  kend = blocks(ii)%kEndor
2241 
2242  i = 2; j = 1; k = 3
2243  i2 = 1; faceid = jmax
2244 
2245  ! Search for neighbors in the subblocks of the given
2246  ! cgns blocks.
2247 
2248  call searchinternalneighbors
2249 
2250  end if
2251 
2252  ! kMin face.
2253 
2254  if (blocks(ii)%kBegor > 1) then
2255 
2256  ! Kmin face is created through splitting. Set some variables
2257  ! for the general treatment.
2258 
2259  indface = blocks(ii)%kBegor
2260  jbeg = blocks(ii)%iBegor
2261  jend = blocks(ii)%iEndor
2262  kbeg = blocks(ii)%jBegor
2263  kend = blocks(ii)%jEndor
2264 
2265  i = 3; j = 1; k = 2
2266  i2 = 2; faceid = kmin
2267 
2268  ! Search for neighbors in the subblocks of the given
2269  ! cgns blocks.
2270 
2271  call searchinternalneighbors
2272 
2273  end if
2274 
2275  ! kMax face.
2276 
2277  if (blocks(ii)%kEndor < cgnsdoms(cgnsid)%kl) then
2278 
2279  ! Kmax face is created through splitting. Set some variables
2280  ! for the general treatment.
2281 
2282  indface = blocks(ii)%kEndor
2283  jbeg = blocks(ii)%iBegor
2284  jend = blocks(ii)%iEndor
2285  kbeg = blocks(ii)%jBegor
2286  kend = blocks(ii)%jEndor
2287 
2288  i = 3; j = 1; k = 2
2289  i2 = 1; faceid = kmax
2290 
2291  ! Search for neighbors in the subblocks of the given
2292  ! cgns blocks.
2293 
2294  call searchinternalneighbors
2295 
2296  end if
2297 
2298  !=================================================================
2299 
2300  contains
2301 
2302  !===============================================================
2303 
2304  subroutine searchinternalneighbors
2305  !
2306  ! searchInternalNeighbors determines block faces created by
2307  ! by the splitting of the original block. The variables set in
2308  ! internalFacesSubblock are used such that a general
2309  ! treatment is possible.
2310  !
2311  implicit none
2312  !
2313  ! Local variables
2314  !
2315  integer(kind=intType) :: mm, jnBeg, jnEnd, knBeg, knEnd
2316 
2317  integer(kind=intType), dimension(3, 2) :: subRange
2318 
2319  integer(kind=intType), dimension(:, :, :), pointer :: ranges
2320 
2321  ! Set the pointer for ranges to make the code more readable.
2322 
2323  ranges => splitinfo%ranges
2324 
2325  ! Loop over the number of blocks into the original block
2326  ! is split.
2327 
2328  do mm = 1, splitinfo%nSubblocks
2329 
2330  ! Check whether the constant index of the face matches
2331  ! the given index of subblock.
2332 
2333  if (ranges(mm, i, i2) == indface) then
2334 
2335  ! Check whether the faces overlap.
2336 
2337  if (jbeg <= ranges(mm, j, 2) .and. &
2338  jend >= ranges(mm, j, 1) .and. &
2339  kbeg <= ranges(mm, k, 2) .and. &
2340  kend >= ranges(mm, k, 1)) then
2341 
2342  ! There is a possible overlap. Determine the nodal
2343  ! range of the subface.
2344 
2345  jnbeg = max(jbeg, ranges(mm, j, 1))
2346  jnend = min(jend, ranges(mm, j, 2))
2347 
2348  knbeg = max(kbeg, ranges(mm, k, 1))
2349  knend = min(kend, ranges(mm, k, 2))
2350 
2351  ! Check whether this is a true subface.
2352 
2353  if (jnend > jnbeg .and. knend > knbeg) then
2354 
2355  ! An overlap occurs. Update the counter jj and check
2356  ! whether enough memory has been allocated.
2357  ! If not, reallocate.
2358 
2359  jj = jj + 1
2360  if (jj > nalloc) call reallocsubfacememory(ii, nalloc)
2361 
2362  ! Set the information of the BCType, BCFaceID,
2363  ! cgnsSubface and the group number. As this face is
2364  ! created internally the latter two variables are
2365  ! set to 0.
2366 
2367  blocks(ii)%BCType(jj) = b2bmatch
2368  blocks(ii)%BCFaceID(jj) = faceid
2369  blocks(ii)%cgnsSubface(jj) = 0
2370  blocks(ii)%groupNum(jj) = 0
2371 
2372  ! Determine the subRange of the subface in the
2373  ! original block.
2374 
2375  subrange(i, 1) = indface; subrange(i, 2) = indface
2376  subrange(j, 1) = jnbeg; subrange(j, 2) = jnend
2377  subrange(k, 1) = knbeg; subrange(k, 2) = knend
2378 
2379  ! Determine the nodal range in the current subblock.
2380 
2381  blocks(ii)%inBeg(jj) = subrange(1, 1) &
2382  - blocks(ii)%iBegor + 1
2383  blocks(ii)%inEnd(jj) = subrange(1, 2) &
2384  - blocks(ii)%iBegor + 1
2385 
2386  blocks(ii)%jnBeg(jj) = subrange(2, 1) &
2387  - blocks(ii)%jBegor + 1
2388  blocks(ii)%jnEnd(jj) = subrange(2, 2) &
2389  - blocks(ii)%jBegor + 1
2390 
2391  blocks(ii)%knBeg(jj) = subrange(3, 1) &
2392  - blocks(ii)%kBegor + 1
2393  blocks(ii)%knEnd(jj) = subrange(3, 2) &
2394  - blocks(ii)%kBegor + 1
2395 
2396  ! Determine the nodal range in the donor block.
2397 
2398  blocks(ii)%dinBeg(jj) = subrange(1, 1) &
2399  - ranges(mm, 1, 1) + 1
2400  blocks(ii)%dinEnd(jj) = subrange(1, 2) &
2401  - ranges(mm, 1, 1) + 1
2402 
2403  blocks(ii)%djnBeg(jj) = subrange(2, 1) &
2404  - ranges(mm, 2, 1) + 1
2405  blocks(ii)%djnEnd(jj) = subrange(2, 2) &
2406  - ranges(mm, 2, 1) + 1
2407 
2408  blocks(ii)%dknBeg(jj) = subrange(3, 1) &
2409  - ranges(mm, 3, 1) + 1
2410  blocks(ii)%dknEnd(jj) = subrange(3, 2) &
2411  - ranges(mm, 3, 1) + 1
2412 
2413  ! Set the neighboring block to mm plus the offset
2414  ! for the current cgns block and set the transformation
2415  ! matrix. The latter is simply 1-2-3, because the
2416  ! orientation of the subblocks is identical to the
2417  ! original block.
2418 
2419  blocks(ii)%neighBlock(jj) = mm + nsubpercgns(cgnsid - 1)
2420 
2421  blocks(ii)%l1(jj) = 1
2422  blocks(ii)%l2(jj) = 2
2423  blocks(ii)%l3(jj) = 3
2424 
2425  end if
2426  end if
2427  end if
2428  end do
2429 
2430  end subroutine searchinternalneighbors
2431 
2432  end subroutine internalfacessubblock
2433  subroutine graphpartitioning(emptyPartitions, commNeglected)
2434  !
2435  ! graphPartitioning partitions the corresponding graph of the
2436  ! computational blocks such that both the number of cells and
2437  ! number of faces is about equal on all processors.
2438  !
2439  use constants
2440  use communication, only: myid, adflow_comm_world, nproc
2441  use partitionmod, only: part, blocks, ubvec, nblocks
2442  use inputparallel, only: loadimbalance
2443  use utils, only: terminate
2444  use sorting, only: qsortintegers, bsearchintegers
2445  implicit none
2446  !
2447  ! Subroutine arguments.
2448  !
2449  logical, intent(out) :: emptyPartitions, commNeglected
2450  !
2451  ! Variables to store the graph in metis format.
2452  !
2453  ! nVertex : Number of vertices in the graph, equals nBlocks.
2454  ! nCon : Number of contraints, 2.
2455  ! xadj(0:nVertex): Number of edges per vertex, cumulative storage
2456  ! format.
2457  ! adjncy(:) : End vertex of the edge; the size of adjncy is
2458  ! xadj(nVertex).
2459  ! vwgt(:,:) : Vertex weights, size equals nCon,nVertex. The
2460  ! vertex weights are stored contiguously.
2461  ! adjwgt(:) : Edge weights, size equals xadj(nVertex). Note
2462  ! that the edge weights of edge i-j can be
2463  ! different from the weight of edge j-i.
2464  ! wgtflag : Whether or not to use weights on edges. Here
2465  ! wgtflag should always be 1 to indicate that
2466  ! edge weights are used.
2467  ! numflag : Flag to indicate the numbering convention,
2468  ! starting from 0 or 1. Here we start from 0.
2469  ! nParts : Number of parts to split the graph. This is
2470  ! nProc.
2471  ! ubvec(2) : Tolerance for the constraints. Stored in the
2472  ! module dpartitionMod.
2473  ! options(5) : Option array; normally the default is used
2474  ! indicated by options(1) = 0.
2475  ! edgecut : On return it contains the edge cut of the
2476  ! distributed graph.
2477  ! part(nVertex) : On return the processor ID for each block.
2478  ! It will be returned in fortran numbering,
2479  ! i.e. starting at 1. Stored in the module
2480  ! distributionMod.
2481 
2482  integer :: nVertex, nCon, wgtflag, numflag, nParts, edgecut
2483  integer, dimension(5) :: options
2484 
2485  integer(kind=intType), dimension(:), allocatable :: xadj, adjncy
2486  integer(kind=intType), dimension(:), allocatable :: adjwgt
2487  integer(kind=intType), dimension(:, :), allocatable :: vwgt
2488  !
2489  ! Local variables.
2490  !
2491  integer :: ierr
2492 
2493  integer(kind=intType) :: i, j
2494  integer(kind=intType) :: nEdges, nEdgesMax, ii, jj, kk
2495 
2496  integer(kind=intType), dimension(0:nProc - 1) :: nBlockPerProc
2497 
2498  integer(kind=intType), dimension(:), allocatable :: tmp
2499 
2500  integer(kind=8) :: nCellsTotal ! 8 byte integers to avoid
2501  integer(kind=8) :: nFacesTotal ! overflow.
2502 
2503  real(kind=4) :: ubvec_temp(2) ! Explict float for ubvec
2504  !
2505  ! Check whether part is allocated from a previous call. If so,
2506  ! release the memory.
2507 
2508  if (allocated(part)) then
2509  deallocate (part, stat=ierr)
2510  if (ierr /= 0) &
2511  call terminate("graphPartitioning", &
2512  "Deallocation failure for part")
2513  end if
2514 
2515  ! Determine the number of edges in the graph and the maximum
2516  ! number for a vertex in the graph.
2517 
2518  nedges = 0
2519  nedgesmax = 0
2520  do i = 1, nblocks
2521  ii = blocks(i)%n1to1
2522  nedges = nedges + ii
2523  nedgesmax = max(nedgesmax, ii)
2524  end do
2525 
2526  ! Initialize some values for the graph.
2527 
2528  nvertex = nblocks
2529  ncon = 2
2530  wgtflag = 1
2531  numflag = 0
2532  nparts = nproc
2533  ubvec(1) = one + loadimbalance
2534  ubvec(2) = one + loadimbalance
2535  options = 0
2536 
2537  ! Allocate the memory to store/build the graph.
2538 
2539  allocate (xadj(0:nvertex), vwgt(ncon, nvertex), adjncy(nedges), &
2540  adjwgt(nedges), part(nvertex), tmp(nedgesmax), stat=ierr)
2541  if (ierr /= 0) &
2542  call terminate("graphPartitioning", &
2543  "Memory allocation failure for graph variables")
2544 
2545  ! Initialize xadj(0) to 0.
2546  ! Furthermore initialize adjwgt to 0, as these values are
2547  ! accumulated due to multiple subfaces between blocks.
2548 
2549  xadj(0) = 0
2550  adjwgt = 0
2551 
2552  ! Loop over the number of blocks to build the graph.
2553 
2554  graphvertex: do i = 1, nblocks
2555 
2556  ! Store the both vertex weights.
2557 
2558  vwgt(1, i) = blocks(i)%nCell
2559  vwgt(2, i) = blocks(i)%nFace
2560 
2561  ! Sort the neighbors in increasing order and neglect the
2562  ! communication to myself, i.e. do not allow an edge to myself.
2563  ! The sorting is necessary, because a block might have several
2564  ! subfaces with another block
2565 
2566  nedges = 0
2567  do j = 1, blocks(i)%n1to1
2568  ii = blocks(i)%nBocos + j
2569  if (blocks(i)%neighBlock(ii) /= i) then
2570  nedges = nedges + 1
2571  tmp(nedges) = blocks(i)%neighBlock(ii)
2572  end if
2573  end do
2574 
2575  ! Sort tmp in increasing order and get rid of the possible
2576  ! multiple entries.
2577 
2578  call qsortintegers(tmp, nedges)
2579 
2580  ii = min(nedges, 1_inttype) ! Be aware of nEdges == 0
2581  do j = 2, nedges
2582  if (tmp(j) /= tmp(ii)) then
2583  ii = ii + 1
2584  tmp(ii) = tmp(j)
2585  end if
2586  end do
2587 
2588  ! Set nEdges to ii and update xadj(i).
2589 
2590  nedges = ii
2591  xadj(i) = xadj(i - 1) + nedges
2592 
2593  ! Repeat the loop over the subfaces, but now store
2594  ! the edge info.
2595 
2596  edges1to1: do j = 1, blocks(i)%n1to1
2597 
2598  ii = blocks(i)%nBocos + j
2599  if (blocks(i)%neighBlock(ii) /= i) then
2600 
2601  ! Search for the block ID and add the offset of xadj(i-1)
2602  ! to obtain the correct index to store the edge info.
2603  ! The -1 to the adjncy is present because C-numbering
2604  ! is used when calling Metis.
2605 
2606  jj = xadj(i - 1) &
2607  + bsearchintegers(blocks(i)%neighBlock(ii), tmp(1:nedges))
2608  adjncy(jj) = blocks(i)%neighBlock(ii) - 1
2609 
2610  ! The weight equals the number of 1st and 2nd level halo
2611  ! cells to be communicated between the blocks. The weights
2612  ! are accumulated, as multiple subfaces between blocks are
2613  ! possible.
2614 
2615  kk = 1
2616  adjwgt(jj) = adjwgt(jj) + 2 * &
2617  (max(abs(blocks(i)%inEnd(ii) - blocks(i)%inBeg(ii)), kk) &
2618  * max(abs(blocks(i)%jnEnd(ii) - blocks(i)%jnBeg(ii)), kk) &
2619  * max(abs(blocks(i)%knEnd(ii) - blocks(i)%knBeg(ii)), kk))
2620  end if
2621 
2622  end do edges1to1
2623 
2624  end do graphvertex
2625 
2626  ! Metis has problems when the total number of cells or faces
2627  ! used in the weights exceeds 2Gb. Therefore the sum of these
2628  ! values is determined and an appropriate weight factor is
2629  ! determined. Note that the type of nCellsTotal and nFacesTotal
2630  ! is integer*8.
2631 
2632  ncellstotal = 0
2633  nfacestotal = 0
2634  do i = 1, nblocks
2635  ncellstotal = ncellstotal + vwgt(1, i)
2636  nfacestotal = nfacestotal + vwgt(2, i)
2637  end do
2638 
2639  if (ncellstotal > 2147483647 .or. nfacestotal > 2147483647) then
2640  ncellstotal = ncellstotal / 2147483647 + 1
2641  nfacestotal = nfacestotal / 2147483647 + 1
2642 
2643  do i = 1, nblocks
2644  vwgt(1, i) = vwgt(1, i) / ncellstotal
2645  vwgt(2, i) = vwgt(2, i) / nfacestotal
2646  end do
2647  end if
2648 
2649  ! Loop over the number of attempts to partition the graph.
2650  ! In the first attempt the communication is taken into account.
2651  ! If not successful, i.e. empty partitions present, the metis
2652  ! routine is called once more, but now with zero adjwgt. This
2653  ! means that the communication cost is neglected and metis
2654  ! normally gives a valid partitioning.
2655  ! Initialize commNeglected to .false. This will change if in
2656  ! the loop below the first call to metis is not successful.
2657 
2658  commneglected = .false.
2659  attemptloop: do ii = 1, 2
2660 
2661  ! Copy ubvec to a float since if you use -r8 flag it gets
2662  ! converted to real
2663  ubvec_temp(1) = real(ubvec(1))
2664  ubvec_temp(2) = real(ubvec(2))
2665 
2666  call metisinterface(nvertex, ncon, xadj, adjncy, vwgt, &
2667  adjwgt, wgtflag, numflag, nparts, &
2668  ubvec_temp, options, edgecut, part)
2669  ! Determine the number of blocks per processor.
2670 
2671  nblockperproc = 0
2672  do i = 1, nblocks
2673  nblockperproc(part(i)) = nblockperproc(part(i)) + 1
2674  end do
2675 
2676  ! Check for empty partitions.
2677 
2678  emptypartitions = .false.
2679  do i = 0, nproc - 1
2680  if (nblockperproc(i) == 0) emptypartitions = .true.
2681  end do
2682 
2683  ! Exit the loop if no empty partitions are present or if
2684  ! this is the second time this loop is executed.
2685 
2686  if (ii == 2 .or. (.not. emptypartitions)) exit attemptloop
2687 
2688  ! The first call to metis resulted in empty partitions.
2689  ! Ignore the communication, i.e. set the number of
2690  ! neighbors to 0, and try again.
2691 
2692  commneglected = .true.
2693  xadj = 0
2694 
2695  end do attemptloop
2696 
2697  ! Deallocate the memory for the graph except part.
2698 
2699  deallocate (xadj, vwgt, adjncy, adjwgt, tmp, stat=ierr)
2700  if (ierr /= 0) &
2701  call terminate("graphPartitioning", &
2702  "Deallocation failure for graph variables")
2703 
2704  end subroutine graphpartitioning
2705 
2706  subroutine checkloadbalance(cellsBalanced, facesBalanced)
2707  !
2708  ! checkLoadBalance determines whether or not the load balance
2709  ! for the cells and faces is met.
2710  !
2711  use constants
2712  use communication, only: adflow_comm_world, myid, nproc
2713  use partitionmod, only: blocks, ubvec, part, nblocks
2714  implicit none
2715  !
2716  ! Subroutine arguments.
2717  !
2718  logical, intent(out) :: cellsBalanced, facesBalanced
2719  !
2720  ! Local variables
2721  !
2722  integer(kind=intType) :: i, j
2723  integer(kind=intType) :: nCellMax, nCellTol
2724  integer(kind=intType) :: nFaceMax, nFaceTol
2725 
2726  integer(kind=intType), dimension(nProc) :: nCell, nFace
2727 
2728  integer(kind=8) :: nCellsEven, nFacesEven ! 8 byte integers to
2729  ! avoid overflow.
2730 
2731  ! Initialize nCell and nFace to 0. These variables will contain
2732  ! the number of cells and faces per partition (== processor)
2733  ! respectively.
2734 
2735  ncell = 0
2736  nface = 0
2737 
2738  ! Determine the number of cells and faces per partition.
2739  ! Note that part(i) is the processor id, which starts at 0.
2740 
2741  do i = 1, nblocks
2742  j = part(i) + 1
2743  ncell(j) = ncell(j) + blocks(i)%nCell
2744  nface(j) = nface(j) + blocks(i)%nFace
2745  end do
2746 
2747  ! Determine the desirable number of cells and faces per processor.
2748 
2749  ncellseven = ncell(1)
2750  nfaceseven = nface(1)
2751 
2752  do i = 2, nproc
2753  ncellseven = ncellseven + ncell(i)
2754  nfaceseven = nfaceseven + nface(i)
2755  end do
2756 
2757  ncellseven = ncellseven / nproc
2758  nfaceseven = nfaceseven / nproc
2759 
2760  ! Determine the maximum value of nCell and nFace
2761  ! and substract the optimal value.
2762 
2763  ncellmax = abs(maxval(ncell) - ncellseven)
2764  nfacemax = abs(maxval(nface) - nfaceseven)
2765 
2766  ! Determine the tolerance for the cells and faces.
2767 
2768  ncelltol = (ubvec(1) - one) * ncellseven
2769  nfacetol = (ubvec(2) - one) * nfaceseven
2770 
2771  ! Check whether the load balance values for the cells and faces
2772  ! are met.
2773 
2774  cellsbalanced = .true.
2775  facesbalanced = .true.
2776 
2777  if (ncellmax > ncelltol) cellsbalanced = .false.
2778  if (nfacemax > nfacetol) facesbalanced = .false.
2779 
2780  ! Determine the load imbalances for the cells and faces
2781  ! and store it in ubvec.
2782 
2783  ubvec(1) = real(ncellmax, realtype) &
2784  / real(ncellseven, realtype)
2785  ubvec(2) = real(nfacemax, realtype) &
2786  / real(nfaceseven, realtype)
2787 
2788  end subroutine checkloadbalance
2789 
2790  subroutine splitblock(compBlock, nSub, nCells, ranges)
2791  !
2792  ! splitBlock tries to split the given computational block into
2793  ! the desired number of subblocks nSub. However it can happen
2794  ! that nSub is a strange number and a different splitting is
2795  ! performed. On return, nSub contains the actual number into the
2796  ! block is split. This number is smaller or equal to nSub on
2797  ! entry. As it is possible that the computational block itself
2798  ! is a subblock of an original cgns block, on return ranges will
2799  ! contain the nodal ranges of the subblocks in the original cgns
2800  ! block. The splitting attempts to keep the needed multigrid
2801  ! capabilities as much as possible.
2802  ! A recursive bisection algorithm is used.
2803  !
2804  use constants
2805  use inputiteration, only: smoother
2806  use partitionmod, onlY: distributionblocktype
2807  implicit none
2808  !
2809  ! Subroutine arguments.
2810  !
2811  type(distributionblocktype), intent(in) :: compBlock
2812  integer(kind=intType), intent(in) :: nCells
2813  integer(kind=intType), intent(inout) :: nSub
2814 
2815  integer(kind=intType), dimension(nSub, 3, 2), intent(out) :: ranges
2816  !
2817  ! Local variables.
2818  !
2819  integer(kind=intType) :: nLevels, level, nn, mm, nTarget
2820  integer(kind=intType) :: nSplit, nSplitNew
2821 
2822  integer(kind=intType), dimension(nSub) :: nSubblocks
2823 
2824  logical, dimension(3) :: viscousDir
2825 
2826  ! Determine the viscous directions of the block.
2827 
2828  viscousdir = .false.
2829  do nn = 1, compblock%nBocos
2830 
2831  ! Check for viscous boundary conditions and if found set the
2832  ! corresponding direction to viscous.
2833 
2834  if (compblock%BCType(nn) == nswalladiabatic .or. &
2835  compblock%BCType(nn) == nswallisothermal) then
2836 
2837  select case (compblock%BCFaceID(nn))
2838  case (imin, imax)
2839  viscousdir(1) = .true.
2840 
2841  case (jmin, jmax)
2842  viscousdir(2) = .true.
2843 
2844  case (kmin, kmax)
2845  viscousdir(3) = .true.
2846  end select
2847  end if
2848  end do
2849 
2850  ! Set viscousDir to .false. if an explicit smoother is used.
2851 
2852  if (smoother == rungekutta) viscousdir = .false.
2853 
2854  ! Determine the number of levels in the bisection.
2855 
2856  nlevels = log(real(nsub, realtype)) / log(two)
2857  if (2**nlevels < nsub) nlevels = nlevels + 1
2858 
2859  ! Initialize the range of the first splittable block to the
2860  ! entire block.
2861 
2862  ranges(1, 1, 1) = 1; ranges(1, 1, 2) = compblock%il
2863  ranges(1, 2, 1) = 1; ranges(1, 2, 2) = compblock%jl
2864  ranges(1, 3, 1) = 1; ranges(1, 3, 2) = compblock%kl
2865 
2866  ! Initialize the number of blocks to be split to 1 and the
2867  ! number of blocks it should be split into to nSub.
2868 
2869  nsplit = 1
2870  nsubblocks(1) = nsub
2871 
2872  ! Loop over the number of levels
2873 
2874  levelloop: do level = 1, nlevels
2875 
2876  ! Initialize nSplitNew to nSplit; it will be used to store the
2877  ! position of the new subblocks. At the end of this loop this
2878  ! value will be set to nSplit for the new level.
2879 
2880  nsplitnew = nsplit
2881 
2882  ! Loop over the number of blocks to be split.
2883 
2884  blocksplitloop: do nn = 1, nsplit
2885 
2886  ! Determine the situation we are having it.
2887 
2888  select case (nsubblocks(nn))
2889 
2890  case (2_inttype, 3_inttype)
2891 
2892  ! Subblock must be split in either two or three. Store
2893  ! the number into the new subblocks must be in nSubblocks.
2894  ! As the routine split2block puts the subblock with the
2895  ! target number of cells in position nSplitNew, this
2896  ! block should not be split any further and thus gets a
2897  ! value of 1. Position nn gets the old value -1, which
2898  ! will be either 1 or 2
2899 
2900  nsplitnew = nsplitnew + 1
2901  nsubblocks(nsplitnew) = 1
2902  nsubblocks(nn) = nsubblocks(nn) - 1
2903 
2904  ! Split the block into two, where a block with nCells
2905  ! should be split off.
2906 
2907  call split2block(nsub, nn, nsplitnew, ncells, ranges, &
2908  viscousdir)
2909 
2910  !===========================================================
2911 
2912  case (4_inttype:)
2913 
2914  ! Subblock must be split in 4 or more. First determine the
2915  ! number of subblocks into the two new subblocks must be
2916  ! split in the next round. Save the old value of
2917  ! nSubblocks(nn) in mm.
2918 
2919  mm = nsubblocks(nn)
2920 
2921  nsplitnew = nsplitnew + 1
2922  nsubblocks(nsplitnew) = nsubblocks(nn) / 2
2923  nsubblocks(nn) = nsubblocks(nn) &
2924  - nsubblocks(nsplitnew)
2925 
2926  ! Determine the number of cells for the new subblock
2927  ! nSplitNew.
2928 
2929  ntarget = (ranges(nn, 1, 2) - ranges(nn, 1, 1)) &
2930  * (ranges(nn, 2, 2) - ranges(nn, 2, 1)) &
2931  * (ranges(nn, 3, 2) - ranges(nn, 3, 1))
2932  ntarget = ntarget * (real(nsubblocks(nsplitnew), realtype) &
2933  / real(mm, realtype))
2934 
2935  ! Split the block into two, where one block should get
2936  ! approximately nTarget cells.
2937 
2938  call split2block(nsub, nn, nsplitnew, ntarget, ranges, &
2939  viscousdir)
2940 
2941  end select
2942 
2943  end do blocksplitloop
2944 
2945  ! Set nSplit to nSplitNew for the next level of bisection.
2946 
2947  nsplit = nsplitnew
2948 
2949  end do levelloop
2950 
2951  ! Get rid of the possible zero cell ranges. In that case
2952  ! nSub will change as well. Add the offset of the computational
2953  ! block, as this might be a subblock itself.
2954 
2955  mm = nsub
2956  nsub = 0
2957  do nn = 1, mm
2958 
2959  ! Test whether this is a non-zero cell subblock.
2960 
2961  if (ranges(nn, 1, 2) > ranges(nn, 1, 1) .and. &
2962  ranges(nn, 2, 2) > ranges(nn, 2, 1) .and. &
2963  ranges(nn, 3, 2) > ranges(nn, 3, 1)) then
2964 
2965  ! This is a valid subblock. Update nSub and copy the range
2966  ! including the offset of the computational block.
2967 
2968  nsub = nsub + 1
2969 
2970  ranges(nsub, 1, 1) = ranges(nn, 1, 1) + compblock%iBegor - 1
2971  ranges(nsub, 2, 1) = ranges(nn, 2, 1) + compblock%jBegor - 1
2972  ranges(nsub, 3, 1) = ranges(nn, 3, 1) + compblock%kBegor - 1
2973 
2974  ranges(nsub, 1, 2) = ranges(nn, 1, 2) + compblock%iBegor - 1
2975  ranges(nsub, 2, 2) = ranges(nn, 2, 2) + compblock%jBegor - 1
2976  ranges(nsub, 3, 2) = ranges(nn, 3, 2) + compblock%kBegor - 1
2977 
2978  end if
2979  end do
2980 
2981  end subroutine splitblock
2982 
2983  !========================================================================
2984 
2985  subroutine split2block(nSub, n1, n2, nTarget, ranges, &
2986  viscousDir)
2987  !
2988  ! split2block splits the block stored in ranges(n1,:,:) into
2989  ! two. The new blocks are stored in ranges(n1,:,:) and
2990  ! ranges(n2,:,:), where the n2 will store the block with the
2991  ! * number of cells closest to nTarget.
2992  !
2993  use inputiteration
2994  implicit none
2995  !
2996  ! Subroutine arguments.
2997  !
2998  integer(kind=intType), intent(in) :: nSub, n1, n2, nTarget
2999 
3000  integer(kind=intType), dimension(nSub, 3, 2), intent(inout) :: ranges
3001 
3002  logical, dimension(3), intent(in) :: viscousDir
3003  !
3004  ! Local variables.
3005  !
3006  integer(kind=intType) :: nMG, mm, ii, jj, i, iiOpt, iOpt
3007  integer(kind=intType) :: nCellOpt
3008 
3009  integer(kind=intType), dimension(3) :: nc, mc, c, nf, prefDir
3010 
3011  ! Determine the maximum of the number of mg levels needed in the
3012  ! cycle and the mg start level. Store this value in nMG.
3013 
3014  nmg = max(nmglevels, mgstartlevel)
3015 
3016  ! Store the number of cells of the original subblock in nc.
3017 
3018  nc(1) = ranges(n1, 1, 2) - ranges(n1, 1, 1)
3019  nc(2) = ranges(n1, 2, 2) - ranges(n1, 2, 1)
3020  nc(3) = ranges(n1, 3, 2) - ranges(n1, 3, 1)
3021 
3022  ! Determine its multigrid capabilities in the three directions.
3023 
3024  loopdir: do mm = 1, 3
3025 
3026  ! Initialize nf(mm), number of mg levels, to 1 and ii to 2.
3027  ! nf is used as a temporary buffer.
3028 
3029  nf(mm) = 1
3030  ii = 2
3031 
3032  ! Loop to determine the number of mg levels.
3033 
3034  do
3035  ! Conditions to exit the loop. The first condition is the
3036  ! truncation to the maximum number of mg levels needed by the
3037  ! iterative algorithm of the solver. The second condition is
3038  ! simply that the number of cells do not allow more standard
3039  ! mg levels.
3040 
3041  if (nf(mm) == nmg) exit
3042  if (mod(nc(mm), ii) /= 0) exit
3043 
3044  ! Update nf(mm) and multiply ii by 2 to test for the
3045  ! next level.
3046 
3047  nf(mm) = nf(mm) + 1
3048  ii = 2 * ii
3049  end do
3050 
3051  ! Determine the constants c(mm) and mc(mm), such that
3052  ! nc(mm) = c(mm)*mc(mm). Mc(mm) is the number of cells in every
3053  ! direction per supercell. A supercell is the smallest possible
3054  ! block able to do multigridding, i.e. mc(mm) = 2**(nf(mm)-1)
3055  ! if the block allows for nMG multigrid levels.
3056 
3057  c(mm) = 2 * nc(mm) / ii
3058  mc(mm) = nc(mm) / c(mm)
3059 
3060  end do loopdir
3061 
3062  ! Test whether the multigrid requirements are not too restrictive
3063  ! for this subblock. If they are, relax them, if possible.
3064 
3065  if (nmg > 1 .and. max(c(1), c(2), c(3)) == 1) then
3066  c(1) = c(1) * 2; c(2) = c(2) * 2; c(3) = c(3) * 2
3067  mc(1) = mc(1) / 2; mc(2) = mc(2) / 2; mc(3) = mc(3) / 2
3068  end if
3069 
3070  ! Determine the number of faces on a slice in the
3071  ! three directions.
3072 
3073  nf(1) = nc(2) * nc(3)
3074  nf(2) = nc(1) * nc(3)
3075  nf(3) = nc(1) * nc(2)
3076 
3077  ! Determine the preferred split direction; this the direction
3078  ! which results in the least amount of additional faces. This is
3079  ! important, because the number of flux evaluations is
3080  ! proportional to the number of faces.
3081 
3082  if (nf(1) < nf(2)) then
3083  if (nf(1) < nf(3)) then
3084  prefdir(1) = 1
3085  if (nf(2) < nf(3)) then
3086  prefdir(2) = 2
3087  prefdir(3) = 3
3088  else
3089  prefdir(2) = 3
3090  prefdir(3) = 2
3091  end if
3092  else
3093  prefdir(1) = 3
3094  prefdir(2) = 1
3095  prefdir(3) = 2
3096  end if
3097  else if (nf(2) < nf(3)) then
3098  prefdir(1) = 2
3099  if (nf(1) < nf(3)) then
3100  prefdir(2) = 1
3101  prefdir(3) = 3
3102  else
3103  prefdir(2) = 3
3104  prefdir(3) = 1
3105  end if
3106  else
3107  prefdir(1) = 3
3108  prefdir(2) = 2
3109  prefdir(3) = 1
3110  end if
3111 
3112  ! Take viscousDir into account to determine the preference
3113  ! direction. For implicit methods it may not be a good idea to
3114  ! split the block parallel to a viscous wall.
3115 
3116  ! Check the 1st preference direction. If it is viscous swap
3117  ! it with the best inviscid direction, if available.
3118 
3119  if (viscousdir(prefdir(1))) then
3120  if (.not. viscousdir(prefdir(2))) then
3121  mm = prefdir(1)
3122  prefdir(1) = prefdir(2)
3123  prefdir(2) = mm
3124  else if (.not. viscousdir(prefdir(3))) then
3125  mm = prefdir(1)
3126  prefdir(1) = prefdir(3)
3127  prefdir(3) = mm
3128  end if
3129  end if
3130 
3131  ! Check the second preference direction. If it is viscous
3132  ! try to swap it whith the third direction.
3133 
3134  if (viscousdir(prefdir(2)) .and. &
3135  .not. viscousdir(prefdir(3))) then
3136  mm = prefdir(2)
3137  prefdir(2) = prefdir(3)
3138  prefdir(3) = mm
3139  end if
3140 
3141  ! Determine the splitting which is best. Initialize nCellOpt
3142  ! to a ridiculously high number.
3143 
3144  ncellopt = 10 * nc(1) * nc(2) * nc(3)
3145 
3146  ! Loop over the three directions.
3147 
3148  do mm = 1, 3
3149 
3150  ! Store the current direction considered in ii and determine
3151  ! the index i, which gives a number of cells as close as
3152  ! possible to the desired number.
3153 
3154  ii = prefdir(mm)
3155  i = nint(real(ntarget, realtype) &
3156  / real(mc(ii) * nf(ii), realtype), inttype)
3157  i = max(i, 1_inttype)
3158 
3159  ! Determine whether the corresponding number of cells is
3160  ! closer to the target number than the currently stored value.
3161  ! If so, store the settings for this splitting.
3162 
3163  jj = i * mc(ii) * nf(ii)
3164  if (abs(jj - ntarget) < abs(ncellopt - ntarget)) then
3165  ncellopt = jj
3166  iiopt = ii
3167  iopt = i
3168  end if
3169 
3170  end do
3171 
3172  ! Copy the range of subblock n1 into n2 as initialization.
3173 
3174  ranges(n2, 1, 1) = ranges(n1, 1, 1); ranges(n2, 1, 2) = ranges(n1, 1, 2)
3175  ranges(n2, 2, 1) = ranges(n1, 2, 1); ranges(n2, 2, 2) = ranges(n1, 2, 2)
3176  ranges(n2, 3, 1) = ranges(n1, 3, 1); ranges(n2, 3, 2) = ranges(n1, 3, 2)
3177 
3178  ! Determine the nodal index where the splitting takes place.
3179 
3180  jj = iopt * mc(iiopt) + ranges(n1, iiopt, 1)
3181 
3182  ! Adapt the corresponding indices of the new subblocks n1 and
3183  ! n2, such that it corresponds to the new situation; n2 should
3184  ! contain the subblock, which contains a number of cells as
3185  ! close as possible to nTarget.
3186 
3187  ranges(n2, iiopt, 2) = jj
3188  ranges(n1, iiopt, 1) = jj
3189 
3190  end subroutine split2block
3191 
3192  subroutine reallocsubfacememory(ii, nAlloc)
3193  !
3194  ! reallocSubfaceMemory reallocates the memory to store the
3195  ! subface information for the given block ii. On entry nAlloc
3196  ! contains the current number of allocated subfaces, on exit
3197  ! this is updated to the new number.
3198  !
3199  use constants
3200  use partitionmod, only: blocks
3201  use utils, only: reallocateinteger
3202  implicit none
3203  !
3204  ! Subroutine arguments.
3205  !
3206  integer(kind=intType), intent(in) :: ii
3207  integer(kind=intType), intent(inout) :: nAlloc
3208  !
3209  ! Local arguments.
3210  !
3211  integer(kind=intType) :: nOld, nNew
3212  !
3213  ! Store the old value of the allocated array and determine the
3214  ! new one. Store this new value in nAlloc.
3215 
3216  nold = nalloc
3217  nnew = nold + 6
3218  nalloc = nnew
3219 
3220  ! Reallocate the memory.
3221 
3222  call reallocateinteger(blocks(ii)%BCType, nnew, nold, .true.)
3223  call reallocateinteger(blocks(ii)%BCFaceID, nnew, nold, .true.)
3224  call reallocateinteger(blocks(ii)%cgnsSubface, nnew, nold, .true.)
3225 
3226  call reallocateinteger(blocks(ii)%inBeg, nnew, nold, .true.)
3227  call reallocateinteger(blocks(ii)%jnBeg, nnew, nold, .true.)
3228  call reallocateinteger(blocks(ii)%knBeg, nnew, nold, .true.)
3229 
3230  call reallocateinteger(blocks(ii)%inEnd, nnew, nold, .true.)
3231  call reallocateinteger(blocks(ii)%jnEnd, nnew, nold, .true.)
3232  call reallocateinteger(blocks(ii)%knEnd, nnew, nold, .true.)
3233 
3234  call reallocateinteger(blocks(ii)%dinBeg, nnew, nold, .true.)
3235  call reallocateinteger(blocks(ii)%djnBeg, nnew, nold, .true.)
3236  call reallocateinteger(blocks(ii)%dknBeg, nnew, nold, .true.)
3237 
3238  call reallocateinteger(blocks(ii)%dinEnd, nnew, nold, .true.)
3239  call reallocateinteger(blocks(ii)%djnEnd, nnew, nold, .true.)
3240  call reallocateinteger(blocks(ii)%dknEnd, nnew, nold, .true.)
3241 
3242  call reallocateinteger(blocks(ii)%neighBlock, nnew, nold, .true.)
3243  call reallocateinteger(blocks(ii)%groupNum, nnew, nold, .true.)
3244 
3245  call reallocateinteger(blocks(ii)%l1, nnew, nold, .true.)
3246  call reallocateinteger(blocks(ii)%l2, nnew, nold, .true.)
3247  call reallocateinteger(blocks(ii)%l3, nnew, nold, .true.)
3248 
3249  end subroutine reallocsubfacememory
3250 
3251 end module loadbalance
splitblocksloadbalance
subroutine splitblocksloadbalance
Definition: loadBalance.F90:881
searchinternalneighbors
subroutine searchinternalneighbors
Definition: loadBalance.F90:2305
splittingisokay
logical function splittingisokay(cgnsID)
Definition: loadBalance.F90:626
splitblockinitialization
subroutine splitblockinitialization(cgnsID)
Definition: loadBalance.F90:678
metisinterface
void metisinterface(int *n, int *ncon, idxtype *xadj, idxtype *adjncy, idxtype *vwgt, idxtype *adjwgt, int *wgtflag, int *numflag, int *nparts, float *ubvec, int *options, int *edgecut, idxtype *part)
Definition: metisInterface.c:87
block
Definition: block.F90:1
block::ndom
integer(kind=inttype) ndom
Definition: block.F90:761
block::flowdoms
type(blocktype), dimension(:, :, :), allocatable, target flowdoms
Definition: block.F90:771
cgnsgrid::cgnsdoms
type(cgnsblockinfotype), dimension(:), allocatable cgnsdoms
Definition: cgnsGrid.F90:495
cgnsgrid::cgnsfamilies
type(cgnsfamilytype), dimension(:), allocatable cgnsfamilies
Definition: cgnsGrid.F90:504
cgnsgrid::cgnsndom
integer(kind=inttype) cgnsndom
Definition: cgnsGrid.F90:491
commonformats::strings
character(len=maxstringlen) strings
Definition: commonFormats.F90:21
communication::adflow_comm_world
integer adflow_comm_world
Definition: communication.F90:175
constants::domaininterfacep
integer(kind=inttype), parameter domaininterfacep
Definition: constants.F90:279
constants::slidinginterface
integer(kind=inttype), parameter slidinginterface
Definition: constants.F90:275
constants::eulerwall
integer(kind=inttype), parameter eulerwall
Definition: constants.F90:262
constants::farfield
integer(kind=inttype), parameter farfield
Definition: constants.F90:263
constants::oversetouterbound
integer(kind=inttype), parameter oversetouterbound
Definition: constants.F90:276
constants::massbleedoutflow
integer(kind=inttype), parameter massbleedoutflow
Definition: constants.F90:269
constants::domaininterfacerhouvw
integer(kind=inttype), parameter domaininterfacerhouvw
Definition: constants.F90:278
constants::domaininterfacerho
integer(kind=inttype), parameter domaininterfacerho
Definition: constants.F90:280
constants::imax
integer(kind=inttype), parameter imax
Definition: constants.F90:293
constants::supersonicinflow
integer(kind=inttype), parameter supersonicinflow
Definition: constants.F90:264
constants::kmin
integer(kind=inttype), parameter kmin
Definition: constants.F90:296
constants::massbleedinflow
integer(kind=inttype), parameter massbleedinflow
Definition: constants.F90:268
constants::jmax
integer(kind=inttype), parameter jmax
Definition: constants.F90:295
constants::supersonicoutflow
integer(kind=inttype), parameter supersonicoutflow
Definition: constants.F90:266
constants::mdot
integer(kind=inttype), parameter mdot
Definition: constants.F90:270
constants::b2bmatch
integer(kind=inttype), parameter b2bmatch
Definition: constants.F90:273
constants::nswalladiabatic
integer(kind=inttype), parameter nswalladiabatic
Definition: constants.F90:260
constants::symm
integer(kind=inttype), parameter symm
Definition: constants.F90:258
constants::timespectral
integer(kind=inttype), parameter timespectral
Definition: constants.F90:115
constants::unsteady
integer(kind=inttype), parameter unsteady
Definition: constants.F90:115
constants::symmpolar
integer(kind=inttype), parameter symmpolar
Definition: constants.F90:259
constants::subsonicoutflow
integer(kind=inttype), parameter subsonicoutflow
Definition: constants.F90:267
constants::rungekutta
integer(kind=inttype), parameter rungekutta
Definition: constants.F90:201
constants::one
real(kind=realtype), parameter one
Definition: constants.F90:72
constants::imin
integer(kind=inttype), parameter imin
Definition: constants.F90:292
constants::domaininterfacetotal
integer(kind=inttype), parameter domaininterfacetotal
Definition: constants.F90:281
constants::two
real(kind=realtype), parameter two
Definition: constants.F90:73
constants::bcthrust
integer(kind=inttype), parameter bcthrust
Definition: constants.F90:271
constants::nswallisothermal
integer(kind=inttype), parameter nswallisothermal
Definition: constants.F90:261
constants::domaininterfaceall
integer(kind=inttype), parameter domaininterfaceall
Definition: constants.F90:277
constants::subsonicinflow
integer(kind=inttype), parameter subsonicinflow
Definition: constants.F90:265
constants::kmax
integer(kind=inttype), parameter kmax
Definition: constants.F90:297
constants::jmin
integer(kind=inttype), parameter jmin
Definition: constants.F90:294
constants::extrap
integer(kind=inttype), parameter extrap
Definition: constants.F90:272
inputiteration::smoother
integer(kind=inttype) smoother
Definition: inputParam.F90:264
inputiteration::mgstartlevel
integer(kind=inttype) mgstartlevel
Definition: inputParam.F90:270
inputiteration::nmglevels
integer(kind=inttype) nmglevels
Definition: inputParam.F90:271
inputmotion::gridmotionspecified
logical gridmotionspecified
Definition: inputParam.F90:481
inputparallel::loadbalanceiter
integer(kind=inttype) loadbalanceiter
Definition: inputParam.F90:502
inputparallel::partitionlikenproc
integer(kind=inttype) partitionlikenproc
Definition: inputParam.F90:502
inputparallel::loadimbalance
real(realtype) loadimbalance
Definition: inputParam.F90:500
inputparallel::splitblocks
logical splitblocks
Definition: inputParam.F90:501
inputphysics::equationmode
integer(kind=inttype) equationmode
Definition: inputParam.F90:583
inputphysics::flowtype
integer(kind=inttype) flowtype
Definition: inputParam.F90:583
inputtimespectral::ntimeintervalsspectral
integer(kind=inttype) ntimeintervalsspectral
Definition: inputParam.F90:645
iteration::deforming_grid
logical deforming_grid
Definition: iteration.f90:69
loadbalance::initflowdoms
subroutine initflowdoms
Definition: loadBalance.F90:1181
loadbalance::checkloadbalance
subroutine checkloadbalance(cellsBalanced, facesBalanced)
Definition: loadBalance.F90:2707
loadbalance::graphpartitioning
subroutine graphpartitioning(emptyPartitions, commNeglected)
Definition: loadBalance.F90:2434
loadbalance::splitblock
subroutine splitblock(compBlock, nSub, nCells, ranges)
Definition: loadBalance.F90:2791
loadbalance::loadbalancegrid
subroutine loadbalancegrid
Definition: loadBalance.F90:6
loadbalance::blockdistribution
subroutine blockdistribution
Definition: loadBalance.F90:410
loadbalance::externalfacessubblock
subroutine externalfacessubblock(cgnsID, ii, jj, nSubPerCGNS, nAlloc, splitInfo)
Definition: loadBalance.F90:1770
loadbalance::bcfacessubblock
subroutine bcfacessubblock(cgnsID, ii, jj)
Definition: loadBalance.F90:1545
loadbalance::reallocsubfacememory
subroutine reallocsubfacememory(ii, nAlloc)
Definition: loadBalance.F90:3193
loadbalance::determinecomputeblocks
subroutine determinecomputeblocks(splitInfo)
Definition: loadBalance.F90:1371
loadbalance::sortsubfaces
subroutine sortsubfaces(oldSubfaceID, blockID)
Definition: loadBalance.F90:1222
loadbalance::internalfacessubblock
subroutine internalfacessubblock(cgnsID, ii, jj, nSubPerCGNS, nAlloc, splitInfo)
Definition: loadBalance.F90:2133
loadbalance::split2block
subroutine split2block(nSub, n1, n2, nTarget, ranges, viscousDir)
Definition: loadBalance.F90:2987
partitionmod::sortrangessplitinfo
subroutine sortrangessplitinfo(splitInfo)
Definition: partitionMod.F90:1117
partitionmod::qsortsubblocksofcgnstype
subroutine qsortsubblocksofcgnstype(arr, nn)
Definition: partitionMod.F90:577
partitionmod::part
integer(kind=inttype), dimension(:), allocatable part
Definition: partitionMod.F90:265
partitionmod::ubvec
real, dimension(2) ubvec
Definition: partitionMod.F90:263
partitionmod::nblocks
integer(kind=inttype) nblocks
Definition: partitionMod.F90:118
partitionmod::blocks
type(distributionblocktype), dimension(:), allocatable blocks
Definition: partitionMod.F90:119
sorting::qsortintegers
subroutine qsortintegers(arr, nn)
Definition: sorting.F90:101
sorting::bsearchintegers
integer(kind=inttype) function bsearchintegers(key, base)
Definition: sorting.F90:18
utils
Definition: utils.F90:1
utils::reallocateinteger
subroutine reallocateinteger(intArray, newSize, oldSize, alwaysFreeMem)
Definition: utils.F90:2783
utils::nullifyflowdompointers
subroutine nullifyflowdompointers(nn, level, sps)
Definition: utils.F90:2580
utils::delta
integer(kind=inttype) function delta(val1, val2)
Definition: utils.F90:2534
utils::terminate
subroutine terminate(routineName, errorMessage)
Definition: utils.F90:502
Generated by doxygen 1.9.1