start adding Jadcobian support --> not working yet

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

12 files changed, 496 insertions, 77 deletions

diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh
index 59dacb8..8bd0225 100644
--- a/src/gr/AHFinderDirect.hh
+++ b/src/gr/AHFinderDirect.hh
@@ -21,10 +21,11 @@ struct	geometry_interpolator_info
 
 //
 // This structure holds all the information we need about the Cactus grid
-// and gridfns outside the top-level driver.  At present we use
-//   CCTK_InterpLocalUniform()
-// to access the Cactus gridfns, and much of this information is specific
-// to that API.
+// and gridfns outside the top-level driver.  At present we use the
+//  CCTK_InterpLocalUniform()  local interpolator to access the Cactus
+// gridfns.  Much of the information in this structure is specific to
+// that API, and (alas) will need changing when we eventually switch to
+// a "global" multiprocessor/distributed interpolator.
 //
 struct	cactus_grid_info
 	{
@@ -67,4 +68,10 @@ extern "C"
 // horizon_function.cc
 void horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
-		      const struct geometry_interpolator_info& gii);
+		      const struct geometry_interpolator_info& gii,
+		      bool Jacobian_flag);
+
+// horizon_Jacobian.cc
+class Jacobian;
+void horizon_Jacobian(const patch_system& ps,
+		      Jacobian& Jac);
diff --git a/src/gr/Jacobian.hh b/src/gr/Jacobian.hh
new file mode 100644
index 0000000..626d080
--- /dev/null
+++ b/src/gr/Jacobian.hh
@@ -0,0 +1,30 @@
+// 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:
+	// 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/auxiliary.maple b/src/gr/auxiliary.maple
index fecfc48..99069e7 100644
--- a/src/gr/auxiliary.maple
+++ b/src/gr/auxiliary.maple
@@ -14,7 +14,7 @@
 # algebraic functions of the basic ADM fields {g_dd,K_dd}, and optionally
 # also generates C code for them (the auxiliar variables).
 #
-# Code files are written to the directory cg/ .
+# Code files are written to the directory ../gr.cg/ .
 #
 # Inputs (Maple):
 #	N
@@ -73,7 +73,7 @@ assert_fnd_exists(g_uu, fnd);
 g_uu__fnd := linalg[inverse](g_dd);
 
 if (cg_flag)
-   then codegen2(g_uu__fnd, 'g_uu', "cg/inverse_metric.c");
+   then codegen2(g_uu__fnd, 'g_uu', "../gr.cg/inverse_metric.c");
 fi;
 
 NULL;
@@ -132,7 +132,7 @@ K__fnd := simplify(
 if (cg_flag)
    then codegen2([K__fnd, K_uu__fnd],
 		 ['K', 'K_uu'],
-		 "cg/extrinsic_curvature_trace_raise.c");
+		 "../gr.cg/extrinsic_curvature_trace_raise.c");
 fi;
 
 NULL;
diff --git a/src/gr/cg.hh b/src/gr/cg.hh
index b6174eb..3095a5a 100644
--- a/src/gr/cg.hh
+++ b/src/gr/cg.hh
@@ -9,6 +9,8 @@
 //	gfn.hh
 //
 
+//******************************************************************************
+
 //
 // The machine-generated code uses the following variables:
 //	patch& p		// current patch
@@ -37,6 +39,8 @@
 #define PARTIAL_SIGMA_SIGMA(ghosted_gridfn_name)	\
   p.partial_sigma_sigma(ghosted_gfns::gfn__ ## ghosted_gridfn_name, irho,isigma)
 
+//******************************************************************************
+
 //
 // ghosted-grid gridfns
 // n.b. since we're always evaluating on the surface, r == h
@@ -44,6 +48,8 @@
 #define h	p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma)
 #define r	h
 
+//******************************************************************************
+
 //
 // nominal-grid gridfns
 //
@@ -81,9 +87,22 @@
 
 #define H	p.gridfn(nominal_gfns::gfn__H, irho,isigma)
 
+#define partial_H_wrt_partial_d_h_1	\
+	p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_1, irho,isigma)
+#define partial_H_wrt_partial_d_h_2	\
+	p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_d_h_2, irho,isigma)
+#define partial_H_wrt_partial_dd_h_11	\
+	p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_11, irho,isigma)
+#define partial_H_wrt_partial_dd_h_12	\
+	p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_12, irho,isigma)
+#define partial_H_wrt_partial_dd_h_22	\
+	p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_22, irho,isigma)
+
+//******************************************************************************
+
 //
-// pseudo-gridfns
-// (used only at a single grid point, not actually stored as gridfns)
+// pseudo-gridfns, i.e. temporaries used only at a single grid point
+// (these are not actually stored as gridfns)
 //
 fp g_uu_11;
 fp g_uu_12;
@@ -121,3 +140,8 @@ fp partial_d_g_uu_313;
 fp partial_d_g_uu_322;
 fp partial_d_g_uu_323;
 fp partial_d_g_uu_333;
+
+fp HA;
+fp HB;
+fp HC;
+fp HD;
diff --git a/src/gr/doit.maple b/src/gr/doit.maple
index 8d5edf1..86cdebe 100644
--- a/src/gr/doit.maple
+++ b/src/gr/doit.maple
@@ -6,10 +6,10 @@ read "setup_gr_gfas.mm";
 setup_gr_gfas();
 
 read "auxiliary.mm";
-read "metric_derivs.mm";
+read "curvature.mm";
 read "horizon.mm";
 
 `saveit/level` := 10;
 auxiliary(true);
-metric_derivs(true);
+curvature(true);
 horizon(true);
diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index fe64dce..e95c783 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -200,7 +200,7 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 //
 if      (STRING_EQUAL(method, "horizon"))
    then {
-	horizon_function(ps, cgi, gii);
+	horizon_function(ps, cgi, gii, false);
 	ps.print_gridfn_with_xyz(nominal_gfns::gfn__H,
 				 true, ghosted_gfns::gfn__h,
 				 "H.dat");
diff --git a/src/gr/gfn.hh b/src/gr/gfn.hh
index 16af52f..24c56fe 100644
--- a/src/gr/gfn.hh
+++ b/src/gr/gfn.hh
@@ -48,6 +48,11 @@ enum	{
 	gfn__partial_d_g_dd_323,
 	gfn__partial_d_g_dd_333,
 	gfn__H,
-	max_gfn = gfn__H
+	gfn__partial_H_wrt_partial_d_h_1,
+	gfn__partial_H_wrt_partial_d_h_2,
+	gfn__partial_H_wrt_partial_dd_h_11,
+	gfn__partial_H_wrt_partial_dd_h_12,
+	gfn__partial_H_wrt_partial_dd_h_22,
+	max_gfn = gfn__partial_H_wrt_partial_dd_h_22	// no comma
 	};
 	  };
diff --git a/src/gr/horizon.maple b/src/gr/horizon.maple
index 6950e56..0707af4 100644
--- a/src/gr/horizon.maple
+++ b/src/gr/horizon.maple
@@ -4,7 +4,8 @@
 # horizon - overall driver
 ## non_unit_normal - compute s_d
 ## non_unit_normal_deriv - compute partial_d_s_d
-## horizon_function - compute HA,HB,HC,HD,H = fn(s_d)
+## horizon_function - compute H[ABCD],H = fn(s_d__fnd) = fn(h)
+## horizon_Jacobian - compute Jacobian coeffs for H[ABCD], H
 #
 
 ################################################################################
@@ -27,6 +28,7 @@ proc(cg_flag::boolean)
 non_unit_normal();
 non_unit_normal_deriv();
 horizon_function(cg_flag);
+horiozn_Jacobian(cg_flag);
 
 NULL;
 end proc;
@@ -34,9 +36,14 @@ end proc;
 ################################################################################
 
 #
-# This function computes the non-unit outward-pointing normal (covector)
-# to the horizon, s_d = n_u(h), using the equation on p.7.1 of my apparent
-# horizon finding notes.  This function does *NOT* generate any C code.
+# This function computes the xyz components of the non-unit
+# outward-pointing normal (covector) to the horizon, s_d = s_d(h),
+# using the equation on p.7.1 of my apparent horizon finding notes.
+#
+# This function does *NOT* generate any C code.
+#
+# Outputs:
+#   s_d = s_d__fnd( x_xyz, X_ud, Diff(h,y_rs) )
 #
 non_unit_normal :=
 proc()
@@ -64,10 +71,18 @@ end proc;
 ################################################################################
 
 #
-# This function computes the 1st derivatives of the non-unit outward-pointing
+# This function computes the xyz derivatives of the non-unit outward-pointing
 # normal (covector) to the horizon, using the equation on p.7.2 of my
-# apparent horizon finding notes.  This function does *NOT* generate any
-# C code.
+# apparent horizon finding notes.
+#
+# This function does *NOT* generate any C code.
+#
+# Inputs:
+#   s_d = s_d__fnd( x_xyz, X_ud, Diff(h,y_rs) )
+#
+# Outputs:
+#   partial_d_s_d = partial_d_s_d__fnd( x_xyz, X_ud, X_udd,
+#					Diff(h,y_rs), Diff(h,y_rs,y_rs) )
 #
 non_unit_normal_deriv :=
 proc()
@@ -113,9 +128,32 @@ end proc;
 # as a function of 1st and 2nd angular partial derivatives of the
 # apparent horizon radius $h$, in the Maple gridfn array H__fnd,
 # using equation (15) in my 1996 apparent horizon finding paper,
-# and using partial_d_g_uu[k,i,j] instead of Diff(g_uu[i,j], x_xyz[k]).
+# but using partial_d_g_uu[k,i,j] instead of Diff(g_uu[i,j], x_xyz[k]).
+#
+# These equation gives H[ABCD] and H as functions of s_d, but here we
+# use s_d__fnd, which we assume is already given in terms of angular
+# derivatives of h.  The result is that we compute H[ABCD] and H directly
+# in terms of 1st and 2nd angular derivatives of h, without s_d ever
+# appearing in our final results.
+#
 # This function also optionally generates C code for this computation.
 #
+# Inputs:
+#   s_d = s_d__fnd( x_xyz, X_ud, Diff(h,y_rs) )
+#   partial_d_s_d = partial_d_s_d__fnd( x_xyz, X_ud, X_udd,
+#					Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+#
+# Outputs:
+#   HA = HA__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
+#		  Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+#   HB = HB__fnd( x_xyz, X_ud, X_udd,
+#		  g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
+#		  Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+#   HC = HC__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
+#   HC = HC__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
+#   ---------------------------------------------------------------
+#   H = H__fnd( HA__fnd, HB__fnd, HC__fnd, HD__fnd )
+#
 horizon_function :=
 proc(cg_flag::boolean)
 global
@@ -163,7 +201,109 @@ HD__fnd := msum('g_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
 H__fnd := HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K;
 
 if (cg_flag)
-   then codegen2(H__fnd, 'H', "cg/horizon_function.c");
+   then codegen2([ HA__fnd, HB__fnd, HC__fnd, HD__fnd, H__fnd],
+		 ['HA',    'HB',    'HC',    'HD',    'H'    ],
+		 "../gr.cg/horizon_function.c");
+fi;
+
+NULL;
+end proc;
+
+################################################################################
+
+#
+# This function computes the Jacobian coefficients for the LHS of the
+# apparent horizon equation with respect to 1st and 2nd angular partial
+# derivatives of the apparent horizon radius $h$.  These coefficients
+# are (or at least should be :) the same as those in equation (A1) in
+# my 1996 apparent horizon finding paper, bue here we compute them in
+# a different manner: we have Maple directly differentiate H__fnd with
+# respect to Diff(h,y_rs[u]) and Diff(h,y_rs[u],y_rs[v]).  We use the
+# Maple frontend() function to do this.
+#
+# This function also optionally generates C code for the Jacobian
+# coefficients.
+#
+# Inputs:
+#   HA = HA__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
+#		  Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+#   HB = HB__fnd( x_xyz, X_ud, X_udd,
+#		  g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
+#		  Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+#   HC = HC__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
+#   HC = HC__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
+#   H = H__fnd( HA__fnd, HB__fnd, HC__fnd, HD__fnd )
+#
+# Outputs:
+#   partial_HA_wrt_partial_d_h		partial_HA_wrt_partial_dd_h
+#   partial_HB_wrt_partial_d_h		partial_HB_wrt_partial_dd_h
+#   partial_HC_wrt_partial_d_h
+#   partial_HD_wrt_partial_d_h
+#   ---------------------------------------------------------------
+#   partial_H_wrt_partial_d_h		partial_H_wrt_partial_dd_h
+#
+horizon_Jacobian :=
+proc(cg_flag::boolean)
+global
+  @include "../maple/coords.minc",
+  @include "../maple/gfa.minc",
+  @include "../gr/gr_gfas.minc";
+local u,v,
+      temp;
+
+printf("%a...\n", procname);
+
+assert_fnd_exists(g_uu);
+assert_fnd_exists(K_uu);
+assert_fnd_exists(partial_d_ln_sqrt_g);
+assert_fnd_exists(partial_d_g_uu);
+assert_fnd_exists(HA)
+assert_fnd_exists(HB)
+assert_fnd_exists(HC)
+assert_fnd_exists(HD)
+
+# Jacobian coefficients of H[ABCD] and H wrt Diff(h,y_rs[u])
+	for u from 1 to N_ang
+	do
+	partial_HA_wrt_partial_d_h__fnd[u]
+		:= frontend('diff', [HA__fnd, Diff(h,y_rs[u])]);
+	partial_HB_wrt_partial_d_h__fnd[u]
+		:= frontend('diff', [HB__fnd, Diff(h,y_rs[u])]);
+	partial_HC_wrt_partial_d_h__fnd[u]
+		:= frontend('diff', [HC__fnd, Diff(h,y_rs[u])]);
+	partial_HD_wrt_partial_d_h__fnd[u]
+		:= frontend('diff', [HD__fnd, Diff(h,y_rs[u])]);
+
+	# equation (A1a) in my 1996 apparent horizon finding paper
+	temp := (3/2)*HA/HD^(5/2) + (1/2)*HB/HD^(3/2) + HC/HD^2;
+	partial_H_wrt_partial_d_h__fnd[u]
+		:= + partial_HA_wrt_partial_d_h__fnd[u] / HD^(3/2)
+		   + partial_HB_wrt_partial_d_h__fnd[u] / HD^(1/2)
+		   + partial_HC_wrt_partial_d_h__fnd[u] / HD
+		   - partial_HD_wrt_partial_d_h__fnd[u] * temp;
+	end do;
+
+# Jacobian coefficients of H[AB] and H wrt Diff(h,y_rs[u],y_rs[v])
+	for u from 1 to N_ang
+	do
+	for v from u to N_ang
+	do
+	partial_HA_wrt_partial_dd_h__fnd[u,v]
+		:= frontend('diff', [HA__fnd, Diff(h,y_rs[u],y_rs[v])]);
+	partial_HB_wrt_partial_dd_h__fnd[u,v]
+		:= frontend('diff', [HB__fnd, Diff(h,y_rs[u],y_rs[v])]);
+
+	# equation (A1b) in my 1996 apparent horizon finding paper
+	partial_H_wrt_partial_dd_h__fnd[u,v]
+		:= + partial_HA_wrt_partial_dd_h__fnd[u,v] / HD^(3/2)
+		   + partial_HB_wrt_partial_dd_h__fnd[u,v] / HD^(1/2)
+	end do;
+
+if (cg_flag)
+   then codegen2([partial_H_wrt_partial_d_h__fnd,
+		  partial_H_wrt_partial_dd_h__fnd],
+		 ['partial_H_wrt_partial_d_h', 'partial_H_wrt_partial_dd_h'],
+		 "../gr.cg/horizon_Jacobian.c");
 fi;
 
 NULL;
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc
new file mode 100644
index 0000000..40fada6
--- /dev/null
+++ b/src/gr/horizon_Jacobian.cc
@@ -0,0 +1,161 @@
+// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h)
+// $Id$
+//
+// <<<prototypes for functions local to this file>>>
+// horizon_Jacobian - 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 "jt/stdc.h"
+#include "jt/util.hh"
+#include "jt/array.hh"
+#include "jt/cpm_map.hh"
+#include "jt/linear_map.hh"
+using jtutil::error_exit;
+
+#include "../config.hh"
+
+#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 "gfn.hh"
+#include "Jacobian.hh"
+#include "AHFinderDirect.hh"
+
+//******************************************************************************
+
+//
+// ***** prototypes for functions local to this file *****
+//
+
+namespace {
+	  };
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function computes the Jacobian matrix of the LHS function H(h)
+// 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
+//	partial_H_wrt_partial_d_h,	# Jacobian coefficients
+//	partial_H_wrt_partial_dd_h,
+//
+// Outputs:
+//	The Jacobian matrix is stored in the Jacobian object Jac.
+void horizon_Jacobian(const patch_system& ps,
+		      Jacobian& Jac);
+{
+CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian");
+
+Jac.zero(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)
+	{
+	//
+	// compute the Jacobian row for this grid point, i.e.
+	//	  partial H(this point x, Jacobian row II)
+	//	---------------------------------------------
+	//	partial h(other points y, Jacobian column JJ)
+	//
+	// FIXME FIXME: we still have to take into account the
+	//		position-dependence of the coefficients,
+	//		cf the difference between J[3H(h)] and J[2H(h)];
+	//		here we're sort of computing the former, but
+	//		Newton's method really wants the latter
+	//
+
+	// Jacobian row index for this point
+	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);
+	const fp Jacobian_coeff_sigma
+	   = xp.gridfn(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);
+	const fp Jacobian_coeff_rho_sigma
+	   = xp.gridfn(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);
+
+	// partial_rho, partial_rho_rho
+	    for (int m_irho = xp.molecule_min_m() ;
+		 m_irho <= xp.molecule_max_m() ;
+		 ++m_irho)
+	    {
+	    const fp Jac_rho = Jacobian_coeff_rho * xp.partial_rho_coeff(m_rho);
+	    const fp Jac_rho_rho = Jacobian_coeff_rho_rho
+				   * xp.partial_rho_rho_coeff(m_rho);
+	    const fp Jac_sum = Jac_rho + Jac_rho_rho;
+	    if (xp.is_in_nominal_grid(x_irho+m_irho, x_isigma))
+	       then {
+		    const int JJ
+			= ps.gpn_of_patch_irho_isigma(xp, x_irho+m_irho,
+							  x_isigma);
+		    Jac.add_to_element(II,JJ, Jac_sum);
+		    }
+	       else {
+		    const patch_edge& xe = (m_irho < 0)
+					   ? xp.min_rho_patch_edge()
+					   : xp.max_rho_patch_edge()
+		    const ghost_zone& gz = xp.ghost_zone_on_edge(xe);
+		    const int x_iperp = xe.iperp_of_irho_isigma(x_irho,
+								x_isigma);
+		    const int x_ipar = xe.ipar_of_irho_isigma(x_irho, x_isigma);
+		    gz.compute_Jacobian(gfn__h, gfn__h);
+		    const patch&      yp = gz.Jacobian_y_patch();
+		    const patch_edge& ye = gz.Jacobian_y_edge();
+		    const int min_y_ipar_m = gz.Jacobian_min_y_ipar_m();
+		    const int max_y_ipar_m = gz.Jacobian_max_y_ipar_m();
+		    const int y_iperp = gz.Jacobian_y_iperp(x_iperp);
+		    const int y_ipar_posn
+			= gz.Jacobian_y_ipar_posn(x_iperp, x_ipar);
+			for (int y_ipar_m = min_y_ipar_m ;
+			     y_ipar_m <= max_y_ipar_m ;
+			     ++y_ipar_m)
+			{
+			const int y_ipar = y_ipar_posn + y_ipar_m;
+			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);
+			const fp gz_Jac = gz.Jacobian(x_iperp,x_ipar, y_ipar_m);
+			Jac.add_to_element(II,JJ, Jac_sum*gz_Jac);
+			}
+			}
+		}
+	}
+	}
+     }
+}
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 6ea0a8d..745760b 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -57,7 +57,9 @@ void compute_H(patch_system& ps);
 //******************************************************************************
 
 //
-// This function computes the LHS function H(h).
+// This function computes the LHS function H(h), and optionally also
+// its Jacobian coefficients (from which the Jacobian matrix may be
+// computed later).
 //
 // Inputs (angular gridfns, on ghosted grid):
 // ... defined on ghosted grid
@@ -86,9 +88,14 @@ void compute_H(patch_system& ps);
 //	partial_d_g_dd_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g_{ij}$
 //	H				# $H = H(h)$
 //
+// Arguments:
+// Jacobian_flag = true to compute the Jacobian coefficients,
+//		   false to skip this.
+//
 void horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
-		      const struct geometry_interpolator_info& ii)
+		      const struct geometry_interpolator_info& ii,
+		      bool Jacobian_flag)
 {
 CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");
 
@@ -101,9 +108,9 @@ setup_xyz_posns(ps);
 // interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$
 interpolate_geometry(ps, cgi, ii);
 
-// compute remaining gridfns --> $H$
+// compute remaining gridfns --> $H$ and optionally Jacobian coefficients
 // by algebraic ops and angular finite differencing
-compute_H(ps);
+compute_H(ps, Jacobian_flag);
 }
 
 //******************************************************************************
@@ -302,8 +309,9 @@ void* const output_arrays[N_OUTPUT_ARRAYS]
 bool first_time = true;
 if (first_time)
    then {
-	// store derivative info in interpolator parameter table
 	first_time = false;
+
+	// store derivative info in interpolator parameter table
 	CCTK_VInfo(CCTK_THORNSTRING,
 		   "         setting up derivative info for interpolator");
 
@@ -404,11 +412,17 @@ if (status < 0)
 //******************************************************************************
 
 //
-// This function computes H(h) by a mixture of algebraic operations
-// and (rho,sigma) finite differencing, on the angular grid.
+// This function computes H(h), and optionally its Jacobian coefficients,
+// (from which the Jacobian matrix may be computed later).  This function
+// uses a mixture of algebraic operations and (rho,sigma) finite differencing.
+// The computation is done (entirely) on the nominal angular grid.
+//
+// Arguments:
+// Jacobian_flag = true to compute the Jacobian coefficients,
+//		   false to skip this.
 //
 namespace {
-void compute_H(patch_system& ps)
+void compute_H(patch_system& ps, bool Jacobian_flag)
 {
 CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
 
@@ -460,22 +474,40 @@ CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
 		// ... each cg/*.c file has a separate set of temp variables,
 		//     and so must be inside its own set of { } braces
 		//
+
+		// gridfn #defins
 		#include "cg.hh"
+
 		  {
-		#include "cg/inverse_metric.c"
+		// g_uu
+		#include "../gr.cg/inverse_metric.c"
 		  }
+
 		  {
-		#include "cg/extrinsic_curvature_trace_raise.c"
+		// K, K_uu
+		#include "../gr.cg/extrinsic_curvature_trace_raise.c"
 		  }
+
 		  {
-		#include "cg/inverse_metric_gradient.c"
+		// partial_d_g_uu
+		#include "../gr.cg/inverse_metric_gradient.c"
 		  }
+
 		  {
-		#include "cg/metric_det_gradient.c"
+		// partial_d_ln_sqrt_g
+		#include "../gr.cg/metric_det_gradient.c"
 		  }
+
 		  {
-		#include "cg/horizon_function.c"
+		// HA, HB, HC, HD, H
+		#include "../gr.cg/horizon_function.c"
 		  }
+
+		if (Jacobian_flag)
+		   then {
+			// partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h
+			#include "../gr.cg/horizon_Jacobian.c"
+			}
 		}
 		}
 	}
diff --git a/src/gr/maple.log b/src/gr/maple.log
index 7732b11..2e5ddd0 100644
--- a/src/gr/maple.log
+++ b/src/gr/maple.log
@@ -634,7 +634,8 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
     assert_fnd_exists(g_dd);
     assert_fnd_exists(g_uu, fnd);
     g_uu__fnd := linalg[inverse](g_dd);
-    if cg_flag then codegen2(g_uu__fnd, 'g_uu', "cg/inverse_metric.c")
+    if cg_flag then
+        codegen2(g_uu__fnd, 'g_uu', "../gr.cg/inverse_metric.c")
     end if;
     NULL
 end proc
@@ -661,7 +662,7 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
         end do
     end do;
     if cg_flag then codegen2([K__fnd, K_uu__fnd], ['K', 'K_uu'],
-        "cg/extrinsic_curvature_trace_raise.c")
+        "../gr.cg/extrinsic_curvature_trace_raise.c")
     end if;
     NULL
 end proc
@@ -692,7 +693,7 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
         end do
     end do;
     if cg_flag then codegen2(partial_d_g_uu__fnd, 'partial_d_g_uu',
-        "cg/inverse_metric_gradient.c")
+        "../gr.cg/inverse_metric_gradient.c")
     end if;
     NULL
 end proc
@@ -714,7 +715,7 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
         'g_uu[i, j]*partial_d_g_dd[k, i, j]', 'i' = 1 .. N, 'j' = 1 .. N))
     end do;
     if cg_flag then codegen2(partial_d_ln_sqrt_g__fnd,
-        'partial_d_ln_sqrt_g', "cg/metric_det_gradient.c")
+        'partial_d_ln_sqrt_g', "../gr.cg/metric_det_gradient.c")
     end if;
     NULL
 end proc
@@ -807,7 +808,8 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
     H__fnd :=
         HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K
         ;
-    if cg_flag then codegen2(H__fnd, 'H', "cg/horizon_function.c") end if;
+    if cg_flag then codegen2(H__fnd, 'H', "../gr.cg/horizon_function.c")
+    end if;
     NULL
 end proc
 
@@ -817,7 +819,7 @@ end proc
 
 > auxiliary(true);
 inverse_metric...
-codegen2(g_uu) --> "cg/inverse_metric.c"
+codegen2(g_uu) --> "../gr.cg/inverse_metric.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
@@ -826,7 +828,7 @@ codegen2(g_uu) --> "cg/inverse_metric.c"
    find temporary variables
       --> `codegen2/temps`
    convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=1000208, alloc=917336, time=0.10
+bytes used=1000148, alloc=917336, time=0.12
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
@@ -834,12 +836,12 @@ bytes used=1000208, alloc=917336, time=0.10
       --> `codegen2/fix_rationals`
    writing C code
 extrinsic_curvature_trace_raise...
-codegen2([K, K_uu]) --> "cg/extrinsic_curvature_trace_raise.c"
+codegen2([K, K_uu]) --> "../gr.cg/extrinsic_curvature_trace_raise.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=2000708, alloc=1441528, time=0.15
+bytes used=2000796, alloc=1441528, time=0.22
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
@@ -847,39 +849,39 @@ bytes used=2000708, alloc=1441528, time=0.15
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
-bytes used=3000940, alloc=1572576, time=0.21
+bytes used=3001776, alloc=1572576, time=0.28
    convert p/q --> RATIONAL(p/q)
       --> `codegen2/fix_rationals`
    writing C code
 > metric_derivs(true);
 inverse_metric_gradient...
-bytes used=4001284, alloc=1638100, time=0.26
-codegen2(partial_d_g_uu) --> "cg/inverse_metric_gradient.c"
+bytes used=4002196, alloc=1638100, time=0.36
+codegen2(partial_d_g_uu) --> "../gr.cg/inverse_metric_gradient.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=5002720, alloc=1638100, time=0.35
+bytes used=5002380, alloc=1638100, time=0.43
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
    convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
       --> `codegen2/fix_Diff`
-bytes used=6002996, alloc=1703624, time=0.41
+bytes used=6004052, alloc=1703624, time=0.50
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
-bytes used=7003240, alloc=1703624, time=0.46
+bytes used=7004228, alloc=1703624, time=0.57
    convert p/q --> RATIONAL(p/q)
       --> `codegen2/fix_rationals`
    writing C code
-bytes used=8003492, alloc=1703624, time=0.49
+bytes used=8004548, alloc=1703624, time=0.63
 metric_det_gradient...
-codegen2(partial_d_ln_sqrt_g) --> "cg/metric_det_gradient.c"
+codegen2(partial_d_ln_sqrt_g) --> "../gr.cg/metric_det_gradient.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=9003660, alloc=1703624, time=0.60
+bytes used=9004804, alloc=1703624, time=0.74
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
@@ -894,40 +896,40 @@ codegen/C/expression:   Unknown function:   RATIONAL   will be left as is.
 > horizon(true);
 non_unit_normal...
 non_unit_normal_deriv...
-bytes used=10003824, alloc=1769148, time=0.66
+bytes used=10004956, alloc=1769148, time=0.85
 horizon_function...
-bytes used=11004072, alloc=1834672, time=0.73
-codegen2(H) --> "cg/horizon_function.c"
+bytes used=11005112, alloc=1834672, time=0.93
+codegen2(H) --> "../gr.cg/horizon_function.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=12004376, alloc=1834672, time=0.79
-bytes used=13004892, alloc=2031244, time=0.87
-bytes used=14005160, alloc=2031244, time=0.96
+bytes used=12005688, alloc=1834672, time=1.01
+bytes used=13005864, alloc=2031244, time=1.12
+bytes used=14006148, alloc=2031244, time=1.21
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
    convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=15005740, alloc=2096768, time=1.01
+bytes used=15006312, alloc=2096768, time=1.28
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
-bytes used=16006176, alloc=2096768, time=1.07
+bytes used=16006476, alloc=2096768, time=1.35
       --> `codegen2/unindex`
-bytes used=17006396, alloc=2096768, time=1.13
+bytes used=17006632, alloc=2096768, time=1.41
    convert p/q --> RATIONAL(p/q)
-bytes used=18006808, alloc=2096768, time=1.18
+bytes used=18006900, alloc=2096768, time=1.48
       --> `codegen2/fix_rationals`
    writing C code
 codegen/C/expression:   Unknown function:   PARTIAL_RHO   will be left as is.
 codegen/C/expression:   Unknown function:   PARTIAL_SIGMA   will be left as is.
-bytes used=19007204, alloc=2096768, time=1.22
+bytes used=19007116, alloc=2096768, time=1.54
 codegen/C/expression:   Unknown function:   PARTIAL_RHO_RHO   
 will be left as is.
 codegen/C/expression:   Unknown function:   PARTIAL_RHO_SIGMA   
 will be left as is.
 codegen/C/expression:   Unknown function:   PARTIAL_SIGMA_SIGMA   
 will be left as is.
-bytes used=20007392, alloc=2096768, time=1.38
+bytes used=20007280, alloc=2096768, time=1.66
 > quit
-bytes used=20239428, alloc=2096768, time=1.40
+bytes used=20238456, alloc=2096768, time=1.70
diff --git a/src/gr/setup_gr_gfas.maple b/src/gr/setup_gr_gfas.maple
index ede82ca..dc079cb 100644
--- a/src/gr/setup_gr_gfas.maple
+++ b/src/gr/setup_gr_gfas.maple
@@ -22,7 +22,7 @@ make_gfa('K_dd', {inert}, [1..N, 1..N], symmetric);
 # index-raised and contracted metric and extrinsic curvature
 make_gfa('g_uu', {inert,fnd}, [1..N, 1..N], symmetric);
 make_gfa('K_uu', {inert,fnd}, [1..N, 1..N], symmetric);
-make_gfa('K', {inert, fnd}, [], none);
+make_gfa('K', {inert,fnd}, [], none);
 
 # xyz partial derivatives of metric
 # ... as far as Maple is concerned, these are indeed separate gridfns;
@@ -37,22 +37,40 @@ make_gfa('partial_d_ln_sqrt_g', {inert,fnd}, [1..N], none);
 make_gfa('partial_d_g_uu', {inert,fnd}, [1..N, 1..N, 1..N], symmetric3_23);
 
 # radius of horizon
-make_gfa('h', {inert, fnd}, [], none);
+make_gfa('h', {inert,fnd}, [], none);
 
 # outward-pointing *non*-unit normal (covector) to horizon
 # and it's xyz-coordinate partial derivatives
 make_gfa('s_d', {inert,fnd}, [1..N], none);
 make_gfa('partial_d_s_d', {inert,fnd}, [1..N, 1..N], none);
 
-# LHS of apparent horizon equation
-make_gfa('H', {inert, fnd}, [], none);
-
 # subexpressions for computing LHS of apparent horizon equation
 # ... these are defined by (15) in my 1996 apparent horizon finding paper
-make_gfa('HA', {inert, fnd}, [], none);
-make_gfa('HB', {inert, fnd}, [], none);
-make_gfa('HC', {inert, fnd}, [], none);
-make_gfa('HD', {inert, fnd}, [], none);
+make_gfa('HA', {inert,fnd}, [], none);
+make_gfa('HB', {inert,fnd}, [], none);
+make_gfa('HC', {inert,fnd}, [], none);
+make_gfa('HD', {inert,fnd}, [], none);
+
+# LHS of apparent horizon equation
+make_gfa('H', {inert,fnd}, [], none);
+
+# Jacobian coefficients for H[ABCD]
+# these are the partial derivatives of H[ABCD]
+# ... wrt Diff(h,x_rs[x])
+make_gfa('partial_HA_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_HB_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_HC_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_HD_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+# ... wrt Diff(h,x_rs[x],x_rs[y])
+make_gfa('partial_HA_wrt_partial_dd_h', {inert,fnd},
+	 [1..N_ang, 1..N_ang], symmetric);
+make_gfa('partial_HB_wrt_partial_dd_h', {inert,fnd},
+	 [1..N_ang, 1..N_ang], symmetric);
+
+# Jacobian coefficients for H itself
+make_gfa('partial_H_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_H_wrt_partial_dd_h', {inert,fnd},
+	 [1..N_ang, 1..N_ang], symmetric);
 
 NULL;
 end proc;
@@ -75,7 +93,7 @@ end proc;
 # var_seq = (in) (varargs) An expression sequence of the variables to
 #			   differentiate with respect to.
 #
-`Diff/gridfn` :=
+`Diff/gridfn2` :=
 proc(operand)			# varargs
 option remember;		# performance optimization
 global