various changes including d/dr terms in Jacobian by numerical perturbation,

tweak I/O parameters,
move printing Jacobian out of Jacobian class into test driver,
drop unused array BLAS routines in jtutil::


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

1 files changed, 188 insertions, 60 deletions

diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index ca978ae..00695e9 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -7,6 +7,8 @@
 /// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords)
 /// setup_ellipsoid - setup up a coordiante ellipsoid in h
 ///
+/// print_Jacobians - print a pair of Jacobians
+///
 
 #include <stdio.h>
 #include <assert.h>
@@ -54,6 +56,10 @@ void setup_Kerr_horizon(patch_system& ps,
 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 print_Jacobians(const patch_system& ps,
+		     const Jacobian& Jac, const Jacobian& NP_Jac,
+		     const char file_name[]);
 	  };
 
 //******************************************************************************
@@ -155,91 +161,124 @@ 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\"",
-		   initial_guess__read_from_file__file_name);
+		   "reading initial guess h from \"%s\"",
+		   h_file_name);
 	ps.read_ghosted_gridfn(ghosted_gfns::gfn__h,
-			       initial_guess__read_from_file__file_name,
+			       h_file_name,
 			       false);		// no ghost zones
 	}
-else if (STRING_EQUAL(initial_guess_method, "ellipsoid"))
-   then setup_ellipsoid(ps,
-			initial_guess__ellipsoid__x_center,
-			initial_guess__ellipsoid__y_center,
-			initial_guess__ellipsoid__z_center,
-			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_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*/
-ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
-				 true, ghosted_gfns::gfn__h,
-				 "h.dat",
-				 false);	// no ghost zones
+   else {
+	if	(STRING_EQUAL(initial_guess_method, "ellipsoid"))
+	   then setup_ellipsoid(ps,
+				initial_guess__ellipsoid__x_center,
+				initial_guess__ellipsoid__y_center,
+				initial_guess__ellipsoid__z_center,
+				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_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*/
+
+	// write initial guess back to this file
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "writing initial guess h to \"%s\"",
+		   h_file_name);
+	ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
+					 true, ghosted_gfns::gfn__h,
+					 h_file_name,
+					 false);	// no ghost zones
+	}
 
 
 //
 // find the apparent horizon
 //
-jtutil::norm<fp> H_norms;
+struct Jacobian_info Jac_info;
+Jac_info.perturbation_amplitude = NP_Jacobian__perturbation_amplitude;
+
+struct solver_info solver_info;
+solver_info.max_Newton_iterations = max_Newton_iterations;
+solver_info.H_norm_for_convergence = H_norm_for_convergence;
+solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
+
 if      (STRING_EQUAL(method, "horizon function"))
    then {
-	horizon_function(ps, cgi, gii, false, H_norms);
+	jtutil::norm<fp> H_norms;
+
+	const int timer_handle = CCTK_TimerCreateI();
+	CCTK_TimerStartI(timer_handle);
+	  horizon_function(ps, cgi, gii, false, &H_norms);
+	CCTK_TimerStopI(timer_handle);
+	printf("timer stats for evaluating H(h):\n");
+	CCTK_TimerPrintDataI(timer_handle, -1);
+
 	CCTK_VInfo(CCTK_THORNSTRING,
-		   "   H(h) rms-norm %.2e, infinity-norm %.2e\n",
+		   "   H(h) rms-norm %.2e, infinity-norm %.2e",
 		   H_norms.rms_norm(), H_norms.infinity_norm());
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "writing H(h) to \"%s\"",
+		   H_of_h_file_name);
 	ps.print_gridfn_with_xyz(nominal_gfns::gfn__H,
 				 true, ghosted_gfns::gfn__h,
-				 "H.dat");
-	}
-else if (STRING_EQUAL(method, "Jacobian"))
-   then {
-	Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
-	horizon_function(ps, cgi, gii, true, H_norms);
-	horizon_Jacobian_SD(ps, Jac);
-	print_Jacobian(Jacobian_file_name, Jac);
+				 H_of_h_file_name);
 	}
 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);
+	jtutil::norm<fp> H_norms;
+	Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
+	horizon_function(ps, cgi, gii, true, &H_norms);
+	horizon_Jacobian(ps, cgi, gii, Jac_info, 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);
+	horizon_function(ps, cgi, gii, true, &H_norms);
+	horizon_Jacobian_NP(ps, cgi, gii, Jac_info, NP_Jac);
 
-	print_Jacobians(Jacobian_file_name, SD_Jac, NP_Jac);
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "writing Jacobians to \"%s\"",
+		   Jacobian_file_name);
+	print_Jacobians(ps, Jac, NP_Jac, Jacobian_file_name);
 	}
-else if (STRING_EQUAL(method, "Newton (NP Jacobian)"))
+else if (STRING_EQUAL(method, "Newton solve"))
    then {
-	Newton_solve(ps, cgi, gii,
-		     Jacobian_type,
-		     NP_Jacobian__perturbation_amplitude,
-		     max_Newton_iterations,
-		     H_norm_for_convergence);
+	Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
+
+	const int timer_handle = CCTK_TimerCreateI();
+	CCTK_TimerStartI(timer_handle);
+	  Newton_solve(ps, cgi, gii, Jac_info, solver_info, Jac);
+	CCTK_TimerStopI(timer_handle);
+	printf("timer stats for finding the apparent horizon:\n");
+	CCTK_TimerPrintDataI(timer_handle, -1);
+
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "writing final horizon shape h to \"%s\"",
+		   h_file_name);
 	ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
 					 true, ghosted_gfns::gfn__h,
-					 "h.dat",
+					 h_file_name,
 					 false);	// no ghost zones
+	CCTK_VInfo(CCTK_THORNSTRING,
+		   "writing H(h) to \"%s\"",
+		   H_of_h_file_name);
+	ps.print_gridfn_with_xyz(nominal_gfns::gfn__H,
+				 true, ghosted_gfns::gfn__h,
+				 H_of_h_file_name);
 	}
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "unknown method=\"%s\"!",
@@ -384,3 +423,92 @@ CCTK_VInfo(CCTK_THORNSTRING,
 	}
 }
 	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function prints a pair of Jacobian matrices (and their difference)
+// to a named output file.
+//
+namespace {
+void print_Jacobians(const patch_system& ps,
+		     const Jacobian& Jac, const Jacobian& NP_Jac,
+		     const char file_name[])
+{
+FILE *fileptr = fopen(file_name, "w");
+if (fileptr == NULL)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "print_Jacobians(): can't open output file \"%s\"!",
+		   file_name);					/*NOTREACHED*/
+
+fprintf(fileptr, "# column 1 = x II\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(II,JJ)\n");
+fprintf(fileptr, "# column 10 = NP_Jac(II,JJ)\n");
+fprintf(fileptr, "# column 11 = abs error = Jac - NP_Jac\n");
+fprintf(fileptr, "# column 12 = rel error = abs error / max(|Jac|,|NP_Jac|)\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& yp = ps.ith_patch(ypn);
+
+		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);
+
+		if (! Jac.is_explicitly_stored(II,JJ))
+		   then continue;			// skip sparse points
+
+		const fp SD = Jac(II,JJ);
+		const fp NP = NP_Jac(II,JJ);
+		const fp abs_error = SD - NP;
+
+		if ((SD == 0.0) && (NP == 0.0))
+		   then continue;		// skip zero values (if == )
+
+		const fp abs_SD = jtutil::abs(SD);
+		const fp abs_NP = jtutil::abs(NP);
+		const fp rel_error = abs_error / jtutil::max(abs_SD, abs_NP);
+
+		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(SD), double(NP),
+			double(abs_error), double(rel_error));
+		}
+		}
+		}
+	}
+	}
+	}
+
+fclose(fileptr);
+}
+	  }