move Jacobian.* to ../elliptic/

update other files to compute more of Jacobian coeffs


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

6 files changed, 177 insertions, 227 deletions

diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh
index 8bd0225..8e5fc8c 100644
--- a/src/gr/AHFinderDirect.hh
+++ b/src/gr/AHFinderDirect.hh
@@ -72,6 +72,7 @@ void horizon_function(patch_system& ps,
 		      bool Jacobian_flag);
 
 // horizon_Jacobian.cc
-class Jacobian;
-void horizon_Jacobian(const patch_system& ps,
+Jacobian& create_Jacobian(const patch_system& ps,
+			  const char Jacobian_type[]);
+void horizon_Jacobian(patch_system& ps,
 		      Jacobian& Jac);
diff --git a/src/gr/Jacobian.hh b/src/gr/Jacobian.hh
deleted file mode 100644
index e49187a..0000000
--- a/src/gr/Jacobian.hh
+++ /dev/null
@@ -1,34 +0,0 @@
-// Jacobian.hh -- data structures for the Jacobian matrix
-// $Id$
-
-//
-// prerequisites:
-//
-
-//
-// A Jacobian object stores the (a) Jacobian matrix of H with respect to h.
-// Jacobian is an abstract base class, from which we derive a separate
-// concrete class for each type of (storage scheme for the) Jacobian matrix.
-//
-// A row/column index of the Jacobian (denoted II/JJ) is a linearization
-// of the patch system's (patch,irho,isigma).
-//
-
-class	Jacobian
-	{
-public:
-	// basic meta-info
-	patch_system& my_patch_system();
-
-public:
-	// access a matrix element
-	virtual const fp& operator()(int II, int JJ) const;	// rvalue
-	virtual       fp& operator()(int II, int JJ) const;	// lvalue
-
-	// zero the whole matrix or a single row
-	virtual void zero();
-	virtual void zero_row(int II);
-
-	// add to an element: Jacobian(II,JJ) += delta
-	virtual void add_to_element(int II, int JJ, fp delta);
-	};
diff --git a/src/gr/cg.hh b/src/gr/cg.hh
index 3095a5a..2688a93 100644
--- a/src/gr/cg.hh
+++ b/src/gr/cg.hh
@@ -3,7 +3,9 @@
 
 //
 // This header file defines the "virtual machine" used by machine-generated
-// code.
+// code.  It is "dangerous" in that it #defines macros for all the gridfns,
+// which will likely produce syntax errors if they are used outside their
+// intentended environment.
 //
 // prerequisites
 //	gfn.hh
diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index 3982be2..0c04dfe 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -3,7 +3,7 @@
 //
 // <<<prototypes for functions local to this file>>>
 // AHFinderDirect_driver - top-level driver
-/// setup_ellipsoid_initial_guess - set up ellipsoid in h gridfn
+/// setup_Kerr_KerrSchild_initial_guess - set up Kerr/Kerr-Schild initial guess
 //
 
 #include <stdio.h>
@@ -32,6 +32,7 @@ using jtutil::error_exit;
 #include "../util/patch_interp.hh"
 #include "../util/ghost_zone.hh"
 #include "../util/patch_system.hh"
+#include "../util/Jacobian.hh"
 
 #include "gfn.hh"
 #include "AHFinderDirect.hh"
@@ -43,10 +44,9 @@ using jtutil::error_exit;
 //
 
 namespace {
-void setup_ellipsoid_initial_guess
-	(patch_system& ps,
-	 fp global_center_x, fp global_center_y, fp global_center_z,
-	 fp radius_x, fp radius_y, fp radius_z);
+void setup_Kerr_KerrSchild_horizon(patch_system& ps,
+				   fp mass, fp spin,
+				   fp xposn, fp yposn, fp zposn);
 	  };
 
 //******************************************************************************
@@ -113,42 +113,22 @@ 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];
-cgi.g_dd_11_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::gxx")
-					  );
-cgi.g_dd_12_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::gxy")
-					  );
-cgi.g_dd_13_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::gxz")
-					  );
-cgi.g_dd_22_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::gyy")
-					  );
-cgi.g_dd_23_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::gyz")
-					  );
-cgi.g_dd_33_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::gzz")
-					  );
-cgi.K_dd_11_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::kxx")
-					  );
-cgi.K_dd_12_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::kxy")
-					  );
-cgi.K_dd_13_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::kxz")
-					  );
-cgi.K_dd_22_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::kyy")
-					  );
-cgi.K_dd_23_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::kyz")
-					  );
-cgi.K_dd_33_data = static_cast<const fp *>(
-			CCTK_VarDataPtr(cctkGH, 0, "einstein::kzz")
-					  );
+// n.b. The  cgi.[gK]_dd_??_data  are actually  const fp *  pointers,
+//	since we won't modify the 3-D gridfn data!  But  static_cast<...>
+//      won't change const modifiers, so we just cast to  fp*  and let
+//	the assignment take care of the const part...
+cgi.g_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxx"));
+cgi.g_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxy"));
+cgi.g_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxz"));
+cgi.g_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gyy"));
+cgi.g_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gyz"));
+cgi.g_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gzz"));
+cgi.K_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxx"));
+cgi.K_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxy"));
+cgi.K_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxz"));
+cgi.K_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kyy"));
+cgi.K_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kyz"));
+cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kzz"));
 
 
 //
@@ -175,15 +155,14 @@ if	(STRING_EQUAL(initial_guess_method, "read from file"))
 			       initial_guess__read_from_file__file_name,
 			       false);		// no ghost zones
 	}
-else if (STRING_EQUAL(initial_guess_method, "ellipsoid"))
+else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild"))
    then {
-	setup_ellipsoid_initial_guess(ps,
-				      initial_guess__ellipsoid__center_global_x,
-				      initial_guess__ellipsoid__center_global_y,
-				      initial_guess__ellipsoid__center_global_z,
-				      initial_guess__ellipsoid__radius_x,
-				      initial_guess__ellipsoid__radius_y,
-				      initial_guess__ellipsoid__radius_z);
+	setup_Kerr_KerrSchild_horizon(ps,
+				      initial_guess__Kerr_KerrSchild__mass,
+				      initial_guess__Kerr_KerrSchild__spin,
+				      initial_guess__Kerr_KerrSchild__xposn,
+				      initial_guess__Kerr_KerrSchild__yposn,
+				      initial_guess__Kerr_KerrSchild__zposn);
 	ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
 					 true, ghosted_gfns::gfn__h,
 					 "h.dat",
@@ -204,6 +183,12 @@ if      (STRING_EQUAL(method, "horizon"))
 				 true, ghosted_gfns::gfn__h,
 				 "H.dat");
 	}
+if      (STRING_EQUAL(method, "horizon Jacobian"))
+   then {
+	Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
+	horizon_function(ps, cgi, gii, true);
+	horizon_Jacobian(ps, Jac);
+	}
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "unknown method=\"%s\"!",
 		   method);					/*NOTREACHED*/
@@ -212,38 +197,30 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 //******************************************************************************
 
 //
-// 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)
+// This function sets up the horizon of a Kerr black hole in Kerr-Schild
+// coordinates, on the nominal grid, in the  h  gridfn.
 //
-// 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
+// Kerr-Schild coordinates are described in MTW Exercise 33.8, page 903,
+// and the horizon is worked out on page 13 of my AHFinderDirect notes.
 //
-// to solve this, we introduce intermediate variables
-//	AU = A - U
-//	BV = B - V
-//	CW = C - W
+// Arguments:
+// (mass,spin) = Describe the Kerr black hole.
+// [xyz]_posn = The position of the Kerr black hole.
 //
 namespace {
-void setup_ellipsoid_initial_guess
-	(patch_system& ps,
-	 fp global_center_x, fp global_center_y, fp global_center_z,
-	 fp radius_x, fp radius_y, fp radius_z)
+void setup_Kerr_KerrSchild_horizon(patch_system& ps,
+				   fp mass, fp spin,
+				   fp xposn, fp yposn, fp zposn)
 {
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "   h = ellipsoid: global_center=(%g,%g,%g)",
-	   global_center_x, global_center_y, global_center_z);
+	   "   setting up Kerr/Kerr-Schild horizon");
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "                  radius=(%g,%g,%g)",
-	   radius_x, radius_y, radius_z);
+	   "      mass=%g, spin=%g, posn=(%g,%g,%g)",
+	   double(mass), double(spin),
+	   double(xposn), double(yposn), double(zposn));
+
+// horizon in Kerr-Schild is coordinate sphere $r = (1 + \sqrt{1-a^2}) m$
+const fp r = (1.0 + sqrt(1.0 - spin*spin)) * mass;
 
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
 	{
@@ -255,44 +232,19 @@ CCTK_VInfo(CCTK_THORNSTRING,
 		     isigma <= p.max_isigma() ;
 		     ++isigma)
 		{
-		const fp rho = p.rho_of_irho(irho);
+		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 = global_center_x - ps.origin_x();
-		const fp BV = global_center_y - ps.origin_y();
-		const fp CW = global_center_z - ps.origin_z();
-		const fp a = radius_x;
-		const fp b = radius_y;
-		const fp c = radius_z;
-
-		// 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"
-"   expected exactly one r>0 solution, 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*/
+		fp theta, phi;
+		p.theta_phi_of_rho_sigma(rho,sigma, theta,phi);
+		fp Kerr_x, Kerr_y, Kerr_z;
+		p.xyz_of_r_rho_sigma(r,rho,sigma, Kerr_x,Kerr_y,Kerr_z);
 
-		// r = horizon radius at this grid point
-		p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r;
+		const fp KS_x = Kerr_x - spin*sin(theta)*sin(phi);
+		const fp KS_y = Kerr_y + spin*sin(theta)*cos(phi);
+		const fp KS_z = Kerr_z;
+		p.ghosted_gridfn(ghosted_gfns::gfn__h,irho,isigma)
+			= jtutil::hypot3(KS_x, KS_y, KS_z);
 		}
 		}
 	}
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc
index 23390ab..c7d397e 100644
--- a/src/gr/horizon_Jacobian.cc
+++ b/src/gr/horizon_Jacobian.cc
@@ -2,7 +2,9 @@
 // $Id$
 //
 // <<<prototypes for functions local to this file>>>
-// horizon_Jacobian - top-level driver
+// create_Jacobian - create a Jacobian matrix of the appropriate type
+// horizon_Jacobian - top-level driver to compute the Jacobian
+/// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian
 //
 
 #include <stdio.h>
@@ -30,9 +32,9 @@ using jtutil::error_exit;
 #include "../util/patch_interp.hh"
 #include "../util/ghost_zone.hh"
 #include "../util/patch_system.hh"
+#include "../util/Jacobian.hh"
 
 #include "gfn.hh"
-#include "Jacobian.hh"
 #include "AHFinderDirect.hh"
 
 //******************************************************************************
@@ -42,14 +44,29 @@ using jtutil::error_exit;
 //
 
 namespace {
-void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp,
-				    int x_irho, int x_isigma,
-				    int xm_irho, int xm_isigma,
-				    fp mol);
+void add_ghost_zone_Jacobian(const patch_system& ps,
+			     const patch& xp, const ghost_zone& xmgz,
+			     int x_II,
+			     int xm_irho, int xm_isigma,
+			     fp mol, Jacobian& Jac);
 	  }
 
 //******************************************************************************
-//******************************************************************************
+
+//
+// This function decodes the  Jacobian_type  parameter and creates the
+// appropriate type of Jacobian matrix object.
+//
+Jacobian& create_Jacobian(const patch_system& ps,
+			  const char Jacobian_type[])
+{
+if	(STRING_EQUAL(Jacobian_type, "dense matrix"))
+   then return *new dense_Jacobian(ps);
+else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "unknown Jacobian_type=\"%s\"!",
+		   Jacobian_type);				/*NOTREACHED*/
+}
+
 //******************************************************************************
 
 //
@@ -65,8 +82,8 @@ void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp,
 //
 // Outputs:
 //	The Jacobian matrix is stored in the Jacobian object Jac.
-void horizon_Jacobian(const patch_system& ps,
-		      Jacobian& Jac);
+//
+void horizon_Jacobian(patch_system& ps, Jacobian& Jac)
 {
 CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian");
 
@@ -80,8 +97,8 @@ Jac.zero();
 	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)
+	x_isigma <= xp.max_isigma() ;
+	++x_isigma)
 	{
 	//
 	// compute the Jacobian row for this grid point, i.e.
@@ -96,17 +113,25 @@ Jac.zero();
 	//		Newton's method really wants the latter
 	//
 
+	// 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(gfn__partial_H_wrt_partial_d_h_1, x_irho,x_isigma);
+	   = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_1,
+		       x_irho, x_isigma);
 	const fp Jacobian_coeff_sigma
-	   = xp.gridfn(gfn__partial_H_wrt_partial_d_h_2, x_irho,x_isigma);
+	   = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_2,
+		       x_irho, x_isigma);
 	const fp Jacobian_coeff_rho_rho
-	   = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_11, x_irho,x_isigma);
+	   = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_11,
+		       x_irho, x_isigma);
 	const fp Jacobian_coeff_rho_sigma
-	   = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_12, x_irho,x_isigma);
+	   = xp.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_12,
+		       x_irho, x_isigma);
 	const fp Jacobian_coeff_sigma_sigma
-	   = xp.gridfn(gfn__partial_H_wrt_partial_dd_h_22, x_irho,x_isigma);
+	   = xp.gridfn(nominal_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() ;
@@ -115,13 +140,16 @@ Jac.zero();
 	    {
 	    const int xm_irho = x_irho + m_irho;
 	    const fp Jac_rho     = Jacobian_coeff_rho
-				   * xp.partial_rho_coeff(m_rho);
+				   * xp.partial_rho_coeff(m_irho);
 	    const fp Jac_rho_rho = Jacobian_coeff_rho_rho
-				   * xp.partial_rho_rho_coeff(m_rho);
-	    add_molecule_point_to_Jacobian(ps, xp,
-					   x_irho, x_isigma,
-					   xm_irho, x_isigma,
-					   Jac_rho + Jac_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, xp, xp.minmax_rho_ghost_zone(m_irho < 0),
+			 II, xm_irho, x_isigma,
+			 Jac_sum, Jac);
 	    }
 
 	// partial_sigma, partial_sigma_sigma
@@ -131,13 +159,16 @@ Jac.zero();
 	    {
 	    const int xm_isigma = x_isigma + m_isigma;
 	    const fp Jac_sigma       = Jacobian_coeff_sigma
-				       * xp.partial_sigma_coeff(m_sigma);
+				       * xp.partial_sigma_coeff(m_isigma);
 	    const fp Jac_sigma_sigma = Jacobian_coeff_sigma_sigma
-				       * xp.partial_sigma_sigma_coeff(m_sigma);
-	    add_molecule_point_to_Jacobian(ps, xp,
-					   x_irho, x_isigma,
-					   x_irho, xm_isigma,
-					   Jac_sigma + Jac_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, xp, xp.minmax_sigma_ghost_zone(m_isigma < 0),
+			 II, x_irho, xm_isigma,
+			 Jac_sum, Jac);
 	    }
 
 	// partial_rho_sigma
@@ -152,12 +183,16 @@ Jac.zero();
 	    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_rho, m_sigma);
-	    add_molecule_point_to_Jacobian(ps, xp,
-					   x_irho, x_isigma,
-					   xm_irho, xm_isigma,
-					   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 add_ghost_zone_Jacobian
+			(ps, xp, 
+		 xp.corner_ghost_zone_containing_point(m_irho < 0, m_isigma < 0,
+						       xm_irho, xm_isigma),
+			 II, xm_irho, xm_isigma,
+			 Jac_rho_sigma, Jac);
 	    }
 	    }
 
@@ -169,59 +204,52 @@ Jac.zero();
 //******************************************************************************
 
 //
-// This function adds all the elements for a single molecule point to
-// a Jacobian matrix, including any contributions from other patches if
-// the specified point lies in a ghost zone.
+// This function adds the ghost-zone Jacobian dependency contributions
+// for a single ghost-zone point, to a Jacobian matrix.
 //
 // Arguments:
 // ps = The patch system.
 // xp = The patch containing the center point of the molecule.
-// x_(irho,isigma) = The coordinates of the center point of the molecule.
-// xm_(irho,isigma) = The coordinates of the x+m 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.
 // mol = The molecule coefficient.
+// Jac = The Jacobian matrix.
 //
 namespace {
-void add_molecule_point_to_Jacobian(const patch_system& ps, const patch& xp,
-				    int x_irho, int x_isigma,
-				    int xm_irho, int xm_isigma,
-				    fp mol)
+void add_ghost_zone_Jacobian(const patch_system& ps,
+			     const patch& xp, const ghost_zone& xmgz,
+			     int x_II,
+			     int xm_irho, int xm_isigma,
+			     fp mol, Jacobian& Jac)
 {
-// Jacobian row
-const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho, x_isigma);
-
-if (xp.is_in_nominal_grid(xm_irho, xm_isigma))
-   then {
-	const int JJ = ps.gpn_of_patch_irho_isigma(xp, xm_irho, xm_isigma);
-	Jac.add_to_element(II,JJ, Jac_coeff);
-	}
-   else {
-	const ghost_zone& xmgz
-		= xp.ghost_zone_containing_point(xm_irho, 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);
-
-	// on what other points ym does this molecule point xm depend
-	// via the patch_system::synchronize() operation?
-	const patch&      ymp = ps.synchronize_Jacobian_y_patch(xmgz);
-	const patch_edge& yme = ps.synchronize_Jacobian_y_edge (xmgz);
-	const int min_ym_ipar_m = ps.synchronize_Jacobian_min_y_ipar_m(xmgz);
-	const int max_ym_ipar_m = ps.synchronize_Jacobian_max_y_ipar_m(xmgz);
-	const int ym_iperp = ps.synchronize_Jacobian_y_iperp(xmgz, xm_iperp);
-	const int ym_ipar_posn
-	   = ps.synchronize_Jacobian_y_ipar_posn(xmgz, xm_iperp, xm_ipar);
-	    for (int ym_ipar_m = min_ym_ipar_m ;
-		 ym_ipar_m <= max_ym_ipar_m ;
-		 ++ym_ipar_m)
-	    {
-	    const int ym_ipar = ym_ipar_posn + ym_ipar_m;
-	    const int ym_irho   = yme.  irho_of_iperp_ipar(ym_iperp,ym_ipar);
-	    const int ym_isigma = yme.isigma_of_iperp_ipar(ym_iperp,ym_ipar);
-	    const int JJ = ps.gpn_of_patch_irho_isigma(ymp, ym_irho, ym_isigma);
-	    const fp sync_Jac = ps.synchronize_Jacobian(xmgz, xm_iperp, xm_ipar,
-							ym_ipar_m);
-	    Jac.add_to_element(II,JJ, mol*sync_Jac);
-	    }
+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);
+
+// on what other points ym does this molecule point xm depend
+// via the patch_system::synchronize() operation?
+const patch&      ymp = ps.synchronize_Jacobian_y_patch(xmgz);
+const patch_edge& yme = ps.synchronize_Jacobian_y_edge (xmgz);
+const int min_ym_ipar_m = ps.synchronize_Jacobian_min_y_ipar_m(xmgz);
+const int max_ym_ipar_m = ps.synchronize_Jacobian_max_y_ipar_m(xmgz);
+const int ym_iperp = ps.synchronize_Jacobian_y_iperp(xmgz, xm_iperp);
+const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn
+				(xmgz, xm_iperp, xm_ipar);
+
+// add the Jacobian contributions from the ym points
+	for (int ym_ipar_m = min_ym_ipar_m ;
+	     ym_ipar_m <= max_ym_ipar_m ;
+	     ++ym_ipar_m)
+	{
+	const int ym_ipar = ym_ipar_posn + ym_ipar_m;
+	const int ym_irho   = yme.  irho_of_iperp_ipar(ym_iperp,ym_ipar);
+	const int ym_isigma = yme.isigma_of_iperp_ipar(ym_iperp,ym_ipar);
+	const int JJ = ps.gpn_of_patch_irho_isigma(ymp, ym_irho, ym_isigma);
+	const fp sync_Jac = ps.synchronize_Jacobian(xmgz, xm_iperp, xm_ipar,
+						    ym_ipar_m);
+	Jac(x_II,JJ) += mol*sync_Jac;
 	}
 }
 	  }
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 50e9254..36d5746 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -33,6 +33,7 @@ using jtutil::error_exit;
 #include "../util/patch_interp.hh"
 #include "../util/ghost_zone.hh"
 #include "../util/patch_system.hh"
+#include "../util/Jacobian.hh"
 
 #include "gfn.hh"
 #include "AHFinderDirect.hh"
@@ -98,8 +99,8 @@ void horizon_function(patch_system& ps,
 {
 CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");
 
-// fill in values of ghosted gridfns in ghost zones
-ps.synchronize(ghosted_gfns::gfn__h, ghosted_gfns::gfn__h);
+// fill in values of all ghosted gridfns in ghost zones
+ps.synchronize();
 
 // set up xyz positions of grid points
 setup_xyz_posns(ps);