Now lo(:) and hi(:) are determined from ng(:), and not vice-versa

README.md

@@ -10,6 +10,7 @@ P. Costa. *A FFT-accelerated multi-block finite-difference solver for massively
 Comput. Phys. Commun. 271 : 108194 (2022) [[DOI:10.1016/j.cpc.2021.108194]](doi.org/10.1016/j.cpc.2021.108194) [[arXiv:2106.03583]](https://arxiv.org/pdf/2106.03583.pdf).
 
 ## News
+*[08/07/2022]* The input files describing the block geometry (under `geo/block.???`) have been simplified. Now, instead of prescribing the lower and upper extents of each block `lo(:)` and `hi(:)`, the number of grid points `ng(:)` is prescribed. This change makes it much easier to refine the grid without, since one does not need to correct extents of neighboring blocks. See the updated [`src/INFO_INPUT.md`](src/INFO_INPUT.md) for more details.
 
 ## Features
 

examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.001

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-1 1 1 ! lo(1:3)
-32 32 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 0. 0. 0. ! lmin(1:3)
 .5 .5 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.002

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-1 33 1 ! lo(1:3)
-32 64 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 0. .5 0. ! lmin(1:3)
 .5 1. 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.003

@@ -1,6 +1,5 @@
 2 1 1 ! dims(1:3)
-33 33 1 ! lo(1:3)
-128 64 64 ! hi(1:3)
+96 32 64 ! ng(1:3)
 .5 .5 0. ! lmin(1:3)
 2.0 1. 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/block.001

@@ -1,8 +1,7 @@
 1 1 1 ! dims(1:3)
-1 1 1 ! lo(1:3)
-320 2 64 ! hi(1:3)
--0.5 0. -0.5 ! lmin(1:3)
-4.5 .03125 0.5 ! lmax(1:3)
+320 2 64 ! ng(1:3)
+-0.5 0. -0.5 ! lmin(1:3)
+4.5 .03125 0.5 ! lmax(1:3)
 0 0 0 ! gt(1:3)
 0. 0. 0. ! gr(1:3)
 D N F F F F ! cbcvel(0:1,1:3,1) [u BC type]

examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.002

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-33 1 1 ! lo(1:3)
-64 32 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 .5 0. 0. ! lmin(1:3)
 1. .5 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.003

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-1 33 1 ! lo(1:3)
-32 64 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 0. .5 0. ! lmin(1:3)
 .5 1. 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.004

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-33 33 1 ! lo(1:3)
-64 64 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 .5 .5 0. ! lmin(1:3)
 1. 1. 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

src/INFO_INPUT.md

@@ -16,7 +16,7 @@ Let us start with block-independent parameters in `dns.in`. Consider the followi
  100 100. 0.1 ! nstep,time_max,tw_max
  T F F ! stop_type(1:3)
  F T ! restart,is_overwrite_save
-10 10 20 5000 10000 2000  ! icheck, iout0d, iout1d, iout2d, iout3d, isave
+ 10 10 20 5000 10000 2000 ! icheck, iout0d, iout1d, iout2d, iout3d, isave
  0. 0. 0. ! bforce(1:3)
  4 ! nthreadsmax
 ~~~
@@ -114,8 +114,7 @@ The geometry, boundary and initial conditions, domain decompositions is set in a
 `geo/block.001`:
 ~~~
 1 1 1 ! dims(1:3)
-1 1 1 ! lo(1:3)
-32 32 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 0. 0. 0. ! lmin(1:3)
 .5 .5 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)
@@ -135,8 +134,7 @@ zer ! inivel
 `geo/block.002`:
 ~~~
 1 1 1 ! dims(1:3)
-1 33 1 ! lo(1:3)
-32 64 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 0. .5 0. ! lmin(1:3)
 .5 1. 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)
@@ -156,8 +154,7 @@ zer ! inivel
 `geo/block.003`:
 ~~~
 2 1 1 ! dims(1:3)
-33 33 1 ! lo(1:3)
-128 64 64 ! hi(1:3)
+96 32 64 ! ng(1:3)
 .5 .5 0. ! lmin(1:3)
 2.0 1. 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)
@@ -180,8 +177,7 @@ zer ! inivel
 
 ~~~
 1 1 1 ! dims(1:3)
-1 1 1 ! lo(1:3)
-32 32 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 0. 0. 0. ! lmin(1:3)
 .5 .5 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)
@@ -192,7 +188,7 @@ These lines set the domain decomposition and computational grid for each block.
 `dims` is the **number of computational** subdomains in each direction.
 
 
-`lo(1:3)` and `hi(1:3)` are the **coordinates of the lower and upper corners** of the block in question, **in index space**.
+`ng(1:3)` is the **number of grid points along each direction** of the block in question.
 
 `lmin(1:3)` and `lmax(1:3)` are the **physical coordinates of the lower and upper corners** of the block in question.
 

src/geo/block.001

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-1 1 1 ! lo(1:3)
-32 32 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 0. 0. 0. ! lmin(1:3)
 .5 .5 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

src/geo/block.002

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-1 33 1 ! lo(1:3)
-32 64 64 ! hi(1:3)
+32 32 64 ! ng(1:3)
 0. .5 0. ! lmin(1:3)
 .5 1. 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

src/geo/block.003

@@ -1,6 +1,5 @@
 2 1 1 ! dims(1:3)
-33 33 1 ! lo(1:3)
-128 64 64 ! hi(1:3)
+96 32 64 ! ng(1:3)
 .5 .5 0. ! lmin(1:3)
 2.0 1. 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

src/initmpi.f90

@@ -6,43 +6,97 @@ module mod_initmpi
  private
  public initmpi
  contains
- subroutine initmpi(my_block,id_first,dims,cbc,bc,periods,lmin,lmax,gt,gr,lo,hi,ng,nb,is_bound,halos)
+ subroutine initmpi(my_block,nblocks,id_first,dims,cbc,bc,periods,lmin,lmax,gt,gr,lo,hi,lo_g,hi_g,ng,nb,is_bound,halos)
  implicit none
- integer , intent(in ) :: my_block,id_first
+ integer , intent(in ) :: my_block,nblocks,id_first
  integer , intent(inout), dimension( 3) :: dims
  character(len=1), intent(in ), dimension(0:1,3) :: cbc
  real(rp) , intent(in ), dimension(0:1,3) :: bc
  integer , intent(in ), dimension( 3) :: periods
  real(rp) , intent(in ), dimension( 3) :: lmin,lmax
  integer , intent(in ), dimension( 3) :: gt
  real(rp) , intent(in ), dimension( 3) :: gr
- integer , intent(inout), dimension( 3) :: lo,hi
+ integer , intent(inout), dimension( 3) :: lo,hi,lo_g,hi_g
  integer , intent(out ), dimension( 3) :: ng
  integer , intent(out ), dimension(0:1,3) :: nb
  logical , intent(out ), dimension(0:1,3) :: is_bound
  type(MPI_DATATYPE), intent(out ), dimension( 3) :: halos
  integer :: nrank
- integer , dimension( 3) :: n,start,coords,lo_g,hi_g
- integer , allocatable, dimension(:,:) :: lo_all,hi_all,gt_all,l,lo_g_all,hi_g_allo_g_all,hi_g_all
+ integer , dimension( 3) :: n,start,coords
+ integer , allocatable, dimension(:,:) :: ng_all,lo_all,hi_all,gt_all,l,lo_g_all,hi_g_allo_g_all,hi_g_all
  real(rp), allocatable, dimension(:,:) :: lmin_all,lmax_all,gr_all
  integer , allocatable, dimension(:) :: blocks_all
+ logical , allocatable, dimension(:) :: is_done,is_unlocked
  real(rp) , allocatable, dimension(:,:,:) :: bc_all
  character(len=1), allocatable, dimension(:,:,:) :: cbc_all
- integer :: i,j,k,idir,iidir,inb,irank
+ integer :: i,j,k,idir,iidir,inb,irank,iblock,ifriend,ioffset,idir_t(2)
  logical :: is_nb,found_friend
  integer(i8) :: ntot,ntot_max,ntot_min,ntot_sum
  integer :: isize
  type(MPI_DATATYPE) :: MPI_INTEGER_I8
+ type(MPI_COMM ) :: comm_leaders
  !
  call MPI_COMM_SIZE(MPI_COMM_WORLD,nrank)
  call MPI_COMM_SPLIT(MPI_COMM_WORLD,my_block,myid,comm_block)
  call MPI_COMM_RANK(comm_block, myid_block)
  !
  ! generate Cartesian topology
  !
- hi_g(:) = hi(:)
- lo_g(:) = lo(:)
- ng(:) = hi_g(:)-lo_g(:)+1
+ ! determine lo_g and hi_g from the prescribed decomposition
+ !
+ call MPI_COMM_SPLIT(MPI_COMM_WORLD,myid_block,myid,comm_leaders)
+ if(myid_block == 0) then
+ !
+ ! gather block size and connectivity information
+ !
+ allocate(lo_all(3,nblocks),ng_all(3,nblocks),cbc_all(0:1,3,nblocks),bc_all(0:1,3,nblocks),lmin_all(3,nblocks))
+ call MPI_ALLGATHER(ng ,3,MPI_INTEGER ,ng_all ,3,MPI_INTEGER ,comm_leaders)
+ call MPI_ALLGATHER(cbc,6,MPI_CHARACTER,cbc_all,6,MPI_CHARACTER,comm_leaders)
+ call MPI_ALLGATHER( bc,6,MPI_REAL_RP , bc_all,6,MPI_REAL_RP ,comm_leaders)
+ !
+ ! determine lower bounds, taking block #1 as reference
+ !
+ allocate(is_done(nblocks),is_unlocked(nblocks))
+ is_done(:) = .false.
+ is_unlocked(:) = .false.
+ iblock = 1
+ lo_all(:,iblock) = 1
+ is_unlocked(iblock) = .true.
+ do while(.not.all(is_done))
+ if(iblock < 1 .or. iblock > nblocks) then
+ iblock = findloc(.not.is_done,.true.,1)
+ write(stderr,*) 'ERROR: invalid connectivity for block ', iblock, '.'
+ write(stderr,*) ''
+ error stop
+ endif
+ do idir=1,3
+ idir_t(:) = pack([1,2,3],[1,2,3] /= idir) ! extract tangential directions
+ do inb=0,1
+ if(cbc_all(inb,idir,iblock) == 'F') then
+ ifriend = nint(bc_all(inb,idir,iblock))
+ if(.not.is_unlocked(ifriend)) then
+ ! periodicity check
+ if(inb == 0 .and. lmin_all(idir,ifriend) > lmin_all(idir,iblock)) cycle
+ if(inb == 1 .and. lmin_all(idir,ifriend) < lmin_all(idir,iblock)) cycle
+ ! normal direction
+ ! inb = 0 -> shift - ng_all(idir,ifriend) | inb = 1 -> shift + ng_all(idir,iblock )
+ ioffset = -(1-inb)*ng_all(idir,ifriend) + inb*ng_all(idir,iblock)
+ lo_all(idir,ifriend) = lo_all(idir,iblock) + ioffset
+ ! tangential directions
+ lo_all(idir_t(:),ifriend) = lo_all(idir_t(:),iblock)
+ is_unlocked(ifriend) = .true.
+ end if
+ end if
+ end do
+ end do
+ is_done(iblock) = .true.
+ iblock = findloc(.not.is_done.and.is_unlocked,.true.,1)
+ end do
+ lo_g(:) = lo_all(:,my_block)
+ deallocate(lo_all,ng_all,cbc_all,bc_all,lmin_all)
+ endif
+ call MPI_BCAST(lo_g,3,MPI_INTEGER,0,comm_block)
+ hi_g(:) = lo_g(:)+ng(:)-1
  !
  ! determine array extents for possibly uneven data
  !

src/main.f90

@@ -214,7 +214,7 @@ program snac
  ! initialize MPI/OpenMP
  !
  !!$call omp_set_num_threads(nthreadsmax) ! ! overwrites the input, disable for now
- call initmpi(my_block,id_first,dims,cbcpre,bcpre,periods,lmin,lmax,gt,gr,lo,hi,ng,nb,is_bound,halos)
+ call initmpi(my_block,nblocks,id_first,dims,cbcpre,bcpre,periods,lmin,lmax,gt,gr,lo,hi,lo_g,hi_g,ng,nb,is_bound,halos)
  lo_1(:) = lo(:) - lo_g(:) + 1 ! lo(:) with 1 as first index in the beginning of each block
  hi_1(:) = hi(:) - lo_g(:) + 1 ! hi(:) with 1 as first index in the beginning of each block
  !

tests/lid_driven_cavity/geo/block.001

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-1 1 1 ! lo(1:3)
-2 32 32 ! hi(1:3)
+2 32 32 ! ng(1:3)
 0. 0. 0. ! lmin(1:3)
 0.3125 .5 .5 ! lmax(1:3)
 0 0 0 ! gt(1:3)

tests/lid_driven_cavity/geo/block.002

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-1 33 1 ! lo(1:3)
-2 64 32 ! hi(1:3)
+2 32 32 ! ng(1:3)
 0. .5 0. ! lmin(1:3)
 0.3125 1. .5 ! lmax(1:3)
 0 0 0 ! gt(1:3)

tests/lid_driven_cavity/geo/block.003

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-1 1 33 ! lo(1:3)
-2 32 64 ! hi(1:3)
+2 32 32 ! ng(1:3)
 0. 0. .5 ! lmin(1:3)
 0.3125 .5 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)

tests/lid_driven_cavity/geo/block.004

@@ -1,6 +1,5 @@
 1 1 1 ! dims(1:3)
-1 33 33 ! lo(1:3)
-2 64 64 ! hi(1:3)
+2 32 32 ! ng(1:3)
 0. .5 .5 ! lmin(1:3)
 0.3125 1. 1. ! lmax(1:3)
 0 0 0 ! gt(1:3)