Carpet: Handle granularity in recursive load balancing

3 files changed, 55 insertions, 18 deletions

diff --git a/Carpet/Carpet/src/LoadBalanceReal/carpet_boxtypes.F90 b/Carpet/Carpet/src/LoadBalanceReal/carpet_boxtypes.F90
index a03d312c6..8358406c8 100644
--- a/Carpet/Carpet/src/LoadBalanceReal/carpet_boxtypes.F90
+++ b/Carpet/Carpet/src/LoadBalanceReal/carpet_boxtypes.F90
@@ -8,8 +8,6 @@ module carpet_boxtypes
   integer :: ghostsize
   real(wp) :: alpha
   logical :: limit_size
-  integer :: granularity
-  integer :: granularity_boundary
   integer :: procid
 
   contains
@@ -437,14 +435,11 @@ module carpet_boxtypes
 !       Find the number of points in each of the 2 chunks while making
 !       sure the chunks are not too small compared to the ghostsize if
 !       limit_size is true.
-!       Take also the granularity into account.
 !       At the outer boundary, take also the boundary size into account,
 !       assuming that the boundary size is equal to the ghost size.
 !       (Boundary points cannot be split, and near a boundary, the
 !        minimum number of interior points is the number of ghost points.)
-        np1 = nint((real(maxcost,wp)*p1)/nprocs/granularity)*granularity
-!       TODO: Take granularity_boundary into account as well!
-!       TODO: Take granularity into account when limiting!
+        np1 = nint((real(maxcost,wp)*p1)/nprocs)
         if (limit_size) then
           lower_is_outer = &
                newregarr(i)%point%outer_boundaries%obound(mydim(1),1)/=0
diff --git a/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.F90 b/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.F90
index 441bf2d78..5667c4e30 100644
--- a/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.F90
+++ b/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.F90
@@ -4,7 +4,7 @@ subroutine splitregions_recursively ( &
      cxx_superregs, nsuperregs, &
      cxx_regs, &
      nprocs, &
-     ghostsize_, alpha_, limit_size_, granularity_, granularity_boundary_, &
+     ghostsize_, alpha_, limit_size_, &
      procid_)
   use carpet_boxtypes
   implicit none
@@ -16,8 +16,6 @@ subroutine splitregions_recursively ( &
   integer,      intent(in) :: ghostsize_
   CCTK_REAL,    intent(in) :: alpha_
   integer,      intent(in) :: limit_size_
-  CCTK_INT,     intent(in) :: granularity_
-  CCTK_INT,     intent(in) :: granularity_boundary_
   integer,      intent(in) :: procid_
   
   type(ptr), allocatable :: sregions(:)
@@ -79,12 +77,10 @@ subroutine splitregions_recursively ( &
   
   
   ! Set global parameters
-  ghostsize            = ghostsize_
-  alpha                = alpha_
-  limit_size           = limit_size_ /= 0
-  granularity          = granularity_
-  granularity_boundary = granularity_boundary_
-  procid               = procid_
+  ghostsize  = ghostsize_
+  alpha      = alpha_
+  limit_size = limit_size_ /= 0
+  procid     = procid_
   
   
   
diff --git a/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.cc b/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.cc
index 893545f32..fd3148864 100644
--- a/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.cc
+++ b/Carpet/Carpet/src/LoadBalanceReal/splitregions_recursively.cc
@@ -162,8 +162,6 @@ namespace Carpet {
                                         int const& ghostsize,
                                         CCTK_REAL const& alpha,
                                         int const& limit_size,
-                                        CCTK_INT const& granularity,
-                                        CCTK_INT const& granularity_boundary,
                                         int const& procid);
   
   void
@@ -230,6 +228,23 @@ namespace Carpet {
       superregs.AT(r).map = r;
     }
     
+    // Handle granularity: artificially reduce the size of the
+    // superregions
+    for (int r=0; r<nsuperregs; ++r) {
+      region_t& superreg = superregs.AT(r);
+      ibbox& ext = superregs.AT(r).extent;
+      ivect lo = ext.lower();
+      ivect hi = ext.upper() + ext.stride();
+      ivect str = ext.stride();
+      // cut off outer boundaries
+      lo += ivect(superreg.outer_boundaries[0]) * granularity_boundary * str;
+      hi -= ivect(superreg.outer_boundaries[1]) * granularity_boundary * str;
+      // multiply stride by granularity
+      str *= granularity;
+      assert(all((hi-lo)%str==0));
+      ext = ibbox(lo, hi-str, str);
+    }
+    
     int const real_nprocs = CCTK_nProcs (cctkGH);
     if (recompose_verbose) cout << "SRMR real_nprocs " << real_nprocs << endl;
     
@@ -268,9 +283,40 @@ namespace Carpet {
     int const procid = CCTK_MyProc(cctkGH);
     CCTK_FNAME(splitregions_recursively)
       (cxx_superregs, nsuperregs, cxx_regs, nprocs,
-       ghostsize, alpha, limit_size, granularity, granularity_boundary, procid);
+       ghostsize, alpha, limit_size, procid);
     int nregs = regs.size();
     
+    // Handle granularity: increase the size of the superregions and
+    // regions again
+    for (int r=0; r<nsuperregs; ++r) {
+      region_t& superreg = superregs.AT(r);
+      ibbox& ext = superregs.AT(r).extent;
+      ivect lo = ext.lower();
+      ivect hi = ext.upper() + ext.stride();
+      ivect str = ext.stride();
+      // divide stride by granularity
+      assert(all(str%granularity==0));
+      str /= granularity;
+      // add outer boundaries again
+      lo -= ivect(superreg.outer_boundaries[0]) * granularity_boundary * str;
+      hi += ivect(superreg.outer_boundaries[1]) * granularity_boundary * str;
+      ext = ibbox(lo, hi-str, str);
+    }
+    for (int r=0; r<nregs; ++r) {
+      region_t& reg = regs.AT(r);
+      ibbox& ext = regs.AT(r).extent;
+      ivect lo = ext.lower();
+      ivect hi = ext.upper() + ext.stride();
+      ivect str = ext.stride();
+      // divide stride by granularity
+      assert(all(str%granularity==0));
+      str /= granularity;
+      // add outer boundaries again
+      lo -= ivect(reg.outer_boundaries[0]) * granularity_boundary * str;
+      hi += ivect(reg.outer_boundaries[1]) * granularity_boundary * str;
+      ext = ibbox(lo, hi-str, str);
+    }
+    
     if (same_number_of_components_on_each_process) {
       // Ensure all processes have the same number of components
       vector<int> ncomps(nprocs, 0);