changes to compute Jacobian coefficients properly

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

7 files changed, 292 insertions, 188 deletions

diff --git a/src/patch/ghost_zone.cc b/src/patch/ghost_zone.cc
index 72d9738..939d345 100644
--- a/src/patch/ghost_zone.cc
+++ b/src/patch/ghost_zone.cc
@@ -585,6 +585,7 @@ void interpatch_ghost_zone::compute_Jacobian
 	 bool want_min_par_corner = true,
 	 bool want_non_corner = true,
 	 bool want_max_par_corner = true)
+	const
 {
 verify_caller_wants_all_of_ghost_zone(want_min_par_corner,
 				      want_non_corner,
diff --git a/src/patch/ghost_zone.hh b/src/patch/ghost_zone.hh
index 408c162..5db2c52 100644
--- a/src/patch/ghost_zone.hh
+++ b/src/patch/ghost_zone.hh
@@ -213,10 +213,8 @@
 //
 
 //
-// N.b. const qualifiers on member functions of ghost_zone and its derived
-// classes refer to the union of the class functions and the underlying
-// gridfn data.  That is, anything which changes the gridfn data is taken
-// to be non-const.
+// N.b. const qualifiers in ghost_zone and its derived classes refer to
+//      the underlying gridfn data.
 // 
 
 // forward declarations
@@ -268,6 +266,7 @@ public:
 				      bool want_min_par_corner = true,
 				      bool want_noncorner = true,
 				      bool want_max_par_corner = true)
+		const
 		= 0;
 
 	//
@@ -331,18 +330,6 @@ public:
 	bool is_min() const { return my_edge().is_min(); }
 	bool is_rho() const { return my_edge().is_rho(); }
 
-	// adjacent ghost zones to our min/max corners
-	const ghost_zone& min_par_adjacent_ghost_zone() const
-		{
-		return my_patch()
-		       .ghost_zone_on_edge( my_edge().min_par_adjacent_edge() );
-		}
-	const ghost_zone& max_par_adjacent_ghost_zone() const
-		{
-		return my_patch()
-		       .ghost_zone_on_edge( my_edge().max_par_adjacent_edge() );
-		}
-
 	// min/max iperp of the ghost zone
 	int min_iperp() const
 		{
@@ -373,14 +360,31 @@ public:
 	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
-	// interpatch ghost zone ==> check consistency of this and the
-	//			     other patch's ghost zones and
-	//			     patch_interp objects
-	virtual void assert_fully_setup() const { }
+	// point membership predicate
+	bool is_in_ghost_zone(int iperp, int ipar)
+		const
+		{
+		// n.b. don't test ipar until we're sure iperp is in range!
+		return (iperp >= min_iperp()) && (iperp <= max_iperp())
+		       && (ipar >= min_ipar(iperp))
+		       && (ipar <= max_ipar(iperp));
+		}
 
+	// adjacent ghost zones to our min/max corners
+	const ghost_zone& min_par_adjacent_ghost_zone() const
+		{
+		return my_patch()
+		       .ghost_zone_on_edge( my_edge().min_par_adjacent_edge() );
+		}
+	const ghost_zone& max_par_adjacent_ghost_zone() const
+		{
+		return my_patch()
+		       .ghost_zone_on_edge( my_edge().max_par_adjacent_edge() );
+		}
+
+	//
+	// ***** constructor, finish setup, destructor *****
+	//
 protected:
 	// ... values for  is_interpatch_in  constructor argument
 	static const bool ghost_zone_is_symmetry = false;
@@ -396,6 +400,14 @@ protected:
 		  is_interpatch_(is_interpatch_in)
 		{ }
 public:
+	// assert() that ghost zone is fully setup:
+	// defined here ==> no-op
+	// symmetry ghost zone ==> unchanged ==> no-op
+	// interpatch ghost zone ==> check consistency of this and the
+	//			     other patch's ghost zones and
+	//			     patch_interp objects
+	virtual void assert_fully_setup() const { }
+
 	// destructor must be virtual to allow destruction
 	// of derived classes via ptr/ref to this class
 	virtual ~ghost_zone() { }
@@ -493,6 +505,7 @@ public:
 			      bool want_min_par_corner = true,
 			      bool want_noncorner = true,
 			      bool want_max_par_corner = true)
+		const
 		{ }
 
 	// to which patch/edge do the y points in this Jacobian row belong?
@@ -544,6 +557,11 @@ public:
 	int min_ipar(int iperp) const { return extreme_min_ipar(); }
 	int max_ipar(int iperp) const { return extreme_max_ipar(); }
 
+	//
+	// ***** constructors, destructor *****
+	//
+public:
+
 	// constructor for mirror-symmetry ghost zone
 	symmetry_ghost_zone(const patch_edge& my_edge_in);
 
@@ -637,7 +655,8 @@ public:
 	void compute_Jacobian(int ghosted_min_gfn, int ghosted_max_gfn,
 			      bool want_min_par_corner = true,
 			      bool want_noncorner = true,
-			      bool want_max_par_corner = true);
+			      bool want_max_par_corner = true)
+		const;
 
 	// to which patch/edge do the y points in this Jacobian row belong?
 	patch&            Jacobian_y_patch() const { return other_patch_; }
@@ -676,10 +695,22 @@ public:
 		return (*Jacobian_buffer_)(y_iperp, x_ipar, y_ipar_m);
 		}
 
+private:
+	// helper function for synchronize() and compute_Jacobian():
+	// verify (no-op if ok, error_exit() if not) that caller wants
+	// to operate on the entire ghost zone, since our implementations
+	// of synchronize() and compute_Jacobian() only support this case
+	void verify_caller_wants_all_of_ghost_zone(bool want_min_par_corner,
+						   bool want_non_corner,
+						   bool want_max_par_corner)
+		const;
+
+
 	//
 	// ***** low-level client interface *****
 	//
 
+public:
 	// basic connectivity info
 	patch& other_patch() const { return other_patch_; }
 	const patch_edge& other_edge() const { return other_edge_; }
@@ -688,29 +719,6 @@ public:
 	// and patch_interp objects
 	void assert_fully_setup() const;
 
-	//
-	// constructor, finish setup, destructor
-	//
-	interpatch_ghost_zone(const patch_edge& my_edge_in,
-			      const patch_edge& other_edge_in,
-			      int N_overlap_points);
-	// finish setup (requires adjacent-side ghost_zone objects
-	// to exist, though not to have finish_setup() called):
-	// - setup par coordinate mapping information
-	// - setup interpolator data pointers & result buffer
-	// - create patch_interp object to interpolate from *other* patch
-	void finish_setup(int interp_handle, int interp_par_table_handle);
-	~interpatch_ghost_zone();
-
-private:
-	// verify (no-op if ok, error_exit() if not) that caller wants
-	// to operate on the entire ghost zone, since our implementations
-	// of synchronize() and compute_Jacobian() only support this case
-	void verify_caller_wants_all_of_ghost_zone(bool want_min_par_corner,
-						   bool want_non_corner,
-						   bool want_max_par_corner)
-		const;
-
 	// min/max ipar of the ghost zone for specified iperp
 	// with possibly "triangular" corners depending on the type
 	// (symmetry vs interpatch) of the adjacent ghost zones
@@ -725,6 +733,23 @@ private:
 		return other_iperp_->map(iperp);
 		}
 
+	//
+	// ***** constructor, finish setup, destructor *****
+	//
+public:
+	interpatch_ghost_zone(const patch_edge& my_edge_in,
+			      const patch_edge& other_edge_in,
+			      int N_overlap_points);
+
+	// finish setup (requires adjacent-side ghost_zone objects
+	// to exist, though not to have finish_setup() called):
+	// - setup par coordinate mapping information
+	// - setup interpolator data pointers & result buffer
+	// - create patch_interp object to interpolate from *other* patch
+	void finish_setup(int interp_handle, int interp_par_table_handle);
+
+	~interpatch_ghost_zone();
+
 private:
         // we forbid copying and passing by value
         // by declaring the copy constructor and assignment operator
@@ -786,13 +811,13 @@ private:
 	//      -------------------
 	//      partial gridfn at y
 	//
-	int Jacobian_min_y_ipar_m_, Jacobian_max_y_ipar_m_;
+	mutable int Jacobian_min_y_ipar_m_, Jacobian_max_y_ipar_m_;
 
 	// other patch's y ipar posn for a Jacobian row
 	// ... subscripts are (oiperp, ipar)
-	jtutil::array2d<CCTK_INT>* Jacobian_y_ipar_posn_;
+	mutable jtutil::array2d<CCTK_INT>* Jacobian_y_ipar_posn_;
 
 	// Jacobian values
 	// ... subscripts are (y_iperp, x_ipar, y_ipar_m)
-	jtutil::array3d<fp>* Jacobian_buffer_;
+	mutable jtutil::array3d<fp>* Jacobian_buffer_;
 	};
diff --git a/src/patch/make.code.defn b/src/patch/make.code.defn
index ad0b9ac..2919722 100644
--- a/src/patch/make.code.defn
+++ b/src/patch/make.code.defn
@@ -10,7 +10,8 @@ SRCS = coords.cc \
        ghost_zone.cc \
        patch_interp.cc \
        patch_info.cc \
-       patch_system.cc
+       patch_system.cc \
+       Jacobian.cc
        # test_patch_system.cc
 
 # Subdirectories containing source files
diff --git a/src/patch/patch.cc b/src/patch/patch.cc
index dfa4369..1002a60 100644
--- a/src/patch/patch.cc
+++ b/src/patch/patch.cc
@@ -8,8 +8,8 @@
 // x_patch::x_patch
 // y_patch::y_patch
 //
-// patch::minmax_ang_ghost_zone
 // patch::ghost_zone_on_edge
+// patch::corner_ghost_zone_containing_point
 // patch::interpatch_ghost_zone_on_edge
 // patch::create_mirror_symmetry_ghost_zone
 // patch::create_periodic_symmetry_ghost_zone
@@ -156,29 +156,58 @@ y_patch::y_patch(patch_system &my_patch_system_in, int patch_number_in,
 //******************************************************************************
 
 //
-// This function returns a reference to the specified ghost zone
-// of this patch.
+// This function returns a reference to the ghost zone on a specified
+// edge, after first assert()ing that the edge belongs to this patch.
 //
-ghost_zone& patch::minmax_ang_ghost_zone(bool want_min, bool want_rho)
+// N.b. This function can't be inline in "patch.hh" because it needs
+//	member functions of class patch_edge, which comes after class patch
+//	in our #include order.
+//
+ghost_zone& patch::ghost_zone_on_edge(const patch_edge& e)
 	const
 {
-return want_min ? (want_rho ? min_rho_ghost_zone()
-			    : min_sigma_ghost_zone())
-		: (want_rho ? max_rho_ghost_zone()
-			    : max_sigma_ghost_zone());
+assert(& e.my_patch() == this);
+return minmax_ang_ghost_zone(e.is_min(), e.is_rho());
 }
 
 //******************************************************************************
 
 //
-// This function returns a reference to the specified ghost zone
-// of this patch.
+// This function determines which of the two adjacent ghost zones meeting
+// at a specified corner, contains a specified point.  If the point isn't
+// in either ghost zone, an abort(0) is done.
 //
-ghost_zone& patch::ghost_zone_on_edge(const patch_edge &edge)
-const
+// Arguments:
+// {rho,sigma}_is_min = Specify the corner (and implicitly the ghost zones).
+// irho,isigma = Specify the point.
+//
+// Results:
+// This function returns (a reference to) the desired ghost zone.
+ghost_zone& patch::corner_ghost_zone_containing_point
+	(bool rho_is_min, bool sigma_is_min,
+	 int irho, int isigma)
+	const
 {
-assert(& edge.my_patch() == this);
-return minmax_ang_ghost_zone(edge.is_min(), edge.is_rho());
+ghost_zone&   rho_gz =   minmax_rho_ghost_zone(  rho_is_min);
+ghost_zone& sigma_gz = minmax_sigma_ghost_zone(sigma_is_min);
+
+const patch_edge&   rho_edge =   rho_gz.my_edge();
+const patch_edge& sigma_edge = sigma_gz.my_edge();
+
+const int   rho_iperp =   rho_edge.iperp_of_irho_isigma(irho, isigma);
+const int   rho_ipar  =   rho_edge. ipar_of_irho_isigma(irho, isigma);
+const int sigma_iperp = sigma_edge.iperp_of_irho_isigma(irho, isigma);
+const int sigma_ipar  = sigma_edge. ipar_of_irho_isigma(irho, isigma);
+
+const bool is_in_rho_ghost_zone
+	=   rho_gz.is_in_ghost_zone(  rho_iperp,   rho_ipar);
+const bool is_in_sigma_ghost_zone
+	= sigma_gz.is_in_ghost_zone(sigma_iperp, sigma_ipar);
+
+// check that point is in exactly one ghost zone
+assert(is_in_rho_ghost_zone ^ is_in_sigma_ghost_zone);
+
+return is_in_rho_ghost_zone ? rho_gz : sigma_gz;
 }
 
 //******************************************************************************
diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index b58882d..6325dfe 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -312,11 +312,41 @@ public:
 		assert(max_sigma_ghost_zone_ != NULL);
 		return *max_sigma_ghost_zone_;
 		}
-	ghost_zone& minmax_ang_ghost_zone(bool want_min, bool want_rho) const;
-	ghost_zone& ghost_zone_on_edge(const patch_edge &edge) const;
 
-	// get ghost zone on specified edge, verify that it is indeed
-	// interpatch, static_cast<> to interpatch_ghost_zone ptr
+	// ... these fns are defined explicitly (unlike for grid:: stuff)
+	//     because they're actually used by Jacobian code in ../gr/
+	ghost_zone& minmax_rho_ghost_zone(bool want_min)
+		const
+		{
+		return want_min ? min_rho_ghost_zone()
+				: max_rho_ghost_zone();
+		}
+	ghost_zone& minmax_sigma_ghost_zone(bool want_min)
+		const
+		{
+		return want_min ? min_sigma_ghost_zone()
+				: max_sigma_ghost_zone();
+		}
+
+	ghost_zone& minmax_ang_ghost_zone(bool want_min, bool want_rho)
+		const
+		{
+		return want_rho ? minmax_rho_ghost_zone(want_min)
+				: minmax_sigma_ghost_zone(want_min);
+		}
+
+	ghost_zone& ghost_zone_on_edge(const patch_edge &e) const;
+
+	// which of the two ghost zones at a specifieid corner,
+	// contains a specified point?
+	ghost_zone& corner_ghost_zone_containing_point
+		(bool rho_is_min, bool sigma_is_min,	// specifies corner
+		 int irho, int isigma)			// specifies point
+		const;
+
+	// get ghost zone on specified edge,
+	// assert() that it is indeed interpatch,
+	// static_cast<> to interpatch_ghost_zone ptr
 	interpatch_ghost_zone& interpatch_ghost_zone_on_edge
 		(const patch_edge &e)
 		const;
diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index ec80133..dcc54c4 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -1151,6 +1151,7 @@ void patch_system::synchronize(int ghosted_min_gfn_to_sync,
 //
 void patch_system::compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync,
 						int ghosted_max_gfn_to_sync)
+	const
 {
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
@@ -1227,89 +1228,89 @@ return xgz.Jacobian(x_iperp, x_ipar, y_ipar_m);
 // specified (unknown-grid) radius gridfn.
 //
 void patch_system::print_unknown_gridfn
-(bool ghosted_flag, int unknown_gfn,
- bool print_xyz_flag, bool radius_is_ghosted_flag,
-		      int unknown_radius_gfn,
- const char output_file_name[], bool want_ghost_zones)
-const
+	(bool ghosted_flag, int unknown_gfn,
+	 bool print_xyz_flag, bool radius_is_ghosted_flag,
+			      int unknown_radius_gfn,
+	 const char output_file_name[], bool want_ghost_zones)
+	const
 {
 if (want_ghost_zones && !ghosted_flag)
-then error_exit(PANIC_EXIT,
+   then error_exit(PANIC_EXIT,
 "***** patch_system::print_unknown_gridfn(unknown_gfn=%d):\n"
 "        can't have want_ghost_zones && !ghosted_flag !\n"
 ,
-	   unknown_gfn);				/*NOTREACHED*/
+		   unknown_gfn);				/*NOTREACHED*/
 if (want_ghost_zones && print_xyz_flag && !radius_is_ghosted_flag)
-then error_exit(PANIC_EXIT,
+   then error_exit(PANIC_EXIT,
 "***** patch_system::print_unknown_gridfn(unknown_gfn=%d):\n"
 "        can't have want_ghost_zones && print_xyz_flag\n"
 "        && !radius_is_ghosted_flag!\n"
 "        unknown_radius_gfn=%d\n"
 ,
-	   unknown_gfn,
-	   unknown_radius_gfn);				/*NOTREACHED*/
+		   unknown_gfn,
+		   unknown_radius_gfn);				/*NOTREACHED*/
 
 FILE *output_fp = std::fopen(output_file_name, "w");
 if (output_fp == NULL)
-then error_exit(ERROR_EXIT,
+   then error_exit(ERROR_EXIT,
 "***** patch_system::print_unknown_gridfn(unknown_gfn=%d):\n"
 "        can't open output file \"%s\"\n!"
 ,
-	   unknown_gfn,
-	   output_file_name);				/*NOTREACHED*/
-
-for (int pn = 0 ; pn < N_patches() ; ++pn)
-{
-const patch& p = ith_patch(pn);
-
-fprintf(output_fp, "### %s patch\n", p.name());
-fprintf(output_fp, "# %s_gfn=%d\n",
-	(ghosted_flag ? "ghosted" : "nominal"), unknown_gfn);
-fprintf(output_fp, "# dpx = %s\n", p.name_of_dpx());
-fprintf(output_fp, "# dpy = %s\n", p.name_of_dpy());
-fprintf(output_fp, "#\n");
-fprintf(output_fp,
-	print_xyz_flag
-	? "# dpx\tdpy\tgridfn\tglobal_x\tglobal_y\tglobal_z\n"
-	: "# dpx\tdpy\tgridfn\n");
+		   unknown_gfn,
+		   output_file_name);				/*NOTREACHED*/
 
-	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)
+	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
-	const fp rho = p.rho_of_irho(irho);
-	const fp sigma = p.sigma_of_isigma(isigma);
-	const fp dpx = p.dpx_of_rho_sigma(rho, sigma);
-	const fp dpy = p.dpy_of_rho_sigma(rho, sigma);
+	const patch& p = ith_patch(pn);
+
+	fprintf(output_fp, "### %s patch\n", p.name());
+	fprintf(output_fp, "# %s_gfn=%d\n",
+		(ghosted_flag ? "ghosted" : "nominal"), unknown_gfn);
+	fprintf(output_fp, "# dpx = %s\n", p.name_of_dpx());
+	fprintf(output_fp, "# dpy = %s\n", p.name_of_dpy());
+	fprintf(output_fp, "#\n");
 	fprintf(output_fp,
-		"%g\t%g\t%#.15g",
-		dpx, dpy, p.unknown_gridfn(ghosted_flag,
-					   unknown_gfn, irho,isigma));
-	if (print_xyz_flag)
-	   then {
-		const fp r = p.unknown_gridfn(radius_is_ghosted_flag,
-					      unknown_radius_gfn,
-					      irho,isigma);
-		fp local_x, local_y, local_z;
-		p.xyz_of_r_rho_sigma(r,rho,sigma,
-				     local_x,local_y,local_z);
-		const fp global_x = global_x_of_local_x(local_x);
-		const fp global_y = global_y_of_local_y(local_y);
-		const fp global_z = global_z_of_local_z(local_z);
+		print_xyz_flag
+		? "# dpx\tdpy\tgridfn\tglobal_x\tglobal_y\tglobal_z\n"
+		: "# dpx\tdpy\tgridfn\n");
+
+		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);
+		const fp dpx = p.dpx_of_rho_sigma(rho, sigma);
+		const fp dpy = p.dpy_of_rho_sigma(rho, sigma);
 		fprintf(output_fp,
-			"\t%#g\t%#g\t%#g",
-			global_x, global_y, global_z);
+			"%g\t%g\t%#.15g",
+			dpx, dpy, p.unknown_gridfn(ghosted_flag,
+						   unknown_gfn, irho,isigma));
+		if (print_xyz_flag)
+		   then {
+			const fp r = p.unknown_gridfn(radius_is_ghosted_flag,
+						      unknown_radius_gfn,
+						      irho,isigma);
+			fp local_x, local_y, local_z;
+			p.xyz_of_r_rho_sigma(r,rho,sigma,
+					     local_x,local_y,local_z);
+			const fp global_x = global_x_of_local_x(local_x);
+			const fp global_y = global_y_of_local_y(local_y);
+			const fp global_z = global_z_of_local_z(local_z);
+			fprintf(output_fp,
+				"\t%#g\t%#g\t%#g",
+				global_x, global_y, global_z);
+			}
+		fprintf(output_fp, "\n");
+		}
+		fprintf(output_fp, "\n");
 		}
 	fprintf(output_fp, "\n");
 	}
-	fprintf(output_fp, "\n");
-	}
-fprintf(output_fp, "\n");
-}
 
 std::fclose(output_fp);
 }
@@ -1329,89 +1330,89 @@ void patch_system::read_unknown_gridfn(bool ghosted_flag, int unknown_gfn,
 			       bool want_ghost_zones)
 {
 if (want_ghost_zones && !ghosted_flag)
-then error_exit(PANIC_EXIT,
+   then error_exit(PANIC_EXIT,
 "***** patch_system::read_unknown_gridfn(unknown_gfn=%d):\n"
 "        can't have want_ghost_zones && !ghosted_flag !\n"
 ,
-	   unknown_gfn);				/*NOTREACHED*/
+		   unknown_gfn);				/*NOTREACHED*/
 
 FILE *input_fp = std::fopen(input_file_name, "r");
 if (input_fp == NULL)
-then error_exit(ERROR_EXIT,
+   then error_exit(ERROR_EXIT,
 "***** patch_system::read_unknown_gridfn(unknown_gfn=%d):\n"
 "        can't open input file \"%s\"\n!"
 ,
-	   unknown_gfn,
-	   input_file_name);				/*NOTREACHED*/
+		   unknown_gfn,
+		   input_file_name);				/*NOTREACHED*/
 
 int line_number = 1;
-for (int pn = 0 ; pn < N_patches() ; ++pn)
-{
-patch& p = ith_patch(pn);
+	for (int pn = 0 ; pn < N_patches() ; ++pn)
+	{
+	patch& p = 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);
-const fp dpx = p.dpx_of_rho_sigma(rho, sigma);
-const fp dpy = p.dpy_of_rho_sigma(rho, sigma);
-
-const int N_buffer = 100;
-char buffer[N_buffer];
-// read/discard comments and blank lines
-	do
+	for (int irho = p.effective_min_irho(want_ghost_zones) ;
+	     irho <= p.effective_max_irho(want_ghost_zones) ;
+	     ++irho)
 	{
-	if (std::fgets(buffer, N_buffer, input_fp) == NULL)
-	   then error_exit(ERROR_EXIT,
+	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);
+	const fp dpx = p.dpx_of_rho_sigma(rho, sigma);
+	const fp dpy = p.dpy_of_rho_sigma(rho, sigma);
+
+	const int N_buffer = 100;
+	char buffer[N_buffer];
+	// read/discard comments and blank lines
+		do
+		{
+		if (std::fgets(buffer, N_buffer, input_fp) == NULL)
+		   then error_exit(ERROR_EXIT,
 "***** patch::read_unknown_gridfn(%s patch, unknown_gfn=%d):\n"
 "        I/O error or unexpected end-of-file on input!\n"
 "        at irho=%d of [%d,%d], isigma=%d of [%d,%d]\n"
 "        dpx=%g dpy=%g\n"
 ,
-			   p.name(), unknown_gfn,
-			   irho, p.effective_min_irho(want_ghost_zones),
-				 p.effective_max_irho(want_ghost_zones),
-			   isigma,
-			      p.effective_min_isigma(want_ghost_zones),
-			      p.effective_max_isigma(want_ghost_zones),
-			   dpx, dpy);			/*NOTREACHED*/
-	++line_number;
-	} while ((buffer[0] == '#') || (buffer[0] == '\n'));
-
-double read_dpx, read_dpy, read_gridfn_value;
-if (sscanf(buffer, "%lf %lf %lf",
-	   &read_dpx, &read_dpy, &read_gridfn_value) != 3)
-   then error_exit(ERROR_EXIT,
+				   p.name(), unknown_gfn,
+				   irho, p.effective_min_irho(want_ghost_zones),
+					 p.effective_max_irho(want_ghost_zones),
+				   isigma,
+				      p.effective_min_isigma(want_ghost_zones),
+				      p.effective_max_isigma(want_ghost_zones),
+				   dpx, dpy);			/*NOTREACHED*/
+		++line_number;
+		} while ((buffer[0] == '#') || (buffer[0] == '\n'));
+
+	double read_dpx, read_dpy, read_gridfn_value;
+	if (sscanf(buffer, "%lf %lf %lf",
+		   &read_dpx, &read_dpy, &read_gridfn_value) != 3)
+	   then error_exit(ERROR_EXIT,
 "***** patch::read_unknown_gridfn(%s patch, unknown_gfn=%d):\n"
 "        bad input data at input line %d!\n"
 ,
-		   p.name(), unknown_gfn,
-		   line_number);			/*NOTREACHED*/
-if (! (    jtutil::fuzzy<fp>::EQ(read_dpx,dpx)
-	&& jtutil::fuzzy<fp>::EQ(read_dpy,dpy)    ) )
-   then error_exit(ERROR_EXIT,
+			   p.name(), unknown_gfn,
+			   line_number);			/*NOTREACHED*/
+	if (! (    jtutil::fuzzy<fp>::EQ(read_dpx,dpx)
+		&& jtutil::fuzzy<fp>::EQ(read_dpy,dpy)    ) )
+	   then error_exit(ERROR_EXIT,
 "***** patch::read_unknown_gridfn(%s patch, unknown_gfn=%d):\n"
 "        wrong (dpx,dpy) at input line %d!\n"
 "        expected (%g,%g)\n"
 "        read     (%g,%g)\n"
 ,
-		   p.name(), unknown_gfn,
-		   line_number,
-		   dpx, dpy,
-		   read_dpx, read_dpy);			/*NOTREACHED*/
+			   p.name(), unknown_gfn,
+			   line_number,
+			   dpx, dpy,
+			   read_dpx, read_dpy);			/*NOTREACHED*/
 
-p.unknown_gridfn(ghosted_flag,
-		 unknown_gfn, irho,isigma) = read_gridfn_value;
-}
-}
+	p.unknown_gridfn(ghosted_flag,
+			 unknown_gfn, irho,isigma) = read_gridfn_value;
+	}
+	}
 
-}
+	}
 
 std::fclose(input_fp);
 }
diff --git a/src/patch/patch_system.hh b/src/patch/patch_system.hh
index a4777f0..a1b4227 100644
--- a/src/patch/patch_system.hh
+++ b/src/patch/patch_system.hh
@@ -27,7 +27,8 @@
 //******************************************************************************
 
 //
-// a  patch_system  object describes a system of interlinked patches.
+// A  patch_system  object describes a system of interlinked patches.
+// Its  const  qualifiers refer to the gridfn data.
 //
 
 class	patch_system
@@ -245,13 +246,26 @@ public:
 	void synchronize(int ghosted_min_gfn_to_sync,
 			 int ghosted_max_gfn_to_sync);
 
+	// ... do this for all ghosted gridfns
+	void synchronize()
+		{ synchronize(ghosted_min_gfn(), ghosted_max_gfn()); }
+
 
 	//
 	// compute Jacobian of  synchronize()  into internal buffers,
 	// taking into account synchronize()'s full 3-phase algorithm
 	//
 	void compute_synchronize_Jacobian(int ghosted_min_gfn_to_sync,
-					  int ghosted_max_gfn_to_sync);
+					  int ghosted_max_gfn_to_sync)
+		const;
+
+	// ... do this for all ghosted gridfns
+	void compute_synchronize_Jacobian()
+		const
+		{
+		compute_synchronize_Jacobian(ghosted_min_gfn(),
+					     ghosted_max_gfn());
+		}
 
 	//
 	// The following functions access the Jacobian computed by
@@ -266,7 +280,8 @@ public:
 		const
 		{ return xgz.Jacobian_y_patch(); }
 	const patch_edge& synchronize_Jacobian_y_edge (const ghost_zone& xgz)
-		const { return xgz.Jacobian_y_edge(); }
+		const
+		{ return xgz.Jacobian_y_edge(); }
 
 	// what is the [min,max] range of m for a Jacobian row?
 	int synchronize_Jacobian_min_y_ipar_m(const ghost_zone& xgz)
@@ -283,7 +298,8 @@ public:
 
 	// what is the  posn  value of the y points in this Jacobian row?
 	int synchronize_Jacobian_y_ipar_posn(const ghost_zone& xgz,
-					     int x_iperp, int x_ipar) const;
+					     int x_iperp, int x_ipar)
+		const;
 
 	// what is the Jacobian
 	//	partial synchronize() gridfn(ghosted_gfn, px, x_iperp, x_ipar)
@@ -294,7 +310,8 @@ public:
 	//	y_posn = Jacobian_y_ipar_posn(x_iperp, x_ipar)
 	// taking into account synchronize()'s full 3-phase algorithm
 	fp synchronize_Jacobian(const ghost_zone& xgz, int x_iperp, int x_ipar,
-				int y_ipar_m) const;
+				int y_ipar_m)
+		const;
 
 
 	//