diff options
Diffstat (limited to 'doc/documentation.tex')
-rw-r--r-- | doc/documentation.tex | 1063 |
1 files changed, 1063 insertions, 0 deletions
diff --git a/doc/documentation.tex b/doc/documentation.tex new file mode 100644 index 0000000..5ac372c --- /dev/null +++ b/doc/documentation.tex @@ -0,0 +1,1063 @@ +%version $Header$ +\documentclass{article} + +\def\eqref#1{$(\ref{#1})$} +\def\cf{cf.\hbox{}} +\def\ie{i.e.\hbox{}} +\def\eg{e.g.\hbox{}} +\def\etal{{\it et~al.\/}} + +\def\defn#1{{\bf #1}} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\begin{document} +\title{LocalInterp} +\author{Jonathan Thornburg, incorporating code from Thomas Radke} +% +% We want CVS to expand the Id keyword on the next line, but we don't +% want TeX to go into math mode to typeset the expansion (because that +% wouldn't look as nice in the output), so we use the "$ $" construct +% to get TeX out of math mode again when typesetting the expansion. +% +\date{$ $Id$ $} +\maketitle + +\abstract{ + This thorn provides processor-local interpolation of + 1-D, 2-D, and 3-D data arrays. It provides several + different interpolators, and supports both the older + \verb|CCTK_InterpLocal()| API, and the newer + \verb|CCTK_InterpLocalUniform()| and + \verb|CCTK_InterpLocalNonUniform()| APIs. + } + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Introduction} + +This thorn provides processor-local interpolation of 1, 2, and 3-D +arrays. At present there are 2~interpolators provided (we hope to +add other interpolators soon): +\begin{description} +\item[Uniform Cartesian] + This is the local interpolator which used to live + in the PUGHInterp thorn. It was written in C by + Thomas Radke in early 2001, drawing on earlier Fortran code + by Paul Walker. It supports the \verb|CCTK_InterpLocal()| API. + It provides Lagrange polynomial interpolation for all + combinations of~1, 2, and 3~dimensions, and orders~1, + 2, and~3. (It would be easy to add additional dimensions + and/or orders if desired.) +\item[Generalized Polynomial] + This interpolator was written in C (plus Maple to generate + the coefficients) by Jonathan Thornburg in winter 2001-2002. + It supports the \verb|CCTK_InterpLocalUniform()| API. + It provides Lagrange polynomial interpolation in + 1~dimension (orders~1-6) and~2 and 3~dimensions (orders~1-4). + (Again, it would be easy to add additional dimensions and/or + orders.) It offers a number of options, described below. +\end{description} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Terminology} + +Within Cactus, each interpolator has a character-string name; +this is mapped to a Cactus \defn{interpolator handle} by +\verb|CCTK_InterpHandle()|. For any given interpolator handle, +there may be a separate interpolator defined for each of the +interpolation APIs (both the processor-local ones provided by this +thorn, and the global ones provided by driver-specific interpolator +thorns such as PUGHInterp). + +Terminology for interpolation seems to differ a bit from one author +to another. Here we refer to the \defn{point-centered} interpolation +of a set of \defn{input arrays} (defining data on a uniformly or +nonuniformly spaced \defn{grid} of \defn{data points}) to a set of +\defn{interpolation points} (specified by a corresponding set of +\defn{interpolation coordinates}), with the results being stored +in a corresponding set of \defn{output arrays}. Alternatively, +we may refer to \defn{cell-centered} interpolation, using a grid +of \defn{data cells} and a set of \defn{interpolation cells}. +(This latter terminology is common in hydrodynamics interpolators.) + +At present all the interpolators do polynomial interpolation, and +allow the interpolation of multiple input arrays (to the same set of +interpolation points) in a single interpolator call, using the basic +algorithm: +\begin{verbatim} +for each interpolation point +{ +choose a finite difference molecule position + somewhere near the interpolation point + for each input array + { + compute an interpolating polynomial which approximates + the input data at the molecule points + output = polynomial(interpolation point) + } +} +\end{verbatim} +In the future, we may add other interpolators where the choice of +molecule is data-dependent (and thus may vary from one input array to +the next), and/or where the entire input grid is used in each interpolation. + +We define the \defn{order} of the interpolation to be the order of +the fitted polynomial. That is, in our terminology, linear interpolation +is order~1, quadratic is order~2, cubic is order~3, etc. An order~$n$ +interpolator thus has $O(\Delta x^{n+1})$ errors for smooth input data. + +However, because the interpolating polynomial generally changes if +the interpolation point moves from one grid cell to another, unless +something special is done the interpolating function isn't smooth, +\ie{} its The 1st~derivative is generically {\em discontinuous\/}, +with $O(\Delta x^n)$ jump discontinuities each time the interpolating +polynomial changes. Correspondingly, the interpolation error is +generically a ``bump function'' which is zero at each grid point +and rises to a local maximum in each grid cell. There are interpolation +algorithms (\eg{} Hermite polynomial and spline interpolation) which +give better smoothness, but none of the present interpolators implement +them. :( + +As given in the function reference section of the Cactus User's Guide, +\verb|interp_coords|, \verb|input_arrays|, and \verb|output_arrays| are +actually all pointers to arrays of \verb|void *| pointers, since we +support a number of different Cactus data types. Internally, the +interpolator casts these \verb|void *| pointers to \verb|CCTK_REAL *| +or whatever the correct Cactus data types are. But when describing +how the interpolator accesses the various arrays, for simplicity we +usually gloss over this casting, \ie{}~we pretend that \verb|interp_coords|, +\verb|input_arrays|, and \verb|output_arrays| are pointers to arrays +of \verb|CCTK_REAL *| pointers. (This may become clearer once you +read the next example.) + +We use \verb|pt|, \verb|in|, and \verb|out| as generic 0-origin integer +subscripts into the arrays of interpolation points, input arrays, and +output arrays respectively. We use \verb|(i,j,k)| as a generic +\verb|N_dims|-vector of integer subscripts into the input array +\verb|input_arrays[in]|. (Thus \defn{{\tt |(i,j,k)|} space} refers to +the grid of data points.) We usually only write array subscripting +expressions for the 3-D case; the restrictions/generalizations to +other dimensions should be obvious. + +For example, for 3-D interpolation, the (x,y,z) coordinates of a typical +interpolation point are given by +\begin{verbatim} +x = interp_coords[0][pt] +y = interp_coords[1][pt] +z = interp_coords[2][pt] +\end{verbatim} +(Notice here that as described above, we've implicitly taken +\verb|interp_coords| to have the C~type +\verb|const CCTK_REAL *interp_coords[]|, glossing over the casting +from its actual C~type of \verb|const void *interp_coords[]|.) + +We take \verb|axis| to be an integer specifying a dimension, +\ie{}~0~for~$x$, 1~for~$y$, 2~for~$z$, \dots. + +When describing Jacobians and domains of dependence, it's useful to +introduce the notion of \defn{molecule position}, a nominal reference +point for the interpolation molecule in \verb|(i,j,k)| coordinate +space. (For example, the molecule position might just be the \verb|(i,j,k)| +coordinates of the molecule's center.) We also introduce +\defn{molecule coordinates} \verb|(mi,mj,jk)|, which are just +\verb|(i,j,k)| coordinates relative to the molecule position. +We use \verb|m| or as a generic molecule coordinate. Thus +(in notation which should be obvious) a generic molecule operation +can be written +\begin{equation} +\verb|output| = \sum_{\tt m} \verb|coeff[posn+m]| \times \verb|input[posn+m]| +\end{equation} +Note that there is no requirement that the output be semantically +located at the position \verb|posn|! (This may become clearer once +you look at the example in +section~\ref{sect-generic-options/Jacobian/fixed-sized-hypercube}.) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Generic Interpolator Options} +\label{sect-generic-options} + +The newer \verb|CCTK_InterpLocalUniform()| and +\verb|CCTK_InterpLocalNonUniform()| APIs specify a \defn{parameter table} +(a Cactus key-value table, manipulated by the \verb|Util_Table*()| APIs) +as a generic mechanism for passing optional input to, and/or returning +optional results from, the interpolator. Different interpolators support +different options; in this section we describe some options which are, +or will plausibly, be common to multiple interpolators. + +Note that except as described in section~\ref{sect-generic-options/caching} +(``Caching''), all interpolator arguments and parameters apply only to +the current interpolator call: there is no visible state kept inside +the interpolator from one call to another. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Interpolation Order} + +The uniform Cartesian interpolator encodes the order in the operator +name, but other interpolators should use a (mandatory) parameter-table +parameter +\begin{verbatim} +const CCTK_INT order; +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Handling of Out-of-Range Interpolation Points} + +By default, interpolators should consider it an error if any interpolation +point lies outside the grid, \ie{} if the ``interpolation'' is actually +an extrapolation. (Some interpolators may apply a small ``fuzz'' to +this test to help avoid floating-point rounding problems.) + +Some interpolators may allow this behavior to be changed by the +optional parameter +\begin{verbatim} +const CCTK_REAL out_of_range_tolerance[N_dims]; +\end{verbatim} +The semantics of this are as follows: For each \verb|axis|, +\begin{description} +\item[\rm If $\hbox{\tt out\_of\_range\_tolerance[axis]} \ge 0.0$,] + then an interpolation point is considered to be ``out of range'' + if and only if the interpolation point is + $> \verb|out_of_range_tolerance[axis]| + \times \verb|coord_delta[axis]|$ + outside the grid in any coordinate. +\item[\rm If $\hbox{\tt out\_of\_range\_tolerance[axis]} = -1.0$,] + then an interpolation point is considered to be ``out of range'' + if and only if a centered molecule (or more generally, + a molecule whose centering is chosen pretending that + the grid is of infinite extent), would require data + from outside the grid. +\end{description} +Other values of \verb|out_of_range_tolerance[axis]| are illegal. + +To provide the ``fuzz'' noted above, \verb|out_of_range_tolerance[]| +should default to having all elements set to a very small positive +value, say $10^{-14}$. (This value should really be scaled with the +floating-point precision used,%%% +\footnote{%%% + If this scaling is done, a reasonable value for + {\tt out\_of\_range\_tolerance[]} would be around + $100 \epsilon$, where the \defn{machine epsilon} + $\epsilon$ is defined to be the smallest positive + floating-point number such that $1.0 + \epsilon$ + compares ``not equal to'' 1.0 in floating-point + arithmetic. Machine epsilon values can be found + in the standard header file {\tt <float.h>}; + for IEEE single and double precision they are + about $1.19{\times}10^{-7}$ and $2.22{\times}10^{-16}$ + respectively. + }%%% +{} but none of the interpolators do this at present.) + +If any interpolation points are out of range, interpolators should +return the error code \verb|CCTK_ERROR_INTERP_POINT_X_RANGE|, and +report further details about the error by setting the following +parameter-table entries: +\begin{verbatim} +/* which interpolation point is out of range? */ +/* (value is pt for out-of-range point) */ +CCTK_INT out_of_range_pt; + +/* in which coordinate axis is the point out of range? */ +/* (value is axis for out-of-range point) */ +CCTK_INT out_of_range_axis; + +/* on which end of the grid is the point out of range? */ +/* (value is -1 for min, +1 for max) */ +CCTK_INT out_of_range_end; +\end{verbatim} + +Note that if multiple points and/or axes are out of range, different +interpolators may vary in which out-of-range error they report. That +is, it is {\em not\/} necessarily the case that the ``first'' such +error will be the one reported.%%% +\footnote{%%% + This is to make life simpler for interpolators + which work in parallel over multiple processors. + }%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Molecule Family} +\label{sect-generic-options/molecule-family} + +An interpolator may support different families/shapes of interpolation +molecules. Hypercube-shaped molecules are the simplest and most common +case, but one could also imagine (say) octagon-shaped molecules in 2-D, +or some generalization of this in higher numbers of dimensions. + +The parameter +\begin{verbatim} +const char molecule_family[]; /* set with Util_TableSetString() */ +\end{verbatim} +may be used to query or change the strategy for choosing the molecules. + +The semantics are that if this key is present with the value \verb|""| +(a 0-length ``empty'' string), this queries the default molecule family: +the interpolator will (re)set the value to a string describing the default +molecule family (\verb|"cube"| for hypercube-shaped molecules%%% +\footnote{%%% + Strictly speaking, these should be called + ``parallelipiped-shaped'' molecules, but this + term is so clumsy (and hard to spell!) that + we just call them hypercube-shaped instead. + }%%% +). +If this key is present with a value which is a non-empty string, this +sets the molecule family/shape based on that string. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Non-Contiguous Input Arrays} +\label{sect-generic-options/non-contiguous-inputs} + +Sometimes one of the input ``arrays'' used by the interpolator might not +be contiguous in memory. For example, we might want to do 2-D interpolation +within a plane of a 3-D grid array, but the plane might or might not be +contiguous in memory. (Another example would be that the input array +might be a member of a compact group.) + +One way to do this would be to use the hyperslab API to explicitly compute +the input data on an appropriate hyperslab, then interpolate within that. +However, some interpolators may support accessing the appropriate hyperslab +of the input grid directly inside the interpolator. If this is supported, +it's probably considerably cheaper than explicitly computing the hyperslab. + +If an interpolator supports this, it should use the following (optional) +parameter-table entries to specify non-contiguous inputs: +\begin{verbatim} +const CCTK_INT input_array_offsets[N_input_arrays]; + +/* the next 3 table entries are shared by all input arrays */ +const CCTK_INT input_array_strides [N_dims]; +const CCTK_INT input_array_min_subscripts[N_dims]; +const CCTK_INT input_array_max_subscripts[N_dims]; +\end{verbatim} + +Then the interpolator would access the input using the generic +subscripting expression +\begin{verbatim} +input_array[in][offset + i*stride_i + j*stride_j + k*stride_k] +\end{verbatim} +where +\begin{verbatim} +offset = input_array_offsets[in] +(stride_i,stride_j,stride_k) = input_array_strides[] +\end{verbatim} +\verb|(i,j,k)| run from \verb|input_array_min_subscripts[]| +to \verb|input_array_max_subscripts[]| inclusive +(n.b.~this is an {\em inclusive\/} range, \ie{} +$\verb|min| \le \verb|(i,j,k)| \le \verb|max|$). + +The defaults are that each input array is contiguous in memory, +\ie{} \verb|input_array_offsets[]| = 0, +\verb|stride| determined from \verb|input_array_dims[]| + in the usual Fortran manner, +\verb|input_array_min_subscripts[]| = all 0s, and +\verb|input_array_max_subscripts[]| = \verb|input_array_dims[]|-1. +If the stride and max subscript are both specified explicitly, then the +explicit \verb|input_array_dims[]| argument to +\verb|CCTK_InterpLocalUniform()| is ignored. + +At present all the interpolators take the output arrays to be +contiguous 1-D arrays. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Derivatives} + +Some interpolators can optionally take derivatives as part of the +interpolation, \ie{} if we view the input data as being samples of +a smooth function, then instead of estimating values of that function +at the interpolation points, the interpolator can instead or additionally +estimate values of various derivatives of that function at the +interpolation points. We don't currently implement it, but one could +also imagine interpolating more general molecule operations such as +Laplacians, running integrals, etc. + +To specify such operations, an interpolator should use the parameter-table +entries +\begin{verbatim} +const CCTK_INT operand_indices[N_output_arrays]; +const CCTK_INT operation_codes[N_output_arrays]; +\end{verbatim} +The semantics here are that +\begin{verbatim} +output_array[out] = op(input_array[operand_indices[out]]) +\end{verbatim} +where \verb|op| is an operator specified by the \verb|operation_codes[out]| +value as described below. + +Note that \verb|operand_indices[]| doesn't directly name the inputs, +but rather gives their (0-origin) subscripts in the list of inputs. +This allows for a more efficient implementation in the (common) case +where some of the input arrays have many different operations applied +to them. + +Negative \verb|operation_codes[out]| values are reserved for future +expansion. An \verb|operation_codes[out]| value which is $\ge 0$ is +taken as a decimal integer encoding a coordinate partial derivative: +each decimal digit means to take the coordinate partial derivative along +that (1-origin) axis; the order of the digits in a number is ignored. +For example: +\begin{flushleft} +\begin{tabular}{cl} +\verb|operation_codes[out]| + & What it means \\ +\hline %----------------------------------------------------------------------- +0 & no derivative, ``just'' interpolate the input array \\ +1 & interpolate $\partial \big/ \partial x^1$ of the input array \\ +2 & interpolate $\partial \big/ \partial x^2$ of the input array \\ +12 or 21 + & interpolate $\partial^2 \big/ \partial x^1 \partial x^2$ + of the input array \\ +33 & interpolate $\partial^2 \big/ \partial (x^3)^2$ of the input array %%%\\ +\end{tabular} +\end{flushleft} + +At present we do {\em not\/} have any \verb|#define|s for the +operation codes in the Cactus header files. However, you can avoid +most of the software-engineering problems of having ``magic numbers'' +for the operation codes, by using the macro +\begin{verbatim} +#define DERIV(op) op +\end{verbatim} +to mark all such \verb|operation_codes[]| values in your code. +There's an example of this in section~\ref{sect-example-derivatives}. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Jacobian and Domain of Dependence} +\label{sect-generic-options/Jacobian} + +Sometimes we want to know the Jacobian of the interpolator, +\begin{equation} +\frac{\partial \hbox{\tt output\_array[out][pt]}} + {\partial \hbox{\tt input\_array[in][(i,j,k)]}} + \label{eqn-Jacobian} +\end{equation} +and/or ``just'' the set of \verb|(i,j,k)| for which this is nonzero +(\ie{} the sparsity structure of the Jacobian, equivalently the domain +of dependence of the output result, or equivalently the interpolation +molecule size and shape) for a given \verb|out|, \verb|in|, and \verb|pt|. + +The complexity of doing this depends (strongly!) on the structure +of the Jacobian, and in particular on the answers to the following +questions: +\begin{itemize} +\item What molecule family is being used?%%% +\footnote{%%% + For something like a spline interpolator the Jacobian + is generally nonzero for all {\tt |(i,j,k)|}, but + exponentially small for most of them; in this case + the Jacobian-query API would probably specify a cutoff + to return an approximate Jacobian with reasonable sparsity. + }%%% +\item Does the interpolator always use the same molecule size + and shape (possibly depending on the interpolation order, + molecule family, or other such parameters), independent + of where the interpolation points are in the grid?%%% +\footnote{%%% + We can always make the interpolation molecules be + the ``same size and shape'' by padding them with + zero entries to bring them all up to the size of + the largest molecule used anywhere in the grid, + but this would be very inefficient if many molecules + were padded in this way. + }%%% +\item Does the interpolator use the same interpolation molecule + size and shape for each output array? (The answer to this + question may depend on whether and/or what derivatives are + being computed.) +\end{itemize} +Because the different cases differ so much in their complexity, +we define several distinct APIs for obtaining the interpolator's +Jacobian and/or domain of dependence. + +%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Querying the Interpolator about the Jacobian's structure} + +To allow generic code to determine which of the different Jacobian-structure +cases applies, (and thus which APIs to use), an interpolator which +supports Jacobian operations should support using the following +parameter-table entries to query the interpolator about the Jacobian +structure: + +The parameter \verb|molecule_family| may may be used to query what +molecule family is being used. This is described in detail in +section~\ref{sect-generic-options/molecule-family}. + +This parameter may be used to query whether the interpolation molecule's +size and/or shape varies with the interpolation coordinates: +\begin{verbatim} +CCTK_INT MSS_is_fn_of_interp_coords; +\end{verbatim} +The semantics of this are that if this key is present%%% +\footnote{%%% + ``MSS'' abbreviates ``{\bf m}olecule {\bf s}ize + and/or {\bf s}hape''. + }%%% +{} (the value doesn't matter), then the interpolator will (re)set +the value to 0~if the molecule size and shape do {\em not\/} vary +with the interpolation coordinates, or 1~if the molecule size and/or +shape {\em do\/} vary with the interpolation coordinates. + +This parameter may be used to query whether the interpolation molecule's +size and/or shape varies from one output array to another: +\begin{verbatim} +CCTK_INT MSS_is_fn_of_which_output; +\end{verbatim} +The semantics of this are that if this key is present (the value +doesn't matter), then the interpolator will (re)set the value to +0~if the molecule size and shape do {\em not\/} vary from one output +array to another (this includes the case where there is only 1~output +array!), or 1~if the molecule size and/or shape {\em do\/} vary from +one output array to another. Note that since the answer to this may +depend on whether and/or what derivatives are being computed, setting +this key may force the interpolator to scan completely through the +derivative specifications. This probably isn't {\em very\/} costly, +but you may want to avoid unnecessarily paying that cost by deleting +the query key (\verb|MSS_is_fn_of_which_output|) from the table +on future interpolator calls. + +%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Fixed-Size Hypercube-Shaped Molecules} +\label{sect-generic-options/Jacobian/fixed-sized-hypercube} + +The simplest case (and the only one for which we have defined an API +at present) is when the molecules are hypercube-shaped and of (small) +fixed size, independent of the interpolation coordinates (though likely +depending on the interpolation order or other such parameters), and +not varying from one output array to another. + +The following parameters may be used to query the molecule size: +\begin{verbatim} +CCTK_INT const interp_molecule_min_m[N_dims]; +CCTK_INT const interp_molecule_max_m[N_dims]; +\end{verbatim} +The semantics of these are that if these keys are present (the values +don't matter), then the interpolator will (re)set the values to give +the (inclusive) minimum and maximum \verb|m|~molecule coordinates. + +The following parameter may be used to query the molecule positions: +\begin{verbatim} +CCTK_INT *const interp_molecule_positions[N_dims]; +\end{verbatim} +The semantics of this is that the caller should set +\verb|interp_molecule_positions[]| to an array of +\verb|N_dims| pointers to (caller-supplied) arrays of +\verb|N_interp_points| \verb|CCTK_INT|s each. If this key exists, +then the interpolator will store the molecule positions in the +pointed-to arrays. + +The following parameters may be used to query the +Jacobian~\eqref{eqn-Jacobian} itself: +\begin{verbatim} +CCTK_REAL *const Jacobian_pointer[N_output_arrays]; +const CCTK_INT Jacobian_offset [N_output_arrays]; + +/* the next 2 table entries are shared by all Jacobians */ +const CCTK_INT Jacobian_interp_point_stride; +const CCTK_INT Jacobian_m_strides[N_dims]; +\end{verbatim} +Then the interpolator would store the Jacobian~\eqref{eqn-Jacobian} in +\begin{verbatim} +Jacobian_pointer[out][offset + pt*Jacobian_interp_point_stride + + mi*stride_i + mj*stride_j + mk*stride_k] +\end{verbatim} +where +\begin{verbatim} +offset = Jacobian_offset[out] +(stride_i,stride_j,stride_k) = Jacobian_m_strides[] +\end{verbatim} + +By appropriately setting the various stride parameters, this allows +a fairly wide variety of possible storage layouts for the Jacobian. + +An example may help to clarify this: Suppose we have a 1-D grid +with 11~grid points, with integer subscripts~0 through~10 inclusive, +and interpolation coordinates given by \verb|coord_origin = 0.0| +and \verb|coord_delta = 0.1|. +Suppose further that we're doing quadratic interpolation, using +3-point (vertex-centered) molecules, and that the interpolator +centers the interpolation molecules as close to the interpolation +point as possible subject to the constraint that the molecules +never require data from outside the grid. +Finally, suppose that we query the Jacobian molecule positions for +the \verb|N_interp_points=14| interpolation points 0.0, 0.04, 0.06, +0.10, 0.14, 0.16, 0.20, 0.80, 0.84, 0.86, 0.90, 0.94, 0.96, 1.00. +Then the queries might return +\begin{verbatim} +interp_molecule_min_m = -1 +interp_molecule_max_m = +1 + /* interp_x molecule */ +interp_molecule_positions[0][ 0] = 1 /* 0.00 [0.0, 0.1, 0.2] */ +interp_molecule_positions[0][ 1] = 1 /* 0.04 [0.0, 0.1, 0.2] */ +interp_molecule_positions[0][ 2] = 1 /* 0.06 [0.0, 0.1, 0.2] */ +interp_molecule_positions[0][ 3] = 1 /* 0.10 [0.0, 0.1, 0.2] */ +interp_molecule_positions[0][ 4] = 1 /* 0.14 [0.0, 0.1, 0.2] */ +interp_molecule_positions[0][ 5] = 2 /* 0.16 [0.1, 0.2, 0.3] */ +interp_molecule_positions[0][ 6] = 2 /* 0.20 [0.1, 0.2, 0.3] */ +interp_molecule_positions[0][ 7] = 8 /* 0.80 [0.7, 0.8, 0.9] */ +interp_molecule_positions[0][ 8] = 8 /* 0.84 [0.7, 0.8, 0.9] */ +interp_molecule_positions[0][ 9] = 9 /* 0.86 [0.8, 0.9, 1.0] */ +interp_molecule_positions[0][10] = 9 /* 0.90 [0.8, 0.9, 1.0] */ +interp_molecule_positions[0][11] = 9 /* 0.94 [0.8, 0.9, 1.0] */ +interp_molecule_positions[0][12] = 9 /* 0.96 [0.8, 0.9, 1.0] */ +interp_molecule_positions[0][13] = 9 /* 1.00 [0.8, 0.9, 1.0] */ +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Caching} +\label{sect-generic-options/caching} + +Some interpolators may support special ``caching'' optimizations to +speed repeated interpolations where some or all of the interpolator +arguments and/or parameters are the same. For example, when interpolating +a tabulated equation of state the number of dimensions, the coordinate +origin and grid spacing, and the input arrays (the tabulated equation +of state data), will probably all be the same from one interpolator +call to another. + +If an interpolator supports caching, the following parameters should +be used to control this: + +\begin{verbatim} +const char cache_type[]; /* set with Util_TableSetString() */ +const char cache_control[]; /* set with Util_TableSetString() */ +CCTK_INT cache_handle; +\end{verbatim} + +There are three basic operations supported: +\begin{description} +\item[Create a Cache] + To set up caching, call the interpolator with \verb|cache_type| + set to describe what arguments and/or parameters will remain + the same in future interpolator calls, \verb|cache_control| + set to the string \verb|"create"|, and \verb|cache_handle| + {\em not\/} in the parameter table. The interpolator will + then do extra (possibly quite time-consuming) work to set + up cached information. The interpolator will delete the + key \verb|cache_control|, and return a handle to the cached + information in \verb|cache_handle|; this allows multiple caches + to be active concurrently. +\item[Use a Cache] + To use a cache (\ie{}~to make an interpolation with the + hoped-for speedup), just call the interpolator with + \verb|cache_handle| set to the value returned when the cache + was created. Note that you still have to provide all the + ``will be the same'' interpolator arguments and/or parameters; + providing a cache handle is essentially just a promise that + these will be the same as in the cache-create interpolator + call. The details of what information is cached, and if/how + the ``will be the same'' arguments are still used, are up to + the interpolator. +\item[Destroy a Cache] + To destroy a cache (\ie{} free any memory allocated when + the cache was created), call the interpolator with + \verb|cache_handle| set to the value returned when the cache + was created, and \verb|cache_control| set to the string + \verb|"destroy"|. The interpolator will delete the keys + \verb|cache_handle| and \verb|cache_control|, and destroy + the cache. +\end{description} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{The Uniform Cartesian Interpolator} + +The uniform Cartesian interpolator is the one which used to live +in the PUGHInterp thorn. It was written in C by Thomas Radke in +early 2001, drawing on earlier Fortran code by Paul Walker. It +supports the \verb|CCTK_InterpLocal()| API, and thus doesn't use +any parameter table at all. It provides Lagrange polynomial +interpolation of orders 1, 2, and 3, registered under the operator +names \verb|"first-order uniform Cartesian"|, +\verb|"second-order uniform Cartesian"|, and +\verb|"third-order uniform Cartesian"| respectively. Each of these +supports 1, 2, and 3-D arrays. The code is hard-wired for +hypercube-shaped interpolation molecules, but it would be easy +to add additional orders and/or dimensions if desired. + +Although the \verb|CCTK_InterpLocal()| API supports both uniform +and nonuniform grids for the input data, the present implementation +assumes a uniform grid (and silently gives wrong results for a +nonuniform grid). + +See the Cactus User's Guide ``Full Description of Functions'' +appendix for a full description of the \verb|CCTK_InterpLocal()| +API, and some examples of how to use it. + +%%%%%%%%%%%%%%%%%%%% + +\subsection{Implementation} + +Internally, this interpolator always does 3-D interpolation, inserting +zero coordinates as appropriate for lower dimensionalities. The +interpolation is done by successive 1-D interpolations along each +axis. See the \verb|README| file in the source code directory +\verb|LocalInterp/src/UniformCartesian/| for further details. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{The Generalized Polynomial Interpolator} + +The generalized polynomial interpolator was written in C +(plus Maple to generate the coefficients) by Jonathan Thornburg in +winter 2001-2002. It provides Lagrange polynomial interpolation of +orders 1-6 for 1-D arrays, and 1-4 for 2- and 3-D +arrays, all registered under the operator name +\verb|"generalized polynomial"|. (Again, it would be easy to add +additional orders and/or dimensions if desired.) The code allows +arbitrarily-shaped interpolation molecules, but at the moment only +hypercube-shaped molecules are implemented. + +This interpolator supports a number of the options described in +section~\ref{sect-generic-options}: +\begin{itemize} +\item interpolation order +\item handling of out-of-range interpolation points + (if there are multiple out-of-range points/axes, the + one reported will be the first, \ie{}~the out-of-range + point with the smallest \verb|pt|, and of that points + out-of-range axes, the one with the smallest \verb|axis|) + +\item non-contiguous input arrays +\item derivatives + (when taking derivatives with this interpolator, + it's most efficient to group all operations on + a given input array together in the \verb|operand_indices| + and \verb|operation_codes| arrays, as in the example in + section~\ref{sect-example-derivatives}) +\end{itemize} +At present this interpolator does not support computing the Jacobian, +but we hope to add this in the near future. This interpolator also +does not support any caching features; at present (Feb 2002) we have +no plans to add these. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Smoothing} + +The way the generalized polynomial interpolator is implemented it's +easy to also do Savitzky-Golay smoothing.%%% +\footnote{%%% + See section 14.8 of Numerical Recipes + 2nd edition for a general discussion of + Savitzky-Golay smoothing. + }%%% +{} This is best described by way of an example: Suppose we're doing +1-D cubic interpolation, \ie{} at each interpolation point we're using +a cubic polynomial fitted to 4~surrounding data points. For +Savitzky-Golay smoothing, we would instead {\em least-squares fit\/} +a cubic polynomial to some {\em larger\/} number of surrounding data +points. This combines interpolation with smoothing, so there's less +amplification of noise in the input data in the interpolation outputs. + +The optional input +\begin{verbatim} +const CCTK_INT smoothing; +\end{verbatim} +specifies how much (how many points) to enlarge the interpolation +molecule for this. The default is 0 (no smoothing). 1 would mean to +enlarge the molecule by 1 point, \eg{} to use a 5-point molecule instead +of the usual 4-point one for cubic interpolation. 2 would mean to +enlarge by 2 points, \eg{} to use a 6-point molecule for cubic +interpolation. Etc etc. + +Note that in $>1$~dimension, the default hypercube-shaped molecules +already use more data points than the number of free parameters in +the interpolating polynomials, \ie{} they already do some Savitzky-Golay +smoothing. For example, in 2-D a generic cubic polynomial +$f(x,y) = \sum_{i+j \leq 3} c_{ij} x^i y^j$ has 10 free parameters +$c_{ij}$, which we least-squares fit to the 16 data points in the +$4 \times 4$ molecule. + +Savitzky-Golay smoothing is basically free apart from the increase in +the molecule size, e.g. a \verb|smoothing|=2 cubic interpolation has +exactly the same cost as any other 6-point--molecule interpolation. + +Alas, at the moment only the trivial case \verb|smoothing|=0 is +implemented, but the framework is all there for more general cases. + +%%%%%%%%%%%%%%%%%%%% + +\subsection{Implementation} + +This interpolator's basic design is to use separate specialized code +for each combination of +$(\verb|N_dims|, \verb|molecule_family|, \verb|order|, \verb|smoothing|)$, +\ie{}~in practice for each distinct choice of interpolation molecule. +Maple is used to generate all the interpolation coefficients. +The C preprocessor is then used to generate all the specialized code +from a single master ``template'' file. The template code uses +\verb|#ifdef|s to handle lower dimensions with no extra overhead, +\eg{}~1-D/2-D interpolation is basically just as efficient as in +a purpose-built 1-D/2-D interpolator, with no extra overhead imposed +by the interpolator also supporting higher-dimensional interpolation. + +The Maple code which generates the interpolation coefficients is quite +general-purpose, and can handle an arbitrary dimensionality and molecule +size/shape. Generating new coefficients can be rather time-consuming, +though, \eg{}~the current coefficients for 3-D for orders~1-4 take about +8~cpu minutes to generate using Maple~7 on a 1.7~GHz~P4. + +See the \verb|README| file in the source code directory +\verb|LocalInterp/src/UniformCartesian/| for further details on the +implementation. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{A Simple Example of {\tt CCTK\_InterpLocalUniform} Usage} + +Here's a simple example of interpolating a \verb|CCTK_REAL| and a +\verb|CCTK_COMPLEX| $10 \times 20$ 2-D array, at 5 interpolation points, +using cubic interpolation. + +\begin{verbatim} +#define N_DIMS 2 +#define N_INTERP_POINTS 5 +#define N_INPUT_ARRAYS 2 +#define N_OUTPUT_ARRAYS 2 + +/* (x,y,z) coordinates of data grid points */ +#define X_ORIGIN ... +#define X_DELTA ... +#define Y_ORIGIN ... +#define Y_DELTA ... +const CCTK_REAL origin[N_DIMS] = { X_ORIGIN, Y_ORIGIN }; +const CCTK_REAL delta [N_DIMS] = { X_DELTA, Y_DELTA }; + +/* (x,y,z) coordinates of interpolation points */ +const CCTK_REAL interp_x[N_INTERP_POINTS]; +const CCTK_REAL interp_y[N_INTERP_POINTS]; +const void *const interp_coords[N_DIMS] + = { (const void *) interp_x, (const void *) interp_y }; + +/* input arrays */ +/* ... note Cactus uses Fortran storage ordering, i.e. X is contiguous */ +#define N_X 10 +#define N_Y 20 +const CCTK_REAL input_real [N_Y][N_X]; +const CCTK_COMPLEX input_complex[N_Y][N_X]; +const CCTK_INT input_array_dims[N_DIMS] = { N_X, N_Y }; +const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS] + = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX }; +const void *const input_arrays[N_INPUT_ARRAYS] + = { (const void *) input_real, (const void *) input_complex }; + +/* output arrays */ +CCTK_REAL output_real [N_INTERP_POINTS]; +CCTK_COMPLEX output_complex[N_INTERP_POINTS]; +const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS] + = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX }; +void *const output_arrays[N_OUTPUT_ARRAYS] + = { (void *) output_real, (void *) output_complex }; + +int operator_handle, param_table_handle; +operator_handle = CCTK_InterpHandle("my interpolation operator"); +if (operator_handle < 0) + CCTK_WARN(-1, "can't get interpolation handle!"); +param_table_handle = Util_TableCreateFromString("order=3"); +if (param_table_handle < 0) + CCTK_WARN(-1, "can't create parameter table!"); + +if (CCTK_InterpLocalUniform(N_DIMS, + operator_handle, param_table_handle, + origin, delta, + N_INTERP_POINTS, + CCTK_VARIABLE_REAL, + interp_coords, + N_INPUT_ARRAYS, + input_array_dims, + input_array_type_codes, + input_arrays, + N_OUTPUT_ARRAYS, + output_array_type_codes, + output_arrays) < 0) + CCTK_WARN(-1, "error return from interpolator!"); +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{An Example of Interplating Derivatives} +\label{sect-example-derivatives} + +Consider the problem described earlier: the computation of all the +2nd~derivatives of the 3-metric at a set of interpolation points on +a 2-sphere. Here's how we could do this in C: +\begin{verbatim} +#define N_DIMS 3 + +/* interpolation points */ +#define N_INTERP_POINTS 1000 +const CCTK_REAL interp_x[N_INTERP_POINTS], + interp_y[N_INTERP_POINTS], + interp_z[N_INTERP_POINTS]; +const void *const interp_coords[N_DIMS] + = { (const void *) interp_x, + (const void *) interp_y, + (const void *) interp_z }; + +/* dimensions of the data grid */ +#define NX 30 +#define NY 40 +#define NZ 50 + +/* input arrays */ +#define N_INPUT_ARRAYS 6 +const CCTK_REAL gxx[NX][NY][NZ], gxy[NX][NY][NZ], gxz[NX][NY][NZ], + gyy[NX][NY][NZ], gyz[NX][NY][NZ], + gzz[NX][NY][NZ]; + +const CCTK_INT input_array_dims[N_DIMS] = {NX, NY, NZ}; +const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS] + = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL }; +const void *const interp_coords[N_INPUT_ARRAYS] + = { (const void *) gxx, (const void *) gxy, (const void *) gxz, + (const void *) gyy, (const void *) gyz, + (const void *) gzz }; + +/* output arrays */ +#define N_OUTPUT_ARRAYS 36 +CCTK_REAL + dxx_gxx[N_INTERP_POINTS], dxy_gxx[N_INTERP_POINTS], dxz_gxx[N_INTERP_POINTS], + dyy_gxx[N_INTERP_POINTS], dyz_gxx[N_INTERP_POINTS], + dzz_gxx[N_INTERP_POINTS], + dxx_gxy[N_INTERP_POINTS], dxy_gxy[N_INTERP_POINTS], dxz_gxy[N_INTERP_POINTS], + dyy_gxy[N_INTERP_POINTS], dyz_gxy[N_INTERP_POINTS], + dzz_gxy[N_INTERP_POINTS], + dxx_gxz[N_INTERP_POINTS], dxy_gxz[N_INTERP_POINTS], dxz_gxz[N_INTERP_POINTS], + dyy_gxz[N_INTERP_POINTS], dyz_gxz[N_INTERP_POINTS], + dzz_gxz[N_INTERP_POINTS], + dxx_gyy[N_INTERP_POINTS], dxy_gyy[N_INTERP_POINTS], dxz_gyy[N_INTERP_POINTS], + dyy_gyy[N_INTERP_POINTS], dyz_gyy[N_INTERP_POINTS], + dzz_gyy[N_INTERP_POINTS], + dxx_gyz[N_INTERP_POINTS], dxy_gyz[N_INTERP_POINTS], dxz_gyz[N_INTERP_POINTS], + dyy_gyz[N_INTERP_POINTS], dyz_gyz[N_INTERP_POINTS], + dzz_gyz[N_INTERP_POINTS], + dxx_gzz[N_INTERP_POINTS], dxy_gzz[N_INTERP_POINTS], dxz_gzz[N_INTERP_POINTS], + dyy_gzz[N_INTERP_POINTS], dyz_gzz[N_INTERP_POINTS], + dzz_gzz[N_INTERP_POINTS]; +const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS] + = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL }; +void *const output_arrays[N_OUTPUT_ARRAYS] + = { (void *) dxx_gxx, (void *) dxy_gxx, (void *) dxz_gxx, + (void *) dyy_gxx, (void *) dyz_gxx, + (void *) dzz_gxx, + (void *) dxx_gxy, (void *) dxy_gxy, (void *) dxz_gxy, + (void *) dyy_gxy, (void *) dyz_gxy, + (void *) dzz_gxy, + (void *) dxx_gxz, (void *) dxy_gxz, (void *) dxz_gxz, + (void *) dyy_gxz, (void *) dyz_gxz, + (void *) dzz_gxz, + (void *) dxx_gyy, (void *) dxy_gyy, (void *) dxz_gyy, + (void *) dyy_gyy, (void *) dyz_gyy, + (void *) dzz_gyy, + (void *) dxx_gyz, (void *) dxy_gyz, (void *) dxz_gyz, + (void *) dyy_gyz, (void *) dyz_gyz, + (void *) dzz_gyz, + (void *) dxx_gzz, (void *) dxy_gzz, (void *) dxz_gzz, + (void *) dyy_gzz, (void *) dyz_gzz, + (void *) dzz_gzz }; + +/* integer codes to specify the derivatives */ +/* (for best efficiency we group all operations on a given input together) */ +const CCTK_INT operand_indices[N_OUTPUT_ARRAYS] + = { 0, 0, 0, 0, 0, 0, + 1, 1, 1, 1, 1, 1, + 2, 2, 2, 2, 2, 2, + 3, 3, 3, 3, 3, 3, + 4, 4, 4, 4, 4, 4, + 5, 5, 5, 5, 5, 5 }; +#define DERIV(x) x +const CCTK_INT operation_codes[N_OUTPUT_ARRAYS] + = { DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33), + DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33), + DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33), + DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33), + DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33), + DERIV(11), DERIV(12), DERIV(13), DERIV(22), DERIV(23), DERIV(33) }; + +int operator_handle, param_table_handle; +operator_handle = CCTK_InterpHandle("my interpolation operator"); +if (operator_handle < 0) + CCTK_WARN(-1, "can't get interpolation handle!"); + +param_table_handle = Util_TableCreate(UTIL_TABLE_FLAGS_DEFAULT); +if (param_table_handle < 0) + CCTK_WARN(-1, "can't create parameter table!"); +if (Util_TableSetInt(param_table_handle, 3, "order") < 0) + CCTK_WARN(-1, "can't set order in parameter table!"); +if (Util_TableSetIntArray(param_table_handle, + N_OUTPUT_ARRAYS, operand_indices, + "operand_indices") < 0) + CCTK_WARN(-1, "can't set operand_indices array in parameter table!"); +if (Util_TableSetIntArray(param_table_handle, + N_OUTPUT_ARRAYS, operation_codes, + "operation_codes") < 0) + CCTK_WARN(-1, "can't set operation_codes array in parameter table!"); + +if (CCTK_InterpLocalUniform(N_DIMS, + operator_handle, param_table_handle, + origin, delta, + N_INTERP_POINTS, + CCTK_VARIABLE_REAL, + interp_coords, + N_INPUT_ARRAYS, + input_array_dims, + input_array_type_codes, + input_arrays, + N_OUTPUT_ARRAYS, + output_array_type_codes, + output_arrays) < 0) + CCTK_WARN(-1, "error return from interpolator!"); +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% +% Automatically created from the ccl files +% Do not worry for now. +% +\include{interface} +\include{param} +\include{schedule} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Acknowledgments} + +Thanks to Ian Hawke for helpful comments on the documentation, +to Tom Goodale and Thomas Radke for many useful design discussions, +and to all the Cactus crew for a great infrastructure! + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\end{document} |