complete interpatch-interpolation Jacobian code

add 1st draft of code in src/util/test_patch_system.cc to test this


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

1 files changed, 188 insertions, 2 deletions

diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc
index 16d89a4..5143f26 100644
--- a/src/patch/test_patch_system.cc
+++ b/src/patch/test_patch_system.cc
@@ -98,12 +98,18 @@ extern "C"
   void test_patch_system(CCTK_ARGUMENTS);
 
 namespace {
+void test_synchronize_Jacobians(patch_system& ps,
+				int test_gfn, int NP_test_gfn,
+				fp perturbation_amplitude,
+				const char Jacobian_file_name[]);
+
 void setup_sym_fn_xyz(patch_system& ps, int ghosted_gfn, bool want_ghost_zones);
 void setup_fn_rho_sigma(patch_system& ps, int ghosted_gfn);
 void finite_diff(patch_system& ps,
 		 int ghosted_gfn_src, int gfn_dst,
 		 int which_derivs);
 void analytic_derivs(patch_system& ps, int gfn_dst, int which_derivs);
+
 void gridfn_minus(patch_system& ps, 
 		  int gfn_x, int gfn_y, int gfn_dst);
 void ghosted_gridfn_minus(patch_system& ps, 
@@ -169,28 +175,40 @@ if      (STRING_EQUAL(which_test, "gridfn"))
 	setup_sym_fn_xyz(ps, test_fn_gfn, true);
 	ps.print_ghosted_gridfn(test_fn_gfn, "test_fn.dat");
 	}
-else if (STRING_EQUAL(which_test, "read-gridfn"))
+
+else if (STRING_EQUAL(which_test, "read gridfn"))
    then {
 	ps.read_ghosted_gridfn(test_fn_gfn, "test_fn.dat");
 	ps.print_ghosted_gridfn(test_fn_gfn, "test_fn2.dat");
 	}
+
 else if (STRING_EQUAL(which_test, "synchronize"))
    then {
 	setup_sym_fn_xyz(ps, test_fn_gfn, false);
 	ps.print_ghosted_gridfn(test_fn_gfn, "test_fn_init.dat");
+
 	setup_sym_fn_xyz(ps, test_fn_copy_gfn, true);
 	ps.print_ghosted_gridfn(test_fn_copy_gfn, "test_fn_copy.dat");
+
 	ps.synchronize_ghost_zones(test_fn_gfn, test_fn_gfn);
 	ps.print_ghosted_gridfn(test_fn_gfn, "test_fn_sync.dat");
+
 	ghosted_gridfn_minus(ps,
 			     test_fn_gfn, test_fn_copy_gfn,
 			     ghosted_error_gfn);
 	ps.print_ghosted_gridfn(ghosted_error_gfn, "ghosted_error.dat");
 	}
+
+else if (STRING_EQUAL(which_test, "synchronize Jacobian"))
+   then test_synchronize_Jacobians(ps,
+				   test_fn_copy_gfn, test_fn_gfn,
+				   NP_Jacobian__perturbation_amplitude,
+				   Jacobian_file_name);
+
 else if (STRING_EQUAL(which_test, "derivatives"))
    then {
 	setup_fn_rho_sigma(ps, test_fn_gfn);
-	ps.print_gridfn(nominal_error_gfn, "test_fn.dat");
+	ps.print_gridfn(test_fn_gfn, "test_fn.dat");
 	finite_diff(ps, test_fn_gfn, FD_derivs_gfn, which_derivs);
 	ps.print_gridfn(FD_derivs_gfn, "FD_derivs.dat");
 	analytic_derivs(ps, analytic_derivs_gfn, which_derivs);
@@ -200,6 +218,7 @@ else if (STRING_EQUAL(which_test, "derivatives"))
 		     nominal_error_gfn);
 	ps.print_gridfn(nominal_error_gfn, "nominal_error.dat");
 	}
+
 else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
 		   "unknown which_test=\"%s\"!",
 		   which_test);					/*NOTREACHED*/
@@ -212,6 +231,173 @@ CCTK_VInfo(CCTK_THORNSTRING, "destroying patch system");
 //******************************************************************************
 
 //
+// This function tests the computation of the Jacobian of the
+// patch_system::synchronize() operation.  In outline, it does the following:
+//
+//	set up a test function in test_gfn on the nominal grid 
+//	synchronize test_gfn
+//	print ths synchronized function
+//	set up the same test function in NP_test_gfn on the nominal grid
+//		for each patch p and ghost zone pgz
+//		{
+//		compute the synchronize() Jacobian for this ghost zone
+//		   (via  pgz.compute_Jacobian()  et al)
+//			for each point x in (p,pgz)
+//			{
+//				for each point y in gz.Jacobian_y_patch()
+//				{
+//				const fp save_y_gridfn = NP_test_gfn at y
+//				NP_test_gfn at y += perturbation_amplitude
+//				synchronize NP_test_gfn
+//				NP_Jacobian = (NP_test_gfn at x - test_gfn at x)
+//					      / perturbation_amplitude
+//				NP_test_gfn at y = save_y_gridfn
+//				Jacobian = pgz.Jacobian(...)
+//				if (m in molecule || the Jacobians differ)
+//				   then print Jacobian, NP_Jacobian,
+//					      Jacobian error
+//				}
+//			}
+//		}
+//
+// Arguments:
+// ps = (in out) THe patch system in/on which to do the computations.
+// test_gfn, NP_test_gfn = The gfns of two ghosted test gridfns.
+// perturbation_amplitude = The perturbation amplitude for the NP Jacobian.
+// Jacobian_file_name = The name of the output file to which both Jacobians
+//			should be written.
+//
+namespace {
+void test_synchronize_Jacobians(patch_system& ps,
+				int test_gfn, int NP_test_gfn,
+				fp perturbation_amplitude,
+				const char Jacobian_file_name[])
+{
+setup_sym_fn_xyz(ps, test_gfn, false);
+ps.synchronize_ghost_zones(test_fn_gfn, test_fn_gfn);
+ps.print_ghosted_gridfn(test_fn_copy_gfn, "test_fn_copy.dat");
+
+setup_sym_fn_xyz(ps, NP_test_gfn, false);
+
+FILE *fileptr = fopen(Jacobian_file_name, "w");
+if (fileptr == NULL)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+"test_synchronize_Jacobians(): can't open output plot file \"%s\"!",
+		   Jacobian_file_name);				/*NOTREACHED*/
+fprintf(fileptr, "# column  1 = x patch number\n");
+fprintf(fileptr, "# column  2 = x ghost_zone is_min()\n");
+fprintf(fileptr, "# column  3 = x ghost_zone is_rho()\n");
+fprintf(fileptr, "# column  4 = x_iperp\n");
+fprintf(fileptr, "# column  5 = x_ipar\n");
+fprintf(fileptr, "# column  6 = x_irho\n");
+fprintf(fileptr, "# column  7 = x_isigma\n");
+fprintf(fileptr, "# column  8 = y patch number\n");
+fprintf(fileptr, "# column  9 = y ghost_zone is_min()\n");
+fprintf(fileptr, "# column 10 = y ghost_zone is_rho()\n");
+fprintf(fileptr, "# column 11 = y_iperp\n");
+fprintf(fileptr, "# column 12 = y_ipar\n");
+fprintf(fileptr, "# column 13 = y_irho\n");
+fprintf(fileptr, "# column 14 = y_isigma\n");
+fprintf(fileptr, "# column 15 = Jacobian\n");
+fprintf(fileptr, "# column 16 = NP_Jacobian\n");
+fprintf(fileptr, "# column 17 = Jacobian error\n");
+
+    //*** for each patch p and ghost zone pgz
+    for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+    {
+    patch& p = ps.ith_patch(pn);
+
+    // n.b. these loops must use _int_ variables for the loop
+    //      to terminate!
+    for (int want_min = false ; want_min <= true ; ++want_min)
+    {
+    for (int want_rho = false ; want_rho <= true ; ++want_rho)
+    {
+    const patch_edge& pe = p.minmax_ang_patch_edge(want_min, want_rho);
+    ghost_zone& pgz      = p.minmax_ang_ghost_zone(want_min, want_rho);
+
+    pgz.compute_Jacobian(test_gfn, test_gfn,
+			 true, true, true);	// want *all* of ghost zone
+    const int Jacobian_y_min_ipar_m = pgz.Jacobian_y_min_ipar_m();
+    const int Jacobian_y_max_ipar_m = pgz.Jacobian_y_max_ipar_m();
+
+	//*** for each point x in (p,pgz)
+	for (int x_iperp = pgz.min_iperp() ;
+	     x_iperp <= pgz.max_iperp() ;
+	     ++x_iperp)
+	{
+	// FIXME FIXME -- this includes corners when it shouldn't
+	for (int x_ipar = pgz.ghost_zone_min_ipar() ;
+	     x_ipar <= pgz.ghost_zone_max_ipar() ;
+	     ++x_ipar)
+	{
+	const int x_irho   = pe.  irho_of_iperp_ipar(x_iperp, x_ipar);
+	const int x_isigma = pe.isigma_of_iperp_ipar(x_iperp, x_ipar);
+
+	patch&            q  = pgz.Jacobian_y_patch();
+	const patch_edge& qe = pgz.Jacobian_y_edge();
+	const int Jacobian_y_iperp = pgz.Jacobian_y_iperp(x_iperp);
+	const int Jacobian_y_ipar_posn
+		= pgz.Jacobian_y_ipar_posn(x_iperp, x_ipar);
+
+	    //*** for each point y in gz.Jacobian_patch()
+	    for (int y_irho = q.min_irho() ; y_irho <= q.max_irho() ; ++y_irho)
+	    {
+	    for (int y_isigma = q.min_isigma() ;
+		 y_isigma <= q.max_isigma() ;
+		 ++y_isigma)
+	    {
+	    // compute the NP Jacobian
+	    const fp save_y_gridfn
+		= q.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma);
+	    q.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma)
+		+= perturbation_amplitude;
+
+	    ps.synchronize_ghost_zones(NP_test_gfn, NP_test_gfn);
+	    const fp NP_Jacobian
+		= (   p.ghosted_gridfn(NP_test_gfn, x_irho,x_isigma)
+		    - p.ghosted_gridfn(   test_gfn, x_irho,x_isigma) )
+		  / perturbation_amplitude;
+	    q.ghosted_gridfn(NP_test_gfn, y_irho,y_isigma) = save_y_gridfn;
+
+	    // compute the query Jacobian
+	    const int y_iperp = qe.iperp_of_irho_isigma(y_irho,y_isigma);
+	    const int y_ipar  = qe. ipar_of_irho_isigma(y_irho,y_isigma);
+	    const int y_ipar_m = y_ipar - Jacobian_y_ipar_posn;
+	    const bool m_in_molecule
+		= (y_iperp == Jacobian_y_iperp)
+		  && (y_ipar_m >= Jacobian_y_min_ipar_m)
+		  && (y_ipar_m <= Jacobian_y_max_ipar_m);
+	    const fp Jacobian = pgz.Jacobian(x_iperp, x_ipar, y_ipar_m);
+
+	    // print the results
+	    if (m_in_molecule || (Jacobian != NP_Jacobian))
+	       then fprintf(fileptr,
+"%d %d %d\t%d %d\t%d %d\t%d %d %d\t%d %d\t%d %d\t%.10g\t%.10g\t%e\n",
+			    p.patch_number(), pe.is_min(), pe.is_rho(),
+			    x_iperp, x_ipar, x_irho, x_isigma,
+			    q.patch_number(), qe.is_min(), qe.is_rho(),
+			    y_iperp, y_ipar, y_irho, y_isigma,
+			    double(Jacobian), double(NP_Jacobian),
+			    double(Jacobian-NP_Jacobian));
+	    }
+	    }
+	}
+	}
+    }
+    }
+
+    }
+
+fclose(fileptr);
+}
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
 // This function sets up the ghosted test function for the gridfn and
 // synchronize tests, symmetrizing the test function to match the
 // symmetry of the patch system.