Add elliptic solver thorn that works with Carpet.

darcs-hash:20030902123557-07bb3-67dd485b689f921257309f4eeb7389821fd5fa7a.gz

4 files changed, 345 insertions, 375 deletions

diff --git a/CarpetDev/CarpetJacobi/par/charge_te_cjac.par b/CarpetDev/CarpetJacobi/par/charge_te_cjac.par
index 8db9c4cf6..3e77742dd 100644
--- a/CarpetDev/CarpetJacobi/par/charge_te_cjac.par
+++ b/CarpetDev/CarpetJacobi/par/charge_te_cjac.par
@@ -1,7 +1,7 @@
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac.par,v 1.4 2004/03/23 11:59:20 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac.par,v 1.1 2003/09/02 14:35:58 schnetter Exp $
 # Erik Schnetter <schnetter@uni-tuebingen.de>
 
-ActiveThorns = "Boundary CartGrid3D CoordBase SymBase CarpetIOASCII IOBasic IOUtil Time Carpet CarpetLib CarpetRegrid CarpetReduce NaNCatcher IDScalarWave WaveToyC IDSWTEsimple TATelliptic CarpetJacobi"
+ActiveThorns = "Boundary CartGrid3D CarpetIOASCII IOBasic IOUtil Time Carpet CarpetLib CarpetRegrid CarpetReduce NaNCatcher IDScalarWave WaveToyC IDSWTEsimple TATelliptic CarpetJacobi"
 
 !DESC "Charged sphere initial data, solved with TATelliptic/CarpetJacobi"
 
@@ -31,7 +31,7 @@ CarpetRegrid::outerbounds       = "[[ [[1,0],[1,0],[1,0]] ]]"
 
 grid::type                      = byspacing
 grid::domain                    = octant
-grid::dxyz                      = 1.2
+grid::dxyz                      = 0.6
 grid::avoid_origin              = no
 
 IDScalarWave::initial_data      = charge-TATelliptic-simple
diff --git a/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par b/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par
index b26fe4d48..6f0131e74 100644
--- a/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par
+++ b/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par
@@ -1,7 +1,7 @@
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par,v 1.4 2004/03/23 11:59:20 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par,v 1.1 2003/09/02 14:35:58 schnetter Exp $
 # Erik Schnetter <schnetter@uni-tuebingen.de>
 
-ActiveThorns = "Boundary CartGrid3D CoordBase SymBase CarpetIOASCII IOBasic IOUtil Time Carpet CarpetLib CarpetRegrid CarpetReduce NaNCatcher IDScalarWave WaveToyC IDSWTEsimple TATelliptic CarpetJacobi"
+ActiveThorns = "Boundary CartGrid3D CarpetIOASCII IOBasic IOUtil Time Carpet CarpetLib CarpetRegrid CarpetReduce NaNCatcher IDScalarWave WaveToyC IDSWTEsimple TATelliptic CarpetJacobi"
 
 !DESC "Charged sphere initial data, solved with TATelliptic/CarpetJacobi"
 
@@ -31,7 +31,7 @@ CarpetRegrid::outerbounds       = "[[ [[1,0],[1,0],[1,0]] ]]"
 
 grid::type                      = byspacing
 grid::domain                    = octant
-grid::dxyz                      = 1.2
+grid::dxyz                      = 0.6
 grid::avoid_origin              = no
 
 IDScalarWave::initial_data      = charge-TATelliptic-simple
diff --git a/CarpetDev/CarpetJacobi/param.ccl b/CarpetDev/CarpetJacobi/param.ccl
index bdeecf41b..ce69dbc91 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.2 2004/01/25 14:57:29 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/param.ccl,v 1.1 2003/09/02 14:35:57 schnetter Exp $
 
 BOOLEAN verbose "Produce log output while running" STEERABLE=always
 {
@@ -18,15 +18,3 @@ 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 c326b66cd..a4686af0a 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.3 2004/01/25 14:57:29 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/Jacobi.cc,v 1.1 2003/09/02 14:35:58 schnetter Exp $
 
 #include <cassert>
 #include <cmath>
@@ -14,18 +14,14 @@
 #include "util_ErrorCodes.h"
 #include "util_Table.h"
 
-#include "TATelliptic.h"
-
 #include "carpet.hh"
-
-#include "th.hh"
+#include "TATelliptic.h"
 
 
 
 namespace Carpet {
   // TODO: fix this
-  void CycleTimeLevels (const cGH* cctkGH);
-  void Restrict (const cGH* cctkGH);
+  void Restrict (const cGH* cgh);
 };
 
 
@@ -52,6 +48,134 @@ 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
@@ -70,6 +194,7 @@ namespace CarpetJacobi {
                           const calcfunc applybnds,
                           void * const userdata)
   {
+    DECLARE_CCTK_ARGUMENTS;
     DECLARE_CCTK_PARAMETERS;
     
     // Error control
@@ -87,19 +212,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);
-      assert (CCTK_GroupDimFromVarI(var[n]) == dim);
+      assert (CCTK_GroupTypeFromVarI(var[n]) == CCTK_GF
+              || CCTK_GroupTypeFromVarI(var[n]) == CCTK_ARRAY);
+      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);
-      assert (CCTK_GroupDimFromVarI(res[n]) == dim);
+      assert (CCTK_GroupTypeFromVarI(res[n]) == CCTK_GF
+              || CCTK_GroupTypeFromVarI(res[n]) == CCTK_ARRAY);
+      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]);
@@ -111,7 +236,47 @@ 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);
     
@@ -145,8 +310,7 @@ 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);
@@ -160,375 +324,193 @@ 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;
+    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];
     
-    assert (is_global_mode() || is_level_mode());
-    
-    // The level to solve on
-    const int solve_level = is_level_mode() ? reflevel : reflevels-1;
+    // 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;
     
-    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;
+    if (verbose || veryverbose) {
+      CCTK_VInfo (CCTK_THORNSTRING, "Solving through adaptive Jacobi iterations");
     }
     
-#warning "TODO: assert that all levels are at the same time"
-    
-    // Switch to global mode
-    BEGIN_GLOBAL_MODE(cctkGH) {
+    int iter;
+    for (iter=0; iter<=maxiters; ++iter) {
       
-      // 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));
+      // 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.");
         }
-      } END_REFLEVEL_LOOP;
+        Util_TableSetInt (options_table, iter, "iters");
+        Util_TableDeleteKey (options_table, "error");
+        goto done;
+      }
       
+      // 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]);
       
-      // 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;
+      // Calculate convergence speed
+      speed = old_norm2 / norm2;
       
+      // Log output
       if (verbose || veryverbose) {
-        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);
-            }
-            * 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);
-            }
-            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;
-            
-          }
-        } END_REFLEVEL_LOOP;
-        
-        if (ierr) {
-          if (verbose || veryverbose) {
-            CCTK_VInfo (CCTK_THORNSTRING,
-                        "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,
-                        "Residual calculation or boundary conditions reported error; aborting solver.");
-          }
-          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)
+        const time_t currenttime = time(0);
+        if (iter == maxiters
 #if HAVE_ISNAN
-              || isnan(norm2)
+            || 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 {
+            || norm2 <= minerror
+            || (iter % outevery == 0 && currenttime >= nexttime)) {
+          if (verbose || veryverbose) {
+            if (did_hit_min) {
               CCTK_VInfo (CCTK_THORNSTRING,
-                          "Iteration %d, time %g, residual %g",
-                          iter, double(currenttime-starttime), double(norm2));
+                          "Convergence factor became too small: artificially kept at minfactor (may lead to instability)");
             }
-            if (outeveryseconds > 0) {
-              while (nexttime <= currenttime) nexttime += outeveryseconds;
+            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 HAVE_ISNAN
-        // Exit if things go bad
-        if (isnan(norm2)) {
-          if (verbose || veryverbose) {
+          if (veryverbose) {
             CCTK_VInfo (CCTK_THORNSTRING,
-                        "Encountered NaN.  Aborting solver.");
+                        "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);
           }
-          break;
-        }
-#endif
-        
-        // Return when desired accuracy is reached
-        if (norm2 <= minerror) {
-          if (verbose || veryverbose) {
-            CCTK_VInfo (CCTK_THORNSTRING, "Done solving.");
+          if (outeveryseconds > 0) {
+            while (nexttime <= currenttime) nexttime += outeveryseconds;
           }
-          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;
+      }
+      
+      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.");
         }
-        
-      } // for iter
+        break;
+      }
+#endif
       
-      // Did not solve
-      if (verbose || veryverbose) {
-        CCTK_VInfo (CCTK_THORNSTRING, "Did not converge.");
+      // 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 = -1;
-    
-    done: ;
       
+      // 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;
+      }
       
+      // Smooth
+      CallLocalFunction ((cGH *)cctkGH, common::smooth);
       
-      // 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);
+      // 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.");
         }
-      } END_REFLEVEL_LOOP;
+        Util_TableSetInt (options_table, iter, "iters");
+        Util_TableDeleteKey (options_table, "error");
+        goto done;
+      }
       
-    } END_GLOBAL_MODE;
+#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
+      
+    } // 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
+    set_reflevel ((cGH *)cctkGH, saved_reflevel);
+    set_mglevel ((cGH *)cctkGH, saved_mglevel);
     
     return ierr;
   }