add patch_system::setup_gridfn_storage()

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

1 files changed, 172 insertions, 6 deletions

diff --git a/src/patch/patch_system.cc b/src/patch/patch_system.cc
index 66d9ef6..9e6e55b 100644
--- a/src/patch/patch_system.cc
+++ b/src/patch/patch_system.cc
@@ -4,6 +4,7 @@
 // patch_system::patch_system
 // patch_system::~patch_system
 // patch_system::construct_patches
+// patch_system::setup_gridfn_storage
 // patch_system::setup_full_sphere_patch_system
 // patch_system::setup_plus_z_hemisphere_patch_system
 // patch_system::setup_plus_xy_quadrant_patch_system
@@ -78,7 +79,9 @@ patch_system::patch_system(fp origin_x_in, fp origin_y_in, fp origin_z_in,
 	: global_coords_(origin_x_in, origin_y_in, origin_z_in),
 	  type_(type_in),
 	  N_patches_(N_patches_of_type(type_in)),
-	  all_patches_(N_patches_)
+	  all_patches_(N_patches_),
+	  gridfn_storage_(NULL),	// set in setup_gridfn_storage()
+	  ghosted_gridfn_storage(NULL)	// set in setup_gridfn_storage()
 {
 const int N_overlap_points = 2*N_extend_points + 1;
 
@@ -157,6 +160,9 @@ default:
 //
 patch_system::~patch_system()
 {
+delete ghosted_gridfn_storage_;
+delete gridfn_storage_;
+
 	for (int pn = N_patches()-1 ; pn >= 0 ; --pn)
 	{
 	delete &ith_patch(pn);
@@ -168,16 +174,21 @@ delete all_patches_;
 //******************************************************************************
 
 //
-// This function constructs a set of patches as described by an array
-// of patch_info structures and associated arguments, and makes these
-// patches members of this patch system.  This function does *NOT*
-// create any of the ghost zones or frontiers, and does *NOT* set up
-// any gridfns.
+// This function is called from the patch_system:: constructor to
+// construct a set of patches as described by an array of patch_info
+// structures and associated arguments, and make these patches members
+// of this patch system.  This function also correctly sets
+//	N_grid_points_
+//	N_ghosted_grid_points_
+// This function does *NOT* create any of the ghost zones or frontiers,
+// and does *NOT* set up any gridfns.
 //
 void patch_system::construct_patches(const struct patch_info patch_info[],
 				     int N_ghost_points, int N_extend_points,
 				     fp delta_drho_dsigma)
 {
+N_grid_points_ = 0;
+ghosted_N_grid_points_ = 0;
 	for (int pn = 0 ; pn < N_patches() ; ++pn)
 	{
 	const struct patch_info& pi = patch_info[pn];
@@ -219,6 +230,161 @@ void patch_system::construct_patches(const struct patch_info patch_info[],
 		}
 
 	all_patches_[pn] = p;
+
+	N_grid_points_+= p->N_grid_points();
+	ghosted_N_grid_points_ += p->ghosted_N_grid_points();
+	}
+}
+
+//******************************************************************************
+
+//
+// This function is called from the patch_system:: constructor to set
+// up the storage for all gridfns in all patches, giving each gridfn a
+// contiguous-across-all-patches storage array.  This function assumes
+// that all the patches have already been constructed.
+//
+// For example, given the patches {x,y,z}, the ghosted gridfns {H,J},
+// and the nominal gridfns {a,b,c}, we might picture the storage like
+// this:
+//
+//	xa xa xa ya ya za za za za
+//	xb xb xb yb yb zb zb zb zb
+//	xc xc xc yc yc zc zc zc zc
+//
+//	xH xH xH xH yH yH yH zH zH zH zH zH
+//	xJ xJ xJ xJ yJ yJ yJ zJ zJ zJ zJ zJ
+//
+// Here the upper/lower blocks are for nominal/ghosted  gridfns.
+// The storage is taken as being contiguous within each row (in fact
+// within each block).  Thus the storage for all the nominal gridfns
+// (or all the ghosted gridfns) in a single patch is *non*-contiguous.
+//
+// The creation of patches is done in several phases: first the patches
+// are constructed with no gridfn storage, then we are called to set up
+// the gridfn storage (taking into account the sizes of the other patches),
+// then finally ghost zones and patch frontiers are constructed and
+// interlinked.
+//
+// FIXME: We should pad the gridfn storage as necessary to avoid cache
+//        conflicts, but we don't do this at present.
+//
+void setup_gridfn_storage(int min_gfn, int max_gfn,
+			  int ghosted_min_gfn, int ghosted_max_gfn)
+{
+const int N_gridfns = jtutil::how_many_in_range(min_gfn, max_gfn);
+const int ghosted_N_gridfns
+	= jtutil::how_many_in_range(ghosted_min_gfn, ghosted_max_gfn);
+
+const int gfn_stride = N_grid_points();
+const int ghosted_gfn_stride = ghosted_N_grid_points();
+
+const int N_storage = gfn_stride * N_gridfns;
+const int ghosted_N_storage = ghosted_gfn_stride * ghosted_N_gridfns;
+
+// storage arrays for all gridfns
+gridfn_storage_         = new fp[N_storage];
+ghosted_gridfn_storage_ = new fp[ghosted_N_storage];
+
+// divide up the storage array among the patches
+int posn = 0;
+int ghosted_posn = 0;
+	for (pn = 0 ; pn < N_patches() ; ++pn)
+	{
+	const struct grid_arrays::gridfn_pars gridfn_pars
+		= {
+		  min_gfn, max_gfn,
+		  & gridfn_storage_[posn],
+		  gfn_stride, 0, 0
+		  };
+	const struct grid_arrays::gridfn_pars ghosted_gridfn_pars
+		= {
+		  ghosted_min_gfn, ghosted_max_gfn,
+		  & ghosted_gridfn_storage_[ghosted_posn],
+		  ghosted_gfn_stride, 0, 0
+		  };
+
+	patch& p = ith_patch(pn);
+	p.setup_gridfn_storage(gridfn_pars, ghosted_gridfn_pars);
+
+	posn += p.N_grid_points();
+	ghosted_posn += p.ghosted_N_grid_points();
+	}
+assert(posn == N_grid_points());
+assert(ghosted_posn == ghosted_N_grid_points());
+
+// check to make sure storage for distinct gridfns
+// forms a partition of the overall storage array
+const patch& pfirst = ith_patch(0);
+const patcn& plast  = ith_patch(N_patches()-1);
+	for (int gfn = min_gfn() ; gfn+1 < max_gfn() ; ++gfn)
+	{
+	// range of storage occupied by gridfns:
+	//	gfn   --> [gfn_first, gfn_last]
+	//	gfn+1 --> [gfn1_first, gfn1_last]
+	const fp* const gfn_last
+		= & plast.gridfn(gfn, plast.max_irho(),
+				      plast.max_isigma());
+	const fp* const gfn1_first
+		= & pfirst.gridfn(gfn+1, pfirst.min_irho(),
+					 pfirst.min_isigma());
+	assert(gfn1_first == gfn_last + 1);
+	}
+
+	for (int gfn = ghosted_min_gfn() ; gfn+1 < ghosted_max_gfn() ; ++gfn)
+	{
+	// range of storage occupied by ghosted gridfns:
+	//	gfn   --> [gfn_first, gfn_last]
+	//	gfn+1 --> [gfn1_first, gfn1_last]
+	const fp* const gfn1_last
+		= & plast.ghosted_gridfn(gfn,
+					 plast.ghosted_max_irho(),
+					 plast.ghosted_max_isigma());
+	const fp* const gfn1_first
+		= & pfirst.ghosted_gridfn(gfn+1,
+					  pfirst.ghosted_min_irho(),
+					  pfirst.ghosted_min_isigma());
+	assert(gnosted_gfn1_first == gfn_last + 1);
+	}
+
+// check to make sure storage for distinct patches
+// forms a partition of the storage for each gridfn
+	for (int gfn = min_gfn() ; gfn < max_gfn() ; ++gfn)
+	{
+		for (int pn = 0 ; pn+1 < N_patches() ; ++pn)
+		{
+		const patch& p = ith_patch(pn);
+		const patch& p1 = ith_patch(pn+1);
+
+		// range of storage occupied by gridfn:
+		//	p  --> [p_first, p_last]
+		//	p1 --> [p1_first, p1_last]
+		const fp* const p_last = & p.gridfn(gfn, p.max_irho(),
+							 p.max_isigma());
+		const fp* const p1_first = & p1.gridfn(gfn, p1.min_irho(),
+							    p1.min_isigma());
+		assert(p1_first == p_last + 1);
+		}
+	}
+
+	for (int gfn = ghosted_min_gfn() ; gfn < ghosted_max_gfn() ; ++gfn)
+	{
+		for (int pn = 0 ; pn+1 < N_patches() ; ++pn)
+		{
+		const patch& p = ith_patch(pn);
+		const patch& p1 = ith_patch(pn+1);
+
+		// range of storage occupied by ghosted gridfn:
+		//	p  --> [p_first, p_last]
+		//	p1 --> [p1_first, p1_last]
+		const fp* const p_last
+			= & p.ghosted_gridfn(gfn, p.ghosted_max_irho(),
+						  p.ghosted_max_isigma());
+		const fp* const p1_first
+			= & p1.ghosted_gridfn(gfn, p1.ghosted_min_irho(),
+						   p1.ghosted_min_isigma());
+		assert(p1_first == p_last + 1);
+		}
 	}
 }