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

25 files changed, 229 insertions, 92 deletions

diff --git a/Carpet/Carpet/param.ccl b/Carpet/Carpet/param.ccl
index b98240da2..d17b4851e 100644
--- a/Carpet/Carpet/param.ccl
+++ b/Carpet/Carpet/param.ccl
@@ -101,6 +101,22 @@ CCTK_INT refinement_factor "Refinement factor"
   1:* :: "must be positive"
 } 2
 
+STRING space_refinement_factors "Spatial refinement factors over the coarsest level"
+{
+  "^$" :: "Use the value of refinement_factor"
+  # V = [SDS(,SDS)*]
+  # L = [SVS(,SVS)*]
+  #   = [S[SDS(,SDS)*]S(,S[SDS(,SDS)*]S)*]
+  "^\[[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*(,[[:space:]]*[[:digit:]]+[[:space:]]*)*\][[:space:]]*(,[[:space:]]*\[[[:space:]]*[[:digit:]]+[[:space:]]*(,[[:space:]]*[[:digit:]]+[[:space:]]*)*\][[:space:]]*)*\]$" :: "[ [<ifact>,<jfact>,<kfact>], ... ]"
+} ""
+
+STRING time_refinement_factors "Temporal refinement factors over the coarsest level"
+{
+  "^$" :: "Use the value of refinement_factor"
+  # L = [SDS(,SDS)*]
+  "^\[[[:space:]]*[[:digit:]]+[[:space:]]*(,[[:space:]]*[[:digit:]]+[[:space:]]*)*\]$" :: "[ <tfact>, ... ]"
+} ""
+
 
 
 CCTK_INT convergence_level "Convergence level"
diff --git a/Carpet/Carpet/src/Evolve.cc b/Carpet/Carpet/src/Evolve.cc
index 875f81417..7fff9df97 100644
--- a/Carpet/Carpet/src/Evolve.cc
+++ b/Carpet/Carpet/src/Evolve.cc
@@ -46,7 +46,8 @@ namespace Carpet {
     bool term;
 
     // Early return for non-active reflevels to save the MPI_Allreduce() below
-    if (iteration % (maxreflevelfact / ipow(reffact, reflevels-1)) != 0) {
+    if (iteration % (maxtimereflevelfact / timereffacts.at(reflevels-1)) != 0)
+    {
 
       return false;
 
@@ -139,7 +140,7 @@ namespace Carpet {
       AdvanceTime( cgh, cctk_initial_time );
 
       if ((cgh->cctk_iteration-1)
-          % (maxreflevelfact / ipow(reffact, reflevels-1)) == 0) {
+          % (maxtimereflevelfact / timereffacts.at(reflevels-1)) == 0) {
         Waypoint ("Evolving iteration %d at t=%g",
                   cgh->cctk_iteration, (double)cgh->cctk_time);
       }
@@ -177,7 +178,7 @@ namespace Carpet {
     ++cgh->cctk_iteration;
     if (! adaptive_stepsize) {
       global_time = initial_time
-        + cgh->cctk_iteration * delta_time / maxreflevelfact;
+        + cgh->cctk_iteration * delta_time / maxtimereflevelfact;
       cgh->cctk_time = global_time;
     } else {
       cgh->cctk_time += delta_time;
@@ -191,7 +192,7 @@ namespace Carpet {
 
     for (int rl=0; rl<reflevels; ++rl) {
       const int ml=0;
-      const int do_every = maxreflevelfact / ipow(reffact, rl);
+      const int do_every = maxtimereflevelfact / timereffacts.at(rl);
       if ((cgh->cctk_iteration-1) % do_every == 0) {
         enter_global_mode (cgh, ml);
         enter_level_mode (cgh, rl);
@@ -239,7 +240,7 @@ namespace Carpet {
 
       for (int rl=0; rl<reflevels; ++rl) {
         const int do_every
-          = ipow(mgfact, ml) * (maxreflevelfact / ipow(reffact, rl));
+          = ipow(mgfact, ml) * (maxtimereflevelfact / timereffacts.at(rl));
         if ((cgh->cctk_iteration-1) % do_every == 0) {
           enter_global_mode (cgh, ml);
           enter_level_mode (cgh, rl);
@@ -255,8 +256,8 @@ namespace Carpet {
             vtt.at(m)->advance_time (reflevel, mglevel);
           }
           cgh->cctk_time = (global_time
-                            - delta_time / maxreflevelfact
-                            + delta_time * mglevelfact / reflevelfact);
+                            - delta_time / maxtimereflevelfact
+                            + delta_time * mglevelfact / timereflevelfact);
           CycleTimeLevels (cgh);
 
           Waypoint ("Evolution I at iteration %d time %g%s%s",
@@ -293,7 +294,7 @@ namespace Carpet {
     for (int ml=mglevels-1; ml>=0; --ml) {
       for (int rl=reflevels-1; rl>=0; --rl) {
         const int do_every
-          = ipow(mgfact, ml) * (maxreflevelfact / ipow(reffact, rl));
+          = ipow(mgfact, ml) * (maxtimereflevelfact / timereffacts.at(rl));
         if (cgh->cctk_iteration % do_every == 0) {
           enter_global_mode (cgh, ml);
           enter_level_mode (cgh, rl);
@@ -319,7 +320,7 @@ namespace Carpet {
 
       for (int rl=0; rl<reflevels; ++rl) {
         const int do_every
-          = ipow(mgfact, ml) * (maxreflevelfact / ipow(reffact, rl));
+          = ipow(mgfact, ml) * (maxtimereflevelfact / timereffacts.at(rl));
         if (cgh->cctk_iteration % do_every == 0) {
           enter_global_mode (cgh, ml);
           enter_level_mode (cgh, rl);
@@ -328,7 +329,8 @@ namespace Carpet {
           {
             for (int rl_=0; rl_<reflevels; ++rl_) {
               const int do_every_
-                = ipow(mgfact, ml) * (maxreflevelfact / ipow(reffact, rl_));
+                = (ipow(mgfact, ml)
+                   * (maxtimereflevelfact / timereffacts.at(rl_)));
               if (cgh->cctk_iteration % do_every_ == 0) {
                 finest_active_reflevel = rl_;
               }
diff --git a/Carpet/Carpet/src/Initialise.cc b/Carpet/Carpet/src/Initialise.cc
index befc70784..574e61b23 100644
--- a/Carpet/Carpet/src/Initialise.cc
+++ b/Carpet/Carpet/src/Initialise.cc
@@ -307,7 +307,7 @@ namespace Carpet {
         vtt.at(m)->advance_time (reflevel, mglevel);
       }
       cgh->cctk_time
-        = global_time - tl * delta_time * mglevelfact / reflevelfact;
+        = global_time - tl * delta_time * mglevelfact / timereflevelfact;
       CycleTimeLevels (cgh);
 
       // Set up the initial data
@@ -439,7 +439,7 @@ namespace Carpet {
           vtt.at(m)->advance_time (reflevel, mglevel);
         }
         cgh->cctk_time
-          = global_time + delta_time * mglevelfact / reflevelfact;
+          = global_time + delta_time * mglevelfact / timereflevelfact;
         CycleTimeLevels (cgh);
 
         Waypoint ("Initialisation 3TL evolution I (a) (forwards) at iteration"
@@ -477,7 +477,7 @@ namespace Carpet {
         FlipTimeLevels (cgh);
 
         cgh->cctk_time
-          = global_time + delta_time * mglevelfact / reflevelfact;
+          = global_time + delta_time * mglevelfact / timereflevelfact;
 
         leave_level_mode (cgh);
       } END_MGLEVEL_LOOP;
@@ -549,7 +549,7 @@ namespace Carpet {
           vtt.at(m)->advance_time (reflevel, mglevel);
         }
         cgh->cctk_time
-          = global_time + 2 * delta_time * mglevelfact / reflevelfact;
+          = global_time + 2 * delta_time * mglevelfact / timereflevelfact;
         CycleTimeLevels (cgh);
 
         Waypoint ("Initialisation 3TL evolution I (c) (backwards) at iteration"
@@ -618,7 +618,7 @@ namespace Carpet {
                   (do_meta_mode ? " (meta)" : ""));
 
         const int do_every =
-          ipow(mgfact, mglevel) * (maxreflevelfact / ipow(reffact, rl));
+          ipow(mgfact, mglevel) * (maxtimereflevelfact / timereffacts.at(rl));
         if (cgh->cctk_iteration % do_every == 0)
         {
           // Checkpoint
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index e81b85884..42c65aadc 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -95,7 +95,8 @@ namespace Carpet {
           // Do allow processors with zero grid points
 //           assert (all(bbsss.at(rl).at(c).at(ml).lower() <= bbsssi.at(rl).at(c).at(ml).upper()));
           // Check strides
-          const int str = ipow(reffact, maxreflevels-rl-1) * ipow(mgfact, ml);
+          const ivect str
+            = (maxspacereflevelfact / spacereffacts.at(rl) * ipow(mgfact, ml));
           assert (all(bbsss.at(ml).at(rl).at(c).stride() == str));
           // Check alignments
           assert (all(bbsss.at(ml).at(rl).at(c).lower() % str == 0));
@@ -183,13 +184,13 @@ namespace Carpet {
       for (int rl=0; rl<hh.reflevels(); ++rl) {
         for (int c=0; c<hh.components(rl); ++c) {
           const int convfact = ipow(mgfact, ml);
-          const int levfact = ipow(reffact, rl);
+          const ivect levfact = spacereffacts.at(rl);
           const ivect lower = hh.extents().at(ml).at(rl).at(c).lower();
           const ivect upper = hh.extents().at(ml).at(rl).at(c).upper();
           const ivect stride = hh.extents().at(ml).at(rl).at(c).stride();
-          assert (all(lower * levfact % maxreflevelfact == 0));
-          assert (all(upper * levfact % maxreflevelfact == 0));
-          assert (all(((upper - lower) * levfact / maxreflevelfact)
+          assert (all(lower * levfact % maxspacereflevelfact == 0));
+          assert (all(upper * levfact % maxspacereflevelfact == 0));
+          assert (all(((upper - lower) * levfact / maxspacereflevelfact)
                       % convfact == 0));
           cout << "   [" << ml << "][" << rl << "][" << m << "][" << c << "]"
                << "   exterior: "
@@ -219,12 +220,12 @@ namespace Carpet {
           const ivect lower = hh.extents().at(ml).at(rl).at(c).lower();
           const ivect upper = hh.extents().at(ml).at(rl).at(c).upper();
           const int convfact = ipow(mgfact, ml);
-          const int levfact = ipow(reffact, rl);
+          const ivect levfact = spacereffacts.at(rl);
           cout << "   [" << ml << "][" << rl << "][" << m << "][" << c << "]"
                << "   exterior: "
-               << origin + delta * lower / maxreflevelfact
+               << origin + delta * lower / maxspacereflevelfact
                << " : "
-               << origin + delta * upper / maxreflevelfact
+               << origin + delta * upper / maxspacereflevelfact
                << " : "
                << delta * convfact / levfact << endl;
         }
@@ -316,7 +317,7 @@ namespace Carpet {
         
         // TODO: processor distribution, average load, std deviation
         
-        coarsevolume = totalvolume * ipow (reflevelfact, dim);
+        coarsevolume = totalvolume * prod (rvect (spacereflevelfact));
       }
     }
     cout.precision (oldprecision);
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index d3dd4d54a..9ea25caa9 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -293,8 +293,11 @@ namespace Carpet {
   static void setup_multigrid_info (cGH* cgh, CCTK_INT convergence_level,
                       CCTK_INT convergence_factor,
                       CCTK_INT num_convergence_levels);
-  static void setup_refinement_level_info (CCTK_INT max_refinement_levels,
-                      CCTK_INT refinement_factor);
+  static void
+  setup_refinement_level_info (CCTK_INT const max_refinement_levels,
+                               CCTK_INT const refinement_factor,
+                               char const * const space_refinement_factors,
+                               char const * const time_refinement_factors);
   static void allocate_map(cGH* cgh, int m,
                       CCTK_INT domain_from_coordbase,
                       CCTK_INT convergence_factor, CCTK_INT refinement_factor,
@@ -371,7 +374,9 @@ namespace Carpet {
     setup_multigrid_info ( cgh, convergence_level, convergence_factor,
                            num_convergence_levels );
     
-    setup_refinement_level_info ( max_refinement_levels, refinement_factor );
+    setup_refinement_level_info (max_refinement_levels, refinement_factor,
+                                 space_refinement_factors,
+                                 time_refinement_factors);
     
     // Map information
     carpetGH.maps = maps = num_maps;
@@ -471,11 +476,72 @@ namespace Carpet {
     cgh->cctk_convfac = mgfact;
   }
 
-  void setup_refinement_level_info (CCTK_INT max_refinement_levels,
-                                      CCTK_INT refinement_factor) { 
+  void
+  setup_refinement_level_info (CCTK_INT const max_refinement_levels,
+                               CCTK_INT const refinement_factor,
+                               char const * const space_refinement_factors,
+                               char const * const time_refinement_factors)
+  {
+    // Set maximum number of refinement levels
     maxreflevels = max_refinement_levels;
+    
+#if 0
+    // Set default refinement factor
     reffact = refinement_factor;
-    maxreflevelfact = ipow(reffact, maxreflevels-1);
+#endif
+    
+    // Set the per-level spatial refinement factors
+    if (strcmp (space_refinement_factors, "") == 0) {
+      // Calculate them from the default refinement factor
+      spacereffacts.resize (maxreflevels);
+      for (int n=0; n<maxreflevels; ++n) {
+        spacereffacts.at(n) = ivect (ipow (refinement_factor, n));
+      }
+    } else {
+      // Read them from the parameter
+      try {
+        istringstream srf (space_refinement_factors);
+        srf >> spacereffacts;
+      } catch (input_error) {
+        CCTK_WARN
+          (0, "Could not parse parameter \"space_refinement_factors\"");
+      }
+    }
+    // TODO: turn these into real error messages
+    assert (spacereffacts.size() >= maxreflevels);
+    assert (all (spacereffacts.front() == 1));
+    for (int n=1; n<maxreflevels; ++n) {
+      assert (all (spacereffacts.at(n) >= spacereffacts.at(n-1)));
+      assert (all (spacereffacts.at(n) % spacereffacts.at(n-1) == 0));
+    }
+    
+    // Set the per-level temporal refinement factors
+    if (strcmp (time_refinement_factors, "") == 0) {
+      // Calculate them from the default refinement factor
+      timereffacts.resize (maxreflevels);
+      for (int n=0; n<maxreflevels; ++n) {
+        timereffacts.at(n) = ipow (refinement_factor, n);
+      }
+    } else {
+      // Read them from the parameter
+      try {
+        istringstream trf (time_refinement_factors);
+        trf >> timereffacts;
+      } catch (input_error) {
+        CCTK_WARN (0, "Could not parse parameter \"time_refinement_factors\"");
+      }
+    }
+    // TODO: turn these into real error messages
+    assert (timereffacts.size() >= maxreflevels);
+    assert (timereffacts.front() == 1);
+    for (int n=1; n<maxreflevels; ++n) {
+      assert (timereffacts.at(n) >= timereffacts.at(n-1));
+      assert (timereffacts.at(n) % timereffacts.at(n-1) == 0);
+    }
+    
+    // Calculate the maximum refinement factors
+    maxtimereflevelfact = timereffacts.at (maxreflevels-1);
+    maxspacereflevelfact = spacereffacts.at (maxreflevels-1);
   }
 
   void allocate_map (cGH* cgh, int m,
@@ -558,14 +624,13 @@ namespace Carpet {
     Sanity_check (npoints);
     
     // Base grid extent
-    const int stride = maxreflevelfact;
-    const ivect str (stride);
+    const ivect str (maxspacereflevelfact);
     const ivect lb (0);
     const ivect ub ((npoints - 1) * str);
     const ibbox baseext (lb, ub, str);
     
     // Allocate grid hierarchy
-    vhh.at(m) = new gh (refinement_factor, vertex_centered,
+    vhh.at(m) = new gh (spacereffacts, vertex_centered,
                         convergence_factor, vertex_centered, baseext);
     
     // Allocate data hierarchy
@@ -573,7 +638,7 @@ namespace Carpet {
                         prolongation_order_space, buffer_width);
     
     // Allocate time hierarchy
-    vtt.at(m) = new th (*vhh.at(m), 1.0);
+    vtt.at(m) = new th (*vhh.at(m), timereffacts, 1.0);
       
     check_time_hierarchy (vdd, m, max_refinement_levels, refinement_factor,
                              prolongation_order_space, lghosts, ughosts);
@@ -806,13 +871,13 @@ namespace Carpet {
         const ivect lower = vhh.at(m)->extents().at(ml).at(rl).at(c).lower();
         const ivect upper = vhh.at(m)->extents().at(ml).at(rl).at(c).upper();
         const int convfact = ipow(mgfact, ml);
-        assert (all(lower % maxreflevelfact == 0));
-        assert (all(upper % maxreflevelfact == 0));
-        assert (all(((upper - lower) / maxreflevelfact) % convfact == 0));
+        assert (all(lower % maxspacereflevelfact == 0));
+        assert (all(upper % maxspacereflevelfact == 0));
+        assert (all(((upper - lower) / maxspacereflevelfact) % convfact == 0));
         cout << "   [" << ml << "][" << rl << "][" << m << "][" << c << "]"
-             << "   exterior extent: " << lower / maxreflevelfact
-             << " : " << upper / maxreflevelfact
-             << "   (" << (upper - lower) / maxreflevelfact / convfact + 1
+             << "   exterior extent: " << lower / maxspacereflevelfact
+             << " : " << upper / maxspacereflevelfact
+             << "   (" << (upper - lower) / maxspacereflevelfact / convfact + 1
              << ")" << endl;
       }
     }
@@ -956,9 +1021,11 @@ namespace Carpet {
       
       assert (all(convpowers == convpowers[0]));
       const int amgfact1 = ipow(mgfact, convpowers[0]);
+      const vector<ivect> aspacereffacts (1, ivect (1));
+      const vector<int> atimereffacts (1, 1);
       
       arrdata.at (group).at (0).hh
-        = new gh (refinement_factor, vertex_centered,
+        = new gh (aspacereffacts, vertex_centered,
                   amgfact1, vertex_centered,
                   abaseext);
       
@@ -966,7 +1033,7 @@ namespace Carpet {
         = new dh (*arrdata.at (group).at (0).hh, alghosts, aughosts, 0, 0);
       
       arrdata.at (group).at (0).tt
-        = new th (*arrdata.at (group).at (0).hh, 1.0);
+        = new th (*arrdata.at (group).at (0).hh, atimereffacts, 1.0);
 
       
       
diff --git a/Carpet/Carpet/src/modes.cc b/Carpet/Carpet/src/modes.cc
index cfa3728b9..7c2a922c3 100644
--- a/Carpet/Carpet/src/modes.cc
+++ b/Carpet/Carpet/src/modes.cc
@@ -236,9 +236,10 @@ namespace Carpet {
     Checkpoint ("Entering level mode");
     
     reflevel = rl;
-    reflevelfact = ipow(reffact, reflevel);
-    ivect::ref(cgh->cctk_levfac) = reflevelfact;
-    cgh->cctk_timefac = reflevelfact;
+    timereflevelfact = timereffacts.at (reflevel);
+    spacereflevelfact = spacereffacts.at (reflevel);
+    ivect::ref(cgh->cctk_levfac) = spacereflevelfact;
+    cgh->cctk_timefac = timereflevelfact;
     
     // Set current time
     assert (mglevel>=0 and mglevel<(int)leveltimes.size());
@@ -272,7 +273,8 @@ namespace Carpet {
     }
     
     reflevel = -1;
-    reflevelfact = -deadbeef;
+    timereflevelfact = -deadbeef;
+    spacereflevelfact = ivect(-deadbeef);
     ivect::ref(cgh->cctk_levfac) = -deadbeef;
     cgh->cctk_timefac = -deadbeef;
     
diff --git a/Carpet/Carpet/src/variables.cc b/Carpet/Carpet/src/variables.cc
index 190878acb..8b26a38e7 100644
--- a/Carpet/Carpet/src/variables.cc
+++ b/Carpet/Carpet/src/variables.cc
@@ -18,11 +18,20 @@ namespace Carpet {
   // Refinement levels
   int reflevels;
   
+#if 0
   // Refinement factor
   int reffact;
+#endif
   
-  // Refinement factor on finest grid
-  int maxreflevelfact;
+  // Temporal refinement factors over the coarsest grid
+  vector<int> timereffacts;
+  
+  // Spatial refinement factors over the coarsest grid
+  vector<vect<int,dim> > spacereffacts;
+  
+  // Maximum refinement factors on finest possible grid
+  int maxtimereflevelfact;
+  vect<int,dim> maxspacereflevelfact;
   
   // Base multigrid level
   int basemglevel;
@@ -47,10 +56,11 @@ namespace Carpet {
   int map;
   int component;
   
-  // refinement factor of current level: ipow(refinement_factor, reflevel)
-  int reflevelfact;
+  // Current refinement factors
+  int timereflevelfact;
+  vect<int,dim> spacereflevelfact;
   
-  // multigrid factor of current level: ipow(multigrid_factor, mglevel)
+  // Current multigrid factor
   int mglevelfact;
   
   
diff --git a/Carpet/Carpet/src/variables.hh b/Carpet/Carpet/src/variables.hh
index 7bab18432..c36410ead 100644
--- a/Carpet/Carpet/src/variables.hh
+++ b/Carpet/Carpet/src/variables.hh
@@ -44,11 +44,20 @@ namespace Carpet {
   // Refinement levels
   extern int reflevels;
   
+#if 0
   // Refinement factor
   extern int reffact;
+#endif
   
-  // Refinement factor on finest possible grid
-  extern int maxreflevelfact;
+  // Temporal refinement factors over the coarsest grid
+  extern vector<int> timereffacts;
+  
+  // Spatial refinement factors over the coarsest grid
+  extern vector<vect<int,dim> > spacereffacts;
+  
+  // Maximum refinement factors on finest possible grid
+  extern int maxtimereflevelfact;
+  extern vect<int,dim> maxspacereflevelfact;
   
   // Base multigrid level
   extern int basemglevel;
@@ -73,8 +82,9 @@ namespace Carpet {
   extern int map;
   extern int component;
   
-  // Current refinement factor
-  extern int reflevelfact;
+  // Current refinement factors
+  extern int timereflevelfact;
+  extern vect<int,dim> spacereflevelfact;
   
   // Current multigrid factor
   extern int mglevelfact;
diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc
index fcf1e2f6b..e40880e03 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.cc
+++ b/Carpet/CarpetIOASCII/src/ioascii.cc
@@ -535,7 +535,8 @@ namespace CarpetIOASCII {
                 if (grouptype == CCTK_GF) {
                   for (int d=0; d<dim; ++d) {
                     global_lower[d] = cgh->cctk_origin_space[d];
-                    coord_delta[d] = cgh->cctk_delta_space[d] / maxreflevelfact;
+                    coord_delta[d]
+                      = cgh->cctk_delta_space[d] / maxspacereflevelfact[d];
                   }
                 } else {
                   for (int d=0; d<dim; ++d) {
diff --git a/Carpet/CarpetIOHDF5/src/Recover.cc b/Carpet/CarpetIOHDF5/src/Recover.cc
index e8bcb2e4b..772110a74 100644
--- a/Carpet/CarpetIOHDF5/src/Recover.cc
+++ b/Carpet/CarpetIOHDF5/src/Recover.cc
@@ -523,7 +523,7 @@ int ReadVar (const cGH* const cctkGH, const int vindex,
         gdata* const data = (*ff) (tl, reflevel, component, mglevel);
 
         // Create temporary data storage on processor 0
-        vect<int,dim> str = vect<int,dim>(maxreflevelfact/reflevelfact);
+        vect<int,dim> str = maxspacereflevelfact/spacereflevelfact;
 
         if(grouptype == CCTK_SCALAR || grouptype == CCTK_ARRAY)
           str = vect<int,dim> (1);
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc
index e7798b9ba..1916ef513 100644
--- a/Carpet/CarpetInterp/src/interp.cc
+++ b/Carpet/CarpetInterp/src/interp.cc
@@ -232,7 +232,7 @@ namespace CarpetInterp {
     }
 #else
     rvect const lower = rvect::ref(cgh->cctk_origin_space);
-    rvect const delta = rvect::ref(cgh->cctk_delta_space) / maxreflevelfact;
+    rvect const delta = rvect::ref(cgh->cctk_delta_space) / maxspacereflevelfact;
     rvect const upper = lower + delta * (vhh.at(m)->baseextent.upper() - vhh.at(m)->baseextent.lower());
 #endif
     
@@ -272,7 +272,7 @@ namespace CarpetInterp {
       if (all(pos>=lower && pos<=upper)) {
         for (int rl=maxrl-1; rl>=minrl; --rl) {
           
-          int const fact = maxreflevelfact / ipow(reffact, rl) * ipow(mgfact, mglevel);
+          ivect const fact = maxspacereflevelfact / spacereffacts.at(rl) * ipow(mgfact, mglevel);
           ivect const ipos = ivect(floor((pos - lower) / (delta * fact) + 0.5)) * fact;
           assert (all(ipos % vhh.at(m)->bases().at(ml).at(rl).stride() == 0));
           
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index fcbbace14..d83b6e0a5 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -150,14 +150,18 @@ ostream& output (ostream& os, const vector<T>& v) {
 
 #include "bbox.hh"
 #include "bboxset.hh"
+#include "vect.hh"
 
 template int ipow (int x, int y);
 template CCTK_REAL ipow (CCTK_REAL x, int y);
+template vect<int,3> ipow (vect<int,3> x, int y);
 
+template istream& input (istream& os, vector<int>& v);
 template istream& input (istream& os, vector<bbox<int,3> >& v);
 template istream& input (istream& os, vector<bbox<CCTK_REAL,3> >& v);
 template istream& input (istream& os, vector<vector<bbox<int,3> > >& v);
 template istream& input (istream& os, vector<vector<bbox<CCTK_REAL,3> > >& v);
+template istream& input (istream& os, vector<vect<int,3> >& v);
 template istream& input (istream& os, vector<vect<vect<bool,2>,3> >& v);
 template istream& input (istream& os, vector<vector<vect<vect<bool,2>,3> > >& v);
 
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index 55c23def9..2a0a631f5 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -1,6 +1,7 @@
 #include <cassert>
 #include <cstdlib>
 #include <iostream>
+#include <vector>
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
@@ -8,6 +9,7 @@
 #include "defs.hh"
 #include "dh.hh"
 #include "th.hh"
+#include "vect.hh"
 
 #include "gh.hh"
 
@@ -16,13 +18,19 @@ using namespace std;
 
 
   // Constructors
-gh::gh (const int reffact_, const centering refcent_,
+gh::gh (const vector<ivect> & reffacts_, const centering refcent_,
         const int mgfact_, const centering mgcent_,
         const ibbox baseextent_)
-  : reffact(reffact_), refcent(refcent_),
+  : reffacts(reffacts_), refcent(refcent_),
     mgfact(mgfact_), mgcent(mgcent_),
     baseextent(baseextent_)
 {
+  assert (reffacts.size() >= 1);
+  assert (all (reffacts.front() == 1));
+  for (size_t n = 1; n < reffacts.size(); ++ n) {
+    assert (all (reffacts.at(n) >= reffacts.at(n-1)));
+    assert (all (reffacts.at(n) % reffacts.at(n-1) == 0));
+  }
 }
 
 // Destructors
@@ -136,7 +144,8 @@ void gh::check_refinement_levels ()
   for (int ml=0; ml<mglevels(); ++ml) {
     for (int rl=1; rl<reflevels(); ++rl) {
       assert (all(extents().at(ml).at(rl-1).at(0).stride()
-                  == ivect(reffact) * extents().at(ml).at(rl).at(0).stride()));
+                  == ((reffacts.at(rl) / reffacts.at(rl-1))
+                      * extents().at(ml).at(rl).at(0).stride())));
       // Check contained-ness:
       // first take all coarse grids ...
       ibset all;
@@ -240,7 +249,7 @@ void gh::do_output_bases (ostream& os) const
 ostream& gh::output (ostream& os) const
 {
   os << "gh:"
-     << "reffactor=" << reffact << ",refcentering=" << refcent << ","
+     << "reffacts=" << reffacts << ",refcentering=" << refcent << ","
      << "mgfactor=" << mgfact << ",mgcentering=" << mgcent << ","
      << "extents=" << extents() << ","
      << "outer_boundaries=" << outer_boundaries() << ","
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index ee7b86844..f46fd5c34 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -46,7 +46,7 @@ public:
 public:				// should be readonly
   
   // Fields
-  const int reffact;		// refinement factor
+  const vector<ivect> reffacts; // refinement factors
   const centering refcent;	// vertex or cell centered
  
   const int mgfact;		// default multigrid factor
@@ -68,7 +68,7 @@ private:
 public:
   
   // Constructors
-  gh (const int reffact, const centering refcent,
+  gh (const vector<ivect> & reffacts, const centering refcent,
       const int mgfact, const centering mgcent,
       const ibbox baseextent);
   
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index 84b896b55..0751a03e1 100644
--- a/Carpet/CarpetLib/src/th.cc
+++ b/Carpet/CarpetLib/src/th.cc
@@ -1,6 +1,7 @@
 #include <cassert>
 #include <cmath>
 #include <iostream>
+#include <vector>
 
 #include "cctk.h"
 
@@ -14,9 +15,15 @@ using namespace std;
 
 
 // Constructors
-th::th (gh& h_, const CCTK_REAL basedelta)
-  : h(h_), delta(basedelta)
+th::th (gh& h_, const vector<int> & reffacts_, const CCTK_REAL basedelta)
+  : h(h_), reffacts(reffacts_), delta(basedelta)
 {
+  assert (reffacts.size() >= 1);
+  assert (reffacts.front() == 1);
+  for (size_t n = 1; n < reffacts.size(); ++ n) {
+    assert (reffacts.at(n) >= reffacts.at(n-1));
+    assert (reffacts.at(n) % reffacts.at(n-1) == 0);
+  }
   h.add(this);
 }
 
@@ -44,10 +51,8 @@ void th::recompose ()
       } else {
         times.at(ml).at(rl) = times.at(ml).at(rl-1);
       }
-      if (ml==0 and rl==0) {
-	deltas.at(ml).at(rl) = delta;
-      } else if (ml==0) {
-	deltas.at(ml).at(rl) = deltas.at(ml).at(rl-1) / h.reffact;
+      if (ml==0) {
+	deltas.at(ml).at(rl) = delta / reffacts.at(rl);
       } else {
 	deltas.at(ml).at(rl) = deltas.at(ml-1).at(rl) * h.mgfact;
       }
diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh
index 42bd6504e..3c24c227c 100644
--- a/Carpet/CarpetLib/src/th.hh
+++ b/Carpet/CarpetLib/src/th.hh
@@ -32,6 +32,8 @@ public:				// should be readonly
   
 private:
   
+  const vector<int> reffacts;
+
   CCTK_REAL delta;		// time step
   vector<vector<CCTK_REAL> > times; // current times [ml][rl]
   vector<vector<CCTK_REAL> > deltas; // time steps [ml][rl]
@@ -39,7 +41,7 @@ private:
 public:
   
   // Constructors
-  th (gh& h, const CCTK_REAL basedelta);
+  th (gh& h, const vector<int> & reffacts, const CCTK_REAL basedelta);
   
   // Destructors
   ~th ();
diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc
index 342c3901f..31250ce4b 100644
--- a/Carpet/CarpetReduce/src/mask_carpet.cc
+++ b/Carpet/CarpetReduce/src/mask_carpet.cc
@@ -31,8 +31,6 @@ namespace CarpetMask {
   {
     DECLARE_CCTK_PARAMETERS;
     
-    assert (reffact == 2);
-    
     if (! is_singlemap_mode()) {
       CCTK_WARN (0, "This routine may only be called in singlemap mode");
     }
@@ -47,6 +45,10 @@ namespace CarpetMask {
       
       ibbox const & base = hh.bases().at(mglevel).at(reflevel);
       
+      ivect const reffact
+        = spacereffacts.at(reflevel) / spacereffacts.at(reflevel-1);
+      assert (all (reffact == 2));
+      
       
       
       // Calculate the union of all refined regions
diff --git a/Carpet/CarpetReduce/src/reduce.cc b/Carpet/CarpetReduce/src/reduce.cc
index f51dda6ca..a4d5cf532 100644
--- a/Carpet/CarpetReduce/src/reduce.cc
+++ b/Carpet/CarpetReduce/src/reduce.cc
@@ -1084,7 +1084,7 @@ namespace CarpetReduce {
               weight = (static_cast<CCTK_REAL const *>
                         (CCTK_VarDataPtr (cgh, 0, "CarpetReduce::weight")));
               assert (weight);
-              levfac = ipow((CCTK_REAL)reflevelfact, -grpdim);
+              levfac = 1.0 / prod (rvect (spacereflevelfact));
             } else {
               weight = NULL;
               levfac = 1.0;
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,