Import the recently announced changes:

Import the recently announced changes:
1. Carpet has now an infrastructure for multiple maps (aka "grid
patches").  Instead of a single grid hierarchy there can now be
several.  This is largely untested, because the remainder of Cactus
cannot handle multiple coordinate systems.

2. The order in which the schedule bins are called has changed.  As Ian
Hawke pointed out, the previous order during time evolution was
inconsistent.  The initial data ordering did not allow for recovering
and was not usable for progressively solving elliptic equations for
initial data.

3. Carpet now supports convergence levels.  The convergence level
specifies by how many factors of two the resolution in the parameter
file should be coarsened (or refined, if negative).  This should make
convergence tests and test runs much easier.  It is, in principle, also
possible to run several convergence levels at once.  This has not been
tested because the remainder of Cactus cannot handle multiple
resolutions.  This will be necessary for a multigrid solver, and also
for having a shadow hierarchy to determine where to refine adaptively.

4. Carpet works together with the new CoordBase domain specification
parameters.  Without these, using convergence levels will lead to very
strange results.

5. The "modes" have changed.  There are now:
        meta mode:      the whole simulation
        global mode:    one convergence level
        level mode:     one refinement level
        singlemap mode: one map on one refinement level
        local mode:     as previously
The whole mode handling has been cleaned up.

6. The regridding thorn has been cleaned up.

7. The kind of prolongation stencil is now determined in Carpet, i.e. at
a fairly hight level, instead of in CarpetLib.

8. The low-order prolongation operators have been made much more
efficient (as have previously the higher-order ones).

9. Assorted smaller changes.

For Carpet users, there should be no major incompatibilities.  The major
improvements are 3 and 4 combined.  Here is an example:

CoordBase::domainsize    = extent
CoordBase::spacing       = gridspacing
CoordBase::zero_origin_x = yes
CoordBase::zero_origin_y = yes
CoordBase::zero_origin_z = yes
CoordBase::xextent       = 20.0
CoordBase::yextent       = 20.0
CoordBase::zextent       = 20.0
CoordBase::dx            = 1.0
CoordBase::dy            = 1.0
CoordBase::dz            = 1.0

CoordBase::boundary_shiftout_x_lower = 1
CoordBase::boundary_shiftout_y_lower = 1
CoordBase::boundary_shiftout_z_lower = 1

Carpet::domain_from_coordbase   = yes
Carpet::convergence_level       = 0

grid::type         = coordbase
grid::domain       = octant
grid::avoid_origin = no

This gives you a grid that extends from the origin ("zero_origin") up to
20.0 with a grid spacing of 1.0.  Symmetry zones and boundary zones are
added automatically.  The "shiftout" says that there is no boundary
point on the origin.  The staggering parameters (not shown) default to
"no".  In order to change the resolution, only the convergence level
has to be adjusted.  Note that the old way of specifying the domain
extent still works.

For Carpet developers, one major change is the new mode handling.  As
described in 5, the looping macros (that loop over all refinement
levels, or all components) have changed.

darcs-hash:20040125125727-07bb3-3368611314b2dcb8c8ae58ab3f501b683d7edb8f.gz

2 files changed, 369 insertions, 343 deletions

diff --git a/CarpetDev/CarpetJacobi/param.ccl b/CarpetDev/CarpetJacobi/param.ccl
index ce69dbc91..bdeecf41b 100644
--- a/CarpetDev/CarpetJacobi/param.ccl
+++ b/CarpetDev/CarpetJacobi/param.ccl
@@ -1,5 +1,5 @@
 # Parameter definitions for thorn CarpetJacobi
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/param.ccl,v 1.1 2003/09/02 14:35:57 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/param.ccl,v 1.2 2004/01/25 14:57:29 schnetter Exp $
 
 BOOLEAN verbose "Produce log output while running" STEERABLE=always
 {
@@ -18,3 +18,15 @@ CCTK_INT outeveryseconds "Produce output every n seconds" STEERABLE=always
 {
   0:* :: ""
 } 10
+
+
+
+CCTK_INT solve_minlevel "First level on which to solve" STEERABLE=always
+{
+  0:* :: ""
+} 0
+
+CCTK_INT reflevelpower "Power of the refinement factor" STEERABLE=always
+{
+  0:* :: ""
+} 2
diff --git a/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc
index ad673b342..c326b66cd 100644
--- a/CarpetDev/CarpetJacobi/src/Jacobi.cc
+++ b/CarpetDev/CarpetJacobi/src/Jacobi.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/Jacobi.cc,v 1.2 2003/11/19 14:05:14 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/Jacobi.cc,v 1.3 2004/01/25 14:57:29 schnetter Exp $
 
 #include <cassert>
 #include <cmath>
@@ -14,14 +14,18 @@
 #include "util_ErrorCodes.h"
 #include "util_Table.h"
 
-#include "carpet.hh"
 #include "TATelliptic.h"
 
+#include "carpet.hh"
+
+#include "th.hh"
+
 
 
 namespace Carpet {
   // TODO: fix this
-  void Restrict (const cGH* cgh);
+  void CycleTimeLevels (const cGH* cctkGH);
+  void Restrict (const cGH* cctkGH);
 };
 
 
@@ -48,134 +52,6 @@ namespace CarpetJacobi {
   
   
   
-  namespace common {
-    
-    const int * var;
-    const int * res;
-    int nvars;
-    int options_table;
-    calcfunc calcres;
-    calcfunc applybnds;
-    void * userdata;
-    
-    CCTK_REAL factor;
-    CCTK_REAL * factors;
-    
-    CCTK_REAL norm_count;
-    CCTK_REAL norm_l2;
-    
-    int ierr;
-    
-    void call_calcres (cGH * const cctkGH)
-    {
-      if (ierr) return;
-      ierr = calcres (cctkGH, options_table, userdata);
-    }
-    
-    void norm (cGH * const cctkGH)
-    {
-      int ierr;
-      
-      cGroup groupdata;
-      ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata);
-      assert (!ierr);
-      cGroupDynamicData groupdyndata;
-      ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]),
-                                    &groupdyndata);
-      assert (!ierr);
-      const int thedim = groupdata.dim;
-      assert (thedim>=0 && thedim<=dim);
-      assert (thedim == groupdyndata.dim);
-      int lsh[dim], nghostzones[dim];
-      for (int d=0; d<thedim; ++d) {
-        lsh[d] = groupdyndata.lsh[d];
-        nghostzones[d] = groupdyndata.nghostzones[d];
-      }
-      for (int d=thedim; d<dim; ++d) {
-        lsh[d] = 1;
-        nghostzones[d] = 0;
-      }
-      for (int d=0; d<dim; ++d) {
-        assert (lsh[d]>=0);
-        assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]);
-      }
-      
-      for (int n=0; n<nvars; ++n) {
-        const CCTK_REAL * restrict resptr
-          = (const CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, res[n]);
-        assert (resptr);
-        const CCTK_REAL fac = factors[n];
-        for (int k=nghostzones[2]; k<lsh[2]-nghostzones[2]; ++k) {
-          for (int j=nghostzones[1]; j<lsh[1]-nghostzones[1]; ++j) {
-            for (int i=nghostzones[0]; i<lsh[0]-nghostzones[0]; ++i) {
-              const int ind = i + lsh[0] * (j + lsh[1] * k);
-              ++norm_count;
-              // TODO: scale the norm by the resolution?
-              norm_l2 += pow(fac * resptr[ind], 2);
-            }
-          }
-        }
-      }
-    }
-    
-    void smooth (cGH * const cctkGH)
-    {
-      int ierr;
-      
-      cGroup groupdata;
-      ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata);
-      assert (!ierr);
-      cGroupDynamicData groupdyndata;
-      ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]),
-                                    &groupdyndata);
-      assert (!ierr);
-      const int thedim = groupdata.dim;
-      assert (thedim>=0 && thedim<=dim);
-      assert (thedim == groupdyndata.dim);
-      int lsh[dim], nghostzones[dim];
-      for (int d=0; d<thedim; ++d) {
-        lsh[d] = groupdyndata.lsh[d];
-        nghostzones[d] = groupdyndata.nghostzones[d];
-      }
-      for (int d=thedim; d<dim; ++d) {
-        lsh[d] = 1;
-        nghostzones[d] = 0;
-      }
-      for (int d=0; d<dim; ++d) {
-        assert (lsh[d]>=0);
-        assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]);
-      }
-      
-      for (int n=0; n<nvars; ++n) {
-        CCTK_REAL * restrict varptr
-          = (CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, var[n]);
-        assert (varptr);
-        const CCTK_REAL * restrict resptr
-          = (const CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, res[n]);
-        assert (resptr);
-        assert (resptr != varptr);
-        const CCTK_REAL fac = factor * factors[n];
-        for (int k=nghostzones[2]; k<lsh[2]-nghostzones[2]; ++k) {
-          for (int j=nghostzones[1]; j<lsh[1]-nghostzones[1]; ++j) {
-            for (int i=nghostzones[0]; i<lsh[0]-nghostzones[0]; ++i) {
-              const int ind = i + lsh[0] * (j + lsh[1] * k);
-              varptr[ind] += fac * resptr[ind];
-            }
-          }
-        }
-      }
-    }
-    
-    void call_applybnds (cGH * const cctkGH)
-    {
-      if (ierr) return;
-      ierr = applybnds (cctkGH, options_table, userdata);
-    }
-    
-  } // namespace common
-  
-  
-  
   void CarpetJacobi_register ()
   {
     int const ierr = TATelliptic_RegisterSolver
@@ -194,7 +70,6 @@ namespace CarpetJacobi {
                           const calcfunc applybnds,
                           void * const userdata)
   {
-    DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
     
     // Error control
@@ -212,19 +87,19 @@ namespace CarpetJacobi {
     assert (res);
     for (int n=0; n<nvars; ++n) {
       assert (var[n] >= 0 && var[n] < CCTK_NumVars());
-      assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF
-              || CCTK_GroupTypeFromVarI(var[n]) == CCTK_ARRAY);
-      assert (CCTK_GroupDimFromVarI(var[n]) <= dim);
+      assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF);
+      assert (CCTK_GroupDimFromVarI(var[n]) == dim);
       assert (CCTK_VarTypeI(var[n]) == CCTK_VARIABLE_REAL);
-      assert (CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(var[n])));
+      assert (CCTK_QueryGroupStorageI(cctkGH,
+                                      CCTK_GroupIndexFromVarI(var[n])));
     }
     for (int n=0; n<nvars; ++n) {
       assert (res[n] >= 0 && res[n] < CCTK_NumVars());
-      assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF
-              || CCTK_GroupTypeFromVarI(res[n]) == CCTK_ARRAY);
-      assert (CCTK_GroupDimFromVarI(res[n]) <= dim);
+      assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF);
+      assert (CCTK_GroupDimFromVarI(res[n]) == dim);
       assert (CCTK_VarTypeI(res[n]) == CCTK_VARIABLE_REAL);
-      assert (CCTK_QueryGroupStorageI(cctkGH, CCTK_GroupIndexFromVarI(res[n])));
+      assert (CCTK_QueryGroupStorageI(cctkGH,
+                                      CCTK_GroupIndexFromVarI(res[n])));
     }
     for (int n=0; n<nvars; ++n) {
       assert (var[n] != res[n]);
@@ -236,47 +111,7 @@ namespace CarpetJacobi {
       }
     }
     
-#if 0
-    // Get domain description
-    cGroup groupdata;
-    ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[0]), &groupdata);
-    assert (!ierr);
-    cGroupDynamicData groupdyndata;
-    ierr = CCTK_GroupDynamicData (cctkGH, CCTK_GroupIndexFromVarI(var[0]),
-                                  &groupdyndata);
-    assert (!ierr);
-    const int thedim = groupdata.dim;
-    assert (thedim>=0 && thedim<=dim);
-    assert (thedim == groupdyndata.dim);
-    int lsh[dim], nghostzones[dim];
-    for (int d=0; d<thedim; ++d) {
-      lsh[d] = groupdyndata.lsh[d];
-      nghostzones[d] = groupdyndata.nghostzones[d];
-    }
-    for (int d=thedim; d<dim; ++d) {
-      lsh[d] = 1;
-      nghostzones[d] = 0;
-    }
-    for (int d=0; d<dim; ++d) {
-      assert (lsh[d]>=0);
-      assert (nghostzones[d] >= 0 && 2*nghostzones[d] <= lsh[d]);
-    }
     
-    // Check all variables
-    for (int n=0; n<nvars; ++n) {
-      ierr = CCTK_GroupData (CCTK_GroupIndexFromVarI(var[n]), &groupdata);
-      assert (!ierr);
-      ierr = CCTK_GroupDynamicData
-        (cctkGH, CCTK_GroupIndexFromVarI(var[n]), &groupdyndata);
-      assert (!ierr);
-      assert (groupdata.dim == thedim);
-      assert (groupdyndata.dim == thedim);
-      for (int d=0; d<thedim; ++d) {
-        assert (groupdyndata.lsh[d] == lsh[d]);
-        assert (groupdyndata.nghostzones[d] == nghostzones[d]);
-      }
-    }
-#endif
     
     assert (options_table >= 0);
     
@@ -310,7 +145,8 @@ namespace CarpetJacobi {
     assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
     ierr = Util_TableGetReal (options_table, &deceleration, "deceleration");
     assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
-    ierr = Util_TableGetRealArray (options_table, nvars, &factors[0], "factors");
+    ierr = Util_TableGetRealArray (options_table,
+                                   nvars, &factors[0], "factors");
     assert (ierr==nvars || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
     
     assert (maxiters >= 0);
@@ -324,197 +160,375 @@ namespace CarpetJacobi {
       assert (factors[n] != 0);
     }
     
-    // Check state
-    // (want global mode, or at least level mode)
-    assert (reflevel==-1 || component==-1);
-    const int saved_reflevel = reflevel;
-    const int saved_mglevel = mglevel;
-    if (reflevel!=-1) {
-      set_mglevel ((cGH *)cctkGH, -1);
-      set_reflevel ((cGH *)cctkGH, -1);
-    }
     
-    // Fill common block
-    common::var = var;
-    common::res = res;
-    common::nvars = nvars;
-    common::options_table = options_table;
-    common::calcres = calcres;
-    common::applybnds = applybnds;
-    common::userdata = userdata;
-    common::factor = factor;
-    common::factors = &factors[0];
     
-    // Init statistics
-    CCTK_REAL norm2 = HUGE_VAL;
-    CCTK_REAL speed = 0.0;
-    CCTK_REAL throttle = 1.0;
-    bool did_hit_min = false;
-    bool did_hit_max = false;
-    const time_t starttime = time(0);
-    time_t nexttime = starttime;
+    assert (is_global_mode() || is_level_mode());
+    
+    // The level to solve on
+    const int solve_level = is_level_mode() ? reflevel : reflevels-1;
     
-    if (verbose || veryverbose) {
-      CCTK_VInfo (CCTK_THORNSTRING, "Solving through adaptive Jacobi iterations");
+    if (solve_level < solve_minlevel) {
+      if (verbose || veryverbose) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "Did not solve because the level is too coarse.");
+      }
+      Util_TableSetInt (options_table, 0, "iters");
+      Util_TableSetReal (options_table, -1.0, "error");
+      
+      // Return as if the solving had not converged
+      return -1;
     }
     
-    int iter;
-    for (iter=0; iter<=maxiters; ++iter) {
+#warning "TODO: assert that all levels are at the same time"
+    
+    // Switch to global mode
+    BEGIN_GLOBAL_MODE(cctkGH) {
       
-      // Calculate residual
-      common::ierr = 0;
-      CallLocalFunction ((cGH *)cctkGH, common::call_calcres);
-      ierr = common::ierr;
-      if (ierr != 0) {
-        if (verbose || veryverbose) {
-          CCTK_VInfo (CCTK_THORNSTRING, "Residual calculation reported error; aborting solver.");
+      // Reset the initial data time
+      global_time = 0;
+      delta_time = 1;
+      * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
+      * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time;
+      BEGIN_REFLEVEL_LOOP(cctkGH) {
+        * 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 / ipow(reflevelfact, reflevelpower));
         }
-        Util_TableSetInt (options_table, iter, "iters");
-        Util_TableDeleteKey (options_table, "error");
-        goto done;
-      }
+      } END_REFLEVEL_LOOP;
       
-      // Save old norm
-      const CCTK_REAL old_norm2 = norm2;
       
-      // Calculate norm
-      common::norm_count = 0;
-      common::norm_l2 = 0;
-      CallLocalFunction ((cGH *)cctkGH, common::norm);
-      const int sum_handle = CCTK_ReductionArrayHandle ("sum");
-      assert (sum_handle >= 0);
-      CCTK_REAL reduce_in[2], reduce_out[2];
-      reduce_in[0] = common::norm_count;
-      reduce_in[1] = common::norm_l2;
-      ierr = CCTK_ReduceLocArrayToArray1D (cctkGH, -1, sum_handle,
-                                           reduce_in, reduce_out, 2,
-                                           CCTK_VARIABLE_REAL);
-      norm2 = sqrt(reduce_out[1] / reduce_out[0]);
       
-      // Calculate convergence speed
-      speed = old_norm2 / norm2;
+      // Init statistics
+      vector<CCTK_REAL> norm_counts (solve_level+1);
+      vector<CCTK_REAL> norm_l2s (solve_level+1);
+      CCTK_REAL norm2 = HUGE_VAL;
+      CCTK_REAL speed = 0.0;
+      CCTK_REAL throttle = 1.0;
+      bool did_hit_min = false;
+      bool did_hit_max = false;
+      const time_t starttime = time(0);
+      time_t nexttime = starttime;
       
-      // Log output
       if (verbose || veryverbose) {
-        const time_t currenttime = time(0);
-        if (iter == maxiters
-#if HAVE_ISNAN
-            || isnan(norm2)
-#endif
-            || norm2 <= minerror
-            || (iter % outevery == 0 && currenttime >= nexttime)) {
-          if (verbose || veryverbose) {
-            if (did_hit_min) {
-              CCTK_VInfo (CCTK_THORNSTRING,
-                          "Convergence factor became too small: artificially kept at minfactor (may lead to instability)");
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "Solving through adaptive Jacobi iterations");
+      }
+      
+      const int icoarsestep
+        = ipow (mglevelfact * maxreflevelfact, reflevelpower);
+      const int istep
+        = ipow (mglevelfact * maxreflevelfact / ipow(reffact, solve_level),
+                reflevelpower);
+      int iter = istep;
+      for (;; iter+=istep) {
+        
+        global_time = 1.0 * iter / ipow(maxreflevelfact, reflevelpower);
+        * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
+//         cout << "CJ global iter " << iter << " time " << global_time << flush << endl;
+        
+        // Calculate residual, smooth, and apply boundary conditions
+        ierr = 0;
+        BEGIN_REFLEVEL_LOOP(cctkGH) {
+          const int do_every
+            = ipow (mglevelfact * (maxreflevelfact/reflevelfact),
+                    reflevelpower);
+          if (reflevel <= solve_level && (iter-istep) % do_every == 0) {
+            
+            // Advance time
+            for (int m=0; m<maps; ++m) {
+              vtt[m]->advance_time (reflevel, mglevel);
             }
-            if (did_hit_max) {
-              CCTK_VInfo (CCTK_THORNSTRING,
-                          "Convergence factor became too large: artificially kept at maxfactor (may lead to slower convergence)");
+            * const_cast<CCTK_REAL *> (& cctkGH->cctk_time)
+              = (1.0 * (iter - istep + do_every)
+                 / ipow(maxreflevelfact, reflevelpower));
+            CycleTimeLevels (cctkGH);
+//             cout << "CJ residual iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl;
+            
+            // Calculate residual
+            BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+              BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+                if (! ierr) {
+                  ierr = calcres (cctkGH, options_table, userdata);
+                }
+              } END_LOCAL_COMPONENT_LOOP;
+            } END_MAP_LOOP;
+            
+            // Smooth and apply boundary conditions
+            CCTK_REAL levfac = 0;
+            for (int d=0; d<dim; ++d) {
+              levfac += 1 /
+                pow (cctkGH->cctk_delta_space[d] / cctkGH->cctk_levfac[d], 2);
             }
-            did_hit_min = false;
-            did_hit_max = false;
+            levfac = 1 / levfac;
+            
+            BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+              BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+                if (! ierr) {
+                  
+                  for (int n=0; n<nvars; ++n) {
+                    CCTK_REAL * restrict varptr
+                      = (CCTK_REAL *)CCTK_VarDataPtrI (cctkGH, 0, var[n]);
+                    assert (varptr);
+                    const CCTK_REAL * restrict resptr
+                      = (const CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, res[n]);
+                    assert (resptr);
+                    assert (resptr != varptr);
+                    const CCTK_REAL fac = factor * factors[n] * levfac;
+                    assert (cctkGH->cctk_dim == 3);
+                    for (int k=cctkGH->cctk_nghostzones[2];
+                         k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2];
+                         ++k) {
+                      for (int j=cctkGH->cctk_nghostzones[1];
+                           j<cctkGH->cctk_lsh[1]-cctkGH->cctk_nghostzones[1];
+                           ++j) {
+                        for (int i=cctkGH->cctk_nghostzones[0];
+                             i<cctkGH->cctk_lsh[0]-cctkGH->cctk_nghostzones[0];
+                             ++i) {
+                          const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                          varptr[ind] += fac * resptr[ind];
+                        }
+                      }
+                    }
+                  } // for n
+                  
+                  // Apply boundary conditions
+                  ierr = applybnds (cctkGH, options_table, userdata);
+                  
+                }
+              } END_LOCAL_COMPONENT_LOOP;
+            } END_MAP_LOOP;
+            
           }
-          if (veryverbose) {
+        } END_REFLEVEL_LOOP;
+        
+        if (ierr) {
+          if (verbose || veryverbose) {
             CCTK_VInfo (CCTK_THORNSTRING,
-                        "Iteration %d, time %g, factor %g, throttle %g, residual %g, speed %g",
-                        iter, (double)(currenttime-starttime), (double)factor, (double)throttle, (double)norm2, (double)speed);
-          } else {
+                        "Residual calculation or boundary conditions reported error; aborting solver.");
+          }
+          Util_TableSetInt (options_table, iter, "iters");
+          Util_TableDeleteKey (options_table, "error");
+          goto done;
+        }
+        
+        // Restrict and calculate norm
+        ierr = 0;
+        BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) {
+          const int do_every
+            = ipow (mglevelfact * (maxreflevelfact/reflevelfact),
+                    reflevelpower);
+          if (reflevel <= solve_level && iter % do_every == 0) {
+            
+//             cout << "CJ restrict iter " << iter << " reflevel " << reflevel << " time " << cctkGH->cctk_time << flush << endl;
+            
+            if (! ierr) {
+              if (reflevel < solve_level) {
+                Restrict (cctkGH);
+              }
+            }
+            
+            BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+              BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+                if (! ierr) {
+                  ierr = applybnds (cctkGH, options_table, userdata);
+                }
+              } END_LOCAL_COMPONENT_LOOP;
+            } END_MAP_LOOP;
+            
+            // Initialise norm
+            CCTK_REAL norm_count = 0;
+            CCTK_REAL norm_l2 = 0;
+            
+            BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+              BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+                if (! ierr) {
+                  
+                  // Calculate norm
+                  for (int n=0; n<nvars; ++n) {
+                    const CCTK_REAL * restrict resptr
+                      = (const CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, res[n]);
+                    assert (resptr);
+                    const CCTK_REAL fac = factors[n];
+                    assert (cctkGH->cctk_dim == 3);
+                    for (int k=cctkGH->cctk_nghostzones[2];
+                         k<cctkGH->cctk_lsh[2]-cctkGH->cctk_nghostzones[2];
+                         ++k) {
+                      for (int j=cctkGH->cctk_nghostzones[1];
+                           j<cctkGH->cctk_lsh[1]-cctkGH->cctk_nghostzones[1];
+                           ++j) {
+                        for (int i=cctkGH->cctk_nghostzones[0];
+                             i<cctkGH->cctk_lsh[0]-cctkGH->cctk_nghostzones[0];
+                             ++i) {
+                          const int ind = CCTK_GFINDEX3D(cctkGH, i, j, k);
+                          ++norm_count;
+                          // TODO: scale the norm by the resolution?
+                          norm_l2 += pow(fac * resptr[ind], 2);
+                        }
+                      }
+                    }
+                  } // for n
+                  
+                }
+              } END_LOCAL_COMPONENT_LOOP;
+            } END_MAP_LOOP;
+            
+            const int sum_handle = CCTK_ReductionArrayHandle ("sum");
+            assert (sum_handle >= 0);
+            CCTK_REAL reduce_in[2], reduce_out[2];
+            reduce_in[0] = norm_count;
+            reduce_in[1] = norm_l2;
+            ierr = CCTK_ReduceLocArrayToArray1D
+              (cctkGH, -1, sum_handle, reduce_in, reduce_out, 2,
+               CCTK_VARIABLE_REAL);
+            norm_counts.at(reflevel) = reduce_out[0];
+            norm_l2s.at(reflevel) = reduce_out[1];
+            
+#warning "TODO"
+#if 0
+            CCTK_OutputVarAs (cctkGH, "WaveToyFO::phi", "phi-ell");
+            CCTK_OutputVarAs (cctkGH, "IDSWTEsimpleFO::residual", "residual-ell");
+#endif
+            
+          }
+        } END_REVERSE_REFLEVEL_LOOP;
+        
+        if (ierr) {
+          if (verbose || veryverbose) {
             CCTK_VInfo (CCTK_THORNSTRING,
-                        "Iteration %d, time %g, residual %g",
-                        iter, (double)(currenttime-starttime), (double)norm2);
+                        "Residual calculation or boundary conditions reported error; aborting solver.");
           }
-          if (outeveryseconds > 0) {
-            while (nexttime <= currenttime) nexttime += outeveryseconds;
+          Util_TableSetInt (options_table, iter, "iters");
+          Util_TableDeleteKey (options_table, "error");
+          goto done;
+        }
+        
+        // Save old norm
+        const CCTK_REAL old_norm2 = norm2;
+        
+        // Calculate new norm
+        CCTK_REAL norm_count = 0;
+        CCTK_REAL norm_l2 = 0;
+        for (int rl=0; rl<=solve_level; ++rl) {
+          norm_count += norm_counts[rl];
+          norm_l2 += norm_l2s[rl];
+        }
+        norm2 = sqrt(norm_l2 / norm_count);
+        
+        // Calculate convergence speed
+        speed = old_norm2 / norm2;
+        
+        // Log output
+        if (verbose || veryverbose) {
+          const time_t currenttime = time(0);
+          if ((iter % icoarsestep == 0 && iter >= maxiters)
+#if HAVE_ISNAN
+              || isnan(norm2)
+#endif
+              || norm2 <= minerror
+              || (iter % outevery == 0 && currenttime >= nexttime)) {
+            if (verbose || veryverbose) {
+              if (did_hit_min) {
+                CCTK_VInfo (CCTK_THORNSTRING,
+                            "Convergence factor became too small: artificially kept at minfactor (may lead to instability)");
+              }
+              if (did_hit_max) {
+                CCTK_VInfo (CCTK_THORNSTRING,
+                            "Convergence factor became too large: artificially kept at maxfactor (may lead to slower convergence)");
+              }
+              did_hit_min = false;
+              did_hit_max = false;
+            }
+            if (veryverbose) {
+              CCTK_VInfo (CCTK_THORNSTRING,
+                          "Iteration %d, time %g, factor %g, throttle %g, residual %g, speed %g",
+                          iter, double(currenttime-starttime),
+                          double(factor), double(throttle), double(norm2),
+                          double(speed));
+            } else {
+              CCTK_VInfo (CCTK_THORNSTRING,
+                          "Iteration %d, time %g, residual %g",
+                          iter, double(currenttime-starttime), double(norm2));
+            }
+            if (outeveryseconds > 0) {
+              while (nexttime <= currenttime) nexttime += outeveryseconds;
+            }
           }
         }
-      }
-      
-      if (iter == maxiters) break;
-      
+        
 #if HAVE_ISNAN
-      // Exit if things go bad
-      if (isnan(norm2)) {
-        if (verbose || veryverbose) {
-          CCTK_VInfo (CCTK_THORNSTRING, "Encountered NaN.  Aborting solver.");
+        // Exit if things go bad
+        if (isnan(norm2)) {
+          if (verbose || veryverbose) {
+            CCTK_VInfo (CCTK_THORNSTRING,
+                        "Encountered NaN.  Aborting solver.");
+          }
+          break;
         }
-        break;
-      }
 #endif
-      
-      // Return when desired accuracy is reached
-      if (norm2 <= minerror) {
-        if (verbose || veryverbose) {
-          CCTK_VInfo (CCTK_THORNSTRING, "Done solving.");
+        
+        // Return when desired accuracy is reached
+        if (norm2 <= minerror) {
+          if (verbose || veryverbose) {
+            CCTK_VInfo (CCTK_THORNSTRING, "Done solving.");
+          }
+          Util_TableSetInt (options_table, iter, "iters");
+          Util_TableSetReal (options_table, norm2, "error");
+          ierr = 0;
+          goto done;
         }
-        Util_TableSetInt (options_table, iter, "iters");
-        Util_TableSetReal (options_table, norm2, "error");
-        ierr = 0;
-        goto done;
-      }
+        
+        // Exit if the maximum number of iterations has been reached
+        if (iter % icoarsestep == 0 && iter >= maxiters) {
+          break;
+        }
+        
+        // Adjust speed
+        if (speed > 1.0) {
+          throttle = acceleration;
+        } else {
+          throttle = deceleration;
+        }
+        factor *= throttle;
+        if (factor < minfactor) {
+          did_hit_min = true;
+          factor = minfactor;
+        }
+        if (factor > maxfactor) {
+          did_hit_max = true;
+          factor = maxfactor;
+        }
+        
+      } // for iter
       
-      // Adjust speed
-      if (speed > 1.0) {
-        throttle = acceleration;
-      } else {
-        throttle = deceleration;
-      }
-      factor *= throttle;
-      if (factor < minfactor) {
-        did_hit_min = true;
-        factor = minfactor;
-      }
-      if (factor > maxfactor) {
-        did_hit_max = true;
-        factor = maxfactor;
+      // Did not solve
+      if (verbose || veryverbose) {
+        CCTK_VInfo (CCTK_THORNSTRING, "Did not converge.");
       }
+      Util_TableSetInt (options_table, iter, "iters");
+      Util_TableSetReal (options_table, norm2, "error");
+      ierr = -1;
+    
+    done: ;
       
-      // Smooth
-      CallLocalFunction ((cGH *)cctkGH, common::smooth);
       
-      // Apply boundary conditions
-      common::ierr = 0;
-      CallLocalFunction ((cGH *)cctkGH, common::call_applybnds);
-      ierr = common::ierr;
-      if (ierr != 0) {
-        if (verbose || veryverbose) {
-          CCTK_VInfo (CCTK_THORNSTRING, "Boundary conditions reported error; aborting solver.");
-        }
-        Util_TableSetInt (options_table, iter, "iters");
-        Util_TableDeleteKey (options_table, "error");
-        goto done;
-      }
       
-#if 0
-      // Restrict
-      BEGIN_REVERSE_REFLEVEL_LOOP(cctkGH) {
-        BEGIN_MGLEVEL_LOOP(cctkGH) {
-          Restrict (cctkGH);
-          // TODO: do something here
-//           CCTK_ScheduleTraverse ("CCTK_POSTRESTRICT", (cGH *)cctkGH, CallFunction);
-          CallLocalFunction ((cGH *)cctkGH, common::call_applybnds);
-          assert (! common::ierr);
-        } END_MGLEVEL_LOOP;
-      } END_REVERSE_REFLEVEL_LOOP;
-#endif
+      // Reset the initial time
+#warning "TODO: reset the initial time a bit more intelligently"
+      global_time = 0;
+      delta_time = 1;
+      * const_cast<CCTK_REAL *> (& cctkGH->cctk_time) = global_time;
+      * const_cast<CCTK_REAL *> (& cctkGH->cctk_delta_time) = delta_time;
+      BEGIN_REFLEVEL_LOOP(cctkGH) {
+        * 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);
+        }
+      } END_REFLEVEL_LOOP;
       
-    } // for iter
-    
-    // Did not solve
-    if (verbose || veryverbose) {
-      CCTK_VInfo (CCTK_THORNSTRING, "Did not converge.");
-    }
-    Util_TableSetInt (options_table, iter, "iters");
-    Util_TableSetReal (options_table, norm2, "error");
-    ierr = -1;
-    
-  done:
-    
-    // Restore state
-    if (reflevel!=saved_reflevel) {
-      set_reflevel ((cGH *)cctkGH, saved_reflevel);
-      set_mglevel ((cGH *)cctkGH, saved_mglevel);
-    }
+    } END_GLOBAL_MODE;
     
     return ierr;
   }