assorted changes to start doing ghost-zone and derivative tests

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

1 files changed, 552 insertions, 6 deletions

diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc
index aefe552..62a35cd 100644
--- a/src/patch/test_patch_system.cc
+++ b/src/patch/test_patch_system.cc
@@ -2,9 +2,27 @@
 // $Id$
 
 //
+// <<<misc global data>>>
 // <<<prototypes>>>
+// <<<gfn definitions>>>
+//
 // driver - Cactus interface
 //
+/// setup_sym_fn_xyz - set up symmetrized test function fn(global_[xyz])
+/// setup_fn_rho_sigma - set up test function fn(rho,sigma)
+/// finite_diff - compute linear combination of finite differences
+/// analytic_derivs - compute linear combination of analytic derivatives
+///
+/// sym_fn_xyz - symmetrized test function fn(x,y,z) + fn(-y,x,z) + ...
+/// fn_xyz - test function fn(x,y,z)
+///
+/// gridfn_minus - compute gridfn x - gridfn y --> gridfn z
+/// ghosted_gridfn_minus - compute [ghosted] gridfn x - gridfn y --> gridfn z
+///
+/// fn_rho_sigma - test function fn(rho,sigma)
+/// finite_diff_fn - finite differences of fn(rho,sigma)
+/// analytic_deriv_fn - analytical derivs of fn(rho,sigma) (via Maple)
+//
 
 #include <stdio.h>
 #include <assert.h>
@@ -36,11 +54,80 @@ using jtutil::error_exit;
 //******************************************************************************
 
 //
+// misc global data
+//
+
+// which test are we going to do?
+static const int which_deriv_fn          = 0x1;
+static const int which_deriv_rho         = 0x2;
+static const int which_deriv_sigma       = 0x4;
+static const int which_deriv_rho_rho     = 0x8;
+static const int which_deriv_rho_sigma   = 0x10;
+static const int which_deriv_sigma_sigma = 0x20;
+
+// when testing multiple derivatives, we combine them together
+// with these "quasi-random" weights (chosen from digits of pi)
+static const fp deriv_weight_fn          = 3.14,
+		deriv_weight_rho         = 1.59,
+		deriv_weight_sigma       = 2.65,
+		deriv_weight_rho_rho     = 3.58,
+		deriv_weight_rho_sigma   = 9.79,
+		deriv_weight_sigma_sigma = 3.23;
+
+//******************************************************************************
+
+//
 // ***** prototypes *****
 //
+
 extern "C"
   void test_patch_system(CCTK_ARGUMENTS);
 
+namespace {
+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, 
+			  int ghosted_gfn_x, int ghosted_gfn_y,
+			  int ghosted_gfn_dst);
+
+fp sym_fn_xyz(enum patch_system::patch_system_type type, fp x, fp y, fp z);
+fp fn_xyz(fp x, fp y, fp z);
+
+fp fn_rho_sigma(fp rho, fp sigma);
+fp finite_diff_fn(const patch& p,
+		  int ghosted_gfn_src, int which_derivs,
+		  int irho, int isigma);
+fp analytic_deriv_fn(fp rho, fp sigma, int which_derivs);
+	  };
+
+//******************************************************************************
+
+//
+// ***** gfn definitions *****
+//
+
+// nominal gridfns
+static const int FD_derivs_gfn = -1;
+static const int analytic_derivs_gfn = -2;
+static const int nominal_error_gfn = -3;
+static const int nominal_min_gfn = -3;
+static const int nominal_max_gfn = -1;
+
+// ghosted gridfns
+static const int test_fn_gfn = 1;
+static const int test_fn_copy_gfn = 2;
+static const int ghosted_error_gfn = 3;
+static const int ghosted_min_gfn = 1;
+static const int ghosted_max_gfn = 3;
+
+//******************************************************************************
+//******************************************************************************
 //******************************************************************************
 
 //
@@ -68,18 +155,477 @@ if (interp_par_table_handle < 0)
 		   "bad interpatch-interpolator parameter(s) \"%s\"!",
 		   interpatch_interpolator_pars);		/*NOTREACHED*/
 
-const int min_gfn = 1;
-const int max_gfn = 3;
-const int ghosted_min_gfn = -2;
-const int ghosted_max_gfn = -1;
-
+//
 // create the patch system
+//
 CCTK_VInfo(CCTK_THORNSTRING, "about to create 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,
-		min_gfn, max_gfn,
+		nominal_min_gfn, nominal_max_gfn,
 		ghosted_min_gfn, ghosted_max_gfn,
 		interp_handle, interp_par_table_handle);
 CCTK_VInfo(CCTK_THORNSTRING, "patch system created ok");
+
+//
+// do the actual tests
+//
+FILE *output_file_ptr = fopen(output_file_name, "w");
+if (output_file_ptr == NULL)
+   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "can't open output file \"%s\"!",
+		   output_file_name);				/*NOTREACHED*/
+if      (STRING_EQUAL(which_test, "gridfn"))
+   then {
+	setup_sym_fn_xyz(ps, test_fn_gfn, true);
+	ps.print_ghosted_gridfn(test_fn_gfn, true, output_file_ptr);
+	}
+else if (STRING_EQUAL(which_test, "synchronize"))
+   then {
+	setup_sym_fn_xyz(ps, test_fn_gfn, false);
+	setup_sym_fn_xyz(ps, test_fn_copy_gfn, true);
+	ps.synchronize_ghost_zones(test_fn_gfn, test_fn_gfn);
+	ghosted_gridfn_minus(ps,
+			     test_fn_gfn, test_fn_copy_gfn,
+			     ghosted_error_gfn);
+	ps.print_ghosted_gridfn(ghosted_error_gfn, true, output_file_ptr);
+	}
+else if (STRING_EQUAL(which_test, "derivatives"))
+   then {
+	setup_fn_rho_sigma(ps, test_fn_gfn);
+	finite_diff(ps, test_fn_gfn, FD_derivs_gfn, which_derivs);
+	analytic_derivs(ps, analytic_derivs_gfn, which_derivs);
+	gridfn_minus(ps,
+		     FD_derivs_gfn, analytic_derivs_gfn,
+		     nominal_error_gfn);
+	ps.print_gridfn(nominal_error_gfn, output_file_ptr);
+	}
+else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "unknown which_test=\"%s\"!",
+		   which_test);					/*NOTREACHED*/
+fclose(output_file_ptr);
+
+CCTK_VInfo(CCTK_THORNSTRING, "destroying patch system");
+}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// 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.
+//
+// Arguments:
+// ps = The patch system.
+// ghosted_gfn = Specifies the gridfn to set up.
+// want_ghost_zones = true ==> Set up on ghosted part of ghosted grid
+//		      false ==> Set up on nominal part of ghosted grid
+//
+namespace {
+void setup_sym_fn_xyz(patch_system& ps, int ghosted_gfn, bool want_ghost_zones)
+{
+CCTK_VInfo(CCTK_THORNSTRING, "setting up ghosted test fn(x,y,z)");
+CCTK_VInfo(CCTK_THORNSTRING, "   on %s part of ghosted grid",
+	   want_ghost_zones ? "ghosted" : "nominal");
+
+	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+	{
+	patch& p = ps.ith_patch(pn);
+
+		for (int irho = p.effective_min_irho(want_ghost_zones) ;
+		     irho <= p.effective_max_irho(want_ghost_zones) ;
+		     ++irho)
+		{
+		for (int isigma = p.effective_min_isigma(want_ghost_zones) ;
+		     isigma <= p.effective_max_isigma(want_ghost_zones) ;
+		     ++isigma)
+		{
+		const fp rho = p.rho_of_irho(irho);
+		const fp sigma = p.sigma_of_isigma(isigma);
+
+		fp local_x, local_y, local_z;
+		p.xyz_of_r_rho_sigma(1.0, rho, sigma,
+				     local_x, local_y, local_z);
+
+		p.ghosted_gridfn(ghosted_gfn, irho,isigma)
+			= sym_fn_xyz(ps.type(), local_x, local_y, local_z);
+		}
+		}
+	}
+}
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function sets up the test function for the finite differencing
+// tests.
+//
+// Arguments:
+// ps = The patch system.
+// ghosted_gfn = Specifies the gridfn to set up.
+//
+namespace {
+void setup_fn_rho_sigma(patch_system& ps, int ghosted_gfn)
+{
+CCTK_VInfo(CCTK_THORNSTRING, "setting up ghosted test fn(rho,sigma)");
+
+	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+	{
+	patch& p = ps.ith_patch(pn);
+
+		for (int irho = p.ghosted_min_irho() ;
+		     irho <= p.ghosted_max_irho() ;
+		     ++irho)
+		{
+		for (int isigma = p.ghosted_min_isigma() ;
+		     isigma <= p.ghosted_max_isigma() ;
+		     ++isigma)
+		{
+		const fp rho = p.rho_of_irho(irho);
+		const fp sigma = p.sigma_of_isigma(isigma);
+
+		p.ghosted_gridfn(ghosted_gfn, irho,isigma)
+			= fn_rho_sigma(rho, sigma);
+		}
+		}
+	}
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function computes (on the nominal grid only) the specified
+// linear combination of finite derivatives of a test function.
+//
+// Arguments:
+// ps = The patch system.
+// gfn_src = Specifies the gridfn to finite difference.
+// gfn_dst = Specifies the gridfn in which to store the result.
+// which = Specifies which finite derivatives to include in the test.
+//
+namespace {
+void finite_diff(patch_system& ps,
+		 int ghosted_gfn_src, int gfn_dst,
+		 int which_derivs)
+{
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "finite differencing ghosted_gfn=%d --> gfn=%d (which_derivs=%d)",
+	   ghosted_gfn_src, gfn_dst, which_derivs);
+
+	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)
+	{
+	p.gridfn(gfn_dst, irho,isigma)
+		= finite_diff_fn(p,
+				 ghosted_gfn_src, which_derivs,
+				 irho,isigma);
+	}
+	}
+
+	}
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function computes the specified linear combination of
+// derivatives of the test function, everywhere on the nominal grid.
+//
+namespace {
+void analytic_derivs(patch_system& ps,
+		     int gfn_dst,
+		     int which_derivs)
+{
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "computing analytic derivatives (which_derivs=%d)",
+	   which_derivs);
+
+	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);
+	p.gridfn(gfn_dst, irho,isigma)
+		= analytic_deriv_fn(rho,sigma, which_derivs);
+	}
+	}
+
+	}
+}
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function computes [nominal]  gridfn_x - gridfn_y --> gridfn_z .
+//
+namespace {
+void gridfn_minus(patch_system& ps, 
+		  int gfn_x, int gfn_y, int gfn_dst)
+{
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "[nominal] gfn=%d - gfn=%d --> gfn=%d\n",
+	   gfn_x, gfn_y, gfn_dst);
+
+	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)
+	{
+	p.gridfn(gfn_dst, irho,isigma)
+		=   p.gridfn(gfn_x, irho,isigma)
+		  - p.gridfn(gfn_y, irho,isigma);
+	}
+	}
+
+	}
+}
+	  }
+
+//******************************************************************************
+
+//
+// This function computes [ghosted]  gridfn_x - gridfn_y --> gridfn_z .
+//
+namespace {
+void ghosted_gridfn_minus(patch_system& ps, 
+			  int ghosted_gfn_x, int ghosted_gfn_y,
+			  int ghosted_gfn_dst)
+{
+CCTK_VInfo(CCTK_THORNSTRING,
+	   "[ghosted] gfn=%d - gfn=%d --> gfn=%d\n",
+	   ghosted_gfn_x, ghosted_gfn_y, ghosted_gfn_dst);
+
+	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
+	{
+	patch& p = ps.ith_patch(pn);
+
+		for (int irho = p.ghosted_min_irho() ;
+		     irho <= p.ghosted_max_irho() ;
+		     ++irho)
+		{
+		for (int isigma = p.ghosted_min_isigma() ;
+		     isigma <= p.ghosted_max_isigma() ;
+		     ++isigma)
+		{
+		p.ghosted_gridfn(ghosted_gfn_dst, irho,isigma)
+			=   p.ghosted_gridfn(ghosted_gfn_x, irho,isigma)
+			  - p.ghosted_gridfn(ghosted_gfn_y, irho,isigma);
+		}
+		}
+
+	}
+}
+	  }
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This function symmetrizes  fn_xyz()  (about the origin) to match
+// the patch system's symmetries.
+//
+// To rotate f(x,y) by 90, 180, or 270 degrees:
+//
+//	      (-y,x) |
+//	             |
+//	             |
+//	             |
+//	             |        (x,y)
+//	             |
+//	-------------+-------------
+//	             |
+//	(-x,-y)      |
+//	             |
+//	             |
+//	             |
+//	             | (y,-x)
+//
+namespace {
+fp sym_fn_xyz(enum patch_system::patch_system_type type, fp x, fp y, fp z)
+{
+switch	(type)
+	{
+case patch_system::full_sphere_patch_system:
+	return fn_xyz(x,y,z);
+	break;
+case patch_system::plus_z_hemisphere_patch_system:
+	return fn_xyz(x,y,+z) + fn_xyz(x,y,-z);
+	break;
+case patch_system::plus_xy_quadrant_patch_system:
+	return   fn_xyz(+x,+y,z)
+	       + fn_xyz(-y,+x,z)
+	       + fn_xyz(-x,-y,z)
+	       + fn_xyz(+y,-x,z);
+	break;
+case patch_system::plus_xz_quadrant_patch_system:
+	return   fn_xyz(+x,+y,+z) + fn_xyz(-x,-y,+z)
+	       + fn_xyz(+x,+y,-z) + fn_xyz(-x,-y,-z);
+	break;
+case patch_system::plus_xyz_octant_patch_system:
+	return   fn_xyz(+x,+y,+z) + fn_xyz(+x,+y,-z)
+	       + fn_xyz(-y,+x,+z) + fn_xyz(-y,+x,-z)
+	       + fn_xyz(-x,-y,+z) + fn_xyz(-x,-y,-z)
+	       + fn_xyz(+y,-x,+z) + fn_xyz(+y,-x,-z);
+	break;
+default:
+	error_exit(PANIC_EXIT,
+"***** sym_fn_xyz(): impossible type=(int)%d!\n",
+		   int(type));					/*NOTREACHED*/
+			}
+}
+	  };
+
+//******************************************************************************
+
+//
+// This is the underlying test function for our function and ghost-zone tests.
+//
+namespace {
+fp fn_xyz(fp x, fp y, fp z)
+{
+return (x*(x+0.238) + 2.417*y*(y-0.917) + 1.38*z*(z-0.472))
+       * tanh(jtutil::pow3(cos(z)));
+}
+	  };
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
+// This is the underlying test function for our finite differencing
+// tests.
+//
+namespace {
+fp fn_rho_sigma(fp rho, fp sigma)
+{
+return exp(sin(1.38*rho)) * tanh(0.17+0.83*jtutil::pow2(sin(sigma)));
+}
+	  };
+
+//******************************************************************************
+
+//
+// This function computes the sum of our linear combination of various
+// finite difference approximations to derivatives of fn_rho_sigma().
+//
+// Arguments:
+// ghosted_gfn_src = Specifies the gridfn to finite difference.
+//
+namespace {
+fp finite_diff_fn(const patch& p,
+		  int ghosted_gfn_src, int which_derivs,
+		  int irho, int isigma)
+{
+fp sum = 0.0;
+
+if (which_derivs & which_deriv_fn)
+   then sum += deriv_weight_fn
+	       * p.gridfn(ghosted_gfn_src, irho,isigma);
+
+if (which_derivs & which_deriv_rho)
+   then sum += deriv_weight_rho
+	       * p.partial_rho(ghosted_gfn_src, irho,isigma);
+if (which_derivs & which_deriv_sigma)
+   then sum += deriv_weight_sigma
+	       * p.partial_sigma(ghosted_gfn_src, irho,isigma);
+
+if (which_derivs & which_deriv_rho_rho)
+   then sum += deriv_weight_rho_rho
+	       * p.partial_rho_rho(ghosted_gfn_src, irho,isigma);
+if (which_derivs & which_deriv_rho_sigma)
+   then sum += deriv_weight_rho_sigma
+	       * p.partial_rho_sigma(ghosted_gfn_src, irho,isigma);
+if (which_derivs & which_deriv_sigma_sigma)
+   then sum += deriv_weight_sigma_sigma
+	       * p.partial_sigma_sigma(ghosted_gfn_src, irho,isigma);
+
+return sum;
+}
+	  };
+
+//******************************************************************************
+
+//
+// This function computes the sum of a specified linear combination of
+// various (analytical) derivatives of fn_rho_sigma().
+//
+// The derivatives were machine-generated via Maple's codegen[C]() function:
+// "deriv_patch_system.maple" is the Maple input
+// "deriv_patch_system.out" is the Maple input; code here is cut-n-pasted
+//			   from the Maple codegen[C]() output there
+//
+namespace {
+fp analytic_deriv_fn(fp rho, fp sigma, int which_derivs)
+{
+fp sum = 0.0;
+
+if (which_derivs & which_deriv_fn)
+   then sum += deriv_weight_fn
+	       * fn_rho_sigma(rho,sigma);
+
+if (which_derivs & which_deriv_rho)
+   then sum += deriv_weight_rho
+	       * 0.138E1*cos(0.138E1*rho)*exp(sin(0.138E1*rho))
+		 *tanh(0.17+0.83*pow(sin(sigma),2.0));
+if (which_derivs & which_deriv_sigma)
+   then sum += deriv_weight_sigma
+	       * 0.166E1*exp(sin(0.138E1*rho))
+		 *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
+		 *sin(sigma)*cos(sigma);
+
+if (which_derivs & which_deriv_rho_rho)
+   then sum += deriv_weight_rho_rho
+	       * (
+		 -0.19044E1*sin(0.138E1*rho)*exp(sin(0.138E1*rho))
+		  *tanh(0.17+0.83*pow(sin(sigma),2.0))
+		 +0.19044E1*pow(cos(0.138E1*rho),2.0)*exp(sin(0.138E1*rho))
+		  *tanh(0.17+0.83*pow(sin(sigma),2.0))
+		 );
+if (which_derivs & which_deriv_rho_sigma)
+   then sum += deriv_weight_rho_sigma
+	       * 0.22908E1*cos(0.138E1*rho)*exp(sin(0.138E1*rho))
+		 *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
+		 *sin(sigma)*cos(sigma);
+if (which_derivs & which_deriv_sigma_sigma)
+   then sum += deriv_weight_sigma_sigma
+	       * (
+		 -0.55112E1*exp(sin(0.138E1*rho))
+		  *tanh(0.17+0.83*pow(sin(sigma),2.0))
+		  *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
+		  *pow(sin(sigma),2.0)*pow(cos(sigma),2.0)
+		 +0.166E1*exp(sin(0.138E1*rho))
+		  *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
+		  *pow(cos(sigma),2.0)
+		 -0.166E1*exp(sin(0.138E1*rho))
+		  *(1.0-pow(tanh(0.17+0.83*pow(sin(sigma),2.0)),2.0))
+		  *pow(sin(sigma),2.0)
+		 );
+
+return sum;
 }
+	  };