move these files here from ../gr/

since they're not really GR files, but rather higher-level driver stuff


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@703 f88db872-0e4f-0410-b76b-b9085cfa78c5

11 files changed, 2133 insertions, 0 deletions

diff --git a/src/driver/AHFinderDirect.hh b/src/driver/AHFinderDirect.hh
new file mode 100644
index 0000000..5435701
--- /dev/null
+++ b/src/driver/AHFinderDirect.hh
@@ -0,0 +1,174 @@
+// AHFinderDirect.hh -- misc global-within-this-thorn stuff
+// $Header$
+
+//******************************************************************************
+
+//
+// this enum holds the (a) decoded  method  parameter, i.e. it specifies
+// our top-level method
+//
+enum	method
+	{
+	method__horizon_function,
+	method__Jacobian_test,
+	method__Jacobian_test_NP_only,
+	method__Newton_solve // no comma
+	};
+
+//
+// this enum holds the (a) decoded  verbose_method  parameter, i.e.
+// it specifies which (how many) informational messages we should print
+//
+enum	verbose_level
+	{
+	verbose_level__physics_highlights,
+	verbose_level__physics_details,
+	verbose_level__algorithm_highlights,
+	verbose_level__algorithm_details // no comma
+	};
+
+//
+// this enum holds the (a) decoded  Jacobian_method  parameter,
+// i.e. it specifies how we compute the (a) Jacobian matrix
+//
+enum	Jacobian_method
+	{
+	Jacobian_method__numerical_perturb,
+	Jacobian_method__symbolic_diff_with_FD_dr,
+	Jacobian_method__symbolic_diff // no comma
+	};
+
+//******************************************************************************
+
+//
+// This struct holds parameters for computing the Jacobian matrix.
+//
+struct	Jacobian_info
+	{
+	enum Jacobian_method Jacobian_method;
+	enum Jacobian_type   Jacobian_storage_method;
+	fp perturbation_amplitude;
+	};
+
+//
+// This struct holds parameters for solving the H(h) = 0 equations.
+//
+struct	solver_info
+	{
+	// stuff for Newton_solve()
+	int max_Newton_iterations;
+	bool output_h_and_H_at_each_Newton_iteration;
+	fp max_Delta_h_over_h;
+	fp H_norm_for_convergence;
+	fp Delta_h_norm_for_convergence;
+	bool final_H_update_if_exit_x_H_small;
+	};
+
+//
+// This struct holds info for I/O
+//
+struct	IO_info
+	{
+	const char* h_base_file_name;
+	const char* H_base_file_name;
+	const char* Jacobian_base_file_name;
+
+	// this is used to choose file names
+	int time_iteration;	// the Cactus time interation number
+				// (cctk_iteration)
+	};
+
+//******************************************************************************
+
+//
+// (A single copy of) this struct holds all of our information about
+// a single apparent horizon.
+//
+struct	AH_info
+	{
+	patch_system* ps_ptr;
+	Jacobian* Jac_ptr;
+
+	bool AH_found;
+	fp centroid_x, centroid_y, centroid_z;
+	fp area, mass;
+	};
+
+//
+// (A single copy of) this struct holds all of our state that's
+// persistent across Cactus scheduler calls.
+//
+struct	state
+	{
+	enum method        method;
+	enum verbose_level verbose_level;
+	int timer_handle;
+
+	struct IO_info IO_info;
+	struct Jacobian_info Jac_info;
+	struct solver_info solver_info;
+	struct cactus_grid_info cgi;
+	struct geometry_info gi;
+
+	int N_horizons;
+
+	// this vector is indexed with a 0-origin integer hn ("horizon number")
+	std::vector<AH_info *> AH_info_ptrs;
+	};
+
+//******************************************************************************
+
+//
+// prototypes for functions visible outside their source files
+//
+
+//**************************************
+
+//
+// *** routines called from the Cactus scheduler ***
+//
+
+// setup.cc
+extern "C"
+  void AHFinderDirect_setup(CCTK_ARGUMENTS);
+
+// initial_guess.cc
+extern "C"
+  void AHFinderDirect_initial_guess(CCTK_ARGUMENTS);
+
+// find_horizons.cc
+extern "C"
+  void AHFinderDirect_find_horizons(CCTK_ARGUMENTS);
+
+//**************************************
+
+// Newton.cc
+// return true for success, false for failure to converge
+bool Newton_solve(patch_system& ps,
+		  Jacobian& Jac,
+		  const struct cactus_grid_info& cgi,
+		  const struct geometry_info& gi,
+		  const struct Jacobian_info& Jacobian_info,
+		  const struct solver_info& solver_info,
+		  struct IO_info& IO_info,
+		  int hn, enum verbose_level verbose_level);
+
+// horizon_Jacobian.cc
+void horizon_Jacobian(patch_system& ps,
+		      Jacobian& Jac,
+		      const struct cactus_grid_info& cgi,
+		      const struct geometry_info& gi,
+		      const struct Jacobian_info& Jacobian_info,
+		      bool print_msg_flag);
+
+// io.cc
+void input_gridfn(patch_system& ps, int unknown_gfn,
+		  struct IO_info& IO_info, const char base_file_name[],
+		  int hn, bool print_msg_flag, int AHF_iteration = 0);
+void output_gridfn(patch_system& ps, int unknown_gfn,
+		   struct IO_info& IO_info, const char base_file_name[],
+		   int hn, bool print_msg_flag, int AHF_iteration = 0);
+void print_Jacobians(const patch_system& ps,
+		     const Jacobian* Jac_NP, const Jacobian* Jac_SD_FDdr,
+		     struct IO_info& IO_info, const char base_file_name[],
+		     int hn, bool print_msg_flag, int AHF_iteration = 0);
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
new file mode 100644
index 0000000..63f6243
--- /dev/null
+++ b/src/driver/Newton.cc
@@ -0,0 +1,203 @@
+// Newton.cc -- solve H(h) = 0 via Newton's method
+// $Id$
+//
+// Newton_solve - driver to solve H(h) = 0 via Newton's method
+//
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <vector>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../util/coords.hh"
+#include "../util/grid.hh"
+#include "../util/fd_grid.hh"
+#include "../util/patch.hh"
+#include "../util/patch_edge.hh"
+#include "../util/patch_interp.hh"
+#include "../util/ghost_zone.hh"
+#include "../util/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "AHFinderDirect.hh"
+
+//******************************************************************************
+
+//
+// This function solves H(h) = 0 for h via Newton's method.
+//
+// Arguments:
+// hn = the horizon number (used only in naming output files)
+//
+// Results:
+// This function returns a success flag: this is true if (and only if)
+// it found an h satisfying H(h) = 0 to within the error tolerances,
+// otherwise it's false.
+//
+bool Newton_solve(patch_system& ps,
+		  Jacobian& Jac,
+		  const struct cactus_grid_info& cgi,
+		  const struct geometry_info& gi,
+		  const struct Jacobian_info& Jacobian_info,
+		  const struct solver_info& solver_info,
+		  struct IO_info& IO_info,
+		  int hn, enum verbose_level verbose_level)
+{
+const bool print_more_msg_flag
+	= (verbose_level >= verbose_level__algorithm_details);
+
+	for (int iteration = 1 ;
+	     iteration <= solver_info.max_Newton_iterations ;
+	     ++iteration)
+	{
+	if (print_more_msg_flag)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "Newton iteration %d", iteration);
+
+	if (solver_info.output_h_and_H_at_each_Newton_iteration)
+	   then output_gridfn(ps, gfns::gfn__h,
+			      IO_info, IO_info.h_base_file_name,
+			      hn, print_more_msg_flag, iteration);
+
+	jtutil::norm<fp> H_norms;
+	if (horizon_function(ps,
+			     cgi, gi,
+			     true,	// yes, we want the Jacobian coeffs
+			     print_more_msg_flag,
+			     &H_norms))
+	   then return false;				// *** ERROR RETURN ***
+	if (print_more_msg_flag)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "   H rms-norm=%.2e, infinity-norm=%.2e",
+			   H_norms.rms_norm(), H_norms.infinity_norm());
+
+	if (solver_info.output_h_and_H_at_each_Newton_iteration)
+	   then output_gridfn(ps, gfns::gfn__H,
+			      IO_info, IO_info.H_base_file_name,
+			      hn, print_more_msg_flag, iteration);
+
+	if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence)
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "==> finished (||H|| < %g)",
+			   double(solver_info.H_norm_for_convergence));
+		return true;				// *** NORMAL RETURN ***
+		}
+
+	// compute the Newton step
+	horizon_Jacobian(ps, Jac,
+			 cgi, gi, Jacobian_info,
+			 print_more_msg_flag);
+	if (print_more_msg_flag)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "solving linear system for Delta_h (%d unknowns)",
+			   Jac.NN());
+	const fp rcond = Jac.solve_linear_system(gfns::gfn__H,
+						 gfns::gfn__Delta_h);
+	if	(rcond == 0.0)
+	   then {
+		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "Newton_solve: Jacobian matrix is numerically singular!");
+		return false;				// *** ERROR RETURN ***
+		}
+	if (print_more_msg_flag)
+	   then {
+		if (rcond > 0.0)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "done solving linear system (condition number %.2e)",
+				   double(rcond));
+		   else CCTK_VInfo(CCTK_THORNSTRING,
+				   "done solving linear system");
+		}
+
+	// if the Newton step is too large, scale it down
+	jtutil::norm<fp> h_norms;
+	ps.ghosted_gridfn_norms(gfns::gfn__h, h_norms);
+	const fp max_allowable_Delta_h
+		= solver_info.max_Delta_h_over_h * h_norms.mean();
+	jtutil::norm<fp> Delta_h_norms;
+	ps.gridfn_norms(gfns::gfn__Delta_h, Delta_h_norms);
+	const fp max_Delta_h = Delta_h_norms.infinity_norm();
+	const fp scale = (max_Delta_h > max_allowable_Delta_h)
+			 ? max_allowable_Delta_h / max_Delta_h
+			 : 1.0;
+
+	// take the Newton step (scaled if need be)
+		for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+		{
+		patch& p = ps.ith_patch(pn);
+
+			for (int irho = p.min_irho() ;
+			     irho <= p.max_irho() ;
+			     ++irho)
+			{
+			for (int isigma = p.min_isigma() ;
+			     isigma <= p.max_isigma() ;
+			     ++isigma)
+			{
+			p.ghosted_gridfn(gfns::gfn__h, irho,isigma)
+			   -= scale * p.gridfn(gfns::gfn__Delta_h, irho,isigma);
+			}
+			}
+		}
+	if (print_more_msg_flag)
+	   then {
+		if (scale == 1.0)
+		   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "h += Delta_h (rms-norm=%.2e, infinity-norm=%.2e)",
+				   Delta_h_norms.rms_norm(),
+				   Delta_h_norms.infinity_norm());
+		   else CCTK_VInfo(CCTK_THORNSTRING,
+			   "h += %g * Delta_h (infinity-norm clamped to %.2g)",
+				   scale,
+				   scale * Delta_h_norms.infinity_norm());
+		}
+
+	if ( scale * Delta_h_norms.infinity_norm()
+	     <= solver_info.Delta_h_norm_for_convergence )
+	   then {
+		CCTK_VInfo(CCTK_THORNSTRING,
+			   "==> finished (||Delta_h|| < %g)",
+			   double(solver_info.Delta_h_norm_for_convergence));
+		if (solver_info.final_H_update_if_exit_x_H_small)
+		   then {
+			if (print_more_msg_flag)
+			   then CCTK_VInfo(CCTK_THORNSTRING,
+					   "    doing final H(h) evaluation");
+			if (horizon_function(ps, cgi, gi, false, true))
+			   then return false;		// *** ERROR RETURN ***
+			}
+		return true;				// *** NORMAL RETURN ***
+		}
+	}
+
+
+CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+	   "Newton_solve: no convergence in %d iterations!",
+	   solver_info.max_Newton_iterations);
+if (solver_info.final_H_update_if_exit_x_H_small)
+   then {
+	if (print_more_msg_flag)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "    doing final H(h) evaluation");
+	if (horizon_function(ps, cgi, gi, false, true))
+	   then return false;				// *** ERROR RETURN ***
+	}
+return false;						// *** ERROR RETURN ***
+}
diff --git a/src/driver/ellipsoid.maple b/src/driver/ellipsoid.maple
new file mode 100644
index 0000000..349cc8d
--- /dev/null
+++ b/src/driver/ellipsoid.maple
@@ -0,0 +1,37 @@
+# ellipsoid.maple -- compute equations for offset ellipsoid setup
+# $Header$
+
+#
+# This program finds the intersection of an ellipsoid with an offset ray.
+#
+# ellipsoid has center (A,B,C), radius (a,b,c)
+# angular coordinate system has center (U,V,W)
+#
+# direction cosines wrt angular coordinate center are (alpha,beta,gamma)
+# but Maple predefines gamma = Euler's constant, so we use (xcos,ycos,zcos)
+# instead, i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
+#
+# then the equation of the ellipsoid is
+#	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*r - C)^2
+#	-----------------  +  ----------------  +  -----------------  =  1
+#	        a^2                  b^2                   c^2
+#
+# to solve this, we introduce intermediate variables
+#	AU = A - U
+#	BV = B - V
+#	CW = C - W
+#
+eqn := (xcos*r - AU)^2/a^2 + (ycos*r - BV)^2/b^2 + (zcos*r - CW)^2/c^2 = 1;
+
+read "../maple/util.mm";
+read "../maple/codegen2.mm";
+
+[solve(eqn, r)];
+map(simplify, %);
+[r_plus = %[1], r_minus = %[2]];
+solnlist := [codegen[optimize](%)];
+
+ftruncate("ellipsoid.c");
+print_name_list_dcl(temps_in_eqnlist(solnlist, [r_plus,r_minus]),
+		    "fp", "ellipsoid.c");
+codegen[C](solnlist, filename="ellipsoid.c");
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
new file mode 100644
index 0000000..e12ecf1
--- /dev/null
+++ b/src/driver/find_horizons.cc
@@ -0,0 +1,237 @@
+// find_horizons.cc -- top level driver for finding apparent horizons
+// $Header$
+//
+// <<<access to persistent data>>>
+// <<<prototypes for functions local to this file>>>
+// AHFinderDirect_find_horizons - top-level driver
+//
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <vector>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../util/coords.hh"
+#include "../util/grid.hh"
+#include "../util/fd_grid.hh"
+#include "../util/patch.hh"
+#include "../util/patch_edge.hh"
+#include "../util/patch_interp.hh"
+#include "../util/ghost_zone.hh"
+#include "../util/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "AHFinderDirect.hh"
+
+//******************************************************************************
+
+//
+// ***** access to persistent data *****
+//
+extern struct state state;
+
+//******************************************************************************
+
+//
+// ***** prototypes for functions local to this file
+//
+namespace {
+bool find_horizon(enum method method,
+		  enum verbose_level verbose_level, int timer_handle,
+		  struct IO_info& IO_info,
+		  struct Jacobian_info& Jac_info,
+		  struct solver_info& solver_info,
+		  struct cactus_grid_info& cgi, struct geometry_info& gi,
+		  patch_system& ps, Jacobian* Jac_ptr,
+		  int hn);
+	  };
+
+//******************************************************************************
+
+//
+// This function is called by the Cactus scheduler to find the apparent
+// horizon or horizons in the current slice.
+//
+extern "C"
+  void AHFinderDirect_find_horizons(CCTK_ARGUMENTS)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+state.IO_info.time_iteration = cctk_iteration;
+
+if (state.timer_handle >= 0)
+   then CCTK_TimerResetI(state.timer_handle);
+
+	for (int hn = 0 ; hn < state.N_horizons ; ++hn)
+	{
+	struct AH_info& AH_info = * state.AH_info_ptrs[hn];
+
+	AH_info.AH_found
+	   = find_horizon(state.method,
+			  state.verbose_level, state.timer_handle,
+			  state.IO_info, state.Jac_info, state.solver_info,
+			  state.cgi, state.gi,
+			  * AH_info.ps_ptr, AH_info.Jac_ptr,
+			  hn);
+
+	// FIXME FIXME: compute the centroid, area, mass, ...
+	}
+
+if (state.timer_handle >= 0)
+   then {
+	printf("timer stats for computation:\n");
+	CCTK_TimerPrintDataI(state.timer_handle, -1);
+	}
+}
+
+//******************************************************************************
+
+//
+// This function finds (or more accurately tries to find) a single
+// apparent horizon.
+//
+// Arguments:
+// timer_handle = a valid Cactus timer handle if we want to time the
+//		  apparent horizon process, or -ve to skip this
+//		  (we only time the computation, not the file I/O)
+// Jac_ptr = may be NULL if no Jacobian is needed (depending on  method)
+// hn = the horizon number (used only in naming output files)
+//
+// Results:
+// This function returns a success flag: this is true if (and only if)
+// it found an h satisfying H(h) = 0 to within the error tolerances,
+// otherwise it's false.
+//
+namespace {
+bool find_horizon(enum method method,
+		  enum verbose_level verbose_level, int timer_handle,
+		  struct IO_info& IO_info,
+		  struct Jacobian_info& Jac_info,
+		  struct solver_info& solver_info,
+		  struct cactus_grid_info& cgi, struct geometry_info& gi,
+		  patch_system& ps, Jacobian* Jac_ptr,
+		  int hn)
+{
+bool success = false;
+
+const bool print_io_msg_flag
+	= verbose_level >= verbose_level__algorithm_highlights;
+
+switch	(method)
+	{
+// just evaluate the horizon function
+case method__horizon_function:
+	  {
+	jtutil::norm<fp> H_norms;
+
+	if (timer_handle >= 0)
+	   then CCTK_TimerStartI(timer_handle);
+	horizon_function(ps, cgi, gi, false, true, &H_norms);
+	if (timer_handle >= 0)
+	   then CCTK_TimerStopI(timer_handle);
+
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   H(h) rms-norm %.2e, infinity-norm %.2e",
+		   H_norms.rms_norm(), H_norms.infinity_norm());
+	output_gridfn(ps, gfns::gfn__H,
+		      IO_info, IO_info.H_base_file_name,
+		      hn, true);
+	break;
+	  }
+
+// just compute/print the NP Jacobian
+case method__Jacobian_test_NP_only:
+	  {
+	Jacobian& Jac_NP = *Jac_ptr;
+	horizon_function(ps, cgi, gi, true);
+	Jac_info.Jacobian_method = Jacobian_method__numerical_perturb;
+	horizon_Jacobian(ps, Jac_NP,
+			 cgi, gi, Jac_info,
+			 true);
+
+	print_Jacobians(ps,
+			& Jac_NP, NULL,
+			IO_info, IO_info.Jacobian_base_file_name,
+			hn, true);
+	break;
+	  }
+
+// compute/print the Jacobian by all possible methods
+case method__Jacobian_test:
+	  {
+	Jacobian& Jac_NP = *Jac_ptr;
+	horizon_function(ps, cgi, gi, true);
+	Jac_info.Jacobian_method = Jacobian_method__numerical_perturb;
+	horizon_Jacobian(ps, Jac_NP,
+			 cgi, gi, Jac_info,
+			 true);
+
+	// symbolic differentiation with finite diff d/dr
+	Jacobian& Jac_SD_FDdr
+		= new_Jacobian(ps, Jac_info.Jacobian_storage_method);
+	horizon_function(ps, cgi, gi, true);
+	Jac_info.Jacobian_method = Jacobian_method__symbolic_diff_with_FD_dr;
+	horizon_Jacobian(ps, Jac_SD_FDdr,
+			 cgi, gi, Jac_info,
+			 true);
+
+	print_Jacobians(ps,
+			& Jac_NP, & Jac_SD_FDdr,
+			IO_info, IO_info.Jacobian_base_file_name,
+			hn, true);
+	break;
+	  }
+
+// find the apparent horizon via the Newton solver
+case method__Newton_solve:
+	  {
+	Jacobian& Jac = *Jac_ptr;
+
+	if (timer_handle >= 0)
+	   then CCTK_TimerStartI(timer_handle);
+	success = Newton_solve(ps, Jac,
+			       cgi, gi,
+			       Jac_info, solver_info, IO_info,
+			       hn, verbose_level);
+	if (timer_handle >= 0)
+	   then CCTK_TimerStopI(timer_handle);
+
+	output_gridfn(ps, gfns::gfn__h,
+		      IO_info, IO_info.h_base_file_name,
+		      hn, print_io_msg_flag);
+	output_gridfn(ps, gfns::gfn__H,
+		      IO_info, IO_info.H_base_file_name,
+		      hn, print_io_msg_flag);
+	break;
+	  }
+
+default:
+	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   find_horizons(): unknown method=(int)%d!\n"
+"                    (this should never happen!)"
+,
+		   (int) method);				/*NOTREACHED*/
+	}
+
+return success;
+}
+	  }
diff --git a/src/driver/horizon_Jacobian.cc b/src/driver/horizon_Jacobian.cc
new file mode 100644
index 0000000..0ef098f
--- /dev/null
+++ b/src/driver/horizon_Jacobian.cc
@@ -0,0 +1,480 @@
+// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h)
+// $Id$
+//
+// <<<prototypes for functions local to this file>>>
+//
+// horizon_Jacobian - top-level driver to compute the Jacobian
+///
+/// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation
+///
+/// horizon_Jacobian_partial_SD - compute the partial-deriv terms: symbolic diff
+/// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian
+///
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <vector>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../util/coords.hh"
+#include "../util/grid.hh"
+#include "../util/fd_grid.hh"
+#include "../util/patch.hh"
+#include "../util/patch_edge.hh"
+#include "../util/patch_interp.hh"
+#include "../util/ghost_zone.hh"
+#include "../util/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "AHFinderDirect.hh"
+
+//******************************************************************************
+
+//
+// ***** prototypes for functions local to this file *****
+//
+
+namespace {
+void horizon_Jacobian_NP(patch_system& ps,
+			 Jacobian& Jac,
+			 const struct cactus_grid_info& cgi,
+			 const struct geometry_info& ii,
+			 const struct Jacobian_info& Jacobian_info,
+			 bool print_msg_flag);
+
+void horizon_Jacobian_partial_SD(patch_system& ps,
+				 Jacobian& Jac,
+				 const struct cactus_grid_info& cgi,
+				 const struct geometry_info& gi,
+				 const struct Jacobian_info& Jacobian_info,
+				 bool print_msg_flag);
+void add_ghost_zone_Jacobian(const patch_system& ps,
+			     Jacobian& Jac,
+			     fp mol,
+			     const patch& xp, const ghost_zone& xmgz,
+			     int x_II,
+			     int xm_irho, int xm_isigma);
+void horizon_Jacobian_dr_FD(patch_system& ps,
+			    Jacobian& Jac,
+			    const struct cactus_grid_info& cgi,
+			    const struct geometry_info& gi,
+			    const struct Jacobian_info& Jacobian_info,
+			    bool print_msg_flag);
+	  }
+
+//******************************************************************************
+
+//
+// This function is a top-level driver to compute the Jacobian matrix
+// J[H(h)] of the LHS function H(h).  It just decodes the Jacobian_method
+// parameter and calls the appropriate subfunction.
+//
+void horizon_Jacobian(patch_system& ps,
+		      Jacobian& Jac,
+		      const struct cactus_grid_info& cgi,
+		      const struct geometry_info& gi,
+		      const struct Jacobian_info& Jacobian_info,
+		      bool print_msg_flag)
+{
+switch	(Jacobian_info.Jacobian_method)
+	{
+case Jacobian_method__numerical_perturb:
+	horizon_Jacobian_NP(ps, Jac,
+			    cgi, gi, Jacobian_info,
+			    print_msg_flag);
+	break;
+case Jacobian_method__symbolic_diff:
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   horizon_Jacobian(): Jacobian_method == \"symbolic differentiation\"\n"
+"                       isn't implemented (yet); reverting to\n"
+"                       \"symbolic differentiation with finite diff d/dr\"");
+	// fall through
+case Jacobian_method__symbolic_diff_with_FD_dr:
+	horizon_Jacobian_partial_SD(ps, Jac,
+				    cgi, gi, Jacobian_info,
+				    print_msg_flag);
+	horizon_Jacobian_dr_FD(ps, Jac,
+			       cgi, gi, Jacobian_info,
+			       print_msg_flag);
+	break;
+default:
+	error_exit(PANIC_EXIT,
+"***** horizon_Jacobian(): unknown Jacobian_info.Jacobian_method=(int)%d!\n"
+"                          (this should never happen!)\n"
+,
+		   int(Jacobian_info.Jacobian_method));		/*NOTREACHED*/
+	}
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function computes the Jacobian matrix of the LHS function H(h)
+// by numerical perturbation of the H(h) function.  The algorithm is as
+// follows:
+//
+// we assume that H = H(h) has already been evaluated
+// save_H = H
+//	for each point (y,JJ)
+//	{
+//	const fp save_h_y = h at y;
+//	h at y += perturbation_amplitude;
+//	evaluate H(h) (silently)
+//		for each point (x,II)
+//		{
+//		Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude;
+//		}
+//	h at y = save_h_y;
+//	}
+// H = save_H
+//
+// Inputs (angular gridfns, on ghosted grid):
+//	h				# shape of trial surface
+//	H				# H(h) assumed to already be computed
+//
+// Outputs:
+//	The Jacobian matrix is stored in the Jacobian object Jac.
+//
+namespace {
+void horizon_Jacobian_NP(patch_system& ps,
+			 Jacobian& Jac,
+			 const struct cactus_grid_info& cgi,
+			 const struct geometry_info& ii,
+			 const struct Jacobian_info& Jacobian_info,
+			 bool print_msg_flag)
+{
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   horizon Jacobian (numerical perturbation)");
+const fp epsilon = Jacobian_info.perturbation_amplitude;
+
+ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
+
+	for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
+	{
+	patch& yp = ps.ith_patch(ypn);
+	if (print_msg_flag)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "      perturbing in %s patch",
+			   yp.name());
+
+	for (int y_irho = yp.min_irho() ; y_irho <= yp.max_irho() ; ++y_irho)
+	{
+	for (int y_isigma = yp.min_isigma() ;
+	     y_isigma <= yp.max_isigma() ;
+	     ++y_isigma)
+	{
+	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);
+
+	const fp save_h_y = yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma);
+	yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) += epsilon;
+	horizon_function(ps, cgi, ii);
+		for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
+		{
+		patch& xp = ps.ith_patch(xpn);
+
+		for (int x_irho = xp.min_irho() ;
+		     x_irho <= xp.max_irho() ;
+		     ++x_irho)
+		{
+		for (int x_isigma = xp.min_isigma() ;
+		     x_isigma <= xp.max_isigma() ;
+		     ++x_isigma)
+		{
+		const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);
+		const fp old_H = xp.gridfn(gfns::gfn__save_H, x_irho,x_isigma);
+		const fp new_H = xp.gridfn(gfns::gfn__H, x_irho,x_isigma);
+		Jac(II,JJ) = (new_H - old_H) / epsilon;
+		}
+		}
+		}
+	yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) = save_h_y;
+	}
+	}
+   	} 
+
+ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H);
+}
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function computes the partial derivative terms in the Jacobian
+// matrix of the LHS function H(h), by symbolic differentiation from the
+// Jacobian coefficient (angular) gridfns.  The computation is done a
+// Jacobian row at a time, using equation (25) of my 1996 apparent horizon
+// finding paper.
+//
+// Inputs (angular gridfns, on ghosted grid):
+//	h				# shape of trial surface
+//	H				# H(h) assumed to already be computed
+//	partial_H_wrt_partial_d_h	# Jacobian coefficients
+//	partial_H_wrt_partial_dd_h	# (also assumed to already be computed)
+//
+// Outputs:
+//	The Jacobian matrix is stored in the Jacobian object Jac.
+//
+namespace {
+void horizon_Jacobian_partial_SD(patch_system& ps,
+				 Jacobian& Jac,
+				 const struct cactus_grid_info& cgi,
+				 const struct geometry_info& gi,
+				 const struct Jacobian_info& Jacobian_info,
+				 bool print_msg_flag)
+{
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   horizon Jacobian: partial-deriv terms (symbolic diff)");
+
+Jac.zero_matrix();
+ps.compute_synchronize_Jacobian();
+
+    for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
+    {
+    patch& xp = ps.ith_patch(xpn);
+
+	for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
+	{
+	for (int x_isigma = xp.min_isigma() ;
+	x_isigma <= xp.max_isigma() ;
+	++x_isigma)
+	{
+	//
+	// compute the main Jacobian terms for this grid point, i.e.
+	//	  partial H(this point x, Jacobian row II)
+	//	---------------------------------------------
+	//	partial h(other points y, Jacobian column JJ)
+	//
+
+	// Jacobian row index
+	const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma);
+
+	// Jacobian coefficients for this point
+	const fp Jacobian_coeff_rho
+	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_1,
+		       x_irho, x_isigma);
+	const fp Jacobian_coeff_sigma
+	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_2,
+		       x_irho, x_isigma);
+	const fp Jacobian_coeff_rho_rho
+	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_11,
+		       x_irho, x_isigma);
+	const fp Jacobian_coeff_rho_sigma
+	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_12,
+		       x_irho, x_isigma);
+	const fp Jacobian_coeff_sigma_sigma
+	   = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_22,
+		       x_irho, x_isigma);
+
+	// partial_rho, partial_rho_rho
+	    for (int m_irho = xp.molecule_min_m() ;
+		 m_irho <= xp.molecule_max_m() ;
+		 ++m_irho)
+	    {
+	    const int xm_irho = x_irho + m_irho;
+	    const fp Jac_rho     = Jacobian_coeff_rho
+				   * xp.partial_rho_coeff(m_irho);
+	    const fp Jac_rho_rho = Jacobian_coeff_rho_rho
+				   * xp.partial_rho_rho_coeff(m_irho);
+	    const fp Jac_sum = Jac_rho + Jac_rho_rho;
+	    if (xp.is_in_nominal_grid(xm_irho, x_isigma))
+	       then Jac(II, xp,xm_irho,x_isigma) += Jac_sum;
+	       else add_ghost_zone_Jacobian
+			(ps, Jac,
+			 Jac_sum,
+			 xp, xp.minmax_rho_ghost_zone(m_irho < 0),
+			 II, xm_irho, x_isigma);
+	    }
+
+	// partial_sigma, partial_sigma_sigma
+	    for (int m_isigma = xp.molecule_min_m() ;
+		 m_isigma <= xp.molecule_max_m() ;
+		 ++m_isigma)
+	    {
+	    const int xm_isigma = x_isigma + m_isigma;
+	    const fp Jac_sigma       = Jacobian_coeff_sigma
+				       * xp.partial_sigma_coeff(m_isigma);
+	    const fp Jac_sigma_sigma = Jacobian_coeff_sigma_sigma
+				       * xp.partial_sigma_sigma_coeff(m_isigma);
+	    const fp Jac_sum = Jac_sigma + Jac_sigma_sigma;
+	    if (xp.is_in_nominal_grid(x_irho, xm_isigma))
+	       then Jac(II, xp,x_irho,xm_isigma) += Jac_sum;
+	       else add_ghost_zone_Jacobian
+			(ps, Jac,
+			 Jac_sum,
+			 xp, xp.minmax_sigma_ghost_zone(m_isigma < 0),
+			 II, x_irho, xm_isigma);
+	    }
+
+	// partial_rho_sigma
+	    for (int m_irho = xp.molecule_min_m() ;
+		 m_irho <= xp.molecule_max_m() ;
+		 ++m_irho)
+	    {
+	    for (int m_isigma = xp.molecule_min_m() ;
+		 m_isigma <= xp.molecule_max_m() ;
+		 ++m_isigma)
+	    {
+	    const int xm_irho   = x_irho   + m_irho;
+	    const int xm_isigma = x_isigma + m_isigma;
+	    const fp Jac_rho_sigma
+	       = Jacobian_coeff_rho_sigma
+		 * xp.partial_rho_sigma_coeff(m_irho, m_isigma);
+	    if (xp.is_in_nominal_grid(xm_irho, xm_isigma))
+	       then Jac(II, xp,xm_irho,xm_isigma) += Jac_rho_sigma;
+	       else {
+		    const ghost_zone& xmgz
+		       = xp.corner_ghost_zone_containing_point
+				(m_irho < 0, m_isigma < 0,
+				 xm_irho, xm_isigma);
+		    add_ghost_zone_Jacobian(ps, Jac,
+					    Jac_rho_sigma,
+					    xp, xmgz,
+					    II, xm_irho, xm_isigma);
+		    }
+	    }
+	    }
+
+	}
+	}
+    }
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function adds the ghost-zone Jacobian dependency contributions
+// for a single ghost-zone point, to a Jacobian matrix.
+//
+// Arguments:
+// ps = The patch system.
+// Jac = (out) The Jacobian matrix.
+// mol = The molecule coefficient.
+// xp = The patch containing the center point of the molecule.
+// xmgz = If the x+m point is in a ghost zone, this must be that ghost zone.
+//	  If the x+m point is not in a ghost zone, this argument is ignored.
+// x_II = The Jacobian row of the x point.
+// xm_(irho,isigma) = The coordinates (in xp) of the x+m point of the molecule.
+//
+namespace {
+void add_ghost_zone_Jacobian(const patch_system& ps,
+			     Jacobian& Jac,
+			     fp mol,
+			     const patch& xp, const ghost_zone& xmgz,
+			     int x_II,
+			     int xm_irho, int xm_isigma)
+{
+const patch_edge& xme = xmgz.my_edge();
+const int xm_iperp = xme.iperp_of_irho_isigma(xm_irho, xm_isigma);
+const int xm_ipar  = xme. ipar_of_irho_isigma(xm_irho, xm_isigma);
+
+// FIXME: this won't change from one call to another
+//        ==> it would be more efficient to reuse the same buffer
+//            across multiple calls on this function
+int global_min_ym, global_max_ym;
+ps.synchronize_Jacobian_global_minmax_ym(global_min_ym, global_max_ym);
+jtutil::array1d<fp> Jacobian_buffer(global_min_ym, global_max_ym);
+
+// on what other points y does this molecule point xm depend
+// via the patch_system::synchronize() operation?
+int y_iperp;
+int y_posn, min_ym, max_ym;
+const patch_edge& ye = ps.synchronize_Jacobian(xmgz,
+					       xm_iperp, xm_ipar,
+					       y_iperp,
+					       y_posn, min_ym, max_ym,
+					       Jacobian_buffer);
+const patch& yp = ye.my_patch();
+
+// add the Jacobian contributions from the ym points
+	for (int ym = min_ym ; ym <= max_ym ; ++ym)
+	{
+	const int y_ipar = y_posn + ym;
+	const int y_irho   = ye.  irho_of_iperp_ipar(y_iperp,y_ipar);
+	const int y_isigma = ye.isigma_of_iperp_ipar(y_iperp,y_ipar);
+	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho, y_isigma);
+	Jac(x_II,JJ) += mol*Jacobian_buffer(ym);
+	}
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function sums the d/dr terms into the Jacobian matrix of the
+// LHS function H(h), computing those terms by finite differencing.
+//
+// The basic algorithm is that
+//	Jac += diag[ (H(h+epsilon) - H(h)) / epsilon ]
+//
+// Inputs (angular gridfns, on ghosted grid):
+//	h				# shape of trial surface
+//	H				# H(h) assumed to already be computed
+//
+// Outputs:
+//	Jac += d/dr terms
+//
+namespace {
+void horizon_Jacobian_dr_FD(patch_system& ps,
+			    Jacobian& Jac,
+			    const struct cactus_grid_info& cgi,
+			    const struct geometry_info& gi,
+			    const struct Jacobian_info& Jacobian_info,
+			    bool print_msg_flag)
+{
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   horizon Jacobian: d/dr terms (finite diff)");
+
+const fp epsilon = Jacobian_info.perturbation_amplitude;
+
+// compute H(h+epsilon)
+ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
+ps.add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
+horizon_function(ps, cgi, gi);
+
+	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+	{
+	patch& p = ps.ith_patch(pn);
+		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
+		{
+		for (int isigma = p.min_isigma() ;
+		     isigma <= p.max_isigma() ;
+		     ++isigma)
+		{
+		const int II = ps.gpn_of_patch_irho_isigma(p, irho,isigma);
+		const fp old_H = p.gridfn(gfns::gfn__save_H, irho,isigma);
+		const fp new_H = p.gridfn(gfns::gfn__H, irho,isigma);
+		Jac(II,II) += (new_H - old_H) / epsilon;
+		}
+		}
+	}
+
+// restore h and H
+ps.add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
+ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H);
+}
+	  }
diff --git a/src/driver/initial_guess.cc b/src/driver/initial_guess.cc
new file mode 100644
index 0000000..7bb62df
--- /dev/null
+++ b/src/driver/initial_guess.cc
@@ -0,0 +1,298 @@
+// initial_guess.cc -- top level driver for setting the initial guess
+// $Header$
+//
+// <<<access to persistent data>>>
+// <<<prototypes for functions local to this file>>>
+// AHFinderDirect_initial_guess - top-level driver
+/// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords)
+/// setup_ellipsoid - setup up a coordinate ellipsoid in h
+///
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <vector>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../util/coords.hh"
+#include "../util/grid.hh"
+#include "../util/fd_grid.hh"
+#include "../util/patch.hh"
+#include "../util/patch_edge.hh"
+#include "../util/patch_interp.hh"
+#include "../util/ghost_zone.hh"
+#include "../util/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "AHFinderDirect.hh"
+
+//******************************************************************************
+
+//
+// ***** access to persistent data *****
+//
+extern struct state state;
+
+//******************************************************************************
+
+//
+// ***** prototypes for functions local to this file *****
+//
+
+namespace {
+void setup_Kerr_horizon(patch_system& ps,
+			fp x_posn, fp y_posn, fp z_posn,
+			fp m, fp a,
+			bool Kerr_Schild_flag,
+			bool my_print_msg_flag, bool ellipsoid_print_msg_flag);
+void setup_ellipsoid(patch_system& ps,
+		     fp x_center, fp y_center, fp z_center,
+		     fp x_radius, fp y_radius, fp z_radius,
+		     bool print_msg_flag);
+	  };
+
+//******************************************************************************
+
+//
+// This function is called by the Cactus scheduler to set our initial guess
+// for the horizon position(s).
+//
+extern "C"
+  void AHFinderDirect_initial_guess(CCTK_ARGUMENTS)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+const bool print_msg_flag
+	= state.verbose_level >= verbose_level__algorithm_highlights;
+const bool print_more_msg_flag
+	= state.verbose_level >= verbose_level__algorithm_details;
+
+struct IO_info& IO_info = state.IO_info;
+IO_info.time_iteration = cctk_iteration;
+
+	for (int hn = 0 ; hn < state.N_horizons ; ++hn)
+	{
+	struct AH_info& AH_info = * state.AH_info_ptrs[hn];
+	patch_system& ps = * AH_info.ps_ptr;
+
+	if	(STRING_EQUAL(initial_guess_method, "read from file"))
+	   then input_gridfn(ps, gfns::gfn__h,
+			     IO_info, IO_info.h_base_file_name,
+			     hn, print_msg_flag);
+	else if (STRING_EQUAL(initial_guess_method, "sphere"))
+	   then setup_ellipsoid(ps,
+				initial_guess__sphere__x_center[hn],
+				initial_guess__sphere__y_center[hn],
+				initial_guess__sphere__z_center[hn],
+				initial_guess__sphere__radius[hn],
+				initial_guess__sphere__radius[hn],
+				initial_guess__sphere__radius[hn],
+				print_msg_flag);
+	else if (STRING_EQUAL(initial_guess_method, "ellipsoid"))
+	   then setup_ellipsoid(ps,
+				initial_guess__ellipsoid__x_center[hn],
+				initial_guess__ellipsoid__y_center[hn],
+				initial_guess__ellipsoid__z_center[hn],
+				initial_guess__ellipsoid__x_radius[hn],
+				initial_guess__ellipsoid__y_radius[hn],
+				initial_guess__ellipsoid__z_radius[hn],
+				print_msg_flag);
+	else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr"))
+	   then setup_Kerr_horizon(ps,
+				   initial_guess__Kerr__x_posn[hn],
+				   initial_guess__Kerr__y_posn[hn],
+				   initial_guess__Kerr__z_posn[hn],
+				   initial_guess__Kerr__mass[hn],
+				   initial_guess__Kerr__spin[hn],
+				   false,	// use Kerr coords
+				   print_msg_flag, print_more_msg_flag);
+	else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild"))
+	   then setup_Kerr_horizon(ps,
+				   initial_guess__Kerr__x_posn[hn],
+				   initial_guess__Kerr__y_posn[hn],
+				   initial_guess__Kerr__z_posn[hn],
+				   initial_guess__Kerr__mass[hn],
+				   initial_guess__Kerr__spin[hn],
+				   true,	// use Kerr-Schild coords
+				   print_msg_flag, print_more_msg_flag);
+	else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "unknown initial_guess_method=\"%s\"!",
+			   initial_guess_method);	/*NOTREACHED*/
+
+	// write initial guess back to the data file?
+	if (output_initial_guess != 0)
+	   then output_gridfn(ps, gfns::gfn__h,
+			      IO_info, IO_info.h_base_file_name,
+			      hn, print_msg_flag);
+	}
+}
+
+//******************************************************************************
+
+//
+// This function sets up the horizon of a Kerr black hole in Kerr or
+// Kerr-Schild coordinates, on the nominal grid, in the  h  gridfn.
+//
+// Kerr-Schild coordinates are described in MTW Exercise 33.8, page 903,
+// and the horizon is worked out on page 13.2 of my AHFinderDirect notes.
+//
+// Arguments:
+// [xyz]_posn = The position of the Kerr black hole.
+// (m,a) = Describe the Kerr black hole.  Note that my convention has
+//	   a=J/m^2 dimensionless, while MTW take a=J/m=m*(my a).
+// Kerr_Schild_flag = false to use Kerr coordinates,
+//		      true to use Kerr-Schild coordinates
+// my_print_msg_flag = true to print messages about the Kerr horizon itself,
+//		       false to skip these
+// ellipsoid_print_msg_flag = true to print messages about the ellipsoid
+//				   shape of the Kerr horizon,
+//			      false to skip these
+//
+namespace {
+void setup_Kerr_horizon(patch_system& ps,
+			fp x_posn, fp y_posn, fp z_posn,
+			fp m, fp a,
+			bool Kerr_Schild_flag,
+			bool my_print_msg_flag, bool ellipsoid_print_msg_flag)
+{
+const char* const name = Kerr_Schild_flag ? "Kerr-Schild" : "Kerr";
+
+if (my_print_msg_flag)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "setting up Kerr/%s initial guess",
+		   name);
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "           posn=(%g,%g,%g) mass=%g spin=J/m^2=%g",
+		   double(x_posn), double(y_posn), double(z_posn),
+		   double(m), double(a));
+	}
+
+// horizon in Kerr coordinates is coordinate sphere
+const fp r = m * (1.0 + sqrt(1.0 - a*a));
+
+// horizon in Kerr-Schild coordinates is coordinate ellipsoid
+const fp  z_radius = r;
+const fp xy_radius = Kerr_Schild_flag ? r * sqrt(1.0 + a*a*m*m/(r*r)) : r;
+
+if (ellipsoid_print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "           horizon is coordinate %s",
+		   Kerr_Schild_flag ? "ellipsoid" : "sphere");
+setup_ellipsoid(ps,
+		x_posn, y_posn, z_posn,
+		xy_radius, xy_radius, z_radius,
+		ellipsoid_print_msg_flag);
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function sets up an ellipsoid in the gridfn h, using the
+// formulas in "ellipsoid.maple" and the Maple-generated C code in
+// "ellipsoid.c":
+//
+// ellipsoid has center (A,B,C), radius (a,b,c)
+// angular coordinate system has center (U,V,W)
+//
+// direction cosines wrt angular coordinate center are (xcos,ycos,zcos)
+// i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
+//
+// then the equation of the ellipsoid is
+//	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*r - C)^2
+//	-----------------  +  ----------------  +  -----------------  =  1
+//	        a^2                  b^2                   c^2
+//
+// to solve this, we introduce intermediate variables
+//	AU = A - U
+//	BV = B - V
+//	CW = C - W
+//
+namespace {
+void setup_ellipsoid(patch_system& ps,
+		     fp x_center, fp y_center, fp z_center,
+		     fp x_radius, fp y_radius, fp z_radius,
+		     bool print_msg_flag)
+{
+if (print_msg_flag)
+   then {
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "setting h = ellipsoid: center=(%g,%g,%g)",
+		   double(x_center), double(y_center), double(z_center));
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "                       radius=(%g,%g,%g)",
+		   double(x_radius), double(y_radius), double(z_radius));
+	}
+
+	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+	{
+	patch& p = ps.ith_patch(pn);
+
+		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
+		{
+		for (int isigma = p.min_isigma() ;
+		     isigma <= p.max_isigma() ;
+		     ++isigma)
+		{
+		const fp rho = p.rho_of_irho(irho);
+		const fp sigma = p.sigma_of_isigma(isigma);
+		fp xcos, ycos, zcos;
+		p.xyzcos_of_rho_sigma(rho,sigma, xcos,ycos,zcos);
+
+		// set up variables used by Maple-generated code
+		const fp AU = x_center - ps.origin_x();
+		const fp BV = y_center - ps.origin_y();
+		const fp CW = z_center - ps.origin_z();
+		const fp a = x_radius;
+		const fp b = y_radius;
+		const fp c = z_radius;
+
+		// compute the solutions r_plus and r_minus
+		fp r_plus, r_minus;
+		#include "ellipsoid.c"
+
+		// exactly one of the solutions (call it r) should be positive
+		fp r;
+		if      ((r_plus > 0.0) && (r_minus < 0.0))
+		   then r = r_plus;
+		else if ((r_plus < 0.0) && (r_minus > 0.0))
+		   then r = r_minus;
+		else    CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+				   "\n"
+"   setup_ellipsoid():\n"
+"        expected exactly one r>0 solution to quadratic, got 0 or 2!\n"
+"        %s patch (irho,isigma)=(%d,%d) ==> (rho,sigma)=(%g,%g)\n"
+"        direction cosines (xcos,ycos,zcos)=(%g,%g,%g)\n"
+"        ==> r_plus=%g r_minus=%g\n"
+				   ,
+				   p.name(), irho, isigma,
+				   double(rho), double(sigma),
+				   double(xcos), double(ycos), double(zcos),
+				   double(r_plus), double(r_minus));
+		   						/*NOTREACHED*/
+
+		// r = horizon radius at this grid point
+		p.ghosted_gridfn(gfns::gfn__h, irho,isigma) = r;
+		}
+		}
+	}
+}
+	  }
diff --git a/src/driver/io.cc b/src/driver/io.cc
new file mode 100644
index 0000000..94eb5bb
--- /dev/null
+++ b/src/driver/io.cc
@@ -0,0 +1,290 @@
+// io.cc -- I/O routines for this thorn
+// $Header$
+//
+// <<<prototypes for functions local to this file>>>
+// input_gridfn
+// output_gridfn
+// output_Jacobians
+//
+/// io_file_name
+//
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <vector>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../util/coords.hh"
+#include "../util/grid.hh"
+#include "../util/fd_grid.hh"
+#include "../util/patch.hh"
+#include "../util/patch_edge.hh"
+#include "../util/patch_interp.hh"
+#include "../util/ghost_zone.hh"
+#include "../util/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "AHFinderDirect.hh"
+
+//******************************************************************************
+
+//
+// ***** prototypes for functions local to this file ***** 
+//
+
+namespace {
+const char* io_file_name(struct IO_info& IO_info, const char base_file_name[],
+			 int hn, int AHF_iteration = 0);
+	  };
+
+//******************************************************************************
+
+//
+// This function inputs a gridfn from a data file.
+//
+// We assume that this gridfn is ghosted, but the ghost zones are *not*
+// present in the data file.
+//
+void input_gridfn(patch_system& ps, int unknown_gfn,
+		  struct IO_info& IO_info, const char base_file_name[],
+		  int hn, bool print_msg_flag, int AHF_iteration = 0)
+{
+const char* file_name = io_file_name(IO_info, base_file_name,
+				     hn, AHF_iteration);
+if (print_msg_flag)
+   then {
+	if ( (unknown_gfn == gfns::gfn__h)
+	     && (IO_info.time_iteration == 0) && (AHF_iteration == 0) )
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "   reading initial guess from \"%s\"", file_name);
+	   else CCTK_VInfo(CCTK_THORNSTRING,
+			   "   reading \"%s\"", file_name);
+	}
+
+ps.read_ghosted_gridfn(unknown_gfn,
+		       file_name,
+		       false);			// no ghost zones in data file
+}
+
+//******************************************************************************
+
+//
+// This function outputs a gridfn from a data file.
+//
+// If the gridfn is h, then we also write out the xyz positions of the
+// horizon surface points.
+//
+// FIXME FIXME: if the gridfn is not h, we assume that it's nominal-grid.
+//
+void output_gridfn(patch_system& ps, int unknown_gfn,
+		   struct IO_info& IO_info, const char base_file_name[],
+		   int hn, bool print_msg_flag, int AHF_iteration = 0)
+{
+const char* file_name = io_file_name(IO_info, base_file_name,
+				     hn, AHF_iteration);
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   writing \"%s\"", file_name);
+
+switch	(unknown_gfn)
+	{
+case gfns::gfn__h:
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   writing h to \"%s\"", file_name);
+	ps.print_ghosted_gridfn_with_xyz(unknown_gfn,
+					 true, gfns::gfn__h,
+					 file_name,
+					 false);	// no ghost zones
+	break;
+case gfns::gfn__H:
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   writing H(h) to \"%s\"", file_name);
+	ps.print_gridfn(unknown_gfn,
+			file_name);
+	break;
+default:
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   writing \"%s\"", file_name);
+	ps.print_gridfn(unknown_gfn,
+			file_name);
+	break;
+	}
+}
+
+//******************************************************************************
+
+//
+// This function prints one or two Jacobian matrices (and their difference
+// in the latter case) to a named output file.
+//
+// Arguments:
+// A null Jacobian pointer means to skip that Jacobian.
+//
+void print_Jacobians(const patch_system& ps,
+		     const Jacobian* Jac_NP, const Jacobian* Jac_SD_FDdr,
+		     struct IO_info& IO_info, const char base_file_name[],
+		     int hn, bool print_msg_flag, int AHF_iteration = 0)
+{
+if (Jac_NP == NULL)
+   then {
+	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "print_Jacobinas(): Jac_NP == NULL is not (yet) supported!");
+	return;						// *** ERROR RETURN ***
+	}
+
+const char* file_name = io_file_name(IO_info, base_file_name,
+				     hn, AHF_iteration);
+if (print_msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "   writing %s to \"%s\"",
+		   ((Jac_SD_FDdr == NULL) ? "NP Jacobian"
+					  : "NP and SD_FDdr Jacobians"),
+		   file_name);
+
+FILE *fileptr = fopen(file_name, "w");
+if (fileptr == NULL)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "print_Jacobians(): can't open output file \"%s\"!",
+		   file_name);					/*NOTREACHED*/
+
+fprintf(fileptr, "# column 1 = x patch name\n");
+fprintf(fileptr, "# column 2 = x patch number\n");
+fprintf(fileptr, "# column 3 = x irho\n");
+fprintf(fileptr, "# column 4 = x isigma\n");
+fprintf(fileptr, "# column 5 = x II\n");
+fprintf(fileptr, "# column 6 = y patch name\n");
+fprintf(fileptr, "# column 7 = y patch number\n");
+fprintf(fileptr, "# column 8 = y irho\n");
+fprintf(fileptr, "# column 9 = y isigma\n");
+fprintf(fileptr, "# column 10 = y JJ\n");
+fprintf(fileptr, "# column 11 = Jac_NP(II,JJ)\n");
+if (Jac_SD_FDdr != NULL)
+   then {
+	fprintf(fileptr, "# column 12 = Jac_SD_FDdr(II,JJ)\n");
+	fprintf(fileptr, "# column 13 = abs error in Jac_SD_FDdr\n");
+	fprintf(fileptr, "# column 14 = rel error in Jac_SD_FDdr\n");
+	}
+
+    for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
+    {
+    patch& xp = ps.ith_patch(xpn);
+
+    for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
+    {
+    for (int x_isigma = xp.min_isigma() ;
+	 x_isigma <= xp.max_isigma() ;
+	 ++x_isigma)
+    {
+    const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);
+
+	for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
+	{
+	patch& yp = ps.ith_patch(ypn);
+
+	for (int y_irho = yp.min_irho() ;
+	     y_irho <= yp.max_irho() ;
+	     ++y_irho)
+	{
+	for (int y_isigma = yp.min_isigma() ;
+	     y_isigma <= yp.max_isigma() ;
+	     ++y_isigma)
+	{
+	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);
+
+	if (! Jac_NP->is_explicitly_stored(II,JJ))
+	   then continue;			// skip sparse points
+
+	const fp NP = (*Jac_NP)(II,JJ);
+	const fp SD_FDdr = (Jac_SD_FDdr == NULL) ? 0.0 : (*Jac_SD_FDdr)(II,JJ);
+
+	if ((NP == 0.0) && (SD_FDdr == 0.0))
+	   then continue;		// skip zero values (if == )
+
+	fprintf(fileptr,
+		"%s %d %d %d %d\t%s %d %d %d %d\t%.10g",
+		xp.name(), xpn, x_irho, x_isigma, II,
+		yp.name(), ypn, y_irho, y_isigma, JJ,
+		double(NP));
+
+	if (Jac_SD_FDdr != NULL)
+	   then {
+		const fp abs_NP      = jtutil::abs(NP     );
+		const fp abs_SD_FDdr = jtutil::abs(SD_FDdr);
+		const fp max_abs = jtutil::max(abs_SD_FDdr, abs_NP);
+		const fp SD_FDdr_abs_error = SD_FDdr - NP;
+		const fp SD_FDdr_rel_error = SD_FDdr_abs_error / max_abs;
+
+		fprintf(fileptr,
+			"\t%.10g\t%e\t%e",
+			double(SD_FDdr),
+			double(SD_FDdr_abs_error), double(SD_FDdr_rel_error));
+		}
+
+	fprintf(fileptr, "\n");
+	}
+	}
+	}
+    }
+    }
+    }
+
+fclose(fileptr);
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function encapsulates our file-naming conventions.
+//
+// Arguments:
+// IO_info = note this is *not* const since we may keep things like
+//	     file positions etc here
+// base_file_name[] = from the parameter file
+// time_iteration = the Cactus time iteration number (cctk_iteration)
+// hn = the horizon number
+// AHF_iteration = the apparent horizon finder's internal iteration
+//		   number (>= 1) if this is an intermediate iterate,
+//		   or the default (0) if this is a final computed
+//		   horizon position
+//
+// Results:
+// This function returns (a pointer to) the file name.  The returned
+// result points into a private static buffer; the usual caveats apply.
+//
+namespace {
+const char* io_file_name(struct IO_info& IO_info, const char base_file_name[],
+			 int hn, int AHF_iteration = 0)
+{
+const int N_file_name_buffer = 200;
+static char file_name_buffer[N_file_name_buffer];
+
+if (AHF_iteration == 0)
+   then snprintf(file_name_buffer, N_file_name_buffer,
+		 "%s.t%d.ah%d.dat",
+		 base_file_name, IO_info.time_iteration, hn);
+   else snprintf(file_name_buffer, N_file_name_buffer,
+		 "%s.t%d.ah%d.it%d.dat",
+		 base_file_name, IO_info.time_iteration, hn, AHF_iteration);
+
+return file_name_buffer;
+}
+	  }
diff --git a/src/driver/make.code.defn b/src/driver/make.code.defn
new file mode 100644
index 0000000..cac07f9
--- /dev/null
+++ b/src/driver/make.code.defn
@@ -0,0 +1,12 @@
+# specification of things to be compiled in this directory
+# $Header$
+
+# Source files in this directory
+SRCS = state.cc \
+       setup.cc \
+       initial_guess.cc \
+       find_horizons.cc Newton.cc horizon_Jacobian.cc \
+       io.cc
+
+# Subdirectories containing source files
+SUBDIRS =
diff --git a/src/driver/makefile b/src/driver/makefile
new file mode 100644
index 0000000..4be4001
--- /dev/null
+++ b/src/driver/makefile
@@ -0,0 +1,22 @@
+# Makefile for GR code
+# $Id: makefile,v 1.1 2002-09-03 15:08:15 jthorn Exp $
+
+#
+# Environment Variables:
+#	MAPLE_VERSION	used via @ifdef for version control in Maple code;
+#			typically set to something like  MAPLE_V_RELEASE_4
+#			(not presently used, but may be needed in the future)
+#
+# Targets:
+#	ellipsoid	==> compute ellipsoid equations
+#
+
+ifneq ($(MAPLE_VERSION),)
+MPP_FLAGS := -D$(MAPLE_VERSION)
+endif
+
+###############################################################################
+
+.PHONY		: ellipsoid
+ellipsoid	:
+	maple <ellipsoid.maple >ellipsoid.log
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
new file mode 100644
index 0000000..62e3adc
--- /dev/null
+++ b/src/driver/setup.cc
@@ -0,0 +1,336 @@
+// setup.cc -- top level driver to setup our persistent data structures
+// $Header$
+//
+// <<<prototypes for functions local to this file>>>
+// <<<access to persistent data>>>
+//
+// AHFinderDirect_driver - top-level driver
+///
+/// decode_method - decode the method parameter
+/// decode_verbose_level - decode the verbose_level parameter
+/// decode_Jacobian_method - decode the Jacobian_method parameter
+/// decode_geometry_method - decode the geometry_method parameter
+///
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <vector>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../util/coords.hh"
+#include "../util/grid.hh"
+#include "../util/fd_grid.hh"
+#include "../util/patch.hh"
+#include "../util/patch_edge.hh"
+#include "../util/patch_interp.hh"
+#include "../util/ghost_zone.hh"
+#include "../util/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "AHFinderDirect.hh"
+
+//******************************************************************************
+
+//
+// ***** prototypes for functions local to this file *****
+//
+
+namespace {
+enum method
+  decode_method(const char method_string[]);
+enum verbose_level
+  decode_verbose_level(const char verbose_level_string[]);
+enum Jacobian_method
+  decode_Jacobian_method(const char Jacobian_method_string[]);
+enum geometry_method
+  decode_geometry_method(const char geometry_method_string[]);
+	  };
+
+//******************************************************************************
+
+//
+// ***** access to persistent data *****
+//
+extern struct state state;
+
+//******************************************************************************
+
+//
+// This function is called by the Cactus scheduler to set up the
+// following data structures:
+//	state.gi
+//	state.cgi
+//
+extern "C"
+  void AHFinderDirect_setup(CCTK_ARGUMENTS)
+{
+DECLARE_CCTK_ARGUMENTS
+DECLARE_CCTK_PARAMETERS
+
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "initializing AHFinderDirect data structures");
+
+state.N_horizons = N_horizons;
+CCTK_VInfo(CCTK_THORNSTRING,
+           "             to search for %d horizon%s",
+	   state.N_horizons, ((state.N_horizons == 1) ? "" : "s"));
+
+
+//
+// decode/copy parameters into structures
+//
+
+state.method = decode_method(method);
+state.verbose_level = decode_verbose_level(verbose_level);
+state.timer_handle = (print_timing_stats != 0)
+		     ? CCTK_TimerCreateI()
+		     : -1;
+
+struct IO_info& IO_info = state.IO_info;
+IO_info.h_base_file_name        = h_base_file_name;
+IO_info.H_base_file_name        = H_of_h_base_file_name;
+IO_info.Jacobian_base_file_name = Jacobian_base_file_name;
+IO_info.time_iteration = 0;
+
+struct Jacobian_info& Jac_info = state.Jac_info;
+Jac_info.Jacobian_method = decode_Jacobian_method(Jacobian_method);
+Jac_info.Jacobian_storage_method
+	= decode_Jacobian_type(Jacobian_storage_method);
+Jac_info.perturbation_amplitude = Jacobian_perturbation_amplitude;
+
+struct solver_info& solver_info = state.solver_info;
+solver_info.max_Newton_iterations = max_Newton_iterations;
+solver_info.output_h_and_H_at_each_Newton_iteration
+	= (output_h_and_H_at_each_Newton_iteration != 0);
+solver_info.max_Delta_h_over_h           = max_Delta_h_over_h;
+solver_info.H_norm_for_convergence       = H_norm_for_convergence;
+solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
+solver_info.final_H_update_if_exit_x_H_small
+	= (final_H_update_if_exit_x_H_small != 0);
+
+
+//
+// set up the Cactus grid info
+//
+if (state.verbose_level >= verbose_level__algorithm_highlights)
+   then CCTK_VInfo(CCTK_THORNSTRING, "   setting up Cactus grid info");
+struct cactus_grid_info& cgi = state.cgi;
+cgi.GH = cctkGH;
+cgi.coord_origin[0] = cctk_origin_space[0];
+cgi.coord_origin[1] = cctk_origin_space[1];
+cgi.coord_origin[2] = cctk_origin_space[2];
+cgi.coord_delta[0] = cctk_delta_space[0];
+cgi.coord_delta[1] = cctk_delta_space[1];
+cgi.coord_delta[2] = cctk_delta_space[2];
+cgi.gridfn_dims[0] = cctk_lsh[0];
+cgi.gridfn_dims[1] = cctk_lsh[1];
+cgi.gridfn_dims[2] = cctk_lsh[2];
+// CCTK_VarDataPtr() returns a  void * , but we need a  const fp *;
+// since  static_cast<>  won't change const-ness, we need a 2-stage cast
+// for this:
+#define DATA_PTR_CAST(void_ptr_)	\
+	static_cast<const fp*>( const_cast<const void *>(void_ptr_) )
+cgi.g_dd_11_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxx"));
+cgi.g_dd_12_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxy"));
+cgi.g_dd_13_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxz"));
+cgi.g_dd_22_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyy"));
+cgi.g_dd_23_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyz"));
+cgi.g_dd_33_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gzz"));
+cgi.K_dd_11_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxx"));
+cgi.K_dd_12_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxy"));
+cgi.K_dd_13_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxz"));
+cgi.K_dd_22_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyy"));
+cgi.K_dd_23_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyz"));
+cgi.K_dd_33_data = DATA_PTR_CAST(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz"));
+
+
+//
+// set up the geometry parameters and the geometry interpolator
+//
+struct geometry_info& gi = state.gi;
+gi.geometry_method = decode_geometry_method(geometry_method);
+
+// parameters for geometry_method = "interpolate from Cactus grid"
+if (state.verbose_level >= verbose_level__algorithm_highlights)
+   then CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
+gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
+if (gi.operator_handle < 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"AHFinderDirect_setup(): couldn't find interpolator \"%s\"!",
+		   geometry_interpolator_name);		/*NOTREACHED*/
+gi.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars);
+if (gi.param_table_handle < 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"AHFinderDirect_setup(): bad geometry-interpolator parameter(s) \"%s\"!",
+		   geometry_interpolator_pars);		/*NOTREACHED*/
+
+// parameters for geometry_method = "Schwarzschild/EF"
+gi.geometry__Schwarzschild_EF__mass     = geometry__Schwarzschild_EF__mass;
+gi.geometry__Schwarzschild_EF__x_posn   = geometry__Schwarzschild_EF__x_posn;
+gi.geometry__Schwarzschild_EF__y_posn   = geometry__Schwarzschild_EF__y_posn;
+gi.geometry__Schwarzschild_EF__z_posn   = geometry__Schwarzschild_EF__z_posn;
+gi.geometry__Schwarzschild_EF__epsilon  = geometry__Schwarzschild_EF__epsilon;
+gi.geometry__Schwarzschild_EF__Delta_xyz= geometry__Schwarzschild_EF__Delta_xyz;
+
+
+//
+// create AH-specific info for each AH
+//
+
+// set up the interpatch interpolator
+if (state.verbose_level >= verbose_level__algorithm_highlights)
+   then CCTK_VInfo(CCTK_THORNSTRING, "   setting up interpatch interpolator");
+const int interp_handle = CCTK_InterpHandle(interpatch_interpolator_name);
+if (interp_handle < 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"AHFinderDirect_setup(): couldn't find interpolator \"%s\"!",
+		   interpatch_interpolator_name);		/*NOTREACHED*/
+const int interp_param_table_handle
+	= Util_TableCreateFromString(interpatch_interpolator_pars);
+if (interp_param_table_handle < 0)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"AHFinderDirect_setup(): bad interpatch-interpolator parameter(s) \"%s\"!",
+		   interpatch_interpolator_pars);		/*NOTREACHED*/
+
+	for (int hn = 0 ; hn < N_horizons ; ++hn)
+	{
+	struct AH_info* AH_ptr = new AH_info;
+
+	// create the patch system
+	patch_system& ps
+	   = *new patch_system
+		     (origin_x[hn], origin_y[hn], origin_z[hn],
+		      patch_system::type_of_name(patch_system_type),
+		      N_ghost_points, N_overlap_points, delta_drho_dsigma,
+		      gfns::nominal_min_gfn, gfns::nominal_max_gfn,
+		      gfns::ghosted_min_gfn, gfns::ghosted_max_gfn,
+		      interp_handle, interp_param_table_handle);
+	AH_ptr->ps_ptr = &ps;
+
+	// create the Jacobian matrix
+	AH_ptr->Jac_ptr = (state.method == method__horizon_function)
+			  ? NULL	// no Jacobian used for this case
+			  : &new_Jacobian(ps, Jac_info.Jacobian_storage_method);
+
+	AH_ptr->AH_found = false;
+	AH_ptr->centroid_x = 0.0;
+	AH_ptr->centroid_y = 0.0;
+	AH_ptr->centroid_z = 0.0;
+	AH_ptr->area = 0.0;
+	AH_ptr->mass = 0.0;
+
+	state.AH_info_ptrs.push_back(AH_ptr);
+	}
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function decodes the  method  parameter (string) into an
+// internal enum for future use.
+//
+namespace {
+enum method
+  decode_method(const char method_string[])
+{
+if	(STRING_EQUAL(method_string, "horizon function"))
+   then return method__horizon_function;
+else if (STRING_EQUAL(method_string, "Jacobian test"))
+   then return method__Jacobian_test;
+else if (STRING_EQUAL(method_string, "Jacobian test (NP only)"))
+   then return method__Jacobian_test_NP_only;
+else if (STRING_EQUAL(method_string, "Newton solve"))
+   then return method__Newton_solve;
+else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "decode_method(): unknown method_string=\"%s\"!",
+		   method_string);				/*NOTREACHED*/
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function decodes the  verbose_level  parameter (string) into an
+// internal enum for future use.
+//
+namespace {
+enum verbose_level
+  decode_verbose_level(const char verbose_level_string[])
+{
+if	(STRING_EQUAL(verbose_level_string, "physics highlights"))
+   then return verbose_level__physics_highlights;
+else if (STRING_EQUAL(verbose_level_string, "physics details"))
+   then return verbose_level__physics_details;
+else if (STRING_EQUAL(verbose_level_string, "algorithm highlights"))
+   then return verbose_level__algorithm_highlights;
+else if (STRING_EQUAL(verbose_level_string, "algorithm details"))
+   then return verbose_level__algorithm_details;
+else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"decode_verbose_level(): unknown verbose_level_string=\"%s\"!",
+		   verbose_level_string);			/*NOTREACHED*/
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function decodes the  Jacobian_method  parameter (string) into
+// an internal enum for future use.
+//
+namespace {
+enum Jacobian_method
+  decode_Jacobian_method(const char Jacobian_method_string[])
+{
+if	(STRING_EQUAL(Jacobian_method_string,
+		      "numerical perturbation"))
+   then return Jacobian_method__numerical_perturb;
+else if (STRING_EQUAL(Jacobian_method_string,
+		      "symbolic differentiation with finite diff d/dr"))
+   then return Jacobian_method__symbolic_diff_with_FD_dr;
+else if (STRING_EQUAL(Jacobian_method_string,
+		      "symbolic differentiation"))
+   then return Jacobian_method__symbolic_diff;
+else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"decode_Jacobian_method(): unknown Jacobian_method_string=\"%s\"!",
+		   Jacobian_method_string);			/*NOTREACHED*/
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function decodes the  geometry_method  parameter (string) into
+// an internal enum for future use.
+//
+namespace {
+enum geometry_method
+  decode_geometry_method(const char geometry_method_string[])
+{
+if	(STRING_EQUAL(geometry_method_string, "interpolate from Cactus grid"))
+   then return geometry__interpolate_from_Cactus_grid;
+else if (STRING_EQUAL(geometry_method_string, "Schwarzschild/EF"))
+   then return geometry__Schwarzschild_EF;
+else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"decode_geometry_method(): unknown geometry_method_string=\"%s\"!",
+		   geometry_method_string);			/*NOTREACHED*/
+}
+	  }
diff --git a/src/driver/state.cc b/src/driver/state.cc
new file mode 100644
index 0000000..52d62c8
--- /dev/null
+++ b/src/driver/state.cc
@@ -0,0 +1,44 @@
+// state.cc -- persistent state information (data) for this thorn
+// $Header$
+
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include <vector>
+
+#include "util_Table.h"
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../util/coords.hh"
+#include "../util/grid.hh"
+#include "../util/fd_grid.hh"
+#include "../util/patch.hh"
+#include "../util/patch_edge.hh"
+#include "../util/patch_interp.hh"
+#include "../util/ghost_zone.hh"
+#include "../util/patch_system.hh"
+
+#include "../elliptic/Jacobian.hh"
+
+#include "../gr/gfns.hh"
+#include "../gr/gr.hh"
+
+#include "AHFinderDirect.hh"
+
+//******************************************************************************
+
+//
+// The following data persists across Cactus scheduler calls.
+// The top-level drivers all share it.
+//
+struct state state;