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

4 files changed, 29 insertions, 19 deletions

diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
index c5e8ec62f..ed01bd34c 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
@@ -261,6 +261,8 @@ namespace CarpetAdaptiveRegrid {
       vector<vector<ibbox> > bbss = bbsss.at(0);
       vector<vector<ibbox> > local_bbss = local_bbsss.at(0);
       
+      ivect const reffact = spacereffacts.at(rl+1) / spacereffacts.at(rl);
+
       bool did_regrid = false;
             
       //
@@ -458,15 +460,15 @@ namespace CarpetAdaptiveRegrid {
                     // Correct for the change in level !
                     CCTK_INT ii = (i + cctkGH->cctk_lbnd[0] 
                                    - cctkGH->cctk_nghostzones[0]
-                                   + child_levoff[0]) / reffact - 
+                                   + child_levoff[0]) / reffact[0] - 
                       imin[0];
                     CCTK_INT jj = (j + cctkGH->cctk_lbnd[1] 
                                    - cctkGH->cctk_nghostzones[1]
-                                   + child_levoff[1]) / reffact - 
+                                   + child_levoff[1]) / reffact[1] - 
                       imin[1];
                     CCTK_INT kk = (k + cctkGH->cctk_lbnd[2] 
                                    - cctkGH->cctk_nghostzones[2]
-                                   + child_levoff[1]) / reffact - 
+                                   + child_levoff[1]) / reffact[2] - 
                       imin[2];
                     // Check that this point actually intersects with 
                     // this box (if this component was actually a
@@ -822,7 +824,7 @@ namespace CarpetAdaptiveRegrid {
           rvect hi = int2pos(cctkGH, hh, ihi, reflevel);
           rvect str = base_spacing * 
             ipow((CCTK_REAL)mgfact, basemglevel) /
-            ipow(reffact, reflevel);
+            spacereffacts.at(reflevel);
           rbbox newbbcoord(lo, hi, str);
           
           if (veryverbose) {
@@ -1057,7 +1059,7 @@ namespace CarpetAdaptiveRegrid {
     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;
     
@@ -1084,7 +1086,7 @@ namespace CarpetAdaptiveRegrid {
     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/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc b/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc
index 523c5cf9b..1768aa6ed 100644
--- a/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/manualcoordinatelist.cc
@@ -80,7 +80,7 @@ namespace CarpetAdaptiveRegrid {
       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],
@@ -94,7 +94,7 @@ namespace CarpetAdaptiveRegrid {
             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();
@@ -124,7 +124,7 @@ namespace CarpetAdaptiveRegrid {
         bbvect const & ob = newobss.at(rl-1).at(c);
         // 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,
@@ -230,7 +230,7 @@ namespace CarpetAdaptiveRegrid {
     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/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc
index 99e0ceb6a..cf5407c3f 100644
--- a/CarpetDev/CarpetJacobi/src/Jacobi.cc
+++ b/CarpetDev/CarpetJacobi/src/Jacobi.cc
@@ -194,7 +194,8 @@ namespace CarpetJacobi {
         for (int m=0; m<maps; ++m) {
           vtt[m]->set_time (reflevel, mglevel, 0);
           vtt[m]->set_delta
-            (reflevel, mglevel, 1.0 / ipow(reflevelfact, reflevelpower));
+            (reflevel, mglevel, 1.0 / ipow (maxval (spacereflevelfact),
+                                            reflevelpower));
         }
       } END_REFLEVEL_LOOP;
       
@@ -217,14 +218,16 @@ namespace CarpetJacobi {
       }
       
       const int icoarsestep
-        = ipow (mglevelfact * maxreflevelfact, reflevelpower);
+        = ipow (mglevelfact * maxval (maxspacereflevelfact), reflevelpower);
       const int istep
-        = ipow (mglevelfact * maxreflevelfact / ipow(reffact, solve_level),
+        = ipow (mglevelfact * maxval (maxspacereflevelfact
+                                      / spacereffacts.at(solve_level)),
                 reflevelpower);
       int iter = istep;
       for (;; iter+=istep) {
         
-        global_time = 1.0 * iter / ipow(maxreflevelfact, reflevelpower);
+        global_time = 1.0 * iter / ipow (maxval (maxspacereflevelfact),
+                                         reflevelpower);
         * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
 //         cout << "CJ global iter " << iter << " time " << global_time << flush << endl;
         
@@ -232,7 +235,8 @@ namespace CarpetJacobi {
         ierr = 0;
         BEGIN_REFLEVEL_LOOP(cctkGH) {
           const int do_every
-            = ipow (mglevelfact * (maxreflevelfact/reflevelfact),
+            = ipow (mglevelfact * maxval (maxspacereflevelfact
+                                          / spacereflevelfact),
                     reflevelpower);
           if (reflevel <= solve_level && (iter-istep) % do_every == 0) {
             
@@ -242,7 +246,7 @@ namespace CarpetJacobi {
             }
             * const_cast<CCTK_REAL *> (& cctkGH->cctk_time)
               = (1.0 * (iter - istep + do_every)
-                 / ipow(maxreflevelfact, reflevelpower));
+                 / ipow (maxval (maxspacereflevelfact), reflevelpower));
             CycleTimeLevels (cctkGH);
 //             cout << "CJ residual iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl;
             
@@ -326,7 +330,8 @@ namespace CarpetJacobi {
         ierr = 0;
         BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) {
           const int do_every
-            = ipow (mglevelfact * (maxreflevelfact/reflevelfact),
+            = ipow (mglevelfact * maxval (maxspacereflevelfact
+                                          / spacereflevelfact),
                     reflevelpower);
           if (reflevel <= solve_level && iter % do_every == 0) {
             
@@ -533,7 +538,7 @@ namespace CarpetJacobi {
         * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
         for (int m=0; m<maps; ++m) {
           vtt[m]->set_time (reflevel, mglevel, 0);
-          vtt[m]->set_delta (reflevel, mglevel, 1.0 / reflevelfact);
+          vtt[m]->set_delta (reflevel, mglevel, 1.0 / timereflevelfact);
         }
       } END_REFLEVEL_LOOP;
       
diff --git a/CarpetDev/CarpetMG/src/mg.cc b/CarpetDev/CarpetMG/src/mg.cc
index 24a7ab8e3..92151e134 100644
--- a/CarpetDev/CarpetMG/src/mg.cc
+++ b/CarpetDev/CarpetMG/src/mg.cc
@@ -590,7 +590,10 @@ namespace CarpetMG {
             }
             
             // Recurse
-            CCTK_REAL const coarse_minerror = error / ipow(reffact, dim);
+            ivect const reffact
+              = (spacereffacts.at (reflevel)
+                 / spacereffacts.at (coarse_reflevel));
+            CCTK_REAL const coarse_minerror = error / prod (reffact);
             int const ierr
               = multigrid (cctkGH,
                            var, res, rhs, sav, wgt, aux,