Initial revision

darcs-hash:20010301114010-f6438-12fb8a9ffcc80e86c0a97e37b5b0dae0dbc59b79.gz

9 files changed, 868 insertions, 0 deletions

diff --git a/CarpetDev/CarpetJacobi/README b/CarpetDev/CarpetJacobi/README
new file mode 100644
index 000000000..74967b13f
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/README
@@ -0,0 +1,11 @@
+CVS info   : $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/README,v 1.1 2003/09/02 14:35:57 schnetter Exp $
+
+Cactus Code Thorn CarpetJacobi
+Thorn Author(s)     : Erik Schnetter <schnetter@uni-tuebingen.de>
+Thorn Maintainer(s) : Erik Schnetter <schnetter@uni-tuebingen.de>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+Solve an elliptic equation using Jacobi iterations with varying step
+size.
diff --git a/CarpetDev/CarpetJacobi/doc/documentation.tex b/CarpetDev/CarpetJacobi/doc/documentation.tex
new file mode 100644
index 000000000..11d9ffd8d
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/doc/documentation.tex,v 1.1 2003/09/02 14:35:57 schnetter 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 (later) 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: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/doc/documentation.tex,v 1.1 2003/09/02 14:35:57 schnetter 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@uni-tuebingen.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{}
+
+% the date your document was last changed, if your document is in CVS, 
+% please use:
+%    \date{$ $Date: 2003/09/02 14:35:57 $ $}
+\date{August 31 2003}
+
+\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}
+
+\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/CarpetJacobi/interface.ccl b/CarpetDev/CarpetJacobi/interface.ccl
new file mode 100644
index 000000000..6b61da31a
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/interface.ccl
@@ -0,0 +1,9 @@
+# Interface definition for thorn CarpetJacobi
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/interface.ccl,v 1.1 2003/09/02 14:35:57 schnetter Exp $
+
+IMPLEMENTS: CarpetJacobi
+
+INHERITS: TATelliptic
+
+USES INCLUDE HEADER: carpet.hh
+USES INCLUDE HEADER: TATelliptic.h
diff --git a/CarpetDev/CarpetJacobi/par/charge_te_cjac.par b/CarpetDev/CarpetJacobi/par/charge_te_cjac.par
new file mode 100644
index 000000000..8db9c4cf6
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/par/charge_te_cjac.par
@@ -0,0 +1,60 @@
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/par/charge_te_cjac.par,v 1.4 2004/03/23 11:59:20 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"
+
+!DESC "Charged sphere initial data, solved with TATelliptic/CarpetJacobi"
+
+cactus::cctk_full_warnings      = yes
+cactus::cctk_timer_output       = full
+
+#Cactus::cctk_itlast             = 10
+Cactus::cctk_itlast             = 0
+ 
+Time::dtfac                     = 0.5
+
+driver::global_nsize            = 21
+#driver::global_nsize           = 41
+
+Carpet::verbose                 = yes
+#Carpet::veryverbose             = yes
+
+Carpet::max_refinement_levels   = 2
+#Carpet::prolongation_order_space = 3
+Carpet::prolongation_order_time = 2
+Carpet::init_3_timelevels       = yes
+
+CarpetRegrid::refinement_levels = 1
+CarpetRegrid::refined_regions   = manual-gridpoint-list
+CarpetRegrid::gridpoints        = "[[ ([1,1,1]:[20,20,20]:[1,1,1]) ]]"
+CarpetRegrid::outerbounds       = "[[ [[1,0],[1,0],[1,0]] ]]"
+
+grid::type                      = byspacing
+grid::domain                    = octant
+grid::dxyz                      = 1.2
+grid::avoid_origin              = no
+
+IDScalarWave::initial_data      = charge-TATelliptic-simple
+IDSWTEsimple::solver            = CarpetJacobi
+#IDSWTEsimple::options           = "maxiters=10000 minerror=1e-4 factor=1e-2"
+#IDSWTEsimple::options           = "maxiters=100000 minerror=1e-8 factor=1e-2"
+IDSWTEsimple::options           = "maxiters=10000 minerror=1e-8"
+IDSWTEsimple::radius            = 5.5
+IDSWTEsimple::charge            = 1.0
+
+WaveToyC::bound                 = radiation
+
+CarpetJacobi::verbose           = yes
+
+IO::out_dir                     = "charge_te_cjac"
+IO::out_fileinfo                = none
+
+IOBasic::outInfo_every          = 1
+IOBasic::outInfo_vars           = "wavetoy::phi IDSWTEsimple::residual"
+
+IOBasic::outScalar_every        = 1
+IOBasic::outScalar_vars         = "wavetoy::phi IDSWTEsimple::residual"
+IOBasic::outScalar_style        = "gnuplot"
+
+CarpetIOASCII::out1D_every      = 1
+CarpetIOASCII::out1D_vars       = "wavetoy::phi IDSWTEsimple::residual"
diff --git a/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par b/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par
new file mode 100644
index 000000000..b26fe4d48
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/par/charge_te_cjac_2.par
@@ -0,0 +1,60 @@
+# $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 $
+# 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"
+
+!DESC "Charged sphere initial data, solved with TATelliptic/CarpetJacobi"
+
+cactus::cctk_full_warnings      = yes
+cactus::cctk_timer_output       = full
+
+#Cactus::cctk_itlast             = 10
+Cactus::cctk_itlast             = 0
+ 
+Time::dtfac                     = 0.5
+
+driver::global_nsize            = 21
+#driver::global_nsize           = 41
+
+Carpet::verbose                 = yes
+#Carpet::veryverbose             = yes
+
+Carpet::max_refinement_levels   = 2
+#Carpet::prolongation_order_space = 3
+Carpet::prolongation_order_time = 2
+Carpet::init_3_timelevels       = yes
+
+CarpetRegrid::refinement_levels = 2
+CarpetRegrid::refined_regions   = manual-gridpoint-list
+CarpetRegrid::gridpoints        = "[[ ([1,1,1]:[20,20,20]:[1,1,1]) ]]"
+CarpetRegrid::outerbounds       = "[[ [[1,0],[1,0],[1,0]] ]]"
+
+grid::type                      = byspacing
+grid::domain                    = octant
+grid::dxyz                      = 1.2
+grid::avoid_origin              = no
+
+IDScalarWave::initial_data      = charge-TATelliptic-simple
+IDSWTEsimple::solver            = CarpetJacobi
+#IDSWTEsimple::options           = "maxiters=10000 minerror=1e-4 factor=1e-2"
+#IDSWTEsimple::options           = "maxiters=100000 minerror=1e-8 factor=1e-2"
+IDSWTEsimple::options           = "maxiters=10000 minerror=1e-8"
+IDSWTEsimple::radius            = 5.5
+IDSWTEsimple::charge            = 1.0
+
+WaveToyC::bound                 = radiation
+
+CarpetJacobi::verbose           = yes
+
+IO::out_dir                     = "charge_te_cjac_2"
+IO::out_fileinfo                = none
+
+IOBasic::outInfo_every          = 1
+IOBasic::outInfo_vars           = "wavetoy::phi IDSWTEsimple::residual"
+
+IOBasic::outScalar_every        = 1
+IOBasic::outScalar_vars         = "wavetoy::phi IDSWTEsimple::residual"
+IOBasic::outScalar_style        = "gnuplot"
+
+CarpetIOASCII::out1D_every      = 1
+CarpetIOASCII::out1D_vars       = "wavetoy::phi IDSWTEsimple::residual"
diff --git a/CarpetDev/CarpetJacobi/param.ccl b/CarpetDev/CarpetJacobi/param.ccl
new file mode 100644
index 000000000..bdeecf41b
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/param.ccl
@@ -0,0 +1,32 @@
+# 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 $
+
+BOOLEAN verbose "Produce log output while running" STEERABLE=always
+{
+} "no"
+
+BOOLEAN veryverbose "Produce much log output while running" STEERABLE=always
+{
+} "no"
+
+CCTK_INT outevery "Produce output every n iterations" STEERABLE=always
+{
+  1:* :: ""
+} 1
+
+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/schedule.ccl b/CarpetDev/CarpetJacobi/schedule.ccl
new file mode 100644
index 000000000..009b40f7a
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/schedule.ccl
@@ -0,0 +1,7 @@
+# Schedule definitions for thorn CarpetJacobi
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/schedule.ccl,v 1.1 2003/09/02 14:35:57 schnetter Exp $
+
+SCHEDULE CarpetJacobi_register AT startup
+{
+  LANG: C
+} "Register the elliptic solver"
diff --git a/CarpetDev/CarpetJacobi/src/Jacobi.cc b/CarpetDev/CarpetJacobi/src/Jacobi.cc
new file mode 100644
index 000000000..c326b66cd
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/src/Jacobi.cc
@@ -0,0 +1,536 @@
+// $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>
+#include <cstdio>
+#include <ctime>
+#include <iostream>
+#include <vector>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "util_ErrorCodes.h"
+#include "util_Table.h"
+
+#include "TATelliptic.h"
+
+#include "carpet.hh"
+
+#include "th.hh"
+
+
+
+namespace Carpet {
+  // TODO: fix this
+  void CycleTimeLevels (const cGH* cctkGH);
+  void Restrict (const cGH* cctkGH);
+};
+
+
+
+namespace CarpetJacobi {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  extern "C" {
+    void CarpetJacobi_register ();
+  }
+  
+  int CarpetJacobi_solve (const cGH * const cctkGH,
+                          const int * const var,
+                          const int * const res,
+                          const int nvars,
+                          const int options_table,
+                          const calcfunc calcres,
+                          const calcfunc applybnds,
+                          void * const userdata);
+  
+  
+  
+  void CarpetJacobi_register ()
+  {
+    int const ierr = TATelliptic_RegisterSolver
+      (CarpetJacobi_solve, "CarpetJacobi");
+    assert (!ierr);
+  }
+  
+  
+  
+  int CarpetJacobi_solve (const cGH * const cctkGH,
+                          const int * const var,
+                          const int * const res,
+                          const int nvars,
+                          const int options_table,
+                          const calcfunc calcres,
+                          const calcfunc applybnds,
+                          void * const userdata)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    // Error control
+    int ierr;
+    
+    if (! CCTK_IsThornActive(CCTK_THORNSTRING)) {
+      CCTK_WARN (0, "Thorn " CCTK_THORNSTRING " has not been activated.  It is therefore not possible to call CarpetJacobi_solve.");
+    }
+    
+    // Check arguments
+    assert (cctkGH);
+    
+    assert (nvars > 0);
+    assert (var);
+    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_VarTypeI(var[n]) == CCTK_VARIABLE_REAL);
+      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_VarTypeI(res[n]) == CCTK_VARIABLE_REAL);
+      assert (CCTK_QueryGroupStorageI(cctkGH,
+                                      CCTK_GroupIndexFromVarI(res[n])));
+    }
+    for (int n=0; n<nvars; ++n) {
+      assert (var[n] != res[n]);
+      for (int nn=0; nn<n; ++nn) {
+        assert (var[nn] != var[n]);
+        assert (var[nn] != res[n]);
+        assert (res[nn] != var[n]);
+        assert (res[nn] != res[n]);
+      }
+    }
+    
+    
+    
+    assert (options_table >= 0);
+    
+    assert (calcres);
+    assert (applybnds);
+    
+    // Get options from options table
+    CCTK_INT maxiters = 10000;
+    CCTK_REAL minerror = 1e-8;
+    CCTK_REAL factor = 1.0e-2;    // slow start
+    CCTK_REAL minfactor = 1.0e-8;
+    CCTK_REAL maxfactor = 1.0;
+    CCTK_REAL acceleration = 1.2; // slow acceleration
+    CCTK_REAL deceleration = 0.1; // emergency brake
+    vector<CCTK_REAL> factors (nvars);
+    for (int n=0; n<nvars; ++n) {
+      factors[n] = 1.0;
+    }
+    
+    ierr = Util_TableGetInt (options_table, &maxiters, "maxiters");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &minerror, "minerror");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &factor, "factor");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &minfactor, "minfactor");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &maxfactor, "maxfactor");
+    assert (ierr==1 || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    ierr = Util_TableGetReal (options_table, &acceleration, "acceleration");
+    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");
+    assert (ierr==nvars || ierr == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    
+    assert (maxiters >= 0);
+    assert (minerror >= 0);
+    assert (factor > 0);
+    assert (minfactor > 0 && minfactor <= factor);
+    assert (maxfactor >= factor);
+    assert (acceleration >= 1);
+    assert (deceleration > 0 && deceleration <= 1);
+    for (int n=0; n<nvars; ++n) {
+      assert (factors[n] != 0);
+    }
+    
+    
+    
+    assert (is_global_mode() || is_level_mode());
+    
+    // The level to solve on
+    const int solve_level = is_level_mode() ? reflevel : reflevels-1;
+    
+    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;
+    }
+    
+#warning "TODO: assert that all levels are at the same time"
+    
+    // Switch to global mode
+    BEGIN_GLOBAL_MODE(cctkGH) {
+      
+      // 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));
+        }
+      } END_REFLEVEL_LOOP;
+      
+      
+      
+      // 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;
+      
+      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)
+#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 HAVE_ISNAN
+        // Exit if things go bad
+        if (isnan(norm2)) {
+          if (verbose || veryverbose) {
+            CCTK_VInfo (CCTK_THORNSTRING,
+                        "Encountered NaN.  Aborting solver.");
+          }
+          break;
+        }
+#endif
+        
+        // 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;
+        }
+        
+        // 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
+      
+      // 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: ;
+      
+      
+      
+      // 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;
+      
+    } END_GLOBAL_MODE;
+    
+    return ierr;
+  }
+
+} // namespace CarpetJacobi
diff --git a/CarpetDev/CarpetJacobi/src/make.code.defn b/CarpetDev/CarpetJacobi/src/make.code.defn
new file mode 100644
index 000000000..efd1ad4c4
--- /dev/null
+++ b/CarpetDev/CarpetJacobi/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn CarpetJacobi
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetDev/CarpetJacobi/src/make.code.defn,v 1.1 2003/09/02 14:35:58 schnetter Exp $
+
+# Source files in this directory
+SRCS = Jacobi.cc
+
+# Subdirectories containing source files
+SUBDIRS = 
+