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, 92 insertions, 68 deletions

diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh
index 70b7f3b..1cf1517 100644
--- a/src/patch/ghost_zone.hh
+++ b/src/patch/ghost_zone.hh
@@ -249,6 +249,11 @@ public:
 	// operation into internal buffers; the following functions
 	// provide access to that Jacobian.
 	//
+	// n.b. terminology is
+	//	partial gridfn at x
+	//	-------------------
+	//	partial gridfn at y
+	//
 	// FIXME: should these be moved out into a separate Jacobian
 	//        object/class?
 	//
@@ -262,46 +267,44 @@ public:
 	// The API in the remaining functions implicitly assumes that
 	// the Jacobian is independent of  ghosted_gfn , and also that
 	// the structure of the Jacobian is "simple", in that the set of
-	// points (with nonzero Jacobian values) in a single row of the
+	// y points (with nonzero Jacobian values) in a single row of the
 	// Jacobian matrix (i.e. the set of points on which a single
 	// ghost-zone point depends),
-	// - lies entirely within a single patch
-	// - has the same number of points for every point in this ghost zone
-	// - has a single oiperp value (= iperp value in their patch)
-	//   (depending on our iperp, of course)
-	// - have a contiguous interval of oipar (= ipar in their patch)
-	//   (depending on our iperp and ipar, of course), whose
-	//   size is
+	// - lies entirely within a single y patch
+	// - has the same number of points for every x point in this ghost zone
+	// - has a single yiperp value (depending on our iperp, of course)
+	// - have a contiguous interval of yipar (depending on our iperp
+	//and ipar, of course), whose size is
 	//	[or can be taken to be without an unreasonable
 	//	amount of zero-padding]
 	//   independent of our iperp and ipar; we parameterize this
-	//   interval as  oipar = posn+m  where  posn  is determined by
+	//   interval as  yipar = posn+m  where  posn  is determined by
 	//   our iperp and ipar, and  m  has a fixed range independent
 	//   of our iperp and ipar
 	//
 
-	// to which patch/edge do the points in this Jacobian row belong?
-	virtual const patch& Jacobian_patch() const = 0;
-	virtual const patch_edge& Jacobian_edge() const = 0;
+	// to which patch/edge do the y points in this Jacobian row belong?
+	virtual patch&            Jacobian_y_patch() const = 0;
+	virtual const patch_edge& Jacobian_y_edge()  const = 0;
 
 	// what is the [min,max] range of m for this ghost zone?
-	virtual int Jacobian_min_oipar_m() const = 0;
-	virtual int Jacobian_max_oipar_m() const = 0;
+	virtual int Jacobian_y_min_ipar_m() const = 0;
+	virtual int Jacobian_y_max_ipar_m() const = 0;
 
-	// what is the oiperp of the Jacobian points (= iperp in their patch)?
-	virtual int Jacobian_oiperp(int iperp) const = 0;
+	// what is the iperp of the Jacobian y points in their (y) patch?
+	virtual int Jacobian_y_iperp(int x_iperp) const = 0;
 
-	// what is the  posn  value of the points in this Jacobian row?
-	virtual int Jacobian_oipar_posn(int iperp, int ipar) const = 0;
+	// what is the  posn  value of the y points in this Jacobian row?
+	virtual int Jacobian_y_ipar_posn(int x_iperp, int x_ipar) const = 0;
 
 	// what is the Jacobian
-	//	    partial synchronize gridfn(ghosted_gfn, iperp, ipar)
-	//	------------------------------------------------------------
-	//	partial other patch gridfn(ghosted_gfn, oiperp, posn+ipar_m)
+	//	partial synchronize() px.gridfn(ghosted_gfn, x_iperp, x_ipar)
+	//	-------------------------------------------------------------
+	//	   partial py.gridfn(ghosted_gfn, y_iperp, y_posn+y_ipar_m)
 	// where
-	//	oiperp = Jacobian_oiperp(iperp)
-	//	posn = Jacobian_oipar_posn(iperp, ipar)
-	virtual fp Jacobian(int iperp, int ipar, int ipar_m) const = 0;
+	//	y_iperp = Jacobian_y_iperp(x_iperp)
+	//	y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar)
+	virtual fp Jacobian(int x_iperp, int x_ipar, int y_ipar_m) const = 0;
 
 
 public:
@@ -466,6 +469,11 @@ public:
 	//
 	// ***** Jacobian of synchronize() *****
 	//
+	// n.b. terminology is
+	//      partial gridfn at x
+	//      -------------------
+	//      partial gridfn at y
+	//
 
 	// allocate internal buffers, compute Jacobian
 	// ... this function is a no-op in this class
@@ -475,30 +483,31 @@ public:
 			      bool want_max_par_corner)
 		{ }
 
-	// to which patch/edge do the points in this Jacobian row belong?
-	const patch& Jacobian_patch() const { return symmetry_patch_; }
-	const patch_edge& Jacobian_edge() const { return symmetry_edge_; }
+	// to which patch/edge do the y points in this Jacobian row belong?
+	patch&            Jacobian_y_patch() const { return symmetry_patch_; }
+	const patch_edge& Jacobian_y_edge()  const { return symmetry_edge_; }
 
 	// what is the [min,max] range of m for this ghost zone?
-	int Jacobian_min_oipar_m() const { return 0; }
-	int Jacobian_max_oipar_m() const { return 0; }
+	int Jacobian_y_min_ipar_m() const { return 0; }
+	int Jacobian_y_max_ipar_m() const { return 0; }
 
 	// what is the oiperp of the Jacobian points (= iperp in their patch)?
-	virtual int Jacobian_oiperp(int iperp) const
-		{ return iperp_map_of_iperp(iperp); }
+	virtual int Jacobian_y_iperp(int x_iperp) const
+		{ return iperp_map_of_iperp(x_iperp); }
 
 	// what is the  posn  value of the points in this Jacobian row?
-	int Jacobian_oipar_posn(int iperp, int ipar) const
-		{ return ipar_map_of_ipar(ipar); }
+	int Jacobian_y_ipar_posn(int x_iperp, int x_ipar) const
+		{ return ipar_map_of_ipar(x_ipar); }
 
 	// what is the Jacobian
-	//	    partial synchronize gridfn(ghosted_gfn, iperp, ipar)
-	//	------------------------------------------------------------
-	//	partial other patch gridfn(ghosted_gfn, oiperp, posn+ipar_m)
+	//	partial synchronize() px.gridfn(ghosted_gfn, x_iperp, x_ipar)
+	//	-------------------------------------------------------------
+	//	   partial py.gridfn(ghosted_gfn, y_iperp, y_posn+y_ipar_m)
 	// where
-	//	oiperp = Jacobian_oiperp(iperp)
-	//	posn = Jacobian_oipar_posn(iperp, ipar)
-	fp Jacobian(int iperp, int ipar, int ipar_m) const { return 1.0; }
+	//	y_iperp = Jacobian_y_iperp(x_iperp)
+	//	y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar)
+	fp Jacobian(int x_iperp, int x_ipar, int y_ipar_m) const
+		{ return (y_ipar_m == 0) ? 1.0 : 0.0; }
 
 
 	//
@@ -506,7 +515,7 @@ public:
 	//
 
 	// symmetry-map coordinates
-	const patch& symmetry_patch() const { return symmetry_patch_; }
+	patch& symmetry_patch() const { return symmetry_patch_; }
 	const patch_edge& symmetry_edge() const { return symmetry_edge_; }
 	int iperp_map_of_iperp(int iperp) const
 		{ return iperp_map_->map(iperp); }
@@ -544,7 +553,7 @@ private:
 
 private:
 	// symmetry mapping --> interior of which patch?  which edge?
-	const patch& symmetry_patch_;
+	patch& symmetry_patch_;
 	const patch_edge& symmetry_edge_;
 
 	// symmetry mappings for (iperp,ipar)
@@ -602,6 +611,11 @@ public:
 	//
 	// ***** Jacobian of synchronize() *****
 	//
+	// n.b. terminology is
+	//      partial gridfn at x
+	//      -------------------
+	//      partial gridfn at y
+	//
 
 	// allocate internal buffers, compute Jacobian
 	//
@@ -613,36 +627,43 @@ public:
 			      bool want_noncorner,
 			      bool want_max_par_corner);
 
-	// to which patch/edge do the points in this Jacobian row belong?
-	const patch& Jacobian_patch() const { return other_patch_; }
-	const patch_edge& Jacobian_edge() const { return other_edge_; }
+	// to which patch/edge do the y points in this Jacobian row belong?
+	patch&            Jacobian_y_patch() const { return other_patch_; }
+	const patch_edge& Jacobian_y_edge()  const { return other_edge_; }
 
 	// what is the [min,max] range of m for this ghost zone?
-	int Jacobian_min_oipar_m() const { return min_oipar_m_; }
-	int Jacobian_max_oipar_m() const { return max_oipar_m_; }
+	int Jacobian_y_min_ipar_m() const { return min_y_ipar_m_; }
+	int Jacobian_y_max_ipar_m() const { return max_y_ipar_m_; }
 
-	// what is the oiperp of the Jacobian points (= iperp in their patch)?
-	int Jacobian_oiperp(int iperp) const { return other_iperp(iperp); }
+	// what is the iperp of the Jacobian y points in their (y) patch?
+	// ... the ipar row of grid points is actually the same, so
+	//     we just have to translate x_iperp to the y patch's coordinates
+	int Jacobian_y_iperp(int x_iperp) const { return other_iperp(x_iperp); }
 
-	// what is the  posn  value of the points in this Jacobian row?
-	int Jacobian_oipar_posn(int iperp, int ipar) const
+	// what is the  posn  value of the y points in this Jacobian row?
+	int Jacobian_y_ipar_posn(int x_iperp, int x_ipar) const
 		{
-		assert(Jacobian_oipar_posn_ != NULL);
-		const int oiperp = other_iperp(iperp);
-		return (*Jacobian_oipar_posn_)(oiperp, ipar);
+		assert(Jacobian_y_ipar_posn_ != NULL);
+		const int y_iperp = Jacobian_y_iperp(x_iperp);
+		return (*Jacobian_y_ipar_posn_)(y_iperp, x_ipar);
 		}
 
 	// what is the Jacobian
-	//	    partial synchronize gridfn(ghosted_gfn, iperp, ipar)
-	//	------------------------------------------------------------
-	//	partial other patch gridfn(ghosted_gfn, oiperp, posn+ipar_m)
+	//	partial synchronize() px.gridfn(ghosted_gfn, x_iperp, x_ipar)
+	//	-------------------------------------------------------------
+	//	   partial py.gridfn(ghosted_gfn, y_iperp, y_posn+y_ipar_m)
 	// where
-	//	posn = Jacobian_oipar_posn(iperp, ipar)
-	fp Jacobian(int iperp, int ipar, int ipar_m) const
+	//	y_iperp = Jacobian_y_iperp(x_iperp)
+	//	y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar)
+	fp Jacobian(int x_iperp, int x_ipar, int y_ipar_m) const
 		{
 		assert(Jacobian_buffer_ != NULL);
-		const int oiperp = other_iperp(iperp);
-		return (*Jacobian_buffer_)(oiperp, ipar, ipar_m);
+		if ((y_ipar_m >= min_y_ipar_m_) && (y_ipar_m <= max_y_ipar_m_))
+		   then {
+			const int y_iperp = Jacobian_y_iperp(x_iperp);
+			return (*Jacobian_buffer_)(y_iperp, x_ipar, y_ipar_m);
+			}
+		   else return 0.0;
 		}
 
 	//
@@ -650,10 +671,8 @@ public:
 	//
 
 	// basic connectivity info
-	const patch& other_patch() const
-		{ return other_patch_; }
-	const patch_edge& other_edge() const
-		{ return other_edge_; }
+	patch& other_patch() const { return other_patch_; }
+	const patch_edge& other_edge() const { return other_edge_; }
 
 	// check consistency of this and the other patch's ghost zones
 	// and patch_interp objects
@@ -754,13 +773,18 @@ private:
 	//
 	// stuff computed by  compute_Jacobian()
 	//
-	int min_oipar_m_, max_oipar_m_;
+	// n.b. terminology is
+	//      partial gridfn at x
+	//      -------------------
+	//      partial gridfn at y
+	//
+	int min_y_ipar_m_, max_y_ipar_m_;
 
-	// other patch's ipar posn for a Jacobian row
+	// other patch's y ipar posn for a Jacobian row
 	// ... subscripts are (oiperp, ipar)
-	jtutil::array2d<CCTK_INT>* Jacobian_oipar_posn_;
+	jtutil::array2d<CCTK_INT>* Jacobian_y_ipar_posn_;
 
 	// Jacobian values
-	// ... subscripts are (oiperp, ipar, oipar_m)
+	// ... subscripts are (y_iperp, x_ipar, y_ipar_m)
 	jtutil::array3d<fp>* Jacobian_buffer_;
 	};