first preprocessor-based version -- *** not tested yet ***

discussion of template version, and overall design, based on
JT personal CVS ~/mpe/util/fd_grid.hh


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

1 files changed, 399 insertions, 0 deletions

diff --git a/src/patch/fd_grid.hh b/src/patch/fd_grid.hh
new file mode 100644
index 0000000..69b20ce
--- /dev/null
+++ b/src/patch/fd_grid.hh
@@ -0,0 +1,399 @@
+// fd_grid.hh -- grid with finite differencing operations
+// $Id$
+//
+// *** Implementation Notes -- Overview ***
+// *** Implementation Notes -- Techniques using C++ Templates ***
+// *** Implementation Notes -- Techniques using the C/C++ Preprocessor ***
+// *** Implementation Notes -- Run-Time Choice of Molecules ***
+// *** finite difference molecules ***
+// ***** fd_grid - grid with finite differencing operations *****
+//
+
+//
+// prerequisites:
+//    <stdio.h>
+//    <assert.h>
+//    <math.h>		// for M_PI (used by degree/radian conversions)
+//    "jt/util++.hh"	// jtutil:: stuff:
+//			//    how_many_in_range(),
+//			//    degrees_of_radians(), radians_of_degrees(),
+//    "jt/array.hh"
+//    "jt/linear_map.hh"
+//    fp.hh
+//    coords.hh
+//    grid.hh
+//
+
+//******************************************************************************
+
+//
+// *** Implementation Notes -- Overview ***
+//
+
+//
+// The key design problem for our finite differencing is how to
+// implement an entire family of 5(9) finite difference operations in
+// 2D(3D)
+//
+//	partial_rho		partial_sigma
+//	partial_{rho,rho}	partial_{rho,sigma}
+//				partial_{sigma,sigma}
+//
+//	partial_x		partial_y		partial_z
+//	partial_xx		partial_xy		partial_xz
+//				partial_yy		partial_yz
+//							partial_zz
+//
+// without having to write out the finite differencing molecules multiple
+// times, and while still preserving maximum inline-function efficiency.
+// In particular, mixed 2nd-order derivative operations like partial_xy
+// should be automatically composed from the two individual 1st derivative
+// operations (partial_x and partial_y).
+//
+
+//
+// Our basic approach is to define each finite difference molecule in
+// a generic 1-dimensional form using an abstract "data(m)" interface.
+// Here we use the terminology that a finite difference molecule is
+// defined as
+//	out[k] = sum(m) c[m] * in[k+m]
+// where c[] is the vector/matrix of molecule coefficients, and m is
+// the (integer) relative grid coordinate within a molecule.
+//
+// That is, for example, we define the usual 2nd order centered 1st
+// derivative operator as
+//	diff = 0.5*inv_delta_x*(data(+1) - data(-1))
+// leaving unspecified just what the date source is.  We then use this
+// with an appropriate data source (indexing along that gridfn array axis)
+// for each directional derivative operation, and we compose two of
+// these, using the first along x as the data source for the second
+// along y, for the mixed 2nd-order derivative operation.
+//
+
+//******************************************************************************
+
+//
+// *** Implementation Notes -- Techniques using C++ Templates ***
+// 
+
+//
+// There are two plausible ways to use C++ templates
+//	[C++ templates are described in detail in chapter 13 of
+//	Stroustrup's "The C++ Programming Language" (3rd Edition),
+//	hereinafter "C++PL", and chapter 15 of Stroustrup's
+//	"The Design and Evolution of C++", hereinafter "D&EC++".]
+// to write the sort of generic-at-compile-time code we want:
+// - Template specializations for each axis, as discussed in D&EC++
+//   section 15.10.3.
+// - Overloaded functions for each axis, with an argument type
+//   (possibly that of an extra unused argument) selecting the
+//   appropriate axis and hence the appropriate function.  This
+//   technique is discussed in D&EC++ section 15.6.3.1.
+//
+// Quoting from D&EC++ (section 15.6.3.1),
+//
+//	The fundamental observation is that every property
+//	of a type or an algorithm can be represented by a
+//	type (possibly defined specificaly to do exactly
+//	that).  That done, such a type can be used to guide
+//	the overload resolution to select a function that
+//	depends on the desired property.  [...]
+//
+//	Please note that thanks to inlining this resolution
+//	is done at compile-time, so the appropriate [...]
+//	function will be called directly without any run-time
+//	overhead.
+//
+// Quoting from C++PL3 (section 13.4),
+//
+//	Passing [...] operations as a template parameter has two
+//	significant benefits compared to alternatives such as
+//	passing pointers to functions.  Several operations can
+//	be passed as a single argument with no run-time cost.
+//	In addition, the [...] operators [passed this way] are
+//	trivial to inline, whereas inlininkg a call through a
+//	pointer to function requires exceptional attention from
+//	a compiler.
+//
+
+//
+// In my opinion the template-specialization design is cleaner, and it
+// clearly has no run-time cost (whereas the overloaded-function design
+// may have a run-time cost for constructing and passing unused objects),
+// so we use it here.
+//
+// There are, however, two (non-fatal) problema with this approach:
+// - Unfortunately, it appears C++ (or at least gcc 2.95.1) forbids
+//   template specialization within a class, so some of the functions
+//   which whould logically be class members, must instead be defined
+//   outside any class.  We use the namespace  fd_stuff::  to hide
+//   these from the outside world.
+// - C++PL3, section C.13.3, states that
+//	Only class templates can be template arguments.
+//   so we have to use dummy classes around some of our template
+//   functions.  To avoid extra constructor/destructor overhead, we
+//   make these template functions static.
+//
+
+//******************************************************************************
+
+//
+// *** Implementation Notes -- Techniques using the C/C++ Preprocessor ***
+//
+
+//
+// The fundamental problem with the template approaches is portability:
+// Although the C++ standard describes powerful template facilities, not
+// all C++ compilers yet fully support these.  As an alternative, we can
+// use the C/C++ preprocessor.  This is ugly and dangerous (global names!),
+// but with care, can provide the same finite differencing functionality
+// and efficiency as the template-based approaches.
+//
+// Because of its greater portability, we use the preprocessor-based
+// approach here.
+//
+
+//******************************************************************************
+
+
+//
+// *** Implementation Notes -- Run-Time Choice of Molecules ***
+//
+// *If* we want to allow the finite differencing scheme to be changed
+// at run-time (e.g. from a parameter file), there are two plausible
+// ways to do this:
+// - Using  switch(molecule_type) , as is standard in C.  This is
+//   simple, and for this particular application quite well-structured
+//   and maintainable (there are only a few different molecule types,
+//   all centralized in this file).
+// - Using virtual functions, with  molecule  a virtual base class
+//   and individual molecules derived from it.  This is elegant, but
+//   may have some performance problems (below).  It also requires some
+//   sort of switch-based "object factory" to interface with with the
+//   molecule-choice parameters.
+//
+// The typical use of these functions will be from within a loop over
+// a whole grid.  In both cases we can expect excellent accuracy from
+// modern hardware branch prediction (and thus minimal performance loss
+// from the branching).  It's reasonable to expect a compiler to fully
+// inline the switch-based code, exposing all the gridfn array subscriptings
+// to strength reduction etc, but this is much trickier for the
+// virtual-function--based code.
+//
+// For this reason, the switch-based design seems superior.
+// However, for the present, we don't even implement this -- we fix
+// the finite differencing scheme at compile time.
+//
+
+//******************************************************************************
+
+//
+// *** finite difference molecules ***
+//
+
+//**************************************
+
+//
+// define the actual molecules
+//
+
+//
+// these  all take the following arguments:
+// inv_delta_x_ = inverse of grid spacing in the finite differencing
+//		  direction
+// data_= a data-fetching function or macro: data_(gfn, irho, isigma)
+//	  is the data to be finite differenced
+// irho_plus_m_ = a function or macro: irho_plus_m_(irho,m) returns the
+//		  rho coordinate to be passed to data_() for the [m]
+//		  molecule coefficient
+// isigma_plus_m_ = same thing, for the sigma coordinate
+//
+// n.b. We grab the variables gfn, irho, and isigma from the calling
+//      environment, and we define assorted local variables as needed!
+//
+
+// 2nd order
+#define FD_GRID__DX_ORDER2(inv_delta_x_, data_,				\
+			   irho_plus_m_, isigma_plus_m_)		\
+	const fp data_p1 = data_(gfn,   irho_plus_m_(irho  ,+1),	\
+				      isigma_plus_m_(isigma,+1));	\
+	const fp data_m1 = data_(gfn,   irho_plus_m_(irho  ,-1),	\
+				      isigma_plus_m_(isigma,-1));	\
+	const fp sum = 0.5 * (data_p1 - data_m1);			\
+	return inv_delta_x_ * sum;					\
+								/* end macro */
+#define FD_GRID__DXX_ORDER2(inv_delta_x_, data_,			\
+			    irho_plus_m_, isigma_plus_m_)		\
+	const fp data_p1 = data_(gfn,   irho_plus_m_(irho  ,+1),	\
+				      isigma_plus_m_(isigma,+1));	\
+	const fp data_0  = data_(gfn,   irho_plus_m_(irho  , 0),	\
+				      isigma_plus_m_(isigma, 0));	\
+	const fp data_m1 = data_(gfn,   irho_plus_m_(irho  ,-1),	\
+				      isigma_plus_m_(isigma,-1));	\
+	const fp sum = data_m1 - 2.0*data_0 + data_p1;			\
+	return jtutil::pow2(inv_delta_x_) * sum;			\
+								/* end macro */
+
+// 4th order
+#define FD_GRID__DX_ORDER4(inv_delta_x_, data_,				\
+			   irho_plus_m_, isigma_plus_m_)		\
+	const fp data_p2 = data_(gfn,   irho_plus_m_(irho  ,+2),	\
+				      isigma_plus_m_(isigma,+2));	\
+	const fp data_p1 = data_(gfn,   irho_plus_m_(irho  ,+1),	\
+				      isigma_plus_m_(isigma,+1));	\
+	const fp data_m1 = data_(gfn,   irho_plus_m_(irho  ,-1),	\
+				      isigma_plus_m_(isigma,-1));	\
+	const fp data_m2 = data_(gfn,   irho_plus_m_(irho  ,-2),	\
+				      isigma_plus_m_(isigma,-2));	\
+	const fp sum = + (8.0/12.0) * (data_p1 - data_m1)		\
+		       + (1.0/12.0) * (data_m2 - data_p2);		\
+	return inv_delta_x_ * sum;					\
+								/* end macro */
+#define FD_GRID__DXX_ORDER4(inv_delta_x_, data_,			\
+			    irho_plus_m_, isigma_plus_m_)		\
+	const fp data_p2 = data_(gfn,   irho_plus_m_(irho  ,+2),	\
+				      isigma_plus_m_(isigma,+2));	\
+	const fp data_p1 = data_(gfn,   irho_plus_m_(irho  ,+1),	\
+				      isigma_plus_m_(isigma,+1));	\
+	const fp data_0  = data_(gfn,   irho_plus_m_(irho  , 0),	\
+				      isigma_plus_m_(isigma, 0));	\
+	const fp data_m1 = data_(gfn,   irho_plus_m_(irho  ,-1),	\
+				      isigma_plus_m_(isigma,-1));	\
+	const fp data_m2 = data_(gfn,   irho_plus_m_(irho  ,-2),	\
+				      isigma_plus_m_(isigma,-2));	\
+	const fp sum = - (30.0/12.0) * data_0				\
+		       + (16.0/12.0) * (data_m1 + data_p1)		\
+		       - ( 1.0/12.0) * (data_m2 + data_p2);		\
+	return jtutil::pow2(inv_delta_x_) * sum;			\
+								/* end macro */
+
+//**************************************
+
+//
+// choose finite differencing order via preprocessor symbol FINITE_DIFF_ORDER
+//
+#ifndef FINITE_DIFF_ORDER
+  #error "must define FINITE_DIFF_ORDER!"
+#endif
+
+#if   (FINITE_DIFF_ORDER == 2)
+  #define FD_GRID__DX	FD_GRID__DX_ORDER2
+  #define FD_GRID__DXX	FD_GRID__DXX_ORDER2
+#elif (FINITE_DIFF_ORDER == 4)
+  #define FD_GRID__DX	FD_GRID__DX_ORDER4
+  #define FD_GRID__DXX	FD_GRID__DXX_ORDER4
+#else
+  #error "unknown value " FINITE_DIFF_ORDER " for FINITE_DIFF_ORDER!"
+#endif
+
+//******************************************************************************
+
+//
+// ***** fd_grid - grid with finite differencing operations *****
+//
+// An  fd_grid  is identical to a  grid  except that it also defines
+// (rho,sigma)-coordinate finite differencing operations on gridfns.
+//
+
+class	fd_grid
+	: public grid
+	{
+private:
+	//
+	// helper functions to compute (irho,isigma) + [m]
+	// along each axis
+	//
+	static
+	  int     rho_axis__irho_plus_m(int irho  , int m) { return irho  +m; }
+	static
+	  int   rho_axis__isigma_plus_m(int isigma, int m) { return isigma  ; }
+	static
+	  int   sigma_axis__irho_plus_m(int irho  , int m) { return irho    ; }
+	static
+	  int sigma_axis__isigma_plus_m(int isigma, int m) { return isigma+m; }
+
+private:
+	//
+	// helper functions for use as DATA_() argument to the
+	// above finite differencing macros
+	//
+	fp gridfn_data__m_is_rho(int gfn, int irho, int isigma,   int m) const
+		{ return gridfn(gfn, irho+m, isigma); }
+	fp gridfn_data__m_is_sigma(int gfn, int irho, int isigma,   int m) const
+		{ return gridfn(gfn, irho, isigma+m); }
+
+
+public:
+	//
+	// ***** user functions for finite differencing *****
+	//
+
+	// 1st derivatives
+	fp partial_rho(int gfn,  int irho, int isigma)
+		const
+		{
+		FD_GRID__DX(inverse_delta_rho(),
+			    gridfn,
+			    rho_axis__irho_plus_m,
+			    rho_axis__isigma_plus_m);
+		}
+	fp partial_sigma(int gfn,  int irho, int isigma)
+		const
+		{
+		FD_GRID__DX(inverse_delta_sigma(),
+			    gridfn,
+			    sigma_axis__irho_plus_m,
+			    sigma_axis__isigma_plus_m);
+		}
+
+	// mixed 2nd partial derivative
+	fp partial_rho_sigma(int gfn,   int irho, int isigma)
+		const
+		{
+		FD_GRID__DX(inverse_delta_rho(),
+			    partial_sigma,
+			    rho_axis__irho_plus_m,
+			    rho_axis__isigma_plus_m);
+		}
+
+	// "pure" 2nd derivatives
+	fp partial_rho_rho(int gfn,  int irho, int isigma)
+		const
+		{
+		FD_GRID__DXX(inverse_delta_rho(),
+			     gridfn,
+			     rho_axis__irho_plus_m,
+			     rho_axis__isigma_plus_m);
+		}
+	fp partial_sigma_sigma(int gfn,  int irho, int isigma)
+		const
+		{
+		FD_GRID__DXX(inverse_delta_sigma(),
+			     gridfn,
+			     sigma_axis__irho_plus_m,
+			     sigma_axis__isigma_plus_m);
+		}
+
+
+public:
+	//
+	// ***** constructor, destructor *****
+	//
+
+	// constructor: pass through to grid:: constructor
+	fd_grid(const grid_arrays::array_pars &grid_array_pars,
+		const grid::grid_pars &grid_pars_in,
+		int N_gridfns_in)
+		: grid(grid_array_pars_in,
+		       grid_pars_in,
+		       N_gridfns_in)
+		{ }
+	// compiler-generated default destructor is ok
+
+private:
+        // we forbid copying and passing by value
+        // by declaring the copy constructor and assignment operator
+        // private, but never defining them
+	fd_grid(const fd_grid &rhs);
+	fd_grid& operator=(const fd_grid &rhs);
+	};