CarpetAdaptiveRegrid: parallel support

Make CAR work on more than one processor (that is, I've tested it on
my laptop using mpirun -np 2...).

Note that at present whilst it runs and gives reasonable answers you
do _not_ get the same answers as you would on 1 proc. This is a bug,
as yet unfound. But, there probably are lots of other bugs in this
code still to be found...

darcs-hash:20050224223702-0ff1f-f766381f5309c747ea438dbd10976a99c86abd48.gz

3 files changed, 105 insertions, 21 deletions

diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
index 331ca504b..d20a3e677 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
@@ -27,6 +27,8 @@ namespace CarpetAdaptiveRegrid {
 
 
   static gh::mexts local_bbsss;
+  static gh::rbnds local_obss;
+  
   static CCTK_INT last_iteration = -1;  
   
   using namespace std;
@@ -75,10 +77,23 @@ namespace CarpetAdaptiveRegrid {
     if (local_bbsss.empty()) { // It's the first call
       // Is this really the right thing to do on 
       // multiprocessors?
-      local_bbsss = bbsss;
+      //      local_bbsss = bbsss;
+      const ibbox& baseext = 
+        vdd.at(Carpet::map)->bases.at(mglevel).at(reflevel).exterior;
+      vector<ibbox> tmp_bbs;
+      tmp_bbs.push_back (baseext);
+      vector<bbvect> tmp_obs;
+      tmp_obs.push_back (bbvect(true));
+      vector<vector<ibbox> > tmp_bbss(1);
+      vector<vector<bbvect> > tmp_obss(1);
+      tmp_bbss.at(0) = tmp_bbs;
+      tmp_obss.at(0) = tmp_obs;
+      MakeMultigridBoxes(cctkGH, tmp_bbss, tmp_obss, local_bbsss);
+      local_obss = tmp_obss;
       last_iteration = cctkGH->cctk_iteration;
       CCTK_INT do_recompose = 
-        ManualCoordinateList (cctkGH, hh, bbsss, obss, pss, local_bbsss);
+        ManualCoordinateList (cctkGH, hh, bbsss, obss, pss, 
+                              local_bbsss, local_obss);
 
       if (verbose) {
         ostringstream buf;
@@ -148,6 +163,8 @@ namespace CarpetAdaptiveRegrid {
     CCTK_INT do_recompose;
     do_recompose = 1;
 
+    CCTK_INT sum_handle = CCTK_ReductionArrayHandle("sum");
+
     CCTK_INT called_on_ml = mglevel;
     CCTK_INT called_on_rl = reflevel;
     CCTK_INT called_on_map = carpetGH.map;
@@ -165,7 +182,8 @@ namespace CarpetAdaptiveRegrid {
 
       if (verbose) {
         ostringstream buf;
-        buf << "Entering level " << rl << " map " << called_on_map;
+        buf << "Entering level " << rl << " (of " << finest_current_rl 
+            << "), map " << called_on_map;
         CCTK_INFO(buf.str().c_str());
       }
       
@@ -193,7 +211,8 @@ namespace CarpetAdaptiveRegrid {
       
       stack<ibbox> final;
       
-      vector<vector<ibbox> > bbss = local_bbsss.at(0);
+      vector<vector<ibbox> > bbss = bbsss.at(0);
+      vector<vector<ibbox> > local_bbss = local_bbsss.at(0);
       
       bool did_regrid = false;
             
@@ -225,7 +244,8 @@ namespace CarpetAdaptiveRegrid {
           CCTK_INFO(buf.str().c_str());
         }
         
-        vector<CCTK_INT> mask(prod(bb.shape()/bb.stride()), 0);
+        vector<CCTK_INT> local_mask(prod(bb.shape()/bb.stride()), 0),
+          mask(prod(bb.shape()/bb.stride()), 0);
         
         if (veryverbose) {
           ostringstream buf;
@@ -292,7 +312,7 @@ namespace CarpetAdaptiveRegrid {
                     CCTK_INT mindex = ii + 
                       (imax[0] - imin[0] + 1)*
                       (jj + (imax[1] - imin[1] + 1) * kk);
-                    mask[mindex] = 1;
+                    local_mask[mindex] = 1;
                     if (veryverbose) {
                       CCTK_VInfo(CCTK_THORNSTRING, "In error at point"
                                  "\n(%g,%g,%g) [%d,%d,%d] [[%d,%d,%d]]",
@@ -311,7 +331,7 @@ namespace CarpetAdaptiveRegrid {
         // Instead check the error on child level, if exists
         // This should fix the "orphaned grandchild" problem
         
-        if (bbss.size() > reflevel+1) {
+        if (local_bbss.size() > reflevel+1) {
                 
           CCTK_INT currentml = mglevel;
           CCTK_INT currentrl = reflevel;
@@ -376,7 +396,7 @@ namespace CarpetAdaptiveRegrid {
                       CCTK_INT mindex = ii + 
                         (imax[0] - imin[0] + 1)*
                         (jj + (imax[1] - imin[1] + 1) * kk);
-                      mask[mindex] = 1;
+                      local_mask[mindex] = 1;
                       if (veryverbose) {
                         CCTK_VInfo(CCTK_THORNSTRING, "In error at point"
                                    "\n(%g,%g,%g) [%d,%d,%d] [[%d,%d,%d]]",
@@ -400,11 +420,38 @@ namespace CarpetAdaptiveRegrid {
 
         }
         
-        // FIXME
-        
         // Reduce errors (MPI sum)
+
+        CCTK_INT ierr = 
+          CCTK_ReduceLocArrayToArray1D (cctkGH,
+                                        -1,
+                                        sum_handle,
+                                        &local_mask.front(),
+                                        &mask.front(),
+                                        local_mask.size(),
+                                        CCTK_VARIABLE_INT );
+        
+        if (ierr) {
+          ostringstream buf;
+          buf << "The reduction on level " << reflevel << " failed "
+              << "with error " << ierr;
+          CCTK_WARN(0, buf.str().c_str());
+        }
+
         // Set errors to 0/1
         
+        for (CCTK_INT k = 0; k < imax[2] - imin[2] + 1; k++) {
+          for (CCTK_INT j = 0; j < imax[1] - imin[1] + 1; j++) {
+            for (CCTK_INT i = 0; i < imax[0] - imin[0] + 1; i++) {
+              CCTK_INT index =  i + 
+                (imax[0] - imin[0] + 1)*(j + (imax[1] - imin[1] + 1) * k);
+              if (mask[index]) {
+                mask[index] = 1;
+              }
+            }
+          }
+        }
+        
         // Pad the errors: stage 1 - buffer points marked as 2.
         
         for (CCTK_INT k = 0; k < imax[2] - imin[2] + 1; k++) {
@@ -792,13 +839,15 @@ namespace CarpetAdaptiveRegrid {
           if (verbose) {
             CCTK_INFO("Adding new refinement level");
           }
+          local_bbss.resize(reflevel+2);
           bbss.resize(reflevel+2);
+          local_obss.resize(reflevel+2);
           obss.resize(reflevel+2);
           pss.resize(reflevel+2);
         }
-        bbss.at(reflevel+1) = bbs;
-        obss.at(reflevel+1) = obs;
-        MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss);
+        local_bbss.at(reflevel+1) = bbs;
+        local_obss.at(reflevel+1) = obs;
+        MakeMultigridBoxes (cctkGH, local_bbss, local_obss, local_bbsss);
         
         // make multiprocessor aware
         gh::cprocs ps;
@@ -811,15 +860,17 @@ namespace CarpetAdaptiveRegrid {
       } // did_regrid?
       else
       {
-        if (bbss.size() > reflevel+1) {
+        if (local_bbss.size() > reflevel+1) {
           if (verbose) {
             CCTK_INFO("Removing refinement level");
           }
         }
+        local_bbss.resize(reflevel+1);
         bbss.resize(reflevel+1);
+        local_obss.resize(reflevel+1);
         obss.resize(reflevel+1);
         // Set local bbsss
-        MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss);
+        MakeMultigridBoxes (cctkGH, local_bbss, local_obss, local_bbsss);
         
         pss.resize(reflevel+1);
         
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh
index ca91b9511..eb6011da4 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh
@@ -45,7 +45,9 @@ namespace CarpetAdaptiveRegrid {
                             gh::mexts  & bbsss,
                             gh::rbnds  & obss,
                             gh::rprocs & pss,
-                            gh::mexts  & local_bbsss);
+                            gh::mexts  & local_bbsss,
+                            gh::rbnds  & local_obss);
+
     
   void ManualCoordinates_OneLevel (const cGH * const cctkGH,
                                    const gh & hh,
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc b/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc
index 20e5dfc0e..523c5cf9b 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc
@@ -26,7 +26,8 @@ namespace CarpetAdaptiveRegrid {
                             gh::mexts  & bbsss,
                             gh::rbnds  & obss,
                             gh::rprocs & pss,
-                            gh::mexts  & local_bbsss)
+                            gh::mexts  & local_bbsss,
+                            gh::rbnds  & local_obss)
   {
     DECLARE_CCTK_PARAMETERS;
     int ierr;
@@ -38,6 +39,7 @@ namespace CarpetAdaptiveRegrid {
     
     assert (bbsss.size() >= 1);
     vector<vector<ibbox> > bbss = bbsss.at(0);
+    vector<vector<ibbox> > local_bbss = local_bbsss.at(0);
     
     jjvect nboundaryzones, is_internal, is_staggered, shiftout;
     ierr = GetBoundarySpecification
@@ -55,7 +57,9 @@ namespace CarpetAdaptiveRegrid {
     assert (!ierr);
     
     bbss.resize (refinement_levels);
+    local_bbss.resize (refinement_levels);
     obss.resize (refinement_levels);
+    local_obss.resize (refinement_levels);
     pss.resize (refinement_levels);
     
     vector<vector<rbbox> > newbbss;
@@ -136,17 +140,26 @@ namespace CarpetAdaptiveRegrid {
            ext.lower(), ext.upper(), ob, bbs, obs);
       }
 
-      bbss.at(rl) = bbs;
-      obss.at(rl) = obs;
+      local_bbss.at(rl) = bbs;
+      local_obss.at(rl) = obs;
+
+      if (verbose) {
+        ostringstream buf;
+        buf << "About to make multigrid aware for first time whilst doing " 
+            << "manual coordinate list, level " << rl 
+            << ". Total list is:"
+            << endl << bbss << endl << "with obss being "<< local_obss;
+        CCTK_INFO(buf.str().c_str());
+      }      
       
       // make multigrid aware
-      MakeMultigridBoxes (cctkGH, bbss, obss, local_bbsss);
+      MakeMultigridBoxes (cctkGH, local_bbss, local_obss, local_bbsss);
 
       if (verbose) {
         ostringstream buf;
         buf << "Doing manual coordinate list, level " << rl 
             << ". Total list is:"
-            << endl << local_bbsss;
+            << endl << local_bbsss << endl << "with obss being "<< local_obss;
         CCTK_INFO(buf.str().c_str());
       }      
       
@@ -157,10 +170,28 @@ namespace CarpetAdaptiveRegrid {
       bbss.at(rl) = bbs;
       obss.at(rl) = obs;
       pss.at(rl) = ps;
+
+      if (verbose) {
+        ostringstream buf;
+        buf << "Doing manual coordinate list, multiprocs, level " << rl 
+            << "." << endl << "Total list is:"
+            << endl << bbss << endl << "with obss being " << obss 
+            << endl << "Sizes are (bb) " << bbss.at(rl).size() 
+            << " and (ob) " << obss.at(rl).size();
+        CCTK_INFO(buf.str().c_str());
+      }      
               
       // make multigrid aware
       MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
 
+      if (verbose) {
+        ostringstream buf;
+        buf << "Did manual coordinate list, multiprocs, level " << rl 
+            << "." << endl << "Total list is:"
+            << endl << bbsss;
+        CCTK_INFO(buf.str().c_str());
+      }      
+
     } // for rl
     
     return 1;