add an option to do the NP Jacobian only for testing

(useful if SD Jacobian core-dumps)


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

1 files changed, 90 insertions, 58 deletions

diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index cf26d0c..5da7b3f 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -65,7 +65,8 @@ void setup_ellipsoid(patch_system& ps,
 		     fp x_radius, fp y_radius, fp z_radius);
 
 void print_Jacobians(const patch_system& ps,
-		     const Jacobian& Jac_NP, const Jacobian& Jac_SD_FDdr,
+		     const Jacobian* Jac_NP,
+		     bool do_non_NP_Jac, const Jacobian* Jac_SD_FDdr,
 		     const char file_name[]);
 	  };
 
@@ -267,24 +268,36 @@ if      (STRING_EQUAL(method, "horizon function"))
 				 true, gfns::gfn__h,
 				 H_of_h_file_name);
 	}
-else if (STRING_EQUAL(method, "Jacobian test"))
+else if (    STRING_EQUAL(method, "Jacobian test")
+	  || STRING_EQUAL(method, "Jacobian test (NP only)")    )
    then {
-	// symbolic differentiation with finite diff d/dr
-	Jacobian& Jac_SD_FDdr = new_Jacobian(ps, Jacobian_storage_method);
-	horizon_function(ps, cgi, gi);
-	Jacobian_info.Jacobian_method = Jacobian__symbolic_diff_with_FD_dr;
-	horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_SD_FDdr);
+	const bool NP_only = STRING_EQUAL(method, "Jacobian test (NP only)");
+	const bool do_non_NP_Jac = ! NP_only;
 
 	// numerical perturbation
-	Jacobian& Jac_NP = new_Jacobian(ps, Jacobian_storage_method);
-	horizon_function(ps, cgi, gi);
+	Jacobian* const Jac_NP = & new_Jacobian(ps, Jacobian_storage_method);
+	horizon_function(ps, cgi, gi, true);
 	Jacobian_info.Jacobian_method = Jacobian__numerical_perturb;
-	horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_NP);
+	horizon_Jacobian(ps, cgi, gi, Jacobian_info, *Jac_NP);
+
+	// symbolic differentiation with finite diff d/dr
+	Jacobian* Jac_SD_FDdr = NULL;
+	if (do_non_NP_Jac)
+	   then {
+		Jac_SD_FDdr = & new_Jacobian(ps, Jacobian_storage_method);
+		horizon_function(ps, cgi, gi, true);
+		Jacobian_info.Jacobian_method
+			= Jacobian__symbolic_diff_with_FD_dr;
+		horizon_Jacobian(ps, cgi, gi, Jacobian_info, *Jac_SD_FDdr);
+		}
 
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "writing Jacobians to \"%s\"",
-		   Jacobian_file_name);
-	print_Jacobians(ps, Jac_NP, Jac_SD_FDdr, Jacobian_file_name);
+		   "writing Jacobian%s to \"%s\"",
+		   (NP_only ? "" : "s"), Jacobian_file_name);
+	print_Jacobians(ps,
+			Jac_NP,
+			do_non_NP_Jac, Jac_SD_FDdr,
+			Jacobian_file_name);
 	}
 else if (STRING_EQUAL(method, "Newton solve"))
    then {
@@ -505,12 +518,17 @@ CCTK_VInfo(CCTK_THORNSTRING,
 //******************************************************************************
 
 //
-// This function prints a pair of Jacobian matrices (and their difference)
-// to a named output file.
+// This function prints one or two Jacobian matrices (and their difference
+// in the latter case) to a named output file.
+//
+// Arguments:
+// do_non_NP_Jac = true if we should also print *Jac_SD_FDdr
+//		   false to skip this (and print only *Jac_NP)
 //
 namespace {
 void print_Jacobians(const patch_system& ps,
-		     const Jacobian& Jac_NP, const Jacobian& Jac_SD_FDdr,
+		     const Jacobian* Jac_NP,
+		     bool do_non_NP_Jac, const Jacobian* Jac_SD_FDdr,
 		     const char file_name[])
 {
 FILE *fileptr = fopen(file_name, "w");
@@ -519,54 +537,67 @@ if (fileptr == NULL)
 		   "print_Jacobians(): can't open output file \"%s\"!",
 		   file_name);					/*NOTREACHED*/
 
-fprintf(fileptr, "# column 1 = x II\n");
+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 = y JJ\n");
-fprintf(fileptr, "# column 6 = y patch number\n");
-fprintf(fileptr, "# column 7 = y irho\n");
-fprintf(fileptr, "# column 8 = y isigma\n");
-fprintf(fileptr, "# column 9 = Jac_NP(II,JJ)\n");
-fprintf(fileptr, "# column 10 = Jac_SD_FDdr(II,JJ)\n");
-fprintf(fileptr, "# column 11 = abs error in Jac_SD_FDdr\n");
-fprintf(fileptr, "# column 12 = rel error in Jac_SD_FDdr\n");
-
-	for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
+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 (do_non_NP_Jac)
+   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& xp = ps.ith_patch(xpn);
+	patch& yp = ps.ith_patch(ypn);
 
-	for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
+	for (int y_irho = yp.min_irho() ;
+	     y_irho <= yp.max_irho() ;
+	     ++y_irho)
 	{
-	for (int x_isigma = xp.min_isigma() ;
-	     x_isigma <= xp.max_isigma() ;
-	     ++x_isigma)
+	for (int y_isigma = yp.min_isigma() ;
+	     y_isigma <= yp.max_isigma() ;
+	     ++y_isigma)
 	{
-	const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);
+	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);
 
-		for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
-		{
-		patch& yp = ps.ith_patch(ypn);
+	if (! Jac_NP->is_explicitly_stored(II,JJ))
+	   then continue;			// skip sparse points
 
-		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 NP = (*Jac_NP)(II,JJ);
+	const fp SD_FDdr = do_non_NP_Jac ? (*Jac_SD_FDdr)(II,JJ) : 0.0;
 
-		if (! Jac_SD_FDdr.is_explicitly_stored(II,JJ))
-		   then continue;			// skip sparse points
+	if ((NP == 0.0) && (SD_FDdr == 0.0))
+	   then continue;		// skip zero values (if == )
 
-		const fp NP      = Jac_NP     (II,JJ);
-		const fp SD_FDdr = 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 (do_non_NP_Jac)
+	   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);
@@ -574,17 +605,18 @@ fprintf(fileptr, "# column 12 = rel error in Jac_SD_FDdr\n");
 		const fp SD_FDdr_rel_error = SD_FDdr_abs_error / max_abs;
 
 		fprintf(fileptr,
-			"%d %d %d %d\t%d %d %d %d\t%.10g\t%.10g\t%e\t%e\n",
-			II, xpn, x_irho, x_isigma,
-			JJ, ypn, y_irho, y_isigma,
-			double(NP), double(SD_FDdr),
+			"\t%.10g\t%e\t%e",
+			double(SD_FDdr),
 			double(SD_FDdr_abs_error), double(SD_FDdr_rel_error));
 		}
-		}
-		}
+
+	fprintf(fileptr, "\n");
 	}
 	}
 	}
+    }
+    }
+    }
 
 fclose(fileptr);
 }