add a whole bunch of changes from working at home

--> AHFinderDirect now finds AH correctly for Kerr/offset!!!


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

7 files changed, 287 insertions, 114 deletions

diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh
index 8e5fc8c..42b77c9 100644
--- a/src/gr/AHFinderDirect.hh
+++ b/src/gr/AHFinderDirect.hh
@@ -69,10 +69,26 @@ extern "C"
 void horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
 		      const struct geometry_interpolator_info& gii,
-		      bool Jacobian_flag);
+		      bool Jacobian_flag,
+		      jtutil::norm<fp>& H_norms,
+		      bool msg_flag = true);
 
 // horizon_Jacobian.cc
-Jacobian& create_Jacobian(const patch_system& ps,
+Jacobian& create_Jacobian(patch_system& ps,
 			  const char Jacobian_type[]);
-void horizon_Jacobian(patch_system& ps,
-		      Jacobian& Jac);
+void horizon_Jacobian_SD(patch_system& ps,
+			 Jacobian& Jac);
+void horizon_Jacobian_NP(patch_system& ps,
+			 const struct cactus_grid_info& cgi,
+			 const struct geometry_interpolator_info& ii,
+			 Jacobian& Jac,
+			 fp perturbation_amplitude);
+
+// Newton.cc
+void Newton_solve(patch_system& ps,
+		  const struct cactus_grid_info& cgi,
+		  const struct geometry_interpolator_info& gii,
+		  const char Jacobian_type[],
+		  fp perturbation_amplitude,
+		  int max_Newton_iterations,
+		  fp H_norm_for_convergence);
diff --git a/src/gr/cg.hh b/src/gr/cg.hh
index 9008159..cc0ab31 100644
--- a/src/gr/cg.hh
+++ b/src/gr/cg.hh
@@ -102,6 +102,8 @@
 #define partial_H_wrt_partial_dd_h_22	\
 	p.gridfn(nominal_gfns::gfn__partial_H_wrt_partial_dd_h_22, irho,isigma)
 
+#define Delta_h	p.gridfn(nominal_gfns::gfn__Delta_h, irho,isigma)
+
 //******************************************************************************
 
 //
diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index 546a488..ca978ae 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -3,9 +3,10 @@
 //
 // <<<prototypes for functions local to this file>>>
 // AHFinderDirect_driver - top-level driver
+///
+/// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords)
 /// setup_ellipsoid - setup up a coordiante ellipsoid in h
-/// setup_Kerr_KerrSchild_horizon - set up Kerr/Kerr-Schild horizon in h
-//
+///
 
 #include <stdio.h>
 #include <assert.h>
@@ -46,12 +47,13 @@ using jtutil::error_exit;
 //
 
 namespace {
+void setup_Kerr_horizon(patch_system& ps,
+			fp x_center, fp y_center, fp z_center,
+			fp m, fp a,
+			bool Kerr_Schild_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);
-void setup_Kerr_KerrSchild_horizon(patch_system& ps,
-				   fp x_center, fp y_center, fp z_center,
-				   fp mass, fp spin);
 	  };
 
 //******************************************************************************
@@ -139,7 +141,6 @@ cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz"));
 //
 // create the patch system and initialize the xyz derivative coefficients
 //
-CCTK_VInfo(CCTK_THORNSTRING, "   creating patch system");
 patch_system ps(origin_x, origin_y, origin_z,
 		patch_system::type_of_name(patch_system_type),
 		N_ghost_points, N_overlap_points, delta_drho_dsigma,
@@ -154,7 +155,7 @@ patch_system ps(origin_x, origin_y, origin_z,
 if	(STRING_EQUAL(initial_guess_method, "read from file"))
    then {
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   reading initial guess from \"%s\"",
+		   "reading initial guess from \"%s\"",
 		   initial_guess__read_from_file__file_name);
 	ps.read_ghosted_gridfn(ghosted_gfns::gfn__h,
 			       initial_guess__read_from_file__file_name,
@@ -168,13 +169,22 @@ else if (STRING_EQUAL(initial_guess_method, "ellipsoid"))
 			initial_guess__ellipsoid__x_radius,
 			initial_guess__ellipsoid__y_radius,
 			initial_guess__ellipsoid__z_radius);
+else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr"))
+   then setup_Kerr_horizon(ps,
+			   initial_guess__Kerr_KerrSchild__x_center,
+			   initial_guess__Kerr_KerrSchild__y_center,
+			   initial_guess__Kerr_KerrSchild__z_center,
+			   initial_guess__Kerr_KerrSchild__mass,
+			   initial_guess__Kerr_KerrSchild__spin,
+			   false);	// use Kerr coords
 else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild"))
-   then setup_Kerr_KerrSchild_horizon(ps,
-				      initial_guess__Kerr_KerrSchild__x_center,
-				      initial_guess__Kerr_KerrSchild__y_center,
-				      initial_guess__Kerr_KerrSchild__z_center,
-				      initial_guess__Kerr_KerrSchild__mass,
-				      initial_guess__Kerr_KerrSchild__spin);
+   then setup_Kerr_horizon(ps,
+			   initial_guess__Kerr_KerrSchild__x_center,
+			   initial_guess__Kerr_KerrSchild__y_center,
+			   initial_guess__Kerr_KerrSchild__z_center,
+			   initial_guess__Kerr_KerrSchild__mass,
+			   initial_guess__Kerr_KerrSchild__spin,
+			   true);	// use Kerr-Schild coords
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "unknown initial_guess_method=\"%s\"!",
 		   initial_guess_method);			/*NOTREACHED*/
@@ -187,18 +197,49 @@ ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
 //
 // find the apparent horizon
 //
-if      (STRING_EQUAL(method, "horizon"))
+jtutil::norm<fp> H_norms;
+if      (STRING_EQUAL(method, "horizon function"))
    then {
-	horizon_function(ps, cgi, gii, false);
+	horizon_function(ps, cgi, gii, false, H_norms);
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "   H(h) rms-norm %.2e, infinity-norm %.2e\n",
+		   H_norms.rms_norm(), H_norms.infinity_norm());
 	ps.print_gridfn_with_xyz(nominal_gfns::gfn__H,
 				 true, ghosted_gfns::gfn__h,
 				 "H.dat");
 	}
-else if (STRING_EQUAL(method, "horizon Jacobian"))
+else if (STRING_EQUAL(method, "Jacobian"))
    then {
 	Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
-	horizon_function(ps, cgi, gii, true);
-	horizon_Jacobian(ps, Jac);
+	horizon_function(ps, cgi, gii, true, H_norms);
+	horizon_Jacobian_SD(ps, Jac);
+	print_Jacobian(Jacobian_file_name, Jac);
+	}
+else if (STRING_EQUAL(method, "Jacobian test"))
+   then {
+	Jacobian& SD_Jac = create_Jacobian(ps, Jacobian_type);
+	horizon_function(ps, cgi, gii, true, H_norms);
+	horizon_Jacobian_SD(ps, SD_Jac);
+
+	Jacobian& NP_Jac = create_Jacobian(ps, Jacobian_type);
+	horizon_function(ps, cgi, gii, true, H_norms);
+	horizon_Jacobian_NP(ps, cgi, gii,
+			    NP_Jac,
+			    NP_Jacobian__perturbation_amplitude);
+
+	print_Jacobians(Jacobian_file_name, SD_Jac, NP_Jac);
+	}
+else if (STRING_EQUAL(method, "Newton (NP Jacobian)"))
+   then {
+	Newton_solve(ps, cgi, gii,
+		     Jacobian_type,
+		     NP_Jacobian__perturbation_amplitude,
+		     max_Newton_iterations,
+		     H_norm_for_convergence);
+	ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
+					 true, ghosted_gfns::gfn__h,
+					 "h.dat",
+					 false);	// no ghost zones
 	}
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "unknown method=\"%s\"!",
@@ -206,6 +247,56 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// 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]_center = 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
+//
+namespace {
+void setup_Kerr_horizon(patch_system& ps,
+			fp x_center, fp y_center, fp z_center,
+			fp m, fp a,
+			bool Kerr_Schild_flag)
+{
+const char* const name = Kerr_Schild_flag ? "Kerr-Schild" : "Kerr";
+
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "setting up horizon for Kerr in %s coords",
+	   name);
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "   mass=%g, spin=J/m^2=%g, posn=(%g,%g,%g)",
+	   double(m), double(a),
+	   double(x_center), double(y_center), double(z_center));
+
+// 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;
+
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "   horizon is coordinate %s",
+	   Kerr_Schild_flag ? "ellipsoid" : "sphere");
+setup_ellipsoid(ps,
+		x_center, y_center, z_center,
+		xy_radius, xy_radius, z_radius);
+}
+	  }
+
+//******************************************************************************
 
 //
 // This function sets up an ellipsoid in the gridfn h, using the
@@ -234,10 +325,10 @@ void setup_ellipsoid(patch_system& ps,
 		     fp x_radius, fp y_radius, fp z_radius)
 {
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "   h = ellipsoid: center=(%g,%g,%g)",
+	   "setting h = ellipsoid: center=(%g,%g,%g)",
 	   double(x_center), double(y_center), double(z_center));
 CCTK_VInfo(CCTK_THORNSTRING,
-	   "                  radius=(%g,%g,%g)",
+	   "                       radius=(%g,%g,%g)",
 	   double(x_radius), double(y_radius), double(z_radius));
 
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
@@ -275,7 +366,7 @@ CCTK_VInfo(CCTK_THORNSTRING,
 		   then r = r_minus;
 		else    CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 				   "\n"
-"   expected exactly one r>0 solution, got 0 or 2!\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"
@@ -293,62 +384,3 @@ CCTK_VInfo(CCTK_THORNSTRING,
 	}
 }
 	  }
-
-//******************************************************************************
-
-//
-// This function sets up the horizon of a Kerr black hole in 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 of my AHFinderDirect notes.
-//
-// Arguments:
-// [xyz]_center = The position of the Kerr black hole.
-// (mass,spin) = Describe the Kerr black hole.
-//
-// FIXME FIXME: at present [xyz}_center are ignored.
-//
-namespace {
-void setup_Kerr_KerrSchild_horizon(patch_system& ps,
-				   fp x_center, fp y_center, fp z_center,
-				   fp mass, fp spin)
-{
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "   setting up Kerr/Kerr-Schild horizon");
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "      mass=%g, spin=%g, posn=(%g,%g,%g)",
-	   double(mass), double(spin),
-	   double(x_center), double(y_center), double(z_center));
-
-// 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)
-	{
-	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 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);
-
-		const fp KS_x = Kerr_x - mass*spin*sin(theta)*sin(phi);
-		const fp KS_y = Kerr_y + mass*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/gfn.hh b/src/gr/gfn.hh
index 24c56fe..9f54b5c 100644
--- a/src/gr/gfn.hh
+++ b/src/gr/gfn.hh
@@ -53,6 +53,7 @@ enum	{
 	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
+	gfn__Delta_h,
+	max_gfn = gfn__Delta_h		// no comma
 	};
 	  };
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc
index 7718796..bcaec78 100644
--- a/src/gr/horizon_Jacobian.cc
+++ b/src/gr/horizon_Jacobian.cc
@@ -3,9 +3,11 @@
 //
 // <<<prototypes for functions local to this file>>>
 // create_Jacobian - create a Jacobian matrix of the appropriate type
-// horizon_Jacobian - top-level driver to compute the Jacobian
+//
+// horizon_Jacobian_SD - compute the Jacobian by symbolic differentiation
 /// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian
 //
+// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation
 
 #include <stdio.h>
 #include <assert.h>
@@ -53,12 +55,18 @@ void add_ghost_zone_Jacobian(const patch_system& ps,
 	  }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
 // This function decodes the  Jacobian_type  parameter and creates the
 // appropriate type of Jacobian matrix object.
 //
-Jacobian& create_Jacobian(const patch_system& ps,
+// FIXME: the patch system shouldn't really have to be non-const, but
+//	  the Jacobian constructors all require this to allow the
+//	  linear solvers to directly update gridfns
+//
+Jacobian& create_Jacobian(patch_system& ps,
 			  const char Jacobian_type[])
 {
 if	(STRING_EQUAL(Jacobian_type, "dense matrix"))
@@ -69,12 +77,14 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
 // 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.
+// 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
@@ -84,9 +94,10 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 // Outputs:
 //	The Jacobian matrix is stored in the Jacobian object Jac.
 //
-void horizon_Jacobian(patch_system& ps, Jacobian& Jac)
+void horizon_Jacobian_SD(patch_system& ps,
+			 Jacobian& Jac)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian");
+CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian_SD");
 
 ps.compute_synchronize_Jacobian();
 Jac.zero();
@@ -189,7 +200,7 @@ Jac.zero();
 	    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, 
+			(ps, xp,
 		 xp.corner_ghost_zone_containing_point(m_irho < 0, m_isigma < 0,
 						       xm_irho, xm_isigma),
 			 II, xm_irho, xm_isigma,
@@ -254,3 +265,91 @@ const int ym_ipar_posn = ps.synchronize_Jacobian_y_ipar_posn
 	}
 }
 	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// 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:
+//
+// evaluate H(h)
+// array-copy H(h) to scratch array save_H
+//	for each point (y,JJ)
+//	{
+//	const fp save_h_y = h at y;
+//	h at y += perturbation_amplitude;
+//	evaluate H(h)
+//		for each point (x,II)
+//		{
+//		Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude;
+//		}
+//	h at y = save_h_y;
+//	}
+// array-copy save_H back to H(h)
+//
+void horizon_Jacobian_NP(patch_system& ps,
+			 const struct cactus_grid_info& cgi,
+			 const struct geometry_interpolator_info& ii,
+			 Jacobian& Jac,
+			 fp perturbation_amplitude)
+{
+CCTK_VInfo(CCTK_THORNSTRING, "   horizon Jacobian_NP");
+
+// evaluate H(h)
+jtutil::norm<fp> H_norms;
+horizon_function(ps, cgi, ii, false, H_norms);
+
+// array-copy H(h) to scratch array save_H
+fp *save_H = new fp[Jac.NN()];
+jtutil::array_copy(Jac.NN(), ps.gridfn_data(nominal_gfns::gfn__H), save_H);
+
+	for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
+	{
+	patch& yp = ps.ith_patch(ypn);
+	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(ghosted_gfns::gfn__h,
+					      y_irho,y_isigma);
+	yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma)
+		+= perturbation_amplitude;
+	H_norms.reset();
+	horizon_function(ps, cgi, ii, false, H_norms, false);
+		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 new_xH = xp.gridfn(nominal_gfns::gfn__H,
+					    x_irho,x_isigma);
+		Jac(II,JJ) = (new_xH - save_H[II]) / perturbation_amplitude;
+		}
+		}
+		}
+	yp.ghosted_gridfn(ghosted_gfns::gfn__h, y_irho,y_isigma)
+		= save_h_y;
+	}
+	}
+   	} 
+
+jtutil::array_copy(Jac.NN(), save_H, ps.gridfn_data(nominal_gfns::gfn__H));
+delete[] save_H;
+}
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 82f5f31..9c6cb17 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -46,11 +46,15 @@ using jtutil::error_exit;
 //
 
 namespace {
-void setup_xyz_posns(patch_system& ps);
+void setup_xyz_posns(patch_system& ps, bool msg_flag);
 void interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
-			  const struct geometry_interpolator_info& ii);
-void compute_H(patch_system& ps, bool Jacobian_flag);
+			  const struct geometry_interpolator_info& ii,
+			  bool msg_flag);
+void compute_H(patch_system& ps,
+	       bool Jacobian_flag,
+	       jtutil::norm<fp>& H_norms,
+	       bool msg_flag);
 	  }
 
 //******************************************************************************
@@ -92,26 +96,30 @@ void compute_H(patch_system& ps, bool Jacobian_flag);
 // Arguments:
 // Jacobian_flag = true to compute the Jacobian coefficients,
 //		   false to skip this.
+// H_norms = (out) This is set to the gridwise norms of the H(h) function.
 //
 void horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
 		      const struct geometry_interpolator_info& ii,
-		      bool Jacobian_flag)
+		      bool Jacobian_flag,
+		      jtutil::norm<fp>& H_norms,
+		      bool msg_flag = true)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");
 
 // fill in values of all ghosted gridfns in ghost zones
 ps.synchronize();
 
 // set up xyz positions of grid points
-setup_xyz_posns(ps);
+setup_xyz_posns(ps, msg_flag);
 
 // interpolate $g_{ij}$, $K_{ij}$ --> $g_{ij}$, $K_{ij}$, $\partial_k g_{ij}$
-interpolate_geometry(ps, cgi, ii);
+interpolate_geometry(ps, cgi, ii, msg_flag);
 
 // compute remaining gridfns --> $H$ and optionally Jacobian coefficients
 // by algebraic ops and angular finite differencing
-compute_H(ps, Jacobian_flag);
+compute_H(ps, Jacobian_flag, H_norms, msg_flag);
 }
 
 //******************************************************************************
@@ -121,9 +129,11 @@ compute_H(ps, Jacobian_flag);
 // in the gridfns global_[xyz].  These will be used by interplate_geometry().
 //
 namespace {
-void setup_xyz_posns(patch_system& ps)
+void setup_xyz_posns(patch_system& ps, bool msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "      xyz positions and derivative coefficients");
 
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
 	{
@@ -188,10 +198,12 @@ CCTK_VInfo(CCTK_THORNSTRING, "      xyz positions and derivative coefficients");
 namespace {
 void interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
-			  const struct geometry_interpolator_info& gii)
+			  const struct geometry_interpolator_info& gii,
+			  bool msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "      interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid");
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "      interpolate $g_{ij}$, $K_{ij}$ from Cactus 3-D grid");
 
 int status;
 
@@ -313,8 +325,9 @@ if (first_time)
 	first_time = false;
 
 	// store derivative info in interpolator parameter table
-	CCTK_VInfo(CCTK_THORNSTRING,
-		   "         setting up derivative info for interpolator");
+	if (msg_flag)
+	   then CCTK_VInfo(CCTK_THORNSTRING,
+			   "         setting up derivative info for interpolator");
 
 	status = Util_TableSetIntArray(gii.param_table_handle,
 				       N_OUTPUT_ARRAYS, operand_indices,
@@ -339,9 +352,10 @@ if (first_time)
 			   status);				/*NOTREACHED*/
 	}
 
-CCTK_VInfo(CCTK_THORNSTRING,
-	   "         calling interpolator (%d points)",
-	   N_interp_points);
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING,
+		   "         calling interpolator (%d points)",
+		   N_interp_points);
 status = CCTK_InterpLocalUniform(N_GRID_DIMS,
 				 gii.operator_handle, gii.param_table_handle,
 				 cgi.coord_origin, cgi.coord_delta,
@@ -421,11 +435,16 @@ if (status < 0)
 // Arguments:
 // Jacobian_flag = true to compute the Jacobian coefficients,
 //		   false to skip this.
+// H_norms = (out) This is set to the gridwise norms of the H(h) function.
 //
 namespace {
-void compute_H(patch_system& ps, bool Jacobian_flag)
+void compute_H(patch_system& ps,
+	       bool Jacobian_flag,
+	       jtutil::norm<fp>& H_norms,
+	       bool msg_flag)
 {
-CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
+if (msg_flag)
+   then CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
 
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
 	{
@@ -504,6 +523,9 @@ CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");
 		#include "../gr.cg/horizon_function.c"
 		  }
 
+		// update running norms of H(h) function
+		H_norms.data(H);
+
 		if (Jacobian_flag)
 		   then {
 			// partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h
diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn
index 7def0d5..8e9745a 100644
--- a/src/gr/make.code.defn
+++ b/src/gr/make.code.defn
@@ -4,7 +4,8 @@
 # Source files in this directory
 SRCS = driver.cc \
        horizon_function.cc \
-       horizon_Jacobian.cc
+       horizon_Jacobian.cc \
+       Newton.cc
 
 # Subdirectories containing source files
 SUBDIRS =