global: Add varying refinement factors

Add support for varying refinement factors.  The spatial refinement
factors can be different in different directions, can be different
from the time refinement factor, and can be different on each level.
(However, the underlying spatial transport operators do currently not
handle any factors except two.)

darcs-hash:20050501205010-891bb-8d3a74abaad55ee6c77ef18d51fca2a2b69740de.gz

7 files changed, 23 insertions, 17 deletions

diff --git a/Carpet/CarpetRegrid/src/automatic.cc b/Carpet/CarpetRegrid/src/automatic.cc
index 987ef133b..82dbd5796 100644
--- a/Carpet/CarpetRegrid/src/automatic.cc
+++ b/Carpet/CarpetRegrid/src/automatic.cc
@@ -111,6 +111,9 @@ namespace CarpetRegrid {
     const int tl = 0;
     const int ml = 0;
     
+    // Determine refinement factor
+    ivect const reffact = spacereffacts.at(rl+1) / spacereffacts.at(rl);
+    
     if (veryverbose) {
       cout << endl << "MRA: Choosing regions to refine on level " << rl << " in " << hh.components(rl) << " components" << endl;
     }
@@ -122,7 +125,7 @@ namespace CarpetRegrid {
       
       const data<CCTK_REAL>& errordata = *errorgf(tl,rl,c,ml);
       
-      Automatic_Recursive (cctkGH, hh, errordata, bbl, region);
+      Automatic_Recursive (cctkGH, hh, errordata, bbl, region, reffact);
     }
     
     if (veryverbose) {
@@ -167,7 +170,8 @@ namespace CarpetRegrid {
                             const gh & hh,
                             const data<CCTK_REAL> & errordata,
                             list<ibbox> & bbl,
-                            const ibbox & region)
+                            const ibbox & region,
+                            const ivect & reffact)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -202,7 +206,7 @@ namespace CarpetRegrid {
     
     if (cnt == 0) {
       // Don't refine
-    } else if (width*reffact < 2*minwidth or fraction >= minfraction) {
+    } else if (any (reffact*width < 2*minwidth) or fraction >= minfraction) {
       // Refine the whole region
       const ivect lo(region.lower());
       const ivect up(region.upper());
@@ -231,8 +235,8 @@ namespace CarpetRegrid {
       assert ((region1 & region2).empty());
       assert (region1 + region2 == region);
       list<ibbox> bbl1, bbl2;
-      Automatic_Recursive (cctkGH, hh, errordata, bbl1, region1);
-      Automatic_Recursive (cctkGH, hh, errordata, bbl2, region2);
+      Automatic_Recursive (cctkGH, hh, errordata, bbl1, region1, reffact);
+      Automatic_Recursive (cctkGH, hh, errordata, bbl2, region2, reffact);
       // Combine regions if possible
       up2 += str-str/reffact;
       up2[dir] = lo2[dir];
diff --git a/Carpet/CarpetRegrid/src/centre.cc b/Carpet/CarpetRegrid/src/centre.cc
index f9b384142..b9a3f918e 100644
--- a/Carpet/CarpetRegrid/src/centre.cc
+++ b/Carpet/CarpetRegrid/src/centre.cc
@@ -53,8 +53,8 @@ namespace CarpetRegrid {
       ivect const oldrub = rub;
       
       // refined boxes have smaller stride
-      assert (all(rstr%hh.reffact == 0));
-      rstr /= hh.reffact;
+      assert (all(rstr%(hh.reffacts.at(rl)/hh.reffacts.at(rl-1)) == 0));
+      rstr /= hh.reffacts.at(rl)/hh.reffacts.at(rl-1);
       
       // calculate new extent
       ivect const quarter = (rub - rlb) / 4 / rstr * rstr;
diff --git a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
index 9e8588bbb..b444549d4 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
@@ -78,7 +78,7 @@ namespace CarpetRegrid {
         for (size_t c=0; c<newobss.at(rl).size(); ++c) {
           for (int d=0; d<dim; ++d) {
             assert (mglevel==0);
-            rvect const spacing = base_spacing * ipow((CCTK_REAL)mgfact, basemglevel) / ipow(reffact, rl+1);
+            rvect const spacing = base_spacing * ipow((CCTK_REAL)mgfact, basemglevel) / spacereffacts.at(rl+1);
             ierr = ConvertFromPhysicalBoundary
               (dim, &physical_min[0], &physical_max[0],
                &interior_min[0], &interior_max[0],
@@ -92,7 +92,7 @@ namespace CarpetRegrid {
               lo[d] = exterior_min[d];
               newbbss.at(rl).at(c) = rbbox(lo, up, str);
             }
-            newobss.at(rl).at(c)[d][1] = abs(newbbss.at(rl).at(c).upper()[d] - physical_max[d]) < 1.0e-6 * base_spacing[d] / ipow(reffact, rl);
+            newobss.at(rl).at(c)[d][1] = abs(newbbss.at(rl).at(c).upper()[d] - physical_max[d]) < 1.0e-6 * base_spacing[d] / spacereffacts.at(rl)[d];
             if (newobss.at(rl).at(c)[d][1]) {
               rvect lo = newbbss.at(rl).at(c).lower();
               rvect up = newbbss.at(rl).at(c).upper();
@@ -156,7 +156,7 @@ namespace CarpetRegrid {
         // assert (domain_from_coordbase);
         // TODO: why can basemglevel not be used here?
         // rvect const spacing = base_spacing * ipow(CCTK_REAL(mgfact), basemglevel) / ipow(reffact, rl);
-        rvect const spacing = base_spacing / ipow(reffact, rl);
+        rvect const spacing = base_spacing / spacereffacts.at(rl);
         if (! all(abs(ext.stride() - spacing) < spacing * 1.0e-10)) {
           assert (dim==3);
           CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
diff --git a/Carpet/CarpetRegrid/src/manualcoordinates.cc b/Carpet/CarpetRegrid/src/manualcoordinates.cc
index e64aea0e4..f5f84c7b3 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinates.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinates.cc
@@ -117,7 +117,7 @@ namespace CarpetRegrid {
     const ivect global_extent (hh.baseextent.upper() - hh.baseextent.lower());
     
     const rvect scale  = rvect(global_extent) / (global_upper - global_lower);
-    const int levfac = ipow(hh.reffact, rl);
+    const ivect levfac = hh.reffacts.at(rl);
     assert (all (hh.baseextent.stride() % levfac == 0));
     const ivect istride = hh.baseextent.stride() / levfac;
     
@@ -147,7 +147,7 @@ namespace CarpetRegrid {
     const ivect global_extent (hh.baseextent.upper() - hh.baseextent.lower());
     
     const rvect scale  = rvect(global_extent) / (global_upper - global_lower);
-    const int levfac = ipow(hh.reffact, rl);
+    const ivect levfac = hh.reffacts.at(rl);
     assert (all (hh.baseextent.stride() % levfac == 0));
     const ivect istride = hh.baseextent.stride() / levfac;
     
diff --git a/Carpet/CarpetRegrid/src/manualgridpoints.cc b/Carpet/CarpetRegrid/src/manualgridpoints.cc
index a85998997..5028d2064 100644
--- a/Carpet/CarpetRegrid/src/manualgridpoints.cc
+++ b/Carpet/CarpetRegrid/src/manualgridpoints.cc
@@ -96,7 +96,7 @@ namespace CarpetRegrid {
     const ivect rlb  = hh.baseextent.lower();
     const ivect rub  = hh.baseextent.upper();
     
-    const int levfac = ipow(hh.reffact, rl);
+    const ivect levfac = hh.reffacts.at(rl);
     assert (all (rstr % levfac == 0));
     const ivect str (rstr / levfac);
     const ivect lb  (ilower);
diff --git a/Carpet/CarpetRegrid/src/moving.cc b/Carpet/CarpetRegrid/src/moving.cc
index 7936865c3..44a5c2a16 100644
--- a/Carpet/CarpetRegrid/src/moving.cc
+++ b/Carpet/CarpetRegrid/src/moving.cc
@@ -49,10 +49,11 @@ namespace CarpetRegrid {
         (moving_centre_x + moving_circle_radius * cos(argument),
          moving_centre_y + moving_circle_radius * sin(argument),
          moving_centre_z);
-      CCTK_REAL const radius = moving_region_radius / ipow(reffact, rl-1);
+      rvect const radius
+        (rvect(moving_region_radius) / spacereffacts.at(rl-1));
       
-      rvect const rlb (symmetric.ifthen (rvect(0),      pos - rvect(radius)));
-      rvect const rub (symmetric.ifthen (rvect(radius), pos + rvect(radius)));
+      rvect const rlb (symmetric.ifthen (rvect(0), pos - radius));
+      rvect const rub (symmetric.ifthen (radius  , pos + radius));
       
       vector<ibbox> bbs;
       gh::cbnds obs;
diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh
index fa0dd59f1..32710087e 100644
--- a/Carpet/CarpetRegrid/src/regrid.hh
+++ b/Carpet/CarpetRegrid/src/regrid.hh
@@ -131,7 +131,8 @@ namespace CarpetRegrid {
                             const gh & hh,
                             const data<CCTK_REAL> & errorvar,
                             list<ibbox> & bbl,
-                            const ibbox & region);
+                            const ibbox & region,
+                            const ivect & reffact);
   
   void Automatic_Recombine (list<ibbox> & bbl1,
                             list<ibbox> & bbl2,