diff options
-rw-r--r-- | src/patch/fd_grid.hh | 399 |
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); + }; |