* redo min_ipar() and max_ipar() functions in ghost_zone::

  --> now these are virtual and take iperp as an argument
      (ignored for  symmetry_ghost_zone)
  --> new functions  extreme_min_ipar()  and  extreme_max_ipar()
      to give the extreme range (may even be a bit conservative),
      i.e. the range including all of the corners
* fix a bug in   patch_interp::molecule_minmax_ipar_m()
  where we mistakenly returned the min ipar m for both the min and max
  returns :(

* finish modifying the   src/util/test_patch_system.cc   test driver
  to test the Jacobian computation


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

4 files changed, 45 insertions, 50 deletions

diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc
index 3197e2e..ab40101 100644
--- a/src/patch/ghost_zone.cc
+++ b/src/patch/ghost_zone.cc
@@ -60,7 +60,7 @@ iperp_map_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 				     my_edge_in.fp_grid_outer_iperp());
 
 // ipar_map_: identity map
-ipar_map_ = new jtutil::cpm_map<fp>(min_ipar(), max_ipar());
+ipar_map_ = new jtutil::cpm_map<fp>(extreme_min_ipar(), extreme_max_ipar());
 }
 
 //******************************************************************************
@@ -96,7 +96,7 @@ iperp_map_ = new jtutil::cpm_map<fp>(min_iperp(), max_iperp(),
 //
 // parallel map
 //
-ipar_map_ = new jtutil::cpm_map<fp>(min_ipar(), max_ipar(),
+ipar_map_ = new jtutil::cpm_map<fp>(extreme_min_ipar(), extreme_max_ipar(),
 				    my_edge_sample_ipar,
 				    symmetry_edge_sample_ipar,
 				    ipar_map_is_plus);
@@ -130,7 +130,7 @@ void symmetry_ghost_zone::synchronize(int ghosted_min_gfn, int ghosted_max_gfn,
 	{
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
-	for (int ipar = min_ipar() ; ipar <= max_ipar() ; ++ipar)
+	for (int ipar = min_ipar(iperp) ; ipar <= max_ipar(iperp) ; ++ipar)
 	{
 	// do we want to do this point?
 	if (! my_edge().ipar_is_in_selected_corner(want_min_par_corner,
@@ -396,24 +396,17 @@ max_other_iperp_ = std::max(other_iperp(min_iperp()), other_iperp(max_iperp()));
 
 
 //
-// compute extreme range of ipar that we'll use,
-// and set up arrays giving actual [min,max] ipar that we'll use
-//     at each other_iperp
+// set up arrays giving actual [min,max] ipar that we'll use
+// at each other_iperp
 // ... we will pass these arrays by reference to the other patch's
 //     patch_interp object, with ipar being parindex there
 //
-extreme_min_ipar_ = INT_MAX;
-extreme_max_ipar_ = INT_MIN;
 min_ipar_used_ = new jtutil::array1d<int>(min_other_iperp_, max_other_iperp_);
 max_ipar_used_ = new jtutil::array1d<int>(min_other_iperp_, max_other_iperp_);
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
-	const int min_ipar_here = min_ipar(iperp);
-	const int max_ipar_here = max_ipar(iperp);
-	extreme_min_ipar_ = std::min(extreme_min_ipar_, min_ipar_here);
-	extreme_max_ipar_ = std::max(extreme_max_ipar_, max_ipar_here);
-	(*min_ipar_used_)(other_iperp(iperp)) = min_ipar_here;
-	(*max_ipar_used_)(other_iperp(iperp)) = max_ipar_here;
+	(*min_ipar_used_)(other_iperp(iperp)) = min_ipar(iperp);
+	(*max_ipar_used_)(other_iperp(iperp)) = max_ipar(iperp);
 	}
 
 
@@ -422,7 +415,7 @@ max_ipar_used_ = new jtutil::array1d<int>(min_other_iperp_, max_other_iperp_);
 //
 
 other_par_ = new jtutil::array2d<fp>(min_other_iperp_, max_other_iperp_,
-				     extreme_min_ipar_, extreme_max_ipar_);
+				     extreme_min_ipar(), extreme_max_ipar());
 
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
@@ -454,7 +447,7 @@ interp_result_buffer_
 	= new jtutil::array3d<fp>(my_patch().ghosted_min_gfn(),
 					  my_patch().ghosted_max_gfn(),
 				  min_other_iperp_, max_other_iperp_,
-				  extreme_min_ipar_, extreme_max_ipar_);
+				  extreme_min_ipar(), extreme_max_ipar());
 
 //
 // construct the patch_interp object to interpolate from the *other* patch
@@ -544,23 +537,21 @@ other_patch_interp_->interpolate(ghosted_min_gfn, ghosted_max_gfn,
 				 *interp_result_buffer_);
 
 // store the results back into our gridfns
-	for (int gfn = ghosted_min_gfn ; gfn <= ghosted_max_gfn ; ++gfn)
-	{
+    for (int gfn = ghosted_min_gfn ; gfn <= ghosted_max_gfn ; ++gfn)
+    {
 	for (int iperp = min_iperp() ; iperp <= max_iperp() ; ++iperp)
 	{
 	const int oiperp = other_iperp(iperp);
 
-		for (int ipar = min_ipar(iperp) ;
-		     ipar <= max_ipar(iperp) ;
-		     ++ipar)
-		{
-		int irho   = my_edge().  irho_of_iperp_ipar(iperp,ipar);
-		int isigma = my_edge().isigma_of_iperp_ipar(iperp,ipar);
-		my_patch().ghosted_gridfn(gfn, irho,isigma)
-			= (*interp_result_buffer_)(gfn, oiperp,ipar);
-		}
+	for (int ipar = min_ipar(iperp) ; ipar <= max_ipar(iperp) ; ++ipar)
+	{
+	int irho   = my_edge().  irho_of_iperp_ipar(iperp,ipar);
+	int isigma = my_edge().isigma_of_iperp_ipar(iperp,ipar);
+	my_patch().ghosted_gridfn(gfn, irho,isigma)
+		= (*interp_result_buffer_)(gfn, oiperp,ipar);
 	}
 	}
+    }
 }
 
 //******************************************************************************
@@ -594,13 +585,13 @@ other_patch_interp_->molecule_minmax_ipar_m(min_y_ipar_m_, max_y_ipar_m_);
 if (Jacobian_y_ipar_posn_ == NULL)
    then Jacobian_y_ipar_posn_ = new jtutil::array2d<CCTK_INT>
 				     (min_other_iperp_, max_other_iperp_,
-				      extreme_min_ipar_, extreme_max_ipar_);
+				      extreme_min_ipar(), extreme_max_ipar());
 other_patch_interp_->molecule_posn(*Jacobian_y_ipar_posn_);
 
 if (Jacobian_buffer_ == NULL)
    then Jacobian_buffer_ = new jtutil::array3d<fp>
 				 (min_other_iperp_, max_other_iperp_,
-				  extreme_min_ipar_, extreme_max_ipar_,
+				  extreme_min_ipar(), extreme_max_ipar(),
 				  min_y_ipar_m_, max_y_ipar_m_);
 other_patch_interp_->Jacobian(*Jacobian_buffer_);
 }
diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh
index 1cf1517..6b81ee2 100644
--- a/src/patch/ghost_zone.hh
+++ b/src/patch/ghost_zone.hh
@@ -199,6 +199,9 @@
 //	or other elliptic stuff (it would matter for stability in
 //	time evolutions using multipatch finite differencing).]
 //
+// In general, then, a ghost zone is trapezoid-shaped: iperp lies in a
+// fixed interval, while ipar lies in an interval which may depend on iperp.
+//
 
 //*****************************************************************************
 
@@ -249,14 +252,14 @@ public:
 	// operation into internal buffers; the following functions
 	// provide access to that Jacobian.
 	//
+	// FIXME: should these be moved out into a separate Jacobian
+	//        object/class?
+	//
 	// n.b. terminology is
 	//	partial gridfn at x
 	//	-------------------
 	//	partial gridfn at y
 	//
-	// FIXME: should these be moved out into a separate Jacobian
-	//        object/class?
-	//
 	virtual void compute_Jacobian(int ghosted_min_gfn, int ghosted_max_gfn,
 				      bool want_min_par_corner,
 				      bool want_noncorner,
@@ -354,13 +357,18 @@ public:
 	int outer_iperp() const
 		{ return is_min() ? min_iperp() : max_iperp(); }
 
-	// min/max ipar that might possibly be part of this ghost zone
+	// extreme min/max ipar that might possibly be part of this ghost zone
 	// (derived classes may actually use a subset of this)
-	int ghost_zone_min_ipar() const
+	int extreme_min_ipar() const
 		{ return my_edge().min_ipar_with_corners(); }
-	int ghost_zone_max_ipar() const
+	int extreme_max_ipar() const
 		{ return my_edge().max_ipar_with_corners(); }
 
+	// actual min/max ipar in the ghost zone at a particular iperp
+	// (only defined in derived classes)
+	virtual int min_ipar(int iperp) const = 0;
+	virtual int max_ipar(int iperp) const = 0;
+
 	// assert() that ghost zone is fully setup:
 	// defined here ==> no-op
 	// symmetry ghost zone ==> unchanged ==> no-op
@@ -526,11 +534,11 @@ public:
 	fp sign_of_ipar_map() const
 		{ return ipar_map_->fp_sign(); }
 
-	// min/max/size ipar of the ghost zone
+	// min/max ipar of the ghost zone
 	// ... we always include the corners
 	//     (cf. the example at the start of this file)
-	int min_ipar() const { return my_edge().min_ipar_with_corners(); }
-	int max_ipar() const { return my_edge().max_ipar_with_corners(); }
+	int min_ipar(int iperp) const { return extreme_min_ipar(); }
+	int max_ipar(int iperp) const { return extreme_max_ipar(); }
 
 	// constructor for mirror-symmetry ghost zone
 	symmetry_ghost_zone(const patch_edge& my_edge_in);
@@ -701,8 +709,8 @@ private:
 						   bool want_max_par_corner)
 		const;
 
-	// min/max/size ipar of the ghost zone for specified iperp
-	// taking into account how we treat the corners
+	// min/max ipar of the ghost zone for specified iperp
+	// taking into account our "triangular" corners
 	// (cf. the example at the start of this file)
 	int min_ipar(int iperp) const;
 	int max_ipar(int iperp) const;
@@ -752,9 +760,6 @@ private:
 	jtutil::array1d<int>* min_ipar_used_;
 	jtutil::array1d<int>* max_ipar_used_;
 
-	// extreme values of [min,max]_ipar used at each other_iperp
-	int extreme_min_ipar_, extreme_max_ipar_;
-
 	// other patch's (fp) parallel coordinates of our ghost zone points
 	// ... we will pass this array by reference
 	//     to the other patch's patch_interp object
diff --git a/src/patch/patch_interp.cc b/src/patch/patch_interp.cc
index e23ce04..0c7181f 100644
--- a/src/patch/patch_interp.cc
+++ b/src/patch/patch_interp.cc
@@ -466,7 +466,7 @@ if ((status1 < 0) || (status2 < 0))
 ,
 		   status1, status2);				/*NOTREACHED*/
 min_ipar_m = molecule_min_m;
-max_ipar_m = molecule_min_m;
+max_ipar_m = molecule_max_m;
 
 //
 // delete Jacobian-sparsity-pattern query entries from the parameter table
diff --git a/src/patch/test_patch_system.cc b/src/patch/test_patch_system.cc
index 5143f26..6371305 100644
--- a/src/patch/test_patch_system.cc
+++ b/src/patch/test_patch_system.cc
@@ -201,14 +201,14 @@ else if (STRING_EQUAL(which_test, "synchronize"))
 
 else if (STRING_EQUAL(which_test, "synchronize Jacobian"))
    then test_synchronize_Jacobians(ps,
-				   test_fn_copy_gfn, test_fn_gfn,
+				   test_fn_gfn, test_fn_copy_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(test_fn_gfn, "test_fn.dat");
+	ps.print_ghosted_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);
@@ -275,7 +275,7 @@ void test_synchronize_Jacobians(patch_system& ps,
 {
 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");
+ps.print_ghosted_gridfn(test_fn_gfn, "test_fn.dat");
 
 setup_sym_fn_xyz(ps, NP_test_gfn, false);
 
@@ -326,9 +326,8 @@ fprintf(fileptr, "# column 17 = Jacobian error\n");
 	     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() ;
+	for (int x_ipar = pgz.min_ipar(x_iperp) ;
+	     x_ipar <= pgz.max_ipar(x_iperp) ;
 	     ++x_ipar)
 	{
 	const int x_irho   = pe.  irho_of_iperp_ipar(x_iperp, x_ipar);