add some code for the d/dr term in the Jacobian -- not finished yet

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

3 files changed, 46 insertions, 19 deletions

diff --git a/src/gr/AHFinderDirect.hh b/src/gr/AHFinderDirect.hh
index 42b77c9..a8305d2 100644
--- a/src/gr/AHFinderDirect.hh
+++ b/src/gr/AHFinderDirect.hh
@@ -55,6 +55,26 @@ struct	cactus_grid_info
 	const fp* K_dd_33_data;
 	};
 
+//
+// This struct holds the various parameters used in computing the Jacobian
+// matrix.
+//
+struct	Jacobian_info
+	{
+	fp perturbation_amplitude;
+	};
+
+//
+// This struct holds the various parameters used in solving the H(h) = 0
+// equations.
+//
+struct	solver_info
+	{
+	// stuff for Newton_solve()
+	int max_Newton_iterations;
+	fp H_norm_for_convergence;
+	};
+
 //******************************************************************************
 
 //
@@ -70,19 +90,17 @@ void horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
 		      const struct geometry_interpolator_info& gii,
 		      bool Jacobian_flag,
-		      jtutil::norm<fp>& H_norms,
+		      jtutil::norm<fp>* H_norms_ptr,
 		      bool msg_flag = true);
 
 // horizon_Jacobian.cc
 Jacobian& create_Jacobian(patch_system& ps,
 			  const char Jacobian_type[]);
-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);
+void horizon_Jacobian(patch_system& ps,
+		      const struct cactus_grid_info& cgi,
+		      const struct geometry_interpolator_info& gii,
+		      const struct Jacobian_info& Jac_info,
+		      Jacobian& Jac);
 
 // Newton.cc
 void Newton_solve(patch_system& ps,
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc
index bcaec78..99d8d44 100644
--- a/src/gr/horizon_Jacobian.cc
+++ b/src/gr/horizon_Jacobian.cc
@@ -4,7 +4,7 @@
 // <<<prototypes for functions local to this file>>>
 // create_Jacobian - create a Jacobian matrix of the appropriate type
 //
-// horizon_Jacobian_SD - compute the Jacobian by symbolic differentiation
+// horizon_Jacobian - compute the Jacobian (symbolic differentiation)
 /// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian
 //
 // horizon_Jacobian_NP - compute the Jacobian by numerical perturbation
@@ -84,12 +84,14 @@ else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 // This function computes 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.
+// equation (25) of my 1996 apparent horizon finding paper, except that
+// the d/dr term is done by numerical finite differencing.
 //
 // 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,
+//	H				# H(h) assumed to already be computed
+//	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.
@@ -211,6 +213,10 @@ Jac.zero();
 	}
 	}
     }
+
+// compute d/dr term by numerical finite differencing
+perturb_h(Jacobian_info.perturbation_amplitude);
+horizon_function(ps, cgi, gii, false, NULL, false);
 }
 
 //******************************************************************************
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 9c6cb17..95bb6b6 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -53,7 +53,7 @@ void interpolate_geometry(patch_system& ps,
 			  bool msg_flag);
 void compute_H(patch_system& ps,
 	       bool Jacobian_flag,
-	       jtutil::norm<fp>& H_norms,
+	       jtutil::norm<fp>& H_norms_ptr,
 	       bool msg_flag);
 	  }
 
@@ -96,13 +96,16 @@ void compute_H(patch_system& ps,
 // 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.
+// H_norms_ptr = (out) If this pointer is non-NULL, the norm object it
+//		       points to is updated with all the H values in the
+//		       grid.  This norm object can then be used to compute
+//		       various (gridwise) norms of H.
 //
 void horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
 		      const struct geometry_interpolator_info& ii,
 		      bool Jacobian_flag,
-		      jtutil::norm<fp>& H_norms,
+		      jtutil::norm<fp>* H_norms_ptr,
 		      bool msg_flag = true)
 {
 if (msg_flag)
@@ -119,7 +122,7 @@ 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, H_norms, msg_flag);
+compute_H(ps, Jacobian_flag, H_norms_ptr, msg_flag);
 }
 
 //******************************************************************************
@@ -435,12 +438,11 @@ 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,
-	       jtutil::norm<fp>& H_norms,
+	       jtutil::norm<fp>* H_norms_ptr,
 	       bool msg_flag)
 {
 if (msg_flag)
@@ -524,7 +526,8 @@ if (msg_flag)
 		  }
 
 		// update running norms of H(h) function
-		H_norms.data(H);
+		if (H_norms_ptr != NULL)
+		   then H_norms_ptr->data(H);
 
 		if (Jacobian_flag)
 		   then {