CarpetMG: A multigrid solver for Carpet

darcs-hash:20041124061259-891bb-566d93667d53c7746838a83c4171f6874dd496f2.gz

12 files changed, 1933 insertions, 0 deletions

diff --git a/CarpetDev/CarpetMG/README b/CarpetDev/CarpetMG/README
new file mode 100644
index 000000000..86232041b
--- /dev/null
+++ b/CarpetDev/CarpetMG/README
@@ -0,0 +1,8 @@
+Cactus Code Thorn CarpetMG
+Thorn Author(s)     : Erik Schnetter <schnetter@aei.mpg.de>
+Thorn Maintainer(s) : Erik Schnetter <schnetter@aei.mpg.de>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+A generic multigrid solver for Carpet.
diff --git a/CarpetDev/CarpetMG/configuration.ccl b/CarpetDev/CarpetMG/configuration.ccl
new file mode 100644
index 000000000..0f47593ed
--- /dev/null
+++ b/CarpetDev/CarpetMG/configuration.ccl
@@ -0,0 +1,3 @@
+# Configuration definition for thorn CarpetMG
+
+REQUIRES Carpet CarpetLib TATelliptic
diff --git a/CarpetDev/CarpetMG/doc/documentation.tex b/CarpetDev/CarpetMG/doc/documentation.tex
new file mode 100644
index 000000000..3d1a4ac7c
--- /dev/null
+++ b/CarpetDev/CarpetMG/doc/documentation.tex
@@ -0,0 +1,143 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header: /cactusdevcvs/Cactus/doc/ThornGuide/template.tex,v 1.12 2004/01/07 20:12:39 rideout Exp $
+%
+%  Thorn documentation in the latex file doc/documentation.tex
+%  will be included in ThornGuides built with the Cactus make system.
+%  The scripts employed by the make system automatically include
+%  pages about variables, parameters and scheduling parsed from the
+%  relevant thorn CCL files.
+%
+%  This template contains guidelines which help to assure that your
+%  documentation will be correctly added to ThornGuides. More
+%  information is available in the Cactus UsersGuide.
+%
+%  Guidelines:
+%   - Do not change anything before the line
+%       % START CACTUS THORNGUIDE",
+%     except for filling in the title, author, date, etc. fields.
+%        - Each of these fields should only be on ONE line.
+%        - Author names should be separated with a \\ or a comma.
+%   - You can define your own macros, but they must appear after
+%     the START CACTUS THORNGUIDE line, and must not redefine standard
+%     latex commands.
+%   - To avoid name clashes with other thorns, 'labels', 'citations',
+%     'references', and 'image' names should conform to the following
+%     convention:
+%       ARRANGEMENT_THORN_LABEL
+%     For example, an image wave.eps in the arrangement CactusWave and
+%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+%   - Graphics should only be included using the graphicx package.
+%     More specifically, with the "\includegraphics" command.  Do
+%     not specify any graphic file extensions in your .tex file. This
+%     will allow us to create a PDF version of the ThornGuide
+%     via pdflatex.
+%   - References should be included with the latex "\bibitem" command.
+%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+%   - Do not use \appendix, instead include any appendices you need as
+%     standard sections.
+%   - For the benefit of our Perl scripts, and for future extensions,
+%     please use simple latex.
+%
+% *======================================================================*
+%
+% Example of including a graphic image:
+%    \begin{figure}[ht]
+% 	\begin{center}
+%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% 	\end{center}
+% 	\caption{Illustration of this and that}
+% 	\label{MyArrangement_MyThorn_MyLabel}
+%    \end{figure}
+%
+% Example of using a label:
+%   \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+%    \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+%   \bibitem{MyArrangement_MyThorn_Author99}
+%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
+%   1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================*
+
+% If you are using CVS use this line to give version information
+% $Header: /cactusdevcvs/Cactus/doc/ThornGuide/template.tex,v 1.12 2004/01/07 20:12:39 rideout Exp $
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a
+%  thorn without the Cactus Flesh download this from the Cactus
+%  homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Erik Schnetter \textless schnetter@aei.mpg.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{CarpetMG}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+\date{November 19 2004}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here
+%   \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+A generic multigrid solver for Carpet.
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Introduction}
+
+\section{Physical System}
+
+\section{Numerical Implementation}
+
+\section{Using This Thorn}
+
+\subsection{Obtaining This Thorn}
+
+\subsection{Basic Usage}
+
+\subsection{Special Behaviour}
+
+\subsection{Interaction With Other Thorns}
+
+\subsection{Examples}
+
+\subsection{Support and Feedback}
+
+\section{History}
+
+\subsection{Thorn Source Code}
+
+\subsection{Thorn Documentation}
+
+\subsection{Acknowledgements}
+
+
+\begin{thebibliography}{9}
+
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/CarpetDev/CarpetMG/interface.ccl b/CarpetDev/CarpetMG/interface.ccl
new file mode 100644
index 000000000..1023297ee
--- /dev/null
+++ b/CarpetDev/CarpetMG/interface.ccl
@@ -0,0 +1,14 @@
+# Interface definition for thorn CarpetMG
+
+IMPLEMENTS: CarpetMG
+
+INHERITS: TATelliptic
+
+USES INCLUDE HEADER: TATelliptic.h
+USES INCLUDE HEADER: carpet.hh
+
+
+
+CCTK_INT FUNCTION \
+    SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH)
+REQUIRES FUNCTION SymmetryTableHandleForGrid
diff --git a/CarpetDev/CarpetMG/par/charge-ml2.par b/CarpetDev/CarpetMG/par/charge-ml2.par
new file mode 100644
index 000000000..5ea07bc51
--- /dev/null
+++ b/CarpetDev/CarpetMG/par/charge-ml2.par
@@ -0,0 +1,134 @@
+!DESC "Charged sphere initial data, solved with CarpetMG"
+
+Cactus::cctk_full_warnings = yes
+Cactus::cctk_timer_output  = full
+
+Cactus::cctk_itlast = 5120
+
+
+
+ActiveThorns = "NaNCatcher"
+
+
+
+ActiveThorns = "IOUtil"
+
+IO::out_dir = $parfile
+
+
+
+ActiveThorns = "Carpet CarpetLib CarpetInterp CarpetReduce CarpetSlab"
+
+Carpet::domain_from_coordbase = yes
+Carpet::max_refinement_levels = 11
+
+Carpet::ghost_size   = 2
+Carpet::buffer_width = 4
+
+Carpet::prolongation_order_space = 3
+Carpet::prolongation_order_time  = 2
+
+Carpet::convergence_level = 0
+
+Carpet::init_3_timelevels = yes
+
+
+
+ActiveThorns = "TATelliptic CarpetMG"
+
+CarpetMG::verbose = yes
+#CarpetMG::veryverbose = yes
+
+
+
+ActiveThorns = "Boundary CartGrid3D CoordBase ReflectionSymmetry SymBase"
+
+CoordBase::domainsize = minmax
+CoordBase::spacing    = numcells
+
+CoordBase::xmin =  0.0
+CoordBase::ymin =  0.0
+CoordBase::zmin =  0.0
+CoordBase::xmax = 20.0
+CoordBase::ymax = 20.0
+CoordBase::zmax = 20.0
+
+CoordBase::ncells_x = 20
+CoordBase::ncells_y = 20
+CoordBase::ncells_z = 20
+
+CoordBase::boundary_size_x_lower = 2
+CoordBase::boundary_size_y_lower = 2
+CoordBase::boundary_size_z_lower = 2
+
+CoordBase::boundary_shiftout_x_lower = 1
+CoordBase::boundary_shiftout_y_lower = 1
+CoordBase::boundary_shiftout_z_lower = 1
+
+CartGrid3D::type         = coordbase
+CartGrid3D::avoid_origin = no
+
+ReflectionSymmetry::reflection_x = yes
+ReflectionSymmetry::reflection_y = yes
+ReflectionSymmetry::reflection_z = yes
+
+ReflectionSymmetry::avoid_origin_x = no
+ReflectionSymmetry::avoid_origin_y = no
+ReflectionSymmetry::avoid_origin_z = no
+
+
+
+ActiveThorns = "CarpetRegrid"
+
+CarpetRegrid::refinement_levels        = 2
+CarpetRegrid::smart_outer_boundaries   = yes
+
+CarpetRegrid::keep_same_grid_structure = yes
+CarpetRegrid::refined_regions          = manual-coordinate-list
+CarpetRegrid::coordinates = "
+        [ [ ([0.0,0.0,0.0]:[20.0,20.0,20.0]:[0.5,0.5,0.5]) ] ]
+"
+
+
+
+ActiveThorns = "WaveToyC"
+
+WaveToyC::bound	= radiation
+
+
+
+ActiveThorns = "IDScalarWave IDSWTEcarpet"
+
+IDScalarWave::initial_data = charge-TATelliptic-carpet
+
+IDSWTEcarpet::solver  = CarpetMG
+IDSWTEcarpet::options = ""
+IDSWTEcarpet::radius  = 5.5
+IDSWTEcarpet::charge  = 1.0
+
+
+
+ActiveThorns = "Time"
+ 
+Time::dtfac = 0.5
+
+
+
+ActiveThorns = "IOBasic"
+
+IOBasic::outInfo_every = 1
+IOBasic::outInfo_vars  = "wavetoy::phi"
+
+
+
+ActiveThorns = "CarpetIOScalar"
+
+IOScalar::outScalar_every = 1
+IOScalar::outScalar_vars  = "wavetoy::phi IDSWTEcarpet::residuals"
+
+
+
+ActiveThorns = "CarpetIOASCII"
+
+IOASCII::out1D_every = 1
+IOASCII::out1D_vars  = "wavetoy::phi IDSWTEcarpet::residuals"
diff --git a/CarpetDev/CarpetMG/par/charge-ml4.par b/CarpetDev/CarpetMG/par/charge-ml4.par
new file mode 100644
index 000000000..78767862b
--- /dev/null
+++ b/CarpetDev/CarpetMG/par/charge-ml4.par
@@ -0,0 +1,136 @@
+!DESC "Charged sphere initial data, solved with CarpetMG"
+
+Cactus::cctk_full_warnings = yes
+Cactus::cctk_timer_output  = full
+
+Cactus::cctk_itlast = 5120
+
+
+
+ActiveThorns = "NaNCatcher"
+
+
+
+ActiveThorns = "IOUtil"
+
+IO::out_dir = $parfile
+
+
+
+ActiveThorns = "Carpet CarpetLib CarpetInterp CarpetReduce CarpetSlab"
+
+Carpet::domain_from_coordbase = yes
+Carpet::max_refinement_levels = 13
+
+Carpet::ghost_size   = 2
+Carpet::buffer_width = 4
+
+Carpet::prolongation_order_space = 3
+Carpet::prolongation_order_time  = 2
+
+Carpet::convergence_level = 0
+
+Carpet::init_3_timelevels = yes
+
+
+
+ActiveThorns = "TATelliptic CarpetMG"
+
+CarpetMG::verbose = yes
+#CarpetMG::veryverbose = yes
+
+
+
+ActiveThorns = "Boundary CartGrid3D CoordBase ReflectionSymmetry SymBase"
+
+CoordBase::domainsize = minmax
+CoordBase::spacing    = numcells
+
+CoordBase::xmin =  0.0
+CoordBase::ymin =  0.0
+CoordBase::zmin =  0.0
+CoordBase::xmax = 20.0
+CoordBase::ymax = 20.0
+CoordBase::zmax = 20.0
+
+CoordBase::ncells_x = 5
+CoordBase::ncells_y = 5
+CoordBase::ncells_z = 5
+
+CoordBase::boundary_size_x_lower = 2
+CoordBase::boundary_size_y_lower = 2
+CoordBase::boundary_size_z_lower = 2
+
+CoordBase::boundary_shiftout_x_lower = 1
+CoordBase::boundary_shiftout_y_lower = 1
+CoordBase::boundary_shiftout_z_lower = 1
+
+CartGrid3D::type         = coordbase
+CartGrid3D::avoid_origin = no
+
+ReflectionSymmetry::reflection_x = yes
+ReflectionSymmetry::reflection_y = yes
+ReflectionSymmetry::reflection_z = yes
+
+ReflectionSymmetry::avoid_origin_x = no
+ReflectionSymmetry::avoid_origin_y = no
+ReflectionSymmetry::avoid_origin_z = no
+
+
+
+ActiveThorns = "CarpetRegrid"
+
+CarpetRegrid::refinement_levels        = 4
+CarpetRegrid::smart_outer_boundaries   = yes
+
+CarpetRegrid::keep_same_grid_structure = yes
+CarpetRegrid::refined_regions          = manual-coordinate-list
+CarpetRegrid::coordinates = "
+        [ [ ([0.0,0.0,0.0]:[20.0,20.0,20.0]:[2.0,2.0,2.0]) ],
+          [ ([0.0,0.0,0.0]:[20.0,20.0,20.0]:[1.0,1.0,1.0]) ],
+          [ ([0.0,0.0,0.0]:[20.0,20.0,20.0]:[0.5,0.5,0.5]) ] ]
+"
+
+
+
+ActiveThorns = "WaveToyC"
+
+WaveToyC::bound	= radiation
+
+
+
+ActiveThorns = "IDScalarWave IDSWTEcarpet"
+
+IDScalarWave::initial_data = charge-TATelliptic-carpet
+
+IDSWTEcarpet::solver  = CarpetMG
+IDSWTEcarpet::options = "mgiters=10"
+IDSWTEcarpet::radius  = 5.5
+IDSWTEcarpet::charge  = 1.0
+
+
+
+ActiveThorns = "Time"
+ 
+Time::dtfac = 0.5
+
+
+
+ActiveThorns = "IOBasic"
+
+IOBasic::outInfo_every = 1
+IOBasic::outInfo_vars  = "wavetoy::phi"
+
+
+
+ActiveThorns = "CarpetIOScalar"
+
+IOScalar::outScalar_every = 1
+IOScalar::outScalar_vars  = "wavetoy::phi IDSWTEcarpet::residuals"
+
+
+
+ActiveThorns = "CarpetIOASCII"
+
+IOASCII::out1D_every = 1
+IOASCII::out1D_vars  = "wavetoy::phi IDSWTEcarpet::residuals"
diff --git a/CarpetDev/CarpetMG/par/charge.par b/CarpetDev/CarpetMG/par/charge.par
new file mode 100644
index 000000000..675b2ecfb
--- /dev/null
+++ b/CarpetDev/CarpetMG/par/charge.par
@@ -0,0 +1,121 @@
+!DESC "Charged sphere initial data, solved with CarpetMG"
+
+Cactus::cctk_full_warnings = yes
+Cactus::cctk_timer_output  = full
+
+Cactus::cctk_itlast = 5120
+
+
+
+ActiveThorns = "NaNCatcher"
+
+
+
+ActiveThorns = "IOUtil"
+
+IO::out_dir = $parfile
+
+
+
+ActiveThorns = "Carpet CarpetLib CarpetInterp CarpetReduce CarpetSlab"
+
+Carpet::domain_from_coordbase = yes
+Carpet::max_refinement_levels = 10
+
+Carpet::ghost_size   = 2
+Carpet::buffer_width = 4
+
+Carpet::prolongation_order_space = 3
+Carpet::prolongation_order_time  = 2
+
+Carpet::convergence_level = 0
+
+Carpet::init_3_timelevels = yes
+
+
+
+ActiveThorns = "TATelliptic CarpetMG"
+
+CarpetMG::verbose = yes
+#CarpetMG::veryverbose = yes
+
+
+
+ActiveThorns = "Boundary CartGrid3D CoordBase ReflectionSymmetry SymBase"
+
+CoordBase::domainsize = minmax
+CoordBase::spacing    = numcells
+
+CoordBase::xmin =  0.0
+CoordBase::ymin =  0.0
+CoordBase::zmin =  0.0
+CoordBase::xmax = 20.0
+CoordBase::ymax = 20.0
+CoordBase::zmax = 20.0
+
+CoordBase::ncells_x = 40
+CoordBase::ncells_y = 40
+CoordBase::ncells_z = 40
+
+CoordBase::boundary_size_x_lower = 2
+CoordBase::boundary_size_y_lower = 2
+CoordBase::boundary_size_z_lower = 2
+
+CoordBase::boundary_shiftout_x_lower = 1
+CoordBase::boundary_shiftout_y_lower = 1
+CoordBase::boundary_shiftout_z_lower = 1
+
+CartGrid3D::type         = coordbase
+CartGrid3D::avoid_origin = no
+
+ReflectionSymmetry::reflection_x = yes
+ReflectionSymmetry::reflection_y = yes
+ReflectionSymmetry::reflection_z = yes
+
+ReflectionSymmetry::avoid_origin_x = no
+ReflectionSymmetry::avoid_origin_y = no
+ReflectionSymmetry::avoid_origin_z = no
+
+
+
+ActiveThorns = "WaveToyC"
+
+WaveToyC::bound	= radiation
+
+
+
+ActiveThorns = "IDScalarWave IDSWTEcarpet"
+
+IDScalarWave::initial_data = charge-TATelliptic-carpet
+
+IDSWTEcarpet::solver  = CarpetMG
+IDSWTEcarpet::options = ""
+IDSWTEcarpet::radius  = 5.5
+IDSWTEcarpet::charge  = 1.0
+
+
+
+ActiveThorns = "Time"
+ 
+Time::dtfac = 0.5
+
+
+
+ActiveThorns = "IOBasic"
+
+IOBasic::outInfo_every = 1
+IOBasic::outInfo_vars  = "wavetoy::phi"
+
+
+
+ActiveThorns = "CarpetIOScalar"
+
+IOScalar::outScalar_every = 1
+IOScalar::outScalar_vars  = "wavetoy::phi"
+
+
+
+ActiveThorns = "CarpetIOASCII"
+
+IOASCII::out1D_every = 1
+IOASCII::out1D_vars  = "wavetoy::phi"
diff --git a/CarpetDev/CarpetMG/param.ccl b/CarpetDev/CarpetMG/param.ccl
new file mode 100644
index 000000000..aba943937
--- /dev/null
+++ b/CarpetDev/CarpetMG/param.ccl
@@ -0,0 +1,16 @@
+# Parameter definitions for thorn CarpetMG
+
+BOOLEAN verbose "Produce log output while running" STEERABLE=always
+{
+} "no"
+
+BOOLEAN veryverbose "Produce much log output while running" STEERABLE=always
+{
+} "no"
+
+
+
+STRING direct_solver "Direct solver for the coarsest grid"
+{
+  ".*" :: "must be a registered solver for TATelliptic"
+} "TATJacobi"
diff --git a/CarpetDev/CarpetMG/schedule.ccl b/CarpetDev/CarpetMG/schedule.ccl
new file mode 100644
index 000000000..849293412
--- /dev/null
+++ b/CarpetDev/CarpetMG/schedule.ccl
@@ -0,0 +1,11 @@
+# Schedule definitions for thorn CarpetMG
+
+SCHEDULE CarpetMG_register AT startup
+{
+  LANG: C
+} "Register the elliptic solver"
+
+SCHEDULE CarpetMG_paramcheck AT paramcheck
+{
+  LANG: C
+} "Check parameters"
diff --git a/CarpetDev/CarpetMG/src/README b/CarpetDev/CarpetMG/src/README
new file mode 100644
index 000000000..dd65e6254
--- /dev/null
+++ b/CarpetDev/CarpetMG/src/README
@@ -0,0 +1,38 @@
+equation: L(varc) = rhsc
+
+
+
+resf <- L(varf)
+resf <- resf - rhsf
+
+varc <- varc // R(varf)
+savc <- varc
+resc <- 0
+resc <- resc // R(resf)
+rhsc <- resc
+resc <- L(varc)
+rhsc <- rhsc - resc
+
+varc <- L(varc) = rhsc
+
+resc <- varc
+resc <- resc - savc
+resf <- 0
+resf <- P(resc)
+varf <- varf + resf
+
+
+
+rhsc <- R(L(varf) - rhsf) - L(R(varf))
+varc <- L(varc) = R(L(varf) - rhsf) - L(R(varf))
+varf <- varf + P(varc - R(varf))
+
+delc <- L(R(varf) + delc) = R(L(varf) - rhsf) - L(R(varf))
+delf = P(delc)
+
+
+
+L(varf       ) - rhcf = resf
+
+
+L(varf + delf) - rhcf = 0
diff --git a/CarpetDev/CarpetMG/src/make.code.defn b/CarpetDev/CarpetMG/src/make.code.defn
new file mode 100644
index 000000000..2d9211a7c
--- /dev/null
+++ b/CarpetDev/CarpetMG/src/make.code.defn
@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn CarpetMG -*-Makefile-*-
+
+# Source files in this directory
+SRCS = mg.cc
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/CarpetDev/CarpetMG/src/mg.cc b/CarpetDev/CarpetMG/src/mg.cc
new file mode 100644
index 000000000..ab1132758
--- /dev/null
+++ b/CarpetDev/CarpetMG/src/mg.cc
@@ -0,0 +1,1302 @@
+#include <algorithm>
+#include <cassert>
+#include <climits>
+#include <cmath>
+#include <vector>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "carpet.hh"
+#include "TATelliptic.h"
+
+
+
+namespace CarpetMG {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  char const * const solver_name = "CarpetMG";
+  
+  
+  
+  struct options_t {
+    // Width of outer boundary
+    CCTK_INT bndwidth;
+    
+    // Symmetry information
+    CCTK_INT symmetry_handle[2*dim];
+    CCTK_INT symmetry_zone_width[2*dim];
+    
+    // Coarsest level on which the full multigrid algorithm should start
+    CCTK_INT firstlevel;
+    
+    // Maximum number of direct solver iterations
+    CCTK_INT maxiters;
+    
+    // Weight factors for the equations
+    CCTK_REAL factor;
+    vector<CCTK_REAL> factors;
+    
+    // Maximum number of multigrid iterations (per level) and smoothing steps
+    CCTK_INT mgiters;
+    CCTK_INT presteps;
+    CCTK_INT poststeps;
+  };
+  
+  
+  
+  extern "C" {
+    
+    int
+    CarpetMG_register ();
+    
+    void
+    CarpetMG_paramcheck (CCTK_ARGUMENTS);
+    
+  }
+  
+  int
+  solve (cGH const * restrict const cctkGH,
+         int const * restrict const var,
+         int const * restrict const res,
+         int const nvar,
+         int const options_table,
+         calcfunc const calcres,
+         calcfunc const applybnds,
+         void * const userdata);
+  int
+  multigrid (cGH const * restrict const cctkGH,
+             vector<CCTK_INT> const & var,
+             vector<CCTK_INT> const & res,
+             vector<CCTK_INT> const & rhs,
+             vector<CCTK_INT> const & sav,
+             vector<CCTK_INT> const & wgt,
+             vector<CCTK_INT> const & aux,
+             int const options_table,
+             calcfunc const calcres,
+             calcfunc const applybnds,
+             void * const userdata,
+             options_t const & options,
+             CCTK_REAL const minerror);
+  
+  int
+  direct_solve (cGH const * restrict const cctkGH,
+                vector<CCTK_INT> const & var,
+                vector<CCTK_INT> const & res,
+                vector<CCTK_INT> const & rhs,
+                vector<CCTK_INT> const & wgt,
+                vector<CCTK_INT> const & aux,
+                int const options_table,
+                calcfunc const calcres,
+                calcfunc const applybnds,
+                void * const userdata,
+                options_t const & options,
+                CCTK_REAL const minerror);
+  
+  
+  
+  void
+  smooth (cGH const * restrict const cctkGH,
+          vector<CCTK_INT> const & var,
+          vector<CCTK_INT> const & res,
+          vector<CCTK_INT> const & rhs,
+          vector<CCTK_INT> const & wgt,
+          options_t const & options,
+          CCTK_REAL const old_error,
+          CCTK_REAL & error);
+  
+  void
+  norm (cGH const * restrict const cctkGH,
+        vector<CCTK_INT> const & res,
+        vector<CCTK_INT> const & rhs,
+        options_t const & options,
+        CCTK_REAL & error);
+  
+  
+  
+  void
+  residual (cGH const * restrict const cctkGH,
+            int const options_table,
+            calcfunc const calcres,
+            void * const userdata)
+    throw (char const *);
+  
+  void
+  boundary (cGH const * restrict const cctkGH,
+            int const options_table,
+            calcfunc const applybnds,
+            void * const userdata)
+    throw (char const *);
+  
+  
+  
+  void
+  restrict_var (cGH const * restrict const cctkGH,
+                vector<CCTK_INT> const & var,
+                options_t const & options);
+  
+  void
+  prolongate_var (cGH const * restrict const cctkGH,
+                  vector<CCTK_INT> const & var,
+                  options_t const & options);
+  
+  
+  
+  void
+  zero (cGH const * restrict const cctkGH,
+        vector<CCTK_INT> const & dst,
+        options_t const & options);
+  
+  void
+  copy (cGH const * restrict const cctkGH,
+        vector<CCTK_INT> const & dst,
+        vector<CCTK_INT> const & src,
+        options_t const & options);
+  
+  void
+  add (cGH const * restrict const cctkGH,
+       vector<CCTK_INT> const & dst,
+       vector<CCTK_INT> const & src,
+       options_t const & options);
+  
+  void
+  subtract (cGH const * restrict const cctkGH,
+            vector<CCTK_INT> const & dst,
+            vector<CCTK_INT> const & src,
+            options_t const & options);
+  
+    
+    
+  void
+  getvars (int const options_table,
+           char const * restrict const name,
+           vector<CCTK_INT> & vars);
+  
+  void
+  checkvar (int const vi,
+            vector<int> & usedvars);
+  
+  void
+  interior_shape (cGH const * restrict const cctkGH,
+                  options_t const & options,
+                  int * restrict const imin,
+                  int * restrict const imax);
+  
+  
+  
+  int
+  indwidth ();
+  
+  
+  
+  // Register this solver
+  int
+  CarpetMG_register ()
+  {
+    int const ierr = TATelliptic_RegisterSolver (solve, solver_name);
+    assert (! ierr);
+    
+    return 0;
+  }
+  
+  
+  
+  // Check parameters
+  void
+  CarpetMG_paramcheck (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    
+    if (CCTK_EQUALS (direct_solver, solver_name)) {
+      CCTK_VParamWarn (CCTK_THORNSTRING,
+                       "It is not possible to use \"%s\" as direct solver for \"%s\" -- this would lead to an infinite recursion",
+                       solver_name, solver_name);
+    }
+  }
+  
+  
+  
+  // Solve the system
+  int
+  solve (cGH const * restrict const cctkGH,
+         int const * restrict const var_,
+         int const * restrict const res_,
+         int const nvar,
+         int const options_table,
+         calcfunc const calcres,
+         calcfunc const applybnds,
+         void * const userdata)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    // Check arguments
+    assert (cctkGH);
+    assert (var_);
+    assert (res_);
+    assert (nvar >= 0);
+    assert (calcres);
+    assert (applybnds);
+    
+    
+    
+    if (verbose) {
+      CCTK_VInfo (CCTK_THORNSTRING,
+                  "%*s[%d] beginning to solve",
+                  indwidth(), "", reflevel);
+    }
+    
+    
+    
+    // Solve the equation on this and some coarser levels, leaving all
+    // finer levels untouched.
+    assert (is_level_mode());
+    
+    
+    
+    // Copy solution variables from arguments
+    vector<CCTK_INT> var (nvar);
+    for (int n=0; n<nvar; ++n) {
+      var.at(n) = var_[n];
+    }
+    
+    // Copy residual variables from arguments
+    vector<CCTK_INT> res (nvar);
+    for (int n=0; n<nvar; ++n) {
+      res.at(n) = res_[n];
+    }
+    
+    // Get RHS variables from options
+    vector<CCTK_INT> rhs (nvar);
+    {
+      int const icnt
+        = Util_TableGetIntArray (options_table, nvar, &rhs.front(), "rhs");
+      assert (icnt == nvar);
+    }
+    
+    // Get save variables from options
+    vector<CCTK_INT> sav (nvar);
+    {
+      int const icnt
+        = Util_TableGetIntArray (options_table, nvar, &sav.front(), "sav");
+      assert (icnt == nvar);
+    }
+    
+    // Get weight variables from options
+    vector<CCTK_INT> wgt;
+    getvars (options_table, "wgt", wgt);
+    int const nwgt = wgt.size();
+    assert (nwgt == 0 or nwgt == nvar);
+    
+    // Get auxiliary variables from options
+    vector<CCTK_INT> aux;
+    getvars (options_table, "aux", aux);
+    int const naux = aux.size();
+    assert (naux >= 0);
+    
+    
+    
+    // Check variables
+    {
+      vector<int> usedvars;
+      
+      for (int n=0; n<nvar; ++n) {
+        checkvar (var.at(n), usedvars);
+        checkvar (res.at(n), usedvars);
+        checkvar (rhs.at(n), usedvars);
+        checkvar (sav.at(n), usedvars);
+      }
+      
+      for (int n=0; n<nwgt; ++n) {
+        checkvar (wgt.at(n), usedvars);
+      }
+      
+      for (int n=0; n<naux; ++n) {
+        checkvar (aux.at(n), usedvars);
+      }
+    }
+    
+    
+    
+    options_t options;
+    
+    // Get outer boundary information
+    {
+      options.bndwidth = 1;
+      int const icnt
+        = Util_TableGetInt (options_table, &options.bndwidth, "bndwidth");
+      assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    }
+    
+    // Get symmetry information
+    {
+      int icnt;
+      int const symtable = SymmetryTableHandleForGrid (cctkGH);
+      assert (symtable >= 0);
+      
+      icnt = (Util_TableGetIntArray
+              (symtable, 6,
+               options.symmetry_handle, "symmetry_handle"));
+      assert (icnt == 6);
+      
+      icnt = (Util_TableGetIntArray
+              (symtable, 6,
+               options.symmetry_zone_width, "symmetry_zone_width"));
+      assert (icnt == 6);
+    }
+    
+    // Get full multigrid information
+    {
+      options.firstlevel = 0;
+      int const icnt
+        = Util_TableGetInt (options_table, &options.firstlevel, "firstlevel");
+      assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    }
+    
+    // Get direct solver information
+    {
+      options.maxiters = 100;
+      int const icnt
+        = Util_TableGetInt (options_table, &options.maxiters, "maxiters");
+      assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    }
+    
+    // Get equation information
+    {
+      int icnt;
+      
+      options.factor = 0.5;
+      icnt = Util_TableGetReal (options_table, &options.factor, "factor");
+      assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+      
+      options.factors.resize (var.size());
+      for (int n=0; n<options.factors.size(); ++n) {
+        options.factors.at(n) = 1.0;
+      }
+      icnt = (Util_TableGetRealArray
+              (options_table, options.factors.size(),
+               &options.factors.front(), "factors"));
+      assert (icnt == options.factors.size()
+              or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    }
+    
+    // Get multigrid information
+    {
+      int icnt;
+      
+      options.mgiters = 100;
+      icnt = Util_TableGetInt (options_table, &options.mgiters, "mgiters");
+      assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+      
+      options.presteps = 2;
+      icnt = Util_TableGetInt (options_table, &options.presteps, "presteps");
+      assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+      
+      options.poststeps = 2;
+      icnt = Util_TableGetInt (options_table, &options.poststeps, "poststeps");
+      assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    }
+    
+    
+    
+    CCTK_REAL minerror = 1.0e-8;
+    {
+      int const icnt
+        = Util_TableGetReal (options_table, &minerror, "minerror");
+      assert (icnt == 1 or icnt == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    }
+    
+    
+    
+    if (reflevel < options.firstlevel) return +1;
+    
+    
+    
+    // Initialise RHS to zero
+    BEGIN_GLOBAL_MODE (cctkGH) {
+      BEGIN_REFLEVEL_LOOP (cctkGH) {
+        zero (cctkGH, rhs, options);
+      } END_REFLEVEL_LOOP;
+    } END_GLOBAL_MODE;
+    
+    
+    
+    CCTK_REAL error;
+    int const ierr
+      = multigrid (cctkGH,
+                   var, res, rhs, sav, wgt, aux,
+                   options_table,
+                   calcres, applybnds, userdata,
+                   options,
+                   minerror);
+    
+    
+    
+    if (verbose) {
+      CCTK_VInfo (CCTK_THORNSTRING,
+                  "%*s[%d] finished solving",
+                  indwidth(), "", reflevel);
+    }
+    
+    return ierr;
+  }
+  
+  
+  
+  // Perform recursive multigrid cycles until converged
+  // Return values:
+  // positive: did not converge, continue
+  // zero:     did converge, exit
+  // negative: error, abort
+  int
+  multigrid (cGH const * restrict const cctkGH,
+             vector<CCTK_INT> const & var,
+             vector<CCTK_INT> const & res,
+             vector<CCTK_INT> const & rhs,
+             vector<CCTK_INT> const & sav,
+             vector<CCTK_INT> const & wgt,
+             vector<CCTK_INT> const & aux,
+             int const options_table,
+             calcfunc const calcres,
+             calcfunc const applybnds,
+             void * const userdata,
+             options_t const & options,
+             CCTK_REAL const minerror)
+  {
+#if 0
+    DECLARE_CCTK_ARGUMENTS;
+#endif
+    DECLARE_CCTK_PARAMETERS;
+    
+    // Solve on this and some coarser levels
+    assert (is_level_mode());
+    
+#if 0
+    // Calculate the domain size
+    int minsize = INT_MAX;
+    int size = 1;
+    int maxsize = INT_MAX;
+    for (d=0; d<dim; ++d) {
+      int const gsize = cctk_gsh[d];
+      assert (gsize >= 0);
+      minsize = min (minsize, gsize);
+      
+      int const lsize = (cctk_lsh[d]
+                         - (cctk_bbox[2*d  ] ? 0 : cctk_nghostzones[d])
+                         - (cctk_bbox[2*d+1] ? 0 : cctk_nghostzones[d]));
+      assert (lsize >= 0);
+      assert (lsize <= maxsize);
+      size *= lsize;
+      maxsize /= lsize;
+    }
+#endif
+    
+    // Solve directly when on the coarsest level
+    if (reflevel == 0) {
+      return direct_solve (cctkGH,
+                           var, res, rhs, wgt, aux,
+                           options_table,
+                           calcres, applybnds, userdata,
+                           options,
+                           minerror);
+    }
+    
+    
+    
+    try {
+      
+      if (verbose) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "%*s[%d] beginning multigrid solve",
+                    indwidth(), "", reflevel);
+      }
+      
+      // Loop until converged
+      CCTK_REAL error = HUGE_VAL;
+      int iter = 0;
+      while (error > minerror) {
+        
+        if (iter >= options.mgiters) {
+          // Did not converge
+          if (verbose) {
+            CCTK_VInfo (CCTK_THORNSTRING,
+                        "%*s[%d] finished multigrid solve after %d iterations, residual %g",
+                        indwidth(), "", reflevel, iter, double(error));
+          }
+          return 1;
+        }
+        ++iter;
+        
+        
+        
+        // Presmooth
+        {
+          int step = 0;
+          while (error > minerror) {
+            
+            ++step;
+            if (step > options.presteps) break;
+            
+            residual (cctkGH, options_table, calcres, userdata);
+            CCTK_REAL const old_error = error;
+            smooth (cctkGH, var, res, rhs, wgt, options, old_error, error);
+            boundary (cctkGH, options_table, applybnds, userdata);
+            
+            if (veryverbose) {
+              CCTK_VInfo (CCTK_THORNSTRING,
+                          "%*s[%d] iteration %d, presmoothing step %d, residual %g",
+                          indwidth(), "", reflevel,
+                          iter, step, double(error));
+            }
+            
+          } // while error > minerror
+        }
+        
+        
+        
+        // Restrict
+        
+        // solve L(u) = r
+        // set f := L(u) - r
+        //
+        // given:  f = L(u ) - r
+        // wanted: 0 = L(u') - r
+        // known:  R(f) - f is small
+        // assume: R(u'-u) is also small
+        // define: fc = R(f)
+        //         uc = R(u)
+        //         fc = L(uc) - rc   (defines rc)
+        // solve:  L(uc') = rc
+        // use:    R(u'-u) = uc'-uc   (defines u')
+        // yields: f' = L(u') - r
+        // where f' is smaller than f
+        //
+        // uc <- R(u)
+        // fc <- R(f)
+        // rc <- fc
+        // fc <- L(uc) - rc
+        // rc <- fc
+        // then solve L(uc) = rc
+        //
+        // the above is identical to:
+        // uc <- R(u)
+        // rc <- R(f)
+        // rc <- L(uc) - rc
+        // then solve L(uc) = rc
+        //
+        // now change the notation:
+        // varc      <- R(var)
+        // resc-rhsc <- R(res-rhs)
+        // rhsc      <- resc-rhsc
+        // resc-rhsc <- L(varc) - rhsc
+        // rhsc      <- resc-rhsc
+        // then solve L(varc) = rhsc
+        //
+        // the above is identical to:
+        // varc <- R(var)
+        // rhsc <- R(res-rhs)
+        // rhsc <- L(varc) - rhsc
+        // then solve L(varc) = rhsc
+        
+        residual (cctkGH, options_table, calcres, userdata);
+        norm (cctkGH, res, rhs, options, error);
+        subtract (cctkGH, res, rhs, options);
+        // TODO: restrict and fixup aux
+        assert (aux.size() == 0);
+        
+        int const coarse_reflevel = reflevel - 1;
+        BEGIN_GLOBAL_MODE (cctkGH) {
+          enter_level_mode (const_cast<cGH *>(cctkGH), coarse_reflevel);
+          try {
+            
+            restrict_var (cctkGH, var, options);
+            boundary (cctkGH, options_table, applybnds, userdata);
+            copy (cctkGH, sav, var, options);
+            zero (cctkGH, res, options);
+            restrict_var (cctkGH, res, options);
+            copy (cctkGH, rhs, res, options);
+            residual (cctkGH, options_table, calcres, userdata);
+            subtract (cctkGH, rhs, res, options);
+            
+            // Recurse
+            CCTK_REAL const coarse_minerror = error / ipow(reflevelfact, dim);
+            int const ierr
+              = multigrid (cctkGH,
+                           var, res, rhs, sav, wgt, aux,
+                           options_table,
+                           calcres, applybnds, userdata,
+                           options,
+                           coarse_minerror);
+            if (ierr < 0) throw "multigrid";
+            
+            // Prolongate
+            copy (cctkGH, res, var, options);
+            subtract (cctkGH, res, sav, options);
+            
+          } catch (char const *) {
+            // TODO
+            assert (0);
+          }
+          
+          leave_level_mode (const_cast<cGH *>(cctkGH));
+        } END_GLOBAL_MODE;
+            
+        zero (cctkGH, res, options);
+        prolongate_var (cctkGH, res, options);
+        add (cctkGH, var, res, options);
+        boundary (cctkGH, options_table, applybnds, userdata);
+        
+        if (veryverbose) {
+          CCTK_VInfo (CCTK_THORNSTRING,
+                      "%*s[%d] iteration %d, after recursion, residual %g",
+                      indwidth(), "", reflevel,
+                      iter, double(error));
+        }
+        
+        
+        
+        // Postsmooth
+        {
+          int step = 0;
+          while (error > minerror) {
+            
+            ++step;
+            if (step > options.poststeps) break;
+            
+            residual (cctkGH, options_table, calcres, userdata);
+            CCTK_REAL const old_error = error;
+            smooth (cctkGH, var, res, rhs, wgt, options, old_error, error);
+            boundary (cctkGH, options_table, applybnds, userdata);
+            
+            if (veryverbose) {
+              CCTK_VInfo (CCTK_THORNSTRING,
+                          "%*s[%d] iteration %d, postsmoothing step %d, residual %g",
+                          indwidth(), "", reflevel,
+                          iter, step, double(error));
+            }
+            
+          } // while error > minerror
+        }
+        
+        
+        
+      } // while error > minerror
+      
+      if (verbose) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "%*s[%d] finished multigrid solve after %d iterations, residual %g",
+                    indwidth(), "", reflevel, iter, double(error));
+      }
+      
+      // Everything went fine
+      return 0;
+      
+    } catch (char const *) {
+      
+      // There was an error
+      return -1;
+      
+    }
+  }
+  
+  
+  
+  // Solve directly
+  // This assumes that the grid covers the whole domain
+  // Return values:
+  // positive: did not converge, continue
+  // zero:     did converge, exit
+  // negative: error, abort
+  int
+  direct_solve (cGH const * restrict const cctkGH,
+                vector<CCTK_INT> const & var,
+                vector<CCTK_INT> const & res,
+                vector<CCTK_INT> const & rhs,
+                vector<CCTK_INT> const & wgt,
+                vector<CCTK_INT> const & aux,
+                int const options_table,
+                calcfunc const calcres,
+                calcfunc const applybnds,
+                void * const userdata,
+                options_t const & options,
+                CCTK_REAL const minerror)
+  {
+    DECLARE_CCTK_PARAMETERS;    
+    
+    if (verbose) {
+      CCTK_VInfo (CCTK_THORNSTRING,
+                  "%*s[%d] beginning direct solve",
+                  indwidth(), "", reflevel);
+    }
+    
+    assert (is_level_mode());
+    
+    try {
+      
+      // Loop until converged
+      CCTK_REAL error = HUGE_VAL;
+      int iter = 0;
+      while (error > minerror) {
+        
+        if (iter >= options.maxiters) {
+          // Did not converge
+          if (verbose) {
+            CCTK_VInfo (CCTK_THORNSTRING,
+                        "%*s[%d] finished direct solve after %d iterations, residual is %g",
+                        indwidth(), "", reflevel, iter, double(error));
+          }
+          return 1;
+        }
+        ++iter;
+        
+        residual (cctkGH, options_table, calcres, userdata);
+        CCTK_REAL const old_error = error;
+        smooth (cctkGH, var, res, rhs, wgt, options, old_error, error);
+        boundary (cctkGH, options_table, applybnds, userdata);
+        
+        if (veryverbose) {
+          CCTK_VInfo (CCTK_THORNSTRING,
+                      "%*s[%d] iteration %d, residual %g",
+                      indwidth(), "", reflevel,
+                      iter, double(error));
+        }
+        
+      } // while error > min_error
+      
+      if (verbose) {
+        CCTK_VInfo (CCTK_THORNSTRING,
+                    "%*s[%d] finished direct solve after %d iterations, residual is %g",
+                    indwidth(), "", reflevel, iter, double(error));
+      }
+      
+      // Everything went fine
+      return 0;
+      
+    } catch (void *) {
+      
+      // There was an error
+      return -1;
+      
+    }
+  }
+  
+  
+  
+  // Smooth
+  void
+  smooth (cGH const * restrict const cctkGH,
+          vector<CCTK_INT> const & var,
+          vector<CCTK_INT> const & res,
+          vector<CCTK_INT> const & rhs,
+          vector<CCTK_INT> const & wgt,
+          options_t const & options,
+          CCTK_REAL const old_error,
+          CCTK_REAL & error)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    
+    // Initialise errors
+    CCTK_REAL count = 0.0;
+    CCTK_REAL error2 = 0.0;
+    
+    
+    
+    // Calculate grid spacings
+    CCTK_REAL dxinv2 = 0;
+    CCTK_REAL dx[dim];
+    for (int d=0; d<dim; ++d) {
+      // TODO: correct this for solving on grid arrays instead of grid
+      // functions
+      dx[d] = CCTK_DELTA_SPACE(d);
+      dxinv2 += 1.0 / pow(dx[d], 2);
+    }
+    CCTK_REAL const mdxinv2inv = 1.0 / (-2.0 * dxinv2);
+    
+    
+    
+    // Smooth and calculate errors
+    BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+      BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+        DECLARE_CCTK_ARGUMENTS;
+        for (int n=0; n<var.size(); ++n) {
+          
+          CCTK_REAL * restrict const varptr
+            = (static_cast<CCTK_REAL *>
+               (CCTK_VarDataPtrI (cctkGH, 0, var.at(n))));
+          assert (varptr);
+          CCTK_REAL const * restrict const resptr
+            = (static_cast<CCTK_REAL const *>
+               (CCTK_VarDataPtrI (cctkGH, 0, res.at(n))));
+          assert (resptr);
+          CCTK_REAL const * restrict const rhsptr
+            = (static_cast<CCTK_REAL const *>
+               (CCTK_VarDataPtrI (cctkGH, 0, rhs.at(n))));
+          assert (rhsptr);
+          CCTK_REAL const * restrict const wgtptr
+            = (wgt.size() > 0 
+               ? (static_cast<CCTK_REAL const *>
+                  (CCTK_VarDataPtrI (cctkGH, 0, wgt.at(n))))
+               : NULL);
+          assert (wgt.size() == 0 or wgtptr);
+          
+          CCTK_REAL const fac
+            = options.factor * (wgtptr ? 1.0 : options.factors.at(n));
+          
+          int imin[3], imax[3];
+          interior_shape (cctkGH, options, imin, imax);
+          
+          for (int k=imin[2]; k<imax[2]; ++k) {
+            for (int j=imin[1]; j<imax[1]; ++j) {
+              for (int i=imin[0]; i<imax[0]; ++i) {
+                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                
+                CCTK_REAL const diff = resptr[ind] - rhsptr[ind];
+                CCTK_REAL const w = wgtptr ? 1.0 / wgtptr[ind] : mdxinv2inv;
+                
+                varptr[ind] -= fac * w * diff;
+                
+                ++ count;
+                error2 += pow(diff, 2);
+                
+              }
+            }
+          }
+        }
+      } END_LOCAL_COMPONENT_LOOP;
+    } END_MAP_LOOP;
+    
+    
+    
+    // Reduce errors
+    int const sum_handle = CCTK_ReductionArrayHandle ("sum");
+    assert (sum_handle >= 0);
+    CCTK_REAL reduce_in[2], reduce_out[2];
+    reduce_in[0] = count;
+    reduce_in[1] = error2;
+    int const ierr = CCTK_ReduceLocArrayToArray1D
+      (cctkGH, -1, sum_handle, reduce_in, reduce_out, 2, CCTK_VARIABLE_REAL);
+    count = reduce_out[0];
+    error2 = reduce_out[1];
+    if (count > 0) {
+      error = sqrt(error2 / count);
+    } else {
+      error = 0.0;
+    }
+    
+    
+    
+    // Sanity check
+    if (error > old_error) {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Residual increased during smoothing at level %d from %g to %g",
+                  reflevel, double(old_error), double(error));
+    }
+  }
+  
+  
+  
+  // Calculate the norm of the residual without smoothing
+  void
+  norm (cGH const * restrict const cctkGH,
+        vector<CCTK_INT> const & res,
+        vector<CCTK_INT> const & rhs,
+        options_t const & options,
+        CCTK_REAL & error)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    
+    // Initialise errors
+    CCTK_REAL count = 0.0;
+    CCTK_REAL error2 = 0.0;
+    
+    
+    
+    // Calculate errors
+    BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
+      BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+        DECLARE_CCTK_ARGUMENTS;
+        for (int n=0; n<res.size(); ++n) {
+          
+          CCTK_REAL const * restrict const resptr
+            = (static_cast<CCTK_REAL const *>
+               (CCTK_VarDataPtrI (cctkGH, 0, res.at(n))));
+          assert (resptr);
+          CCTK_REAL const * restrict const rhsptr
+            = (static_cast<CCTK_REAL const *>
+               (CCTK_VarDataPtrI (cctkGH, 0, rhs.at(n))));
+          assert (rhsptr);
+          
+          int imin[3], imax[3];
+          interior_shape (cctkGH, options, imin, imax);
+          
+          for (int k=imin[2]; k<imax[2]; ++k) {
+            for (int j=imin[1]; j<imax[1]; ++j) {
+              for (int i=imin[0]; i<imax[0]; ++i) {
+                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                
+                CCTK_REAL const diff = resptr[ind] - rhsptr[ind];
+                
+                ++ count;
+                error2 += pow(diff, 2);
+                
+              }
+            }
+          }
+        }
+      } END_LOCAL_COMPONENT_LOOP;
+    } END_MAP_LOOP;
+    
+    
+    
+    // Reduce errors
+    int const sum_handle = CCTK_ReductionArrayHandle ("sum");
+    assert (sum_handle >= 0);
+    CCTK_REAL reduce_in[2], reduce_out[2];
+    reduce_in[0] = count;
+    reduce_in[1] = error2;
+    int const ierr = CCTK_ReduceLocArrayToArray1D
+      (cctkGH, -1, sum_handle, reduce_in, reduce_out, 2, CCTK_VARIABLE_REAL);
+    count = reduce_out[0];
+    error2 = reduce_out[1];
+    if (count > 0) {
+      error = sqrt(error2 / count);
+    } else {
+      error = 0.0;
+    }
+  }
+  
+  
+  
+  // Calculate the residual
+  void
+  residual (cGH const * restrict const cctkGH,
+            int const options_table,
+            calcfunc const calcres,
+            void * const userdata)
+    throw (char const *)
+  {
+    int ierr = 0;
+    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;
+    if (ierr) throw "calcres";
+  }
+  
+  
+  
+  // Apply the boundary condition
+  void
+  boundary (cGH const * restrict const cctkGH,
+            int const options_table,
+            calcfunc const applybnds,
+            void * const userdata)
+    throw (char const *)
+  {
+    int ierr = 0;
+    ierr = applybnds (cctkGH, options_table, userdata);
+    if (ierr) throw "applybnds";
+    ierr = CallScheduleGroup (const_cast<cGH *>(cctkGH), "ApplyBCs");
+    assert (! ierr);
+  }
+  
+  
+  
+  // Restrict the solution variables from the next finer level
+  void
+  restrict_var (cGH const * restrict const cctkGH,
+                vector<CCTK_INT> const & var,
+                options_t const & options)
+  {
+    assert (reflevel < reflevels - 1);
+    int const tl = 0;
+    for (comm_state<dim> state; !state.done(); state.step()) {
+      for (int m=0; m<maps; ++m) {
+        const CCTK_REAL time = vtt.at(m)->time (tl, reflevel, mglevel);
+        for (int c=0; c<vhh.at(m)->components(reflevel); ++c) {
+          for (int n=0; n<var.size(); ++n) {
+            int const vi = var.at(n);
+            assert (vi >= 0);
+            int const gi = CCTK_GroupIndexFromVarI (vi);
+            assert (gi >= 0);
+            int const v0 = CCTK_FirstVarIndexI (gi);
+            assert (v0 >= 0);
+            arrdata.at(gi).at(m).data.at(vi-v0)->ref_restrict
+              (state, tl, reflevel, c, mglevel, time);
+          } // for n
+        } // for c
+      } // for m
+    } // for state
+  }
+  
+  
+  
+  // Prolongate the solution variables from the next coarser level
+  void
+  prolongate_var (cGH const * restrict const cctkGH,
+                  vector<CCTK_INT> const & var,
+                  options_t const & options)
+  {
+    assert (reflevel > 0);
+    int const tl = 0;
+    for (comm_state<dim> state; !state.done(); state.step()) {
+      for (int m=0; m<maps; ++m) {
+        const CCTK_REAL time = vtt.at(m)->time (tl, reflevel, mglevel);
+        for (int c=0; c<vhh.at(m)->components(reflevel); ++c) {
+          for (int n=0; n<var.size(); ++n) {
+            int const vi = var.at(n);
+            assert (vi >= 0);
+            int const gi = CCTK_GroupIndexFromVarI (vi);
+            assert (gi >= 0);
+            int const v0 = CCTK_FirstVarIndexI (gi);
+            assert (v0 >= 0);
+            arrdata.at(gi).at(m).data.at(vi-v0)->ref_prolongate
+              (state, tl, reflevel, c, mglevel, time);
+          } // for n
+        } // for c
+      } // for m
+    } // for state
+  }
+  
+  
+  
+  void
+  zero (cGH const * restrict const cctkGH,
+        vector<CCTK_INT> const & dst,
+        options_t const & options)
+  {
+    BEGIN_MAP_LOOP (cctkGH, CCTK_GF) {
+      BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
+        DECLARE_CCTK_ARGUMENTS;
+        for (int n=0; n<dst.size(); ++n) {
+          CCTK_REAL * restrict const dstptr
+            = (static_cast<CCTK_REAL *>
+               (CCTK_VarDataPtrI (cctkGH, 0, dst.at(n))));
+          assert (dstptr);
+          for (int k=0; k<cctk_lsh[2]; ++k) {
+            for (int j=0; j<cctk_lsh[1]; ++j) {
+              for (int i=0; i<cctk_lsh[0]; ++i) {
+                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                dstptr[ind] = 0.0;
+              }
+            }
+          }
+        }
+      } END_LOCAL_COMPONENT_LOOP;
+    } END_MAP_LOOP;
+  }
+  
+  
+  
+  void
+  copy (cGH const * restrict const cctkGH,
+        vector<CCTK_INT> const & dst,
+        vector<CCTK_INT> const & src,
+        options_t const & options)
+  {
+    assert (dst.size() == src.size());
+    BEGIN_MAP_LOOP (cctkGH, CCTK_GF) {
+      BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
+        DECLARE_CCTK_ARGUMENTS;
+        for (int n=0; n<dst.size(); ++n) {
+          CCTK_REAL * restrict const dstptr
+            = (static_cast<CCTK_REAL *>
+               (CCTK_VarDataPtrI (cctkGH, 0, dst.at(n))));
+          assert (dstptr);
+          CCTK_REAL const * restrict const srcptr
+            = (static_cast<CCTK_REAL *>
+               (CCTK_VarDataPtrI (cctkGH, 0, src.at(n))));
+          assert (srcptr);
+          for (int k=0; k<cctk_lsh[2]; ++k) {
+            for (int j=0; j<cctk_lsh[1]; ++j) {
+              for (int i=0; i<cctk_lsh[0]; ++i) {
+                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                dstptr[ind] = srcptr[ind];
+              }
+            }
+          }
+        }
+      } END_LOCAL_COMPONENT_LOOP;
+    } END_MAP_LOOP;
+  }
+  
+  
+  
+  void
+  add (cGH const * restrict const cctkGH,
+       vector<CCTK_INT> const & dst,
+       vector<CCTK_INT> const & src,
+       options_t const & options)
+  {
+    assert (dst.size() == src.size());
+    BEGIN_MAP_LOOP (cctkGH, CCTK_GF) {
+      BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
+        DECLARE_CCTK_ARGUMENTS;
+        for (int n=0; n<dst.size(); ++n) {
+          CCTK_REAL * restrict const dstptr
+            = (static_cast<CCTK_REAL *>
+               (CCTK_VarDataPtrI (cctkGH, 0, dst.at(n))));
+          assert (dstptr);
+          CCTK_REAL const * restrict const srcptr
+            = (static_cast<CCTK_REAL *>
+               (CCTK_VarDataPtrI (cctkGH, 0, src.at(n))));
+          assert (srcptr);
+          for (int k=0; k<cctk_lsh[2]; ++k) {
+            for (int j=0; j<cctk_lsh[1]; ++j) {
+              for (int i=0; i<cctk_lsh[0]; ++i) {
+                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                dstptr[ind] += srcptr[ind];
+              }
+            }
+          }
+        }
+      } END_LOCAL_COMPONENT_LOOP;
+    } END_MAP_LOOP;
+  }
+  
+  
+  
+  void
+  subtract (cGH const * restrict const cctkGH,
+            vector<CCTK_INT> const & dst,
+            vector<CCTK_INT> const & src,
+            options_t const & options)
+  {
+    assert (dst.size() == src.size());
+    BEGIN_MAP_LOOP (cctkGH, CCTK_GF) {
+      BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
+        DECLARE_CCTK_ARGUMENTS;
+        for (int n=0; n<dst.size(); ++n) {
+          CCTK_REAL * restrict const dstptr
+            = (static_cast<CCTK_REAL *>
+               (CCTK_VarDataPtrI (cctkGH, 0, dst.at(n))));
+          assert (dstptr);
+          CCTK_REAL const * restrict const srcptr
+            = (static_cast<CCTK_REAL *>
+               (CCTK_VarDataPtrI (cctkGH, 0, src.at(n))));
+          assert (srcptr);
+          for (int k=0; k<cctk_lsh[2]; ++k) {
+            for (int j=0; j<cctk_lsh[1]; ++j) {
+              for (int i=0; i<cctk_lsh[0]; ++i) {
+                int const ind = CCTK_GFINDEX3D (cctkGH, i, j, k);
+                dstptr[ind] -= srcptr[ind];
+              }
+            }
+          }
+        }
+      } END_LOCAL_COMPONENT_LOOP;
+    } END_MAP_LOOP;
+  }
+  
+  
+  
+  // Get a variable list from an options table
+  void
+  getvars (int const options_table,
+           char const * restrict const name,
+           vector<CCTK_INT> & vars)
+  {
+    int const nvars = Util_TableGetIntArray (options_table, 0, NULL, name);
+    if (nvars == UTIL_ERROR_TABLE_NO_SUCH_KEY) {
+      vars.resize (0);
+    } else if (nvars >= 0) {
+      assert (nvars >= 0);
+      vars.resize (nvars);
+      int const icnt
+        = Util_TableGetIntArray (options_table, nvars, &vars.front(), name);
+      assert (icnt == nvars);
+    } else {
+      assert (0);
+    }
+  }
+  
+  
+  
+  // Check a grid variable index
+  void
+  checkvar (int const vi,
+            vector<int> & usedvars)
+  {
+    int ierr;
+    
+    assert (vi>=0 and vi<CCTK_NumVars());
+    
+    int const gi = CCTK_GroupIndexFromVarI (vi);
+    assert (gi>=0);
+    
+    cGroup group;
+    ierr = CCTK_GroupData (gi, &group);
+    assert (! ierr);
+    
+    assert (group.grouptype == CCTK_GF);
+    assert (group.vartype == CCTK_VARIABLE_REAL);
+    assert (group.dim == 3);
+    
+    for (int n=0; n<usedvars.size(); ++n) {
+      assert (vi != usedvars.at(n));
+    }
+    usedvars.push_back (vi);
+  }
+  
+  
+  
+  // Calculate the shape of the interior
+  void
+  interior_shape (cGH const * restrict const cctkGH,
+                  options_t const & options,
+                  int * restrict const imin,
+                  int * restrict const imax)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    
+    for (int d=0; d<dim; ++d) {
+      
+      imin[d] = 0;
+      if (cctk_bbox[2*d]) {
+        if (options.symmetry_handle[2*d] >= 0) {
+          imin[d] += options.symmetry_zone_width[2*d];
+        } else {
+          imin[d] += options.bndwidth;
+        }
+      } else {
+        imin[d] += cctk_nghostzones[d];
+      }
+      
+      imax[d] = cctk_lsh[d];
+      if (cctk_bbox[2*d+1]) {
+        if (options.symmetry_handle[2*d+1] >= 0) {
+          imax[d] -= options.symmetry_zone_width[2*d+1];
+        } else {
+          imax[d] -= options.bndwidth;
+        }
+      } else {
+        imax[d] -= cctk_nghostzones[d];
+      }
+      
+      assert (imin[d] <= imax[d]);
+      
+    } // for d
+  }
+  
+  
+  
+  // Determine indentation
+  int
+  indwidth ()
+  {
+    if (reflevel == -1) return 0;
+    return 3 * (reflevels - reflevel - 1);
+  }
+  
+} // namespace CarpetMG