From 19d9dad1a5eb1ad45243464df92ef86302ad5810 Mon Sep 17 00:00:00 2001
From: Pedro Costa
Date: Fri, 8 Jul 2022 13:08:45 +0000
Subject: [PATCH] Now `lo(:)` and `hi(:)` are determined from `ng(:)`, and not
vice-versa
---
README.md | 1 +
.../L_channel/block.001 | 3 +-
.../L_channel/block.002 | 3 +-
.../L_channel/block.003 | 3 +-
.../kovasznay_flow/block.001 | 7 +-
.../lid_driven_cavity/block.002 | 3 +-
.../lid_driven_cavity/block.003 | 3 +-
.../lid_driven_cavity/block.004 | 3 +-
src/INFO_INPUT.md | 16 ++---
src/geo/block.001 | 3 +-
src/geo/block.002 | 3 +-
src/geo/block.003 | 3 +-
src/initmpi.f90 | 72 ++++++++++++++++---
src/main.f90 | 2 +-
tests/lid_driven_cavity/geo/block.001 | 3 +-
tests/lid_driven_cavity/geo/block.002 | 3 +-
tests/lid_driven_cavity/geo/block.003 | 3 +-
tests/lid_driven_cavity/geo/block.004 | 3 +-
18 files changed, 87 insertions(+), 50 deletions(-)
diff --git a/README.md b/README.md
index ef9f036..011934c 100644
--- a/README.md
+++ b/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
diff --git a/examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.001 b/examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.001
index c2600fc..f900d07 100644
--- a/examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.001
+++ b/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)
diff --git a/examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.002 b/examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.002
index b5ea3ac..29ab78a 100644
--- a/examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.002
+++ b/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)
diff --git a/examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.003 b/examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.003
index 7b94038..e19c329 100644
--- a/examples/__MULTI_BLOCK_GEOMETRY/L_channel/block.003
+++ b/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)
diff --git a/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/block.001 b/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/block.001
index b5bfc0e..168b29d 100644
--- a/examples/__MULTI_BLOCK_GEOMETRY/kovasznay_flow/block.001
+++ b/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]
diff --git a/examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.002 b/examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.002
index 0bf7966..4f557ea 100644
--- a/examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.002
+++ b/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)
diff --git a/examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.003 b/examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.003
index cf33e7a..6414cf1 100644
--- a/examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.003
+++ b/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)
diff --git a/examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.004 b/examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.004
index 2b7b948..4b3ea3b 100644
--- a/examples/__MULTI_BLOCK_GEOMETRY/lid_driven_cavity/block.004
+++ b/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)
diff --git a/src/INFO_INPUT.md b/src/INFO_INPUT.md
index 2379b66..9428da0 100644
--- a/src/INFO_INPUT.md
+++ b/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.
diff --git a/src/geo/block.001 b/src/geo/block.001
index 50ddad3..1910bb2 100644
--- a/src/geo/block.001
+++ b/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)
diff --git a/src/geo/block.002 b/src/geo/block.002
index b5ea3ac..29ab78a 100644
--- a/src/geo/block.002
+++ b/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)
diff --git a/src/geo/block.003 b/src/geo/block.003
index 7b94038..e19c329 100644
--- a/src/geo/block.003
+++ b/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)
diff --git a/src/initmpi.f90 b/src/initmpi.f90
index ece0f99..a652744 100644
--- a/src/initmpi.f90
+++ b/src/initmpi.f90
@@ -6,9 +6,9 @@ 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
@@ -16,23 +16,25 @@ subroutine initmpi(my_block,id_first,dims,cbc,bc,periods,lmin,lmax,gt,gr,lo,hi,n
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)
@@ -40,9 +42,61 @@ subroutine initmpi(my_block,id_first,dims,cbc,bc,periods,lmin,lmax,gt,gr,lo,hi,n
!
! 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
!
diff --git a/src/main.f90 b/src/main.f90
index c05e6aa..ea368dc 100644
--- a/src/main.f90
+++ b/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
!
diff --git a/tests/lid_driven_cavity/geo/block.001 b/tests/lid_driven_cavity/geo/block.001
index d2e9f6f..ab6b264 100644
--- a/tests/lid_driven_cavity/geo/block.001
+++ b/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)
diff --git a/tests/lid_driven_cavity/geo/block.002 b/tests/lid_driven_cavity/geo/block.002
index fe2716f..4ed9ee7 100644
--- a/tests/lid_driven_cavity/geo/block.002
+++ b/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)
diff --git a/tests/lid_driven_cavity/geo/block.003 b/tests/lid_driven_cavity/geo/block.003
index 3d9ae78..6d121bb 100644
--- a/tests/lid_driven_cavity/geo/block.003
+++ b/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)
diff --git a/tests/lid_driven_cavity/geo/block.004 b/tests/lid_driven_cavity/geo/block.004
index 23ca032..b7845f5 100644
--- a/tests/lid_driven_cavity/geo/block.004
+++ b/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)