CarpetRegrid2: Make it possible to regrid all maps at once

darcs-hash:20060804133836-dae7b-10f0524d418ba823dbd5c9db27146347c58ba7c8.gz

1 files changed, 331 insertions, 215 deletions

diff --git a/Carpet/CarpetRegrid2/src/regrid.cc b/Carpet/CarpetRegrid2/src/regrid.cc
index e1a9174a7..f133d043b 100644
--- a/Carpet/CarpetRegrid2/src/regrid.cc
+++ b/Carpet/CarpetRegrid2/src/regrid.cc
@@ -157,6 +157,238 @@ namespace CarpetRegrid2 {
                           CCTK_POINTER          const obss_,
                           CCTK_POINTER          const pss_,
                           CCTK_INT              const force);
+    
+    CCTK_INT
+    CarpetRegrid2_RegridMaps (CCTK_POINTER_TO_CONST const cctkGH_,
+                              CCTK_POINTER          const bbssss_,
+                              CCTK_POINTER          const obsss_,
+                              CCTK_POINTER          const psss_,
+                              CCTK_INT              const force);
+  }
+  
+  
+  
+  void
+  Regrid (cGH const * const cctkGH,
+          gh::rexts & bbss,
+          gh::rbnds & obss)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    assert (Carpet::is_singlemap_mode());
+    gh const & hh = * Carpet::vhh.at (Carpet::map);
+    dh const & dd = * Carpet::vdd.at (Carpet::map);
+    
+    //
+    // Find extent of domain
+    //
+    
+    // This requires that CoordBase is used
+#warning "TODO: check this (for Carpet, and maybe also for CartGrid3D)"
+#warning "TODO: (the check that these two are consistent should be in Carpet)"
+    
+    typedef vect<vect<CCTK_INT,2>,3> jjvect;
+    jjvect nboundaryzones;
+    jjvect is_internal;
+    jjvect is_staggered;
+    jjvect shiftout;
+    {
+      CCTK_INT const ierr = GetBoundarySpecification
+        (2*dim,
+         & nboundaryzones[0][0],
+         & is_internal[0][0],
+         & is_staggered[0][0],
+         & shiftout[0][0]);
+      assert (not ierr);
+    }
+    
+    rvect physical_lower, physical_upper;
+    rvect interior_lower, interior_upper;
+    rvect exterior_lower, exterior_upper;
+    rvect spacing;
+    {
+      CCTK_INT const ierr = GetDomainSpecification
+        (dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    }
+    
+    // Adapt spacing for convergence level
+    spacing *= ipow ((CCTK_REAL) Carpet::mgfact, Carpet::basemglevel);
+    
+    {
+      CCTK_INT const ierr =
+        ConvertFromPhysicalBoundary
+        (dim,
+         & physical_lower[0], & physical_upper[0],
+         & interior_lower[0], & interior_upper[0],
+         & exterior_lower[0], & exterior_upper[0],
+         & spacing[0]);
+      assert (not ierr);
+    }
+    
+    //
+    // Calculate the union of the bounding boxes for all levels
+    //
+    
+    rvect const origin (exterior_lower);
+    rvect const scale (rvect (hh.baseextent.stride()) / spacing);
+    
+    ivect const physical_ilower =
+      rpos2ipos (physical_lower, origin, scale, hh, 0);
+    ivect const physical_iupper =
+      rpos2ipos1 (physical_upper, origin, scale, hh, 0);
+    
+    // The set of refined regions
+    vector <ibboxset> regions (1);
+    
+    // Loop over all centres
+    for (int n = 0; n < centre_description::num_centres(); ++ n) {
+      centre_description centre (n);
+      
+      // Loop over all levels for this centre
+      for (int rl = 1; rl < centre.num_levels; ++ rl) {
+        
+        // Calculate a bbox for this region
+        rvect const rmin = centre.position - centre.radius.at(rl);
+        rvect const rmax = centre.position + centre.radius.at(rl);
+        
+        // Convert to an integer bbox
+        ivect const imin =
+          rpos2ipos (rmin, origin, scale, hh, rl);
+        ivect const imax =
+          rpos2ipos1 (rmax, origin, scale, hh, rl);
+        
+        ivect const levfac = hh.reffacts.at(rl);
+        assert (all (hh.baseextent.stride() % levfac == 0));
+        ivect const istride = hh.baseextent.stride() / levfac;
+        
+        ibbox const region (imin, imax, istride);
+        
+        // Add this region to the list of regions
+        if (static_cast <int> (regions.size()) < rl+1) regions.resize (rl+1);
+        regions.at(rl) |= region;
+        
+      } // for rl
+      
+    } // for n
+    
+      // Ensure that the coarser grids contain the finer ones
+    for (size_t rl = regions.size() - 1; rl >= 2; -- rl) {
+      ibbox coarse = * regions.at(rl-1).begin();
+      
+      regions.at(rl).normalize();
+      ibboxset coarsified;
+      for (ibboxset::const_iterator ibb = regions.at(rl).begin();
+           ibb != regions.at(rl).end();
+           ++ ibb)
+      {
+        ivect const distance (min_distance);
+        ibbox const fbb = * ibb;
+        ibbox const cbb = fbb.expanded_for(coarse);
+        ibbox const ebb = cbb.expand (distance, distance);
+        coarsified |= ebb;
+      }
+      
+      regions.at(rl-1) |= coarsified;
+    }
+    
+    
+    
+    //
+    // Clip at the outer boundary, and convert to (bbss, obss) pair
+    //
+    
+    bbss.resize (regions.size());
+    obss.resize (regions.size());
+    
+    for (size_t rl = 1; rl < regions.size(); ++ rl) {
+      
+      // Find the location of the outer boundary
+      rvect const level_physical_lower = physical_lower;
+      rvect const level_physical_upper = physical_upper;
+      rvect const level_spacing = spacing / rvect (hh.reffacts.at(rl));
+      rvect level_interior_lower, level_interior_upper;
+      rvect level_exterior_lower, level_exterior_upper;
+      {
+        CCTK_INT const ierr =
+          ConvertFromPhysicalBoundary
+          (dim,
+           & level_physical_lower[0], & level_physical_upper[0],
+           & level_interior_lower[0], & level_interior_upper[0],
+           & level_exterior_lower[0], & level_exterior_upper[0],
+           & level_spacing[0]);
+        assert (not ierr);
+      }
+      
+      ivect const level_exterior_ilower =
+        rpos2ipos (level_exterior_lower, origin, scale, hh, rl);
+      ivect const level_exterior_iupper =
+        rpos2ipos1 (level_exterior_upper, origin, scale, hh, rl);
+      
+      // Find the minimum necessary distance to the outer boundary
+      // due to buffer zones.  This is in terms of grid points.
+      i2vect const min_bnd_dist = dd.buffers;
+      
+      // Clip at the outer boundary
+      regions.at(rl).normalize();
+      ibboxset clipped;
+      for (ibboxset::const_iterator ibb = regions.at(rl).begin();
+           ibb != regions.at(rl).end();
+           ++ ibb)
+      {
+        ibbox const bb = * ibb;
+        
+        // Clip boxes that extend outside the boundary.  Enlarge boxes
+        // that are inside but too close to the outer boundary.
+        bvect const lower_is_outer =
+          bb.lower() - min_bnd_dist[0] <= physical_ilower;
+        bvect const upper_is_outer =
+          bb.upper() + min_bnd_dist[1] >= physical_iupper;
+        
+        ibbox const clipped_bb
+          (either (lower_is_outer, level_exterior_ilower, bb.lower()),
+           either (upper_is_outer, level_exterior_iupper, bb.upper()),
+           bb.stride());
+        
+        clipped += clipped_bb;
+      }
+      
+      regions.at(rl) = clipped;
+      
+      // Simplify the set of refined regions
+      regions.at(rl).normalize();
+      
+      // Create a vector of bboxes for this level
+      gh::cexts bbs;
+      gh::cbnds obs;
+      bbs.reserve (regions.at(rl).setsize());
+      obs.reserve (regions.at(rl).setsize());
+      for (ibboxset::const_iterator ibb = regions.at(rl).begin();
+           ibb != regions.at(rl).end();
+           ++ ibb)
+      {
+        ibbox bb = * ibb;
+        
+        bvect const lower_is_outer = bb.lower() <= physical_ilower;
+        bvect const upper_is_outer = bb.upper() >= physical_iupper;
+        
+        bbvect ob;
+        for (int d = 0; d < dim; ++ d) {
+          ob[d][0] = lower_is_outer[d];
+          ob[d][1] = upper_is_outer[d];
+        }
+        bbs.push_back (bb);
+        obs.push_back (ob);
+      }
+      
+      bbss.at(rl) = bbs;
+      obss.at(rl) = obs;
+      
+    } // for rl
   }
   
   
@@ -171,13 +403,8 @@ namespace CarpetRegrid2 {
     DECLARE_CCTK_PARAMETERS;
     
     cGH const * const cctkGH = static_cast <cGH const *> (cctkGH_);
-    gh::mexts  & bbsss = * static_cast <gh::mexts  *> (bbsss_);
-    gh::rbnds  & obss  = * static_cast <gh::rbnds  *> (obss_ );
-    gh::rprocs & pss   = * static_cast <gh::rprocs *> (pss_  );
     
     assert (Carpet::is_singlemap_mode());
-    gh const & hh = * Carpet::vhh.at (Carpet::map);
-    dh const & dd = * Carpet::vdd.at (Carpet::map);
     
     // Decide whether to change the grid hierarchy
     // (We always do)
@@ -191,238 +418,127 @@ namespace CarpetRegrid2 {
         do_recompose = cctkGH->cctk_iteration == 0;
       } else {
         do_recompose =
-          cctkGH->cctk_iteration == 0 or
-          (cctkGH->cctk_iteration > 0 and
-           (cctkGH->cctk_iteration - 1) % regrid_every == 0);
+          reflevel == 0 and
+          (cctkGH->cctk_iteration == 0 or
+           (cctkGH->cctk_iteration > 0 and
+            (cctkGH->cctk_iteration - 1) % regrid_every == 0));
       }
     }
     
     if (do_recompose) {
       
-      //
-      // Find extent of domain
-      //
+      gh::mexts  & bbsss = * static_cast <gh::mexts  *> (bbsss_);
+      gh::rbnds  & obss  = * static_cast <gh::rbnds  *> (obss_ );
+      gh::rprocs & pss   = * static_cast <gh::rprocs *> (pss_  );
       
-      // This requires that CoordBase is used (but this is not
-      // checked)
+      // Make multigrid unaware
+      vector <vector <ibbox> > bbss = bbsss.at(0);
       
-      typedef vect<vect<CCTK_INT,2>,3> jjvect;
-      jjvect nboundaryzones;
-      jjvect is_internal;
-      jjvect is_staggered;
-      jjvect shiftout;
-      {
-        CCTK_INT const ierr = GetBoundarySpecification
-          (2*dim,
-           & nboundaryzones[0][0],
-           & is_internal[0][0],
-           & is_staggered[0][0],
-           & shiftout[0][0]);
-        assert (not ierr);
-      }
+      Regrid (cctkGH, bbss, obss);
       
-      rvect physical_lower, physical_upper;
-      rvect interior_lower, interior_upper;
-      rvect exterior_lower, exterior_upper;
-      rvect spacing;
-      {
-        CCTK_INT const ierr = GetDomainSpecification
-          (dim,
-           & physical_lower[0], & physical_upper[0],
-           & interior_lower[0], & interior_upper[0],
-           & exterior_lower[0], & exterior_upper[0],
-           & spacing[0]);
-        assert (not ierr);
-      }
+      // Make multiprocessor aware
+      pss.resize (bbss.size());
+      for (size_t rl = 0; rl < pss.size(); ++ rl) {
+        Carpet::SplitRegions (cctkGH, bbss.at(rl), obss.at(rl), pss.at(rl));
+      } // for rl
       
-      // Adapt spacing for convergence level
-      spacing *= ipow ((CCTK_REAL) Carpet::mgfact, Carpet::basemglevel);
+      // Make multigrid aware
+      Carpet::MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
       
-      {
-        CCTK_INT const ierr =
-          ConvertFromPhysicalBoundary
-          (dim,
-           & physical_lower[0], & physical_upper[0],
-           & interior_lower[0], & interior_upper[0],
-           & exterior_lower[0], & exterior_upper[0],
-           & spacing[0]);
-        assert (not ierr);
+    } // if do_recompose
+    
+    return do_recompose;
+  }
+  
+  
+  
+  CCTK_INT
+  CarpetRegrid2_RegridMaps (CCTK_POINTER_TO_CONST const cctkGH_,
+                            CCTK_POINTER          const bbssss_,
+                            CCTK_POINTER          const obsss_,
+                            CCTK_POINTER          const psss_,
+                            CCTK_INT              const force)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    cGH const * const cctkGH = static_cast <cGH const *> (cctkGH_);
+    
+    assert (Carpet::is_level_mode());
+    
+    // Decide whether to change the grid hierarchy
+    // (We always do)
+    bool do_recompose;
+    if (force) {
+      do_recompose = true;
+    } else {
+      if (regrid_every == -1) {
+        do_recompose = false;
+      } else if (regrid_every == 0) {
+        do_recompose = cctkGH->cctk_iteration == 0;
+      } else {
+        do_recompose =
+          reflevel == 0 and
+          (cctkGH->cctk_iteration == 0 or
+           (cctkGH->cctk_iteration > 0 and
+            (cctkGH->cctk_iteration - 1) % regrid_every == 0));
       }
+    }
+    
+    if (do_recompose) {
+      vector <gh::mexts>  & bbssss =
+        * static_cast <vector <gh::mexts>  *> (bbssss_);
+      vector <gh::rbnds>  & obsss  =
+        * static_cast <vector <gh::rbnds>  *> (obsss_ );
+      vector <gh::rprocs> & psss   =
+        * static_cast <vector <gh::rprocs> *> (psss_  );
       
-      //
-      // Calculate the union of the bounding boxes for all levels
-      //
-      
-      rvect const origin (exterior_lower);
-      rvect const scale (rvect (hh.baseextent.stride()) / spacing);
-      
-      ivect const physical_ilower =
-        rpos2ipos (physical_lower, origin, scale, hh, 0);
-      ivect const physical_iupper =
-        rpos2ipos1 (physical_upper, origin, scale, hh, 0);
-      
-      // The set of refined regions
-      vector <ibboxset> regions (1);
-      
-      // Loop over all centres
-      for (int n = 0; n < centre_description::num_centres(); ++ n) {
-        centre_description centre (n);
-        
-        // Loop over all levels for this centre
-        for (int rl = 1; rl < centre.num_levels; ++ rl) {
-          
-          // Calculate a bbox for this region
-          rvect const rmin = centre.position - centre.radius.at(rl);
-          rvect const rmax = centre.position + centre.radius.at(rl);
-          
-          // Convert to an integer bbox
-          ivect const imin =
-            rpos2ipos (rmin, origin, scale, hh, rl);
-          ivect const imax =
-            rpos2ipos1 (rmax, origin, scale, hh, rl);
-          
-          ivect const levfac = hh.reffacts.at(rl);
-          assert (all (hh.baseextent.stride() % levfac == 0));
-          ivect const istride = hh.baseextent.stride() / levfac;
-          
-          ibbox const region (imin, imax, istride);
-          
-          // Add this region to the list of regions
-          if (static_cast <int> (regions.size()) < rl+1) regions.resize (rl+1);
-          regions.at(rl) |= region;
-          
-        } // for rl
-        
-      } // for n
-      
-      // Ensure that the coarser grids contain the finer ones
-      for (size_t rl = regions.size() - 1; rl >= 2; -- rl) {
-        ibbox coarse = * regions.at(rl-1).begin();
-        
-        regions.at(rl).normalize();
-        ibboxset coarsified;
-        for (ibboxset::const_iterator ibb = regions.at(rl).begin();
-             ibb != regions.at(rl).end();
-             ++ ibb)
-        {
-          ivect const distance (min_distance);
-          ibbox const fbb = * ibb;
-          ibbox const cbb = fbb.expanded_for(coarse);
-          ibbox const ebb = cbb.expand (distance, distance);
-          coarsified |= ebb;
-        }
-        
-        regions.at(rl-1) |= coarsified;
+      // Make multigrid unaware
+      vector< vector <vector <ibbox> > > bbsss (maps);
+      for (int m = 0; m < maps; ++ m) {
+        bbsss.at(m) = bbssss.at(m).at(0);
       }
       
+      for (int m = 0; m < maps; ++ m) {
+        Regrid (cctkGH, bbsss.at(m), obsss.at(m));
+      }
       
+      // Count levels
+      vector <int> rls (maps);
+      for (int m = 0; m < maps; ++ m) {
+        rls.at(m) = bbsss.at(m).size();
+      }
+      int maxrl = 0;
+      for (int m = 0; m < maps; ++ m) {
+        maxrl = max (maxrl, rls.at(m));
+      }
       
-      //
-      // Clip at the outer boundary, and convert to (bbsss, obss, pss)
-      // triplet
-      //
-      
-      vector <vector <ibbox> > bbss;
-      
-      bbss.resize (regions.size());
-      obss.resize (regions.size());
-      pss .resize (regions.size());
-      
-      bbss.at(0) = bbsss.at(0).at(0);
-      
-      for (size_t rl = 1; rl < regions.size(); ++ rl) {
-        
-        // Find the location of the outer boundary
-        rvect const level_physical_lower = physical_lower;
-        rvect const level_physical_upper = physical_upper;
-        rvect const level_spacing = spacing / rvect (hh.reffacts.at(rl));
-        rvect level_interior_lower, level_interior_upper;
-        rvect level_exterior_lower, level_exterior_upper;
-        {
-          CCTK_INT const ierr =
-            ConvertFromPhysicalBoundary
-            (dim,
-             & level_physical_lower[0], & level_physical_upper[0],
-             & level_interior_lower[0], & level_interior_upper[0],
-             & level_exterior_lower[0], & level_exterior_upper[0],
-             & level_spacing[0]);
-          assert (not ierr);
-        }
-        
-        ivect const level_exterior_ilower =
-          rpos2ipos (level_exterior_lower, origin, scale, hh, rl);
-        ivect const level_exterior_iupper =
-          rpos2ipos1 (level_exterior_upper, origin, scale, hh, rl);
-        
-        // Find the minimum necessary distance to the outer boundary
-        // due to buffer zones.  This is in terms of grid points.
-        i2vect const min_bnd_dist = dd.buffers;
-        
-        // Clip at the outer boundary
-        regions.at(rl).normalize();
-        ibboxset clipped;
-        for (ibboxset::const_iterator ibb = regions.at(rl).begin();
-             ibb != regions.at(rl).end();
-             ++ ibb)
-        {
-          ibbox const bb = * ibb;
-          
-          // Clip boxes that extend outside the boundary.  Enlarge
-          // boxes that are inside but too close to the outer
-          // boundary.
-          bvect const lower_is_outer =
-            bb.lower() - min_bnd_dist[0] <= physical_ilower;
-          bvect const upper_is_outer =
-            bb.upper() + min_bnd_dist[1] >= physical_iupper;
-          
-          ibbox const clipped_bb
-            (either (lower_is_outer, level_exterior_ilower, bb.lower()),
-             either (upper_is_outer, level_exterior_iupper, bb.upper()),
-             bb.stride());
-          
-          clipped += clipped_bb;
+      // Make multiprocessor aware
+      for (int m = 0; m < maps; ++ m) {
+        psss.at(m).resize (bbsss.at(m).size());
+      }
+      for (int rl = 0; rl < maxrl; ++ rl) {
+        vector <vector <ibbox > > bbss (maps);
+        vector <vector <bbvect> > obss (maps);
+        vector <vector <int   > > pss  (maps);
+        for (int m = 0; m < maps; ++ m) {
+          if (rl < rls.at(m)) {
+            bbss.at(m) = bbsss.at(m).at(rl);
+            obss.at(m) = obsss.at(m).at(rl);
+            pss .at(m) = psss .at(m).at(rl);
+          }
         }
-        
-        regions.at(rl) = clipped;
-        
-        // Simplify the set of refined regions
-        regions.at(rl).normalize();
-        
-        // Create a vector of bboxes for this level
-        gh::cexts bbs;
-        gh::cbnds obs;
-        bbs.reserve (regions.at(rl).setsize());
-        obs.reserve (regions.at(rl).setsize());
-        for (ibboxset::const_iterator ibb = regions.at(rl).begin();
-             ibb != regions.at(rl).end();
-             ++ ibb)
-        {
-          ibbox bb = * ibb;
-          
-          bvect const lower_is_outer = bb.lower() <= physical_ilower;
-          bvect const upper_is_outer = bb.upper() >= physical_iupper;
-          
-          bbvect ob;
-          for (int d = 0; d < dim; ++ d) {
-            ob[d][0] = lower_is_outer[d];
-            ob[d][1] = upper_is_outer[d];
+        Carpet::SplitRegionsMaps (cctkGH, bbss, obss, pss);
+        for (int m = 0; m < maps; ++ m) {
+          if (rl < rls.at(m)) {
+            bbsss.at(m).at(rl) = bbss.at(m);
+            obsss.at(m).at(rl) = obss.at(m);
+            psss .at(m).at(rl) = pss .at(m);
           }
-          bbs.push_back (bb);
-          obs.push_back (ob);
         }
-        
-        // Make multiprocessor aware
-        gh::cprocs ps;
-        Carpet::SplitRegions (cctkGH, bbs, obs, ps);
-        
-        bbss.at(rl) = bbs;
-        obss.at(rl) = obs;
-        pss .at(rl) = ps;
-        
       } // for rl
       
       // Make multigrid aware
-      Carpet::MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+      Carpet::MakeMultigridBoxesMaps (cctkGH, bbsss, obsss, bbssss);
       
     } // if do_recompose