diff options
Diffstat (limited to 'doc')
-rw-r--r-- | doc/TODO | 16 | ||||
-rw-r--r-- | doc/documentation.tex | 1735 | ||||
-rw-r--r-- | doc/references | 15 |
3 files changed, 1766 insertions, 0 deletions
diff --git a/doc/TODO b/doc/TODO new file mode 100644 index 0000000..4424b69 --- /dev/null +++ b/doc/TODO @@ -0,0 +1,16 @@ +$Header$ + +Things to do on this thorn: +- make casts from datatype to CCTK_REAL explicit in + GeneralizedPolynomial-Uniform fetch routines +- clean up the horrible inefficiency of the Hermite interpolator + (cf Erik Schnetter's E-mails and CactusBase/1366) + +- allow scalar boundary_{off_centering,extrapolation}_tolerance + values in the parameter table (to mean using this value for each ibndry) + +- maybe change the defaults so off-centering is forbidden unless + the user explicitly asks for it to be allowed? + +- implement a "silent" parameter like in PUGHInterp? + alternatively, print warning msgs for some error returns? diff --git a/doc/documentation.tex b/doc/documentation.tex new file mode 100644 index 0000000..4a4699e --- /dev/null +++ b/doc/documentation.tex @@ -0,0 +1,1735 @@ +% *======================================================================* +% Cactus Thorn template for ThornGuide documentation +% Author: Ian Kelley +% Date: Sun Jun 02, 2002 +% $Header$ +% +% Thorn documentation in the latex file doc/documentation.tex +% will be included in ThornGuides built with the Cactus make system. +% The scripts employed by the make system automatically include +% pages about variables, parameters and scheduling parsed from the +% relevent thorn CCL files. +% +% This template contains guidelines which help to assure that your +% documentation will be correctly added to ThornGuides. More +% information is available in the Cactus UsersGuide. +% +% Guidelines: +% - Do not change anything before the line +% % START CACTUS THORNGUIDE", +% except for filling in the title, author, date etc. fields. +% - Each of these fields should only be on ONE line. +% - Author names should be sparated with a \\ or a comma +% - You can define your own macros, but they must appear after +% the START CACTUS THORNGUIDE line, and must not redefine standard +% latex commands. +% - To avoid name clashes with other thorns, 'labels', 'citations', +% 'references', and 'image' names should conform to the following +% convention: +% ARRANGEMENT_THORN_LABEL +% For example, an image wave.eps in the arrangement CactusWave and +% thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps +% - Graphics should only be included using the graphix package. +% More specifically, with the "includegraphics" command. Do +% not specify any graphic file extensions in your .tex file. This +% will allow us (later) to create a PDF version of the ThornGuide +% via pdflatex. | +% - References should be included with the latex "bibitem" command. +% - Use \begin{abstract}...\end{abstract} instead of \abstract{...} +% - Do not use \appendix, instead include any appendices you need as +% standard sections. +% - For the benefit of our Perl scripts, and for future extensions, +% please use simple latex. +% +% *======================================================================* +% +% Example of including a graphic image: +% \begin{figure}[ht] +% \begin{center} +% \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure} +% \end{center} +% \caption{Illustration of this and that} +% \label{MyArrangement_MyThorn_MyLabel} +% \end{figure} +% +% Example of using a label: +% \label{MyArrangement_MyThorn_MyLabel} +% +% Example of a citation: +% \cite{MyArrangement_MyThorn_Author99} +% +% Example of including a reference +% \bibitem{MyArrangement_MyThorn_Author99} +% {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999), +% 1--16. {\tt http://www.nowhere.com/}} +% +% *======================================================================* + +% If you are using CVS use this line to give version information +% $Header$ + +\documentclass{article} + +% Use the Cactus ThornGuide style file +% (Automatically used from Cactus distribution, if you have a +% thorn without the Cactus Flesh download this from the Cactus +% homepage at www.cactuscode.org) +\usepackage{../../../../doc/ThornGuide/cactus} + +\begin{document} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% The author of the documentation +\author{Jonathan Thornburg, incorporating code from Thomas Radke} + +% The title of the document (not necessarily the name of the Thorn) +\title{Thorn Guide for the {\bf LocalInterp} Thorn} + +% the date your document was last changed, if your document is in CVS, +% please us: +% \date{$ $Date$ $} +\date{$ $Date$ $} + +\maketitle + +% Do not delete next line +% START CACTUS THORNGUIDE + +% Add all definitions used in this documentation here +% \def\mydef etc +\setlength{\unitlength}{1mm} % length unit for latex picture environment +\def\csmash#1{\hbox to 0em{\hss{#1}\hss}} + +\def\defn#1{{\bf #1}} +\def\thorn#1{{\bf #1}} +\def\arrangement#1{{\bf #1}} + +\def\Cplusplus{\hbox{C\hspace{-.05em}\raisebox{.4ex}{\tiny\bf ++}}} +\def\cf{cf.\hbox{}} +\def\Ie{I.e.\hbox{}} +\def\ie{i.e.\hbox{}} +\def\eg{e.g.\hbox{}} + +% Add an abstract for this thorn's documentation +\begin{abstract} +This thorn provides the \verb|CCTK_InterpLocalUniform()| API for +processor-local interpolation of N-dimensional data arrays. The data +arrays (in general there may be many of them) must be defined on a +uniformly spaced grid. + +This thorn provides several variants of Lagrange and Hermite +interpolation. At present the Lagrange interpolation operators support +orders~1 through~6 for 1-D interpolation, and~1 through~4 for 2-D +and 3-D interpolation. The Hermite interpolation operators support +orders~2, 3, and~4 for 1-D interpolation, and~2 and~3 for 2-D +and 3-D interpolation. + +This thorn supports a number of interpolation options, including +parameter-table entries to control the handling of grid boundaries +and out-of-range points, non-contiguous input arrays, taking +derivatives as part of the interpolation, and querying the +Jacobian of the interpolation operation. +\end{abstract} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Introduction} + +This interpolator provides: +\begin{itemize} +\item Lagrange polynomial interpolation of orders 1--6 for 1-D + interpolation and 1--4 for 2-D and 3-D interpolation. + This interpolator actually provides two slightly different + flavors of Lagrange interpolation, \defn{tensor product} + and \defn{maximum degree}, registered under the operator names + \begin{verbatim} + "Lagrange polynomial interpolation (tensor product)" + "Lagrange polynomial interpolation (maximum degree)" + \end{verbatim} + Section~\ref{AEIThorns/AEILocalInterp/sect-multi-dim-interp} + of this thorn guide explains the two variants in detail, but + for most purposes the tensor-product variant is + preferable. For convenience this is also registered under + the operator name + \begin{verbatim} + "Lagrange polynomial interpolation" + \end{verbatim} + and this is probably the best interpolation operator for + general-purpose use. For backwards compatability this interpolator + is also registered under still another operator name, + \begin{verbatim} + "generalized polynomial interpolation" + \end{verbatim} +\item Hermite polynomial interpolation of orders 2-4 for 1-D + interpolation, and 2--3 for 2-D and 3-D interpolation, + registered under the operator name + \begin{verbatim} + "Hermite polynomial interpolation" + \end{verbatim} +\end{itemize} +The code allows arbitrarily-shaped interpolation molecules, +but at the moment only hypercube-shaped molecules are implemented. +It would be easy to add additional orders and/or dimensions if desired. + +This interpolator supports a number of options specified by a +\defn{parameter table} (a Cactus key-value table, manipulated by the +\verb|Util_Table*()| APIs). Note that all the interpolator options +described here apply only to the current interpolator call: there is +no visible state kept inside the interpolator from one call to another. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{History} + +This interpolator was written by Jonathan Thornburg in winter 2001--2002. +Between then and July~2003 it lived in the +\thorn{LocalInterp} thorn in the \arrangement{CactusBase} arrangement, +but in July 2003 it was moved to the (new) +\thorn{AEILocalInterp} thorn in the \arrangement{AEIThorns} arrangement +so it could stay GPL +(Cactus policies forbid GPL code in the CactusBase arrangement). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Basic Terminology} +\label{AEIThorns/AEILocalInterp/sect-basic-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 \thorn{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 an interpolation molecule position + somewhere near the interpolation point + for each output 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\big((\Delta x)^{n+1})\big)$ interpolation errors +for generic smooth input data. +Section~\ref{AEIThorns/AEILocalInterp/sect-multi-dim-interp} explains how the +interpolating polynomial is defined for multi-dimensional interpolation. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{The Non-Smoothness of Interpolation Errors} + +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 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. This is the case, +for example, for tensor-product Lagrange polynomial interpolation. + +[For maximum-degree Lagrange polynomial interpolation, the +interpolation error may not be zero at the grid points, and the +interpolant itself may not be continuous there. +Section~\ref{AEIThorns/AEILocalInterp/sect-multi-dim-interp} explains this +in detail.] + +This thorn also provides Hermite polynomial interpolation, which +guarantees a smooth ($C^1$) interpolant and (for smooth input data) +interpolation error. Unfortunately, this comes at the cost of a +larger molecule size than the same-order Lagrange interplator, and +a much larger interpolation error if the interpolation molecules are +off-centered. (By default, the Hermite interpolator doesn't allow +off-centered molecules, though this can be changed via the +\verb|boundary_off_centering_tolerance[]| and +\verb|excision_off_centering_tolerance[]| parameter-table entries.) + +For most purposes Lagrange polynomial interpolation is better; +only use Hermite polynomial interpolation if you need a smooth +interpolant. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{More Terminology} + +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 a 0-origin integer specifying a coordinate +axis (dimension), \ie{} 0~for~$x$, 1~for~$y$, 2~for~$z$, \dots. +We take \verb|ibndry| to be a 0-origin integer specifying the +combination of a coordinate axis (dimension) and a minimum/maximum +``end'' of the grid; these are enumerated in the order +$x_{\min}$, $x_{\max}$, $y_{\min}$, $y_{\max}$, $z_{\min}$, +$z_{\max}$, \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{AEIThorns/AEILocalInterp/sect-Jacobian/fixed-sized-hypercube}.) +However, by convention we do assume that $\verb|m|=0$ is always a valid +\verb|m| coordinate; this simplifies pointer arithmetic. + +When describing various entries in the parameter table in, +we use \verb|const| qualifiers on table entries to indicate that the +interpolator treats them as \verb|const| variables/arrays, \ie{} it +promises not to change them. In contrast, table entries which are +not shown as \verb|const| are considered read-write by the interpolator; +it typically uses them to return outputs back to the caller. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Interpolator Parameters} + +The \verb|CCTK_InterpLocalUniform()| interpolation API +(described in detail in the ``Function Reference'' appendix of +the Cactus Users' Guide) includes a \defn{parameter table} +(a Cactus key-value table, manipulated by the \verb|Util_Table*()| APIs) +whose table handle is passed to the interpolator. This can be used +to specify a number of options to the interpolator; it can also be +used by the interpolator to return additional information back to +the caller. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Interpolation Order} + +The only mandatory parameter for this interpolator is the interpolation +order: +\begin{verbatim} +const CCTK_INT order; +\end{verbatim} +As noted in section~\ref{AEIThorns/AEILocalInterp/sect-basic-terminology}, +in our terminology linear interpolation is order~1, quadratic is +order~2, cubic is order~3, etc. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Molecule Size and Centering} +\label{AEIThorns/AEILocalInterp/sect-molecule-size+centering} + +If no grid boundaries or excised points are nearby, the interpolator +centers the molecules around the interpolation point as much as possible. +Table~\ref{AEIThorns/AEILocalInterp/tab-molecule-size+centering} gives the molecule +size and details of the centering policy for each interpolation operator +and order. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\begin{table} + +\def\Linterval{\lower0.60ex\hbox{\hskip-0.20em$\big[$}} +\def\Rinterval{\lower0.60ex\hbox{\hskip-0.25em$\big)$}} + +\def\TwoPointMolecule{%%% +\raise0.5ex\hbox{%%% +\begin{picture}(0,0) +\put(-4,0){\circle*{1.5}} +\put(-4,0){\Linterval} +\put(4,0){\Rinterval} +\put(4,0){\circle*{1.5}} +\end{picture}%%% +}}%%% + +\def\ThreePointMolecule{%%% +\raise0.5ex\hbox{%%% +\begin{picture}(0,0) +\put(-8,0){\circle*{1.5}} +\put(-4,0){\Linterval} +\put(0,0){\circle*{1.5}} +\put(4,0){\Rinterval} +\put(8,0){\circle*{1.5}} +\end{picture}%%% +}}%%% + +\def\FourPointMolecule{%%% +\raise0.5ex\hbox{%%% +\begin{picture}(0,0) +\put(-12,0){\circle*{1.5}} +\put(-4,0){\circle*{1.5}} +\put(-4,0){\Linterval} +\put(4,0){\Rinterval} +\put(4,0){\circle*{1.5}} +\put(12,0){\circle*{1.5}} +\end{picture}%%% +}}%%% + +\def\FivePointMolecule{%%% +\raise0.5ex\hbox{%%% +\begin{picture}(0,0) +\put(-16,0){\circle*{1.5}} +\put(-8,0){\circle*{1.5}} +\put(-4,0){\Linterval} +\put(0,0){\circle*{1.5}} +\put(4,0){\Rinterval} +\put(8,0){\circle*{1.5}} +\put(16,0){\circle*{1.5}} +\end{picture}%%% +}}%%% + +\def\SixPointMolecule{%%% +\raise0.5ex\hbox{%%% +\begin{picture}(0,0) +\put(-20,0){\circle*{1.5}} +\put(-12,0){\circle*{1.5}} +\put(-4,0){\circle*{1.5}} +\put(-4,0){\Linterval} +\put(4,0){\Rinterval} +\put(4,0){\circle*{1.5}} +\put(12,0){\circle*{1.5}} +\put(20,0){\circle*{1.5}} +\end{picture}%%% +}}%%% + +\def\SevenPointMolecule{%%% +\raise0.5ex\hbox{%%% +\begin{picture}(0,0) +\put(-24,0){\circle*{1.5}} +\put(-16,0){\circle*{1.5}} +\put(-8,0){\circle*{1.5}} +\put(-4,0){\Linterval} +\put(0,0){\circle*{1.5}} +\put(4,0){\Rinterval} +\put(8,0){\circle*{1.5}} +\put(16,0){\circle*{1.5}} +\put(24,0){\circle*{1.5}} +\end{picture}%%% +}}%%% + +\renewcommand{\arraystretch}{1.5} +\begin{flushleft} +\begin{tabular}{@{}lcc@{\quad\hspace{25mm}}c@{\hspace{25mm}}} + & & Molecule + & \\[-1ex] +Interpolation Operator & Order & Size + & \csmash{Molecule Centering Policy} \\ +\hline %------------------------------------------------------- +\verb|"Lagrange polynomial interpolation"| & 1 & 2-point + & \TwoPointMolecule \\ +\verb|"Lagrange polynomial interpolation"| & 2 & 3-point + & \ThreePointMolecule \\ +\verb|"Lagrange polynomial interpolation"| & 3 & 4-point + & \FourPointMolecule \\ +\verb|"Lagrange polynomial interpolation"| & 4 & 5-point + & \FivePointMolecule \\ +\verb|"Lagrange polynomial interpolation"| & 5 & 6-point + & \SixPointMolecule \\ +\verb|"Lagrange polynomial interpolation"| & 6 & 7-point + & \SevenPointMolecule \\ +\hline %------------------------------------------------------- +\verb|"Hermite polynomial interpolation"| & 2 & 4-point + & \FourPointMolecule \\ +\verb|"Hermite polynomial interpolation"| & 3 & 6-point + & \SixPointMolecule \\ +\verb|"Hermite polynomial interpolation"| & 4 & 6-point + & \SixPointMolecule \\ +\hline %------------------------------------------------------- +\end{tabular} +\end{flushleft} +\caption[Molecule Sizes] + { + This table gives the molecule size and centering for each + interpolation operator and order. $\bullet$~shows the input + grid points, and $[ \quad )$~shows the interval of interpolation + points (relative to the grid) for which the interpolator + will use the molecule shown. + For example, with grid points at integer coordinates, if + the the interpolator is using 5-point molecules, it will + use a molecule containing the grid points $\{ -2, -1, 0, 1, 2 \}$ + for all interpolation points in the range $[-0.5, 0.5)$. + } +\label{AEIThorns/AEILocalInterp/tab-molecule-size+centering} +\end{table} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Handling of Grid Boundaries} + +Near grid boundaries and/or excised points the interpolator can either +off-center the interpolation molecules, or refuse to interpolate +(returning a \verb|CCTK_INTERP_ERROR_POINT_OUTSIDE|\,%%% +\footnote{%%% + For backwards compatability the error code + {\tt CCTK\_INTERP\_ERROR\_POINT\_X\_RANGE} + is also defined; it's a synonym for + {\tt CCTK\_INTERP\_ERROR\_POINT\_OUTSIDE}. + }%%% +{} or \verb|CCTK_INTERP_ERROR_POINT_EXCISED| error code). This behavior +is controlled by the (optional) parameter-table entries +\begin{verbatim} +const CCTK_INT N_boundary_points_to_omit[2*N_dims] +const CCTK_REAL boundary_off_centering_tolerance[2*N_dims] +const CCTK_REAL boundary_extrapolation_tolerance[2*N_dims] +const CCTK_REAL excision_off_centering_tolerance[2*N_dims] # not implemented yet +const CCTK_REAL excision_extrapolation_tolerance[2*N_dims] # not implemented yet +\end{verbatim} +The elements of these arrays are matched up with the grid boundaries +in the order\\ +$[x_{\min}, x_{\max}, y_{\min}, y_{\max}, z_{\min}, z_{\max}, \dots]$. +(For \verb|excision_*_tolerance[]| the minimum/maximum are interpreted +as the corresponding coordinate decreasing/increasing when moving +out of the active region of the grid into the excised region.) +We use \verb|ibndry| as a generic 0-origin index into these arrays. + +%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Omitting Boundary Points} + +For any given grid boundary \verb|ibndry|, the interpolator doesn't use +input data from the\\ +\verb|N_boundary_points_to_omit[ibndry]| grid points closest to the +boundary. In other words, these points are effectively omitted from +the input grid. For example: +\begin{itemize} +\item Setting this parameter to an array of all 0s means the + interpolator may use data from all grid points. + This is the default if this parameter isn't specified. +\item Setting this parameter to an array of all 1s means the + interpolator may not use input data from the outermost + row/plane of grid points on any face. +\end{itemize} + +More generally, this option may be useful for controlling the extent +to which the interpolator uses input data from ghost and/or symmetry +zones. + +%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Off-Centering Molecules near Grid Boundaries + and/or Excision Regions} + +The (optional) parameter-table entries +\begin{verbatim} +const CCTK_REAL boundary_off_centering_tolerance[2*N_dims] +const CCTK_REAL boundary_extrapolation_tolerance[2*N_dims] +const CCTK_REAL excision_off_centering_tolerance[2*N_dims] # not implemented yet +const CCTK_REAL excision_extrapolation_tolerance[2*N_dims] # not implemented yet +\end{verbatim} +control the local interpolator's willingness to off-center +(change away from the default centering) its interpolation molecules +near boundaries and excised points. + +To describe how these parameters work, it's useful to define two +regions in the data coordinate space: +\begin{description} +\item[The \defn{valid-data region}] + is the entire grid minus any ommited-on-boundary + or excised points; to make this a region we take the + Lego-block bounding box of the non-excised grid points. +\item[The \defn{default-centering region}] + is that region in interpolation-point space where the + interpolator can use the default molecule centering + (described in table~\ref{AEIThorns/AEILocalInterp/tab-molecule-size+centering}), + \ie{} where the default-centering molecules require data + only from the data-valid region. +\end{description} +Figure~\ref{AEIThorns/AEILocalInterp/fig-valid-data+default-centering} shows an +example of these regions. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\begin{figure} + +% "X" to mark excised point +\def\ExcisedPoint{%%% +\begin{picture}(0,0) +\put(-2,-2){\line(1,1){4}} +\put(2,-2){\line(-1,1){4}} +\end{picture}%%% +}%%% + +\def\Linterval{\lower0.60ex\hbox{\hskip-0.20em$\big[$}} + +\def\xmin{\begin{picture}(0,0)\put(4,0){\vector(-1,0){4}}\end{picture}} +\def\xmax{\begin{picture}(0,0)\put(-4,0){\vector(1,0){4}}\end{picture}} +\def\ymin{\begin{picture}(0,0)\put(0,4){\vector(0,-1){4}}\end{picture}} +\def\ymax{\begin{picture}(0,0)\put(0,-4){\vector(0,1){4}}\end{picture}} + +\begin{center} +\begin{picture}(150,150) +%%% grid +\multiput(0,0)(10,0){16}{\multiput(0,0)(0,10){16}{\circle*{1.5}}} +%%% excised points +\multiput(70,60)(10,0){2}{\ExcisedPoint} +\multiput(60,70)(10,0){4}{\ExcisedPoint} +\multiput(60,80)(10,0){4}{\ExcisedPoint} +\multiput(70,90)(10,0){2}{\ExcisedPoint} +%%% outer-boundary valid-data region +\multiput(0,0)(2,0){75}{\line(1,0){1}} +\multiput(150,0)(0,2){75}{\line(0,1){1}} +\multiput(150,150)(-2,0){75}{\line(-1,0){1}} +\multiput(0,150)(0,-2){75}{\line(0,-1){1}} +% ... and movement directions +\put(75,10){\ymin} +\put(140,75){\xmax} +\put(75,140){\ymax} +\put(10,75){\xmin} +%%% outer-boundary default-centering regions... +% ... for 3-point molecules +\put(5,5){\line(1,0){140}} +\put(145,5){\line(0,1){140}} +\put(145,145){\line(-1,0){140}} +\put(5,145){\line(0,-1){140}} +% ... for 4-point molecules +\put(10,10){\line(1,0){130}} +\put(140,10){\line(0,1){130}} +\put(140,140){\line(-1,0){130}} +\put(10,140){\line(0,-1){130}} +%%% excised-point valid-data region +\multiput(60,60)(0,-2){5}{\line(0,-1){1}} +\multiput(60,50)(2,0){15}{\line(1,0){1}} +\multiput(90,50)(0,2){5}{\line(0,1){1}} +\multiput(90,60)(2,0){5}{\line(1,0){1}} +\multiput(100,60)(0,2){15}{\line(0,1){1}} +\multiput(100,90)(-2,0){5}{\line(-1,0){1}} +\multiput(90,90)(0,2){5}{\line(0,1){1}} +\multiput(90,100)(-2,0){15}{\line(-1,0){1}} +\multiput(60,100)(0,-2){5}{\line(0,-1){1}} +\multiput(60,90)(-2,0){5}{\line(-1,0){1}} +\multiput(50,90)(0,-2){15}{\line(0,-1){1}} +\multiput(50,60)(2,0){5}{\line(1,0){1}} +%%% excised-point default-centering regions... +% ... for 3-point molecules +\put(55,55){\line(0,-1){10}} +\put(55,45){\line(1,0){40}} +\put(95,45){\line(0,1){10}} +\put(95,55){\line(1,0){10}} +\put(105,55){\line(0,1){40}} +\put(105,95){\line(-1,0){10}} +\put(95,95){\line(0,1){10}} +\put(95,105){\line(-1,0){40}} +\put(55,105){\line(0,-1){10}} +\put(55,95){\line(-1,0){10}} +\put(45,95){\line(0,-1){40}} +\put(45,55){\line(1,0){10}} +% ... for 4-point molecules +\put(50,50){\line(0,-1){10}} +\put(50,40){\line(1,0){50}} +\put(100,40){\line(0,1){10}} +\put(100,50){\line(1,0){10}} +\put(110,50){\line(0,1){50}} +\put(110,100){\line(-1,0){10}} +\put(100,100){\line(0,1){10}} +\put(100,110){\line(-1,0){50}} +\put(50,110){\line(0,-1){10}} +\put(50,100){\line(-1,0){10}} +\put(40,100){\line(0,-1){50}} +\put(40,50){\line(1,0){10}} +% ... and movement directions for 4-point molecules +\put(50,45){\xmax} +\put(75,40){\ymax} +\put(100,45){\xmin} +\put(105,50){\ymax} +\put(110,75){\xmin} +\put(105,100){\ymin} +\put(100,105){\xmin} +\put(75,110){\ymin} +\put(50,105){\xmax} +\put(45,100){\ymin} +\put(40,75){\xmax} +\put(45,50){\ymax} +\end{picture} +\end{center} + +\caption[Example of valid-data and default-centering regions] + { + This figure shows an example of the valid-data and + default-centering regions for a 2-D grid. + $\times$~marks the excised points, and the dashed lines + show the boundaries of the valid-data region. + The solid lines show the boundaries of the default-centering + regions for 3-point and 4-point molecules. The arrows on the + default-centering region for 4-point molecules show which + elements of the \hbox{\tt boundary\_*\_tolerance[ibndry]} + and \hbox{\tt excision\_*\_tolerance[ibndry]} arrays apply + to each line segment of the region boundary:\\ +\centerline{%%% +\renewcommand{\arraystretch}{1.25} +\begin{tabular}{c@{\qquad}r} +$\raise0.50ex\hbox{\xmin \hskip4mm}$ & $x_{\min}$ ({\tt ibndry = 0}) \\ +$\raise0.50ex\hbox{\hskip4mm \xmax}$ & $x_{\max}$ ({\tt ibndry = 1}) \\ +$\lower1mm\hbox{\,\ymin\,}$ & $y_{\min}$ ({\tt ibndry = 2}) \\ +$\raise3mm\hbox{\,\ymax\,}$ & $y_{\max}$ ({\tt ibndry = 3}) %%%\\ +\end{tabular}%%% +} + } +\label{AEIThorns/AEILocalInterp/fig-valid-data+default-centering} +\end{figure} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +Near a boundary or excised regions, the interpolator will off-center +its molecules (if and) only if the interpolation point is both +\begin{itemize} +\item within ($\le$) \verb|*_off_centering_tolerance[ibndry]| + grid spacings of the default-centering region, {\bf and} +\item within ($\le$) \verb|*_extrapolation_tolerance[ibndry]| + grid spacings of the valid-data region. +\end{itemize} +where we use ``\verb|*|'' to denote \verb|boundary| or \verb|excision| +as appropriate. + +There are four cases for these parameters: +\begin{description} +\item[\mathversion{bold} + $\text{\tt *\!\_off\_centering\_tolerance[ibdnry]} = 0.0$ and + $\text{\tt *\!\_extrapolation\_tolerance[ibndry]} = 0.0$:]\mbox{}\\ + No off-centering is allowed: the interpolator reports an + error (\verb|CCTK_ERROR_INTERP_POINT_OUTSIDE| or + \verb|CCTK_ERROR_INTERP_POINT_EXCISED| return code + as appropriate) if any interpolation point is outside + the default-centering region. +\item[\mathversion{bold} + $\text{\tt *\!\_off\_centering\_tolerance[ibdnry]} > 0.0$ and + $\text{\tt *\!\_extrapolation\_tolerance[ibndry]} = 0.0$:]\mbox{}\\ + The interpolator allows interpolation points to be up to ($\le$) + \verb|*_off_centering_tolerance[ibndry]| grid spacings outside + the default-centering region in this direction. + If an interpolation point is beyond this limit in any + direction, then the interpolator reports an error. +\item[\mathversion{bold} + $\text{\tt *\!\_off\_centering\_tolerance[ibdnry]} = \infty$%%% +\footnotemark{}%%% +{} and + $\text{\tt *\!\_extrapolation\_tolerance[ibndry]} > 0.0$:]\mbox{}\\ +\footnotetext{%%% + In practice $999.0$ is a conventional value here. + The actual value doesn't matter so long as it's + larger than any possible molecule radius. + }%%% + The interpolator allows interpolation points to be up to ($\le$) + \verb|*_extrapolation_tolerance[ibndry]| grid spacings outside + the valid-data region in this direction. + If an interpolation point is beyond this limit in any + direction, then the interpolator reports an error. +\item[\mathversion{bold} + $\text{\tt *\!\_off\_centering\_tolerance[ibdnry]} = 0.0$ and + $\text{\tt *\!\_extrapolation\_tolerance[ibndry]} > 0.0$:]\mbox{}\\ + In practice the default-centering region is always a + (normally proper) subset of the valid-data region, so this + case is nonsensical: the positive value for + \verb|*_extrapolation_tolerance[ibndry]| has no effect + because of the $\verb|*_off_centering_tolerance[ibndry]| = 0.0$ + setting. The interpolator gives a warning for this case. +\end{description} + +If any of these table entries aren't specified, the defaults are +\begin{verbatim} +boundary_off_centering_tolerance[ibndry] = 999.0 # effectively infinity +boundary_extrapolation_tolerance[ibndry] = 1.0e-10 +\end{verbatim} +In other words, the interpolation points may be anywhere within the +valid-data region or up to $10^{-10}$ grid spacing outside it. (The +$10^{-10}$ ``fudge factor'' helps to avoid spurious errors in case +floating-point roundoff moves an interpolation point which was supposed +to be just on the boundary of the valid-data region, slightly outside it.) + +[In the near future these defaults may be changed: For Lagrange +polynomial interpolation the new defaults would be +\begin{verbatim} +boundary_off_centering_tolerance[ibndry] = 999.0 # effectively infinity +boundary_extrapolation_tolerance[ibndry] = 1.0e-10 +\end{verbatim} +while for Hermite polynomial interpolation they would be +\begin{verbatim} +boundary_off_centering_tolerance[ibndry] = 1.0e-10 +boundary_extrapolation_tolerance[ibndry] = 0.0 +\end{verbatim} +This would leave Lagrange interpolation unchanged, while for Hermite +interpolation the defaults would forbid any significant off-centering +of the interpolation molecules.] + +%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Out-of-Range Interpolation Point Error Handling} + +If the interpolator finds an interpolation point which is ``out of range'' +as described above, it sets the following parameter-table entries to +give more information about the out-of-range point: +\begin{verbatim} +/* which interpolation point is it? */ +/* (value gives 0-origin pt for the offending point) */ +CCTK_INT error_pt; + +/* in which coordinate axis and direction is the point out of range? */ +/* (value gives 0-origin ibndry for the offending point) */ +CCTK_INT error_ibndry; + +/* in which coordinate axis is the point outside the grid or excised? */ +/* (value gives 0-origin axis for the offending point) */ +CCTK_INT error_axis; + +/* in which direction is the point out of range? */ +/* (value is -1 for min, +1 for max) */ +CCTK_INT error_direction; +\end{verbatim} + +The interpolator then returns a \verb|CCTK_ERROR_INTERP_POINT_OUTSIDE| +error code. + +Note that if the point is out of range in multiple axes, it's undefined +which of them will be reported. Also, if there are multiple out-of-range +points, user code shouldn't make any assumptions about which of them +will be reported -- the only safe assumption is that {\em some\/} +out-of-range point will be reported. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Multi-Dimensional Interpolation} +\label{AEIThorns/AEILocalInterp/sect-multi-dim-interp} + +In multiple dimensions, there are two plausible definitions for the +generic interpolating polynomial of a given order (degree) $n$. For +convenience I'll describe these for the 2-D case, but the generalization +to any number of dimensions should be obvious: + +\begin{itemize} +\item The generic ``\defn{tensor product}'' polynomial of degree $n$ is + \begin{equation} + f(x,y) = \sum_{{0 \le i \le n} \atop {0 \le j \le n}} a_{ij} x^i y^j + \label{AEIThorns/AEILocalInterp/eqn-polynomial-2d-TP} + \end{equation} +\item The generic ``\defn{maximum degree}'' polynomial of degree $n$ is + \begin{equation} + g(x,y) = \sum_{0 \le i+j \le n} a_{ij} x^i y^j + \label{AEIThorns/AEILocalInterp/eqn-polynomial-2d-MD} + \end{equation} +\end{itemize} + +Because it has $(n+1)^2$ coefficients, the tensor-product polynomial $f$ +can (and does) pass through all the $(n+1)^2$ input data points in a +square molecule of size $n+1$. This implies that the interpolation +error vanishes at the input grid points, and that the overall +interpolating function is continuous (up to floating-point roundoff +errors). However, $f$ does have the slightly peculiar property of +having terms up to $x^n y^n$ despite being of ``degree'' $n$. +(For example, the ``linear'' interpolant for $n=1$ would have $xy$ terms, +even though those are formally quadratic in the independent variables +$x$ and $y$.) + +In contrast, the maximum-degree polynomial $g$ has only +$\frac{1}{2} (n+1)(n+2)$ coefficients, so for generic input data +(\ie{} input data which isn't actually sampled from a polynomial of +the form~$(\ref{AEIThorns/AEILocalInterp/eqn-polynomial-2d-MD})$) $g$ can't pass +through all the $(n+1)^2$ input data points in a square molecule of +size $n+1$. The interpolator actually does a least-squares fit of +the polymomial $g$ to the input data in the molecule, so in general +the molecule won't pass through {\em any\/} of the data points! +Moreover, each time the interpolation point crosses a grid square +(for odd~$n$) or the center lines of a grid square (for even~$n$), +the set of points used in the molecule changes +(\cf{}~section~\ref{AEIThorns/AEILocalInterp/sect-molecule-size+centering}), +so the interpolant generally has a jump discontinuity! +For these reasons, the tensor-product +choice~$(\ref{AEIThorns/AEILocalInterp/eqn-polynomial-2d-TP})$ is generally +preferable. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Molecule Family} +\label{AEIThorns/AEILocalInterp/sect-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 (optional) parameter +\begin{verbatim} +/* set with Util_TableSetString() */ +const char molecule_family[]; +\end{verbatim} +may be used to set or query the molecule family. + +If this key is present in the parameter table, the interpolator sets +the molecule family/shape based on the value specified. +If this key {\em isn't\/} present in the parameter table, then the +interpolator sets it to the molecule family being used. + +At present only hypercubed-shaped molecules are implemented; these +are specified by \verb|molecule_family| set to \verb|"cube"|. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Non-Contiguous Input Arrays} +\label{AEIThorns/AEILocalInterp/sect-non-contiguous-inputs} + +Sometimes 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 +arrays might be members of a compact group.) + +The following (optional) parameter-table entries specify non-contiguous +input arrays: +\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} + +In general, the interpolator accesses 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} +and where \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. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Derivatives} +\label{AEIThorns/AEILocalInterp/sect-derivatives} + +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.) + +The following (optional) parameter-table entries are used to specify +this: +\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. (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{AEIThorns/AEILocalInterp/sect-example-derivatives}.) + +Negative \verb|operation_codes[out]| values are reserved for future +use. 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. +Table~\ref{AEIThorns/AEILocalInterp/tab-derivative-codes} summarizes these resulting +derivative codes. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\begin{table} +\begin{center} +\begin{tabular}{cl} +\verb|operation_codes[out]| + & What it means \\ +\hline %--------------------------------------------------------------- +0 & interpolate the input array itself (no derivative) \\[1ex] +1 & interpolate $\partial \big/ \partial x^1$ of the input array \\ +2 & interpolate $\partial \big/ \partial x^2$ of the input array \\ +3 & interpolate $\partial \big/ \partial x^3$ of the input array \\[1ex] +12 or 21& interpolate $\partial^2 \big/ \partial x^1 \partial x^2$ + of the input array \\ +13 or 31& interpolate $\partial^2 \big/ \partial x^1 \partial x^3$ + of the input array \\ +23 or 32& interpolate $\partial^2 \big/ \partial x^2 \partial x^3$ + of the input array \\ +11 & interpolate $\partial^2 \big/ \partial (x^1)^2$ of the input array \\ +22 & interpolate $\partial^2 \big/ \partial (x^2)^2$ of the input array \\ +33 & interpolate $\partial^2 \big/ \partial (x^3)^2$ of the input array \\ +\hline %--------------------------------------------------------------- +\end{tabular} +\end{center} +\caption[Derivative Codes] + { + This table gives the codes in {\tt operation\_codes[out]} + for each possible 1st or 2nd~derivative in 3-D; for 1-D or + 2-D codes referring to nonexistent coordinates are invalid. + } +\label{AEIThorns/AEILocalInterp/tab-derivative-codes} +\end{table} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +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{AEIThorns/AEILocalInterp/sect-example-derivatives}. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Jacobian and Domain of Dependence} +\label{AEIThorns/AEILocalInterp/sect-Jacobian} + +The Jacobian of the interpolator is defined as +\begin{equation} +\frac{\partial \hbox{\tt output\_array[out][pt]}} + {\partial \hbox{\tt input\_array[in][(i,j,k)]}} + \label{AEIThorns/AEILocalInterp/eqn-Jacobian} +\end{equation} +We may want to know the Jacobian itself, and/or ``just'' the set of +\verb|(i,j,k)| for which this is nonzero (\ie{} the Jacobian's sparsity +structure, or 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|.%%% +\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. + }%%% + +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? +\item Does the interpolation molecule size and/or shape depend + on 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 If this interpolator supports computing derivatives + as described in section~\ref{AEIThorns/AEILocalInterp/sect-derivatives}, + does the interpolation molecule size and/or shape depend + on \verb|operation_codes[]|? +\item Does the interpolation molecule size and/or shape depend + on the actual floating-point values being interpolated? + (Examples of this might include ENO (essentially nonoscillatory) + and/or TVD (total-variation diminishing) interpolators for + hydrodynamics calculations.) +\item Do the actual floating-point values of the Jacobian depend + on the actual floating-point values being interpolated? + Equivalently, is the interpolation nonlinear? +\end{itemize} + +Because the different cases differ so much in the amount of information +required to describe the Jacobian, it's hard to define a single API +which covers all cases without burdening the simpler cases with +excessive complexity. Instead, the interpolator defines a +Jacobian-structure--query API to determine which case applies, +together with (conceptually) several different APIs for the +different cases. (At the moment only a single case is implemented.) + +%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Determining the Jacobian's structure} +\label{AEIThorns/AEILocalInterp/sect-Jacobian/structure} + +The following parameter-table entries may be used to query which +of the different Jacobian-structure cases applies: + +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{AEIThorns/AEILocalInterp/sect-molecule-family}. + +If the interpolation molecule size and/or shape vary with the +interpolation coordinates, the interpolator sets the parameter +\begin{verbatim} +CCTK_INT MSS_is_fn_of_interp_coords; +\end{verbatim} +to 1. Otherwise (\ie{} if the interpolation molecule size and shape +are independent of the interpolation coordinates) it should set this +parameter to 0. + +If the interpolator supports computing derivatives as described in +section~\ref{AEIThorns/AEILocalInterp/sect-derivatives}, {\em and\/} +if the interpolation molecule's size and/or shape varies with +\verb|operation_codes[]|, the interpolator sets the parameter +\begin{verbatim} +CCTK_INT MSS_is_fn_of_which_operation; +\end{verbatim} +to 1. Otherwise (\ie{} if the interpolator doesn't support computing +derivatives, or if the interpolator does support computing derivatives +but the interpolation molecule size and shape are independent of the +\verb|operation_code[]| values), the interpolator sets this parameter +to 0. Note that this query tests whether the molecule size and/or shape +depend on \verb|operation_codes[]| in general, independent of whether +there are in fact any distinct values (or even any values at all) passed +in \verb|operation_codes[]| in this particular interpolator call. In +other words, this is a query about the basic design of the interpolator, +not about this particular call. + +If the interpolation molecule's size and/or shape varies with the +actual floating-point values of the input arrays, the interpolator +sets the parameter +\begin{verbatim} +CCTK_INT MSS_is_fn_of_input_array_values; +\end{verbatim} +to 1. Otherwise (\ie{} if the interpolation molecule size and shape +are independent of the input array values; this is a necessary, but +not sufficient, condition for the interpolation to be linear), the +interpolator sets this parameter to 0. + +If the actual floating-point values of the +Jacobian~$(\ref{AEIThorns/AEILocalInterp/eqn-Jacobian})$ +(for a given \verb|out|, \verb|in|, and \verb|pt|) depend on the actual +floating-point values of the input arrays (\ie{} if the interpolation +is nonlinear), the interpolator sets the parameter +\begin{verbatim} +CCTK_INT Jacobian_is_fn_of_input_array_values; +\end{verbatim} +to 1. Otherwise (\ie{} if the interpolation is linear) the interpolator +sets this parameter to 0. + +The current implementation always sets +\begin{verbatim} +MSS_is_fn_of_interp_coords = 0 +MSS_is_fn_of_which_operation = 0 +MSS_is_fn_of_input_array_values = 0 +Jacobian_is_fn_of_input_array_values = 0 +\end{verbatim} + +%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Fixed-Size Hypercube-Shaped Molecules} +\label{AEIThorns/AEILocalInterp/sect-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 +(typically small) fixed size, independent of the interpolation +coordinates and the actual floating-point values in the input arrays +(though presumably depending on the interpolation order and on +\verb|operation_code|). In other words, this case applies if +(and only if) the Jacobian structure information described in +section~\ref{AEIThorns/AEILocalInterp/sect-Jacobian/structure} +returns +\begin{verbatim} +MSS_is_fn_of_interp_coords = 0 +MSS_is_fn_of_which_operation = 0 +MSS_is_fn_of_input_array_values = 0 +Jacobian_is_fn_of_input_array_values = 0 +\end{verbatim} +(These are precisely the values set by the current implementation.) +In the rest of this section we describe the query API for this case. + +The following parameters may be used to query the molecule size: +\begin{verbatim} +CCTK_INT const molecule_min_m[N_dims]; +CCTK_INT const molecule_max_m[N_dims]; +\end{verbatim} +The semantics of these are that if both of these keys are present +(the values and datatypes don't matter), then the interpolator will +(re)set the values to give the (inclusive) minimum and maximum +\verb|m|~molecule coordinates. (Note that either both of these +keys should be present, or neither should be present. This +simplifies the logic in the interpolator slightly.) + +The following parameter may be used to query the molecule positions: +\begin{verbatim} +CCTK_INT* const molecule_positions[N_dims]; +\end{verbatim} +The semantics of this is that the caller should set +\verb|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~$(\ref{AEIThorns/AEILocalInterp/eqn-Jacobian})$ itself: +\begin{verbatim} +CCTK_REAL* const Jacobian_pointer[N_output_arrays]; +const CCTK_INT Jacobian_offset [N_output_arrays]; # optional, default=all 0 + +/* the next 3 table entries are shared by all Jacobians */ +const CCTK_INT Jacobian_interp_point_stride; +const CCTK_INT Jacobian_m_strides[N_dims]; +const CCTK_INT Jacobian_part_stride; # optional, default=1 +\end{verbatim} +If \verb|Jacobian_pointer| is present in the table, then +\verb|Jacobian_interp_point_stride| and \verb|Jacobian_m_strides| +must also be present. For each \verb|out| where +\verb|Jacobian_pointer[out] != NULL|,%%% +\footnote{%%% + That is, setting {\tt Jacobian\_pointer[out] = NULL} + supresses the Jacobian query. If (as is often the + case) the Jacobian is independent of {\tt out}, then + it's a good idea to avoid unnecessary computations by + supressing the queries in this manner for any remaining + output arrays. + }%%% +{} the interpolator would then store the +Jacobian~$(\ref{AEIThorns/AEILocalInterp/eqn-Jacobian})$ in +\begin{verbatim} +Jacobian_pointer[out][offset + + pt*Jacobian_interp_point_stride + + mi*stride_i + mj*stride_j + mk*stride_k + + part*Jacobian_part_stride] +\end{verbatim} +where +\begin{verbatim} +offset = Jacobian_offset[out] +(stride_i,stride_j,stride_k) = Jacobian_m_strides[] +\end{verbatim} +and where \verb|part| is 0 for real values and the real parts of complex +values, and 1 for the imaginary parts of complex values. + +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 Lagrange polynomial interpolation, +with \verb|order = 2| and hence +(by table~\ref{AEIThorns/AEILocalInterp/tab-molecule-size+centering}) +using 3-point molecules. +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{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, using +(by table~\ref{AEIThorns/AEILocalInterp/tab-molecule-size+centering}) +4-point molecules. In other words, at each interpolation point we +use a cubic interpolation 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 maximum-degree Lagrange interpolation +already uses more data points than the number of free parameters in +the interpolating polynomials, \ie{} it already does some Savitzky-Golay +smoothing. For example, in 2-D a generic cubic polynomial +$f(x,y) = \sum_{i+j \le 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. + +The current implementation has all the framework for smoothing, but +only the \verb|smoothing=0| case is implemented at the moment. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Supressing Interpolation} + +Sometimes you may want to call the interpolator, but not actually +do any interpolation for some or all of the arrays: +\begin{itemize} +\item you may be doing a query +\item you may be just checking values in the parameter table +\item you may want to interpolate varying subsets of some set + of arrays +\item you may have no interpolation points, yet still want to + do the interpolator call for other reasons +\end{itemize} + +In such cases there are several possible ways you can supress some +or all of the interpolations: +\begin{description} +\item[\mathversion{bold} + Set $\text{\tt N\_interp\_points} = 0$:]\mbox{}\\ + This supresses all interpolation. In this case + \verb|interp_coords| may also be \verb|NULL| if desired. +\item[\mathversion{bold} + Set $\text{\tt N\_input\_arrays} = 0$ and + $\text{\tt N\_output\_arrays} = 0$:]\mbox{}\\ + This suppresses all interpolation. However, note that + some parameter-table entries like \verb|operand_indices| + and \verb|operation_codes| are specified as being arrays + whose size depends on \verb|N_output_arrays|, and it's an + error for these table entries to be present with the wrong + size (in this case, with any nonzero size). +\item[\mathversion{bold} + Set $\text{\tt output\_arrays[out]} = \text{\tt NULL}$:]\mbox{}\\ + This supresses the interpolation for this particular + output array. This is probably the cleanest way of + selectively turning individual arrays on and off. + If $\verb|output_arrays[out]| = \verb|NULL|$ and + the input arrays aren't needed for any queries you're doing, + then it's useful to also set the corresponding + $\verb|input_arrays[in]| = \verb|NULL|$, to supress the + interpolator fetching data (which in this case will never be used) + from the input arrays. +\end{description} + +Note that the combination +$\verb|input_arrays[in]| = \verb|NULL|$ and +$\verb|output_arrays[out]| \ne \verb|NULL|$ is an error: conceptually +you're asking for an interpolation (the output array pointer is +non-\verb|NULL|) but not supplying any input data (the input array +pointer is \verb|NULL|). + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Implementation} + +This interpolator's basic design is to use separate specialized code +for each combination of +\begin{verbatim} + (N_dims, molecule_family, order, smoothing) +\end{verbatim} +\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. + +Note that when compiling the code in the directory +\verb|src/GeneralizedPolynomial-Uniform| of this thorn, you may get +compiler warnings about casts discarding \verb|const| qualifiers from +pointers in (the \verb|#include|-ed file) \verb|template.c|. Don't +worry -- the code is actually ok.%%% +\footnote{%%% + The warnings are really a design bug in C's + const-qualifier rules when pointers point to + other pointers. See question~11.10 in the + (excellent!) {\tt comp.lang.c} online FAQ + ({\tt http://www.eskimo.com/~scs/C-faq/top.html}) + for further discussion. gcc~2.95.3 gives the + spurious warnings, but gcc~3.2.1 doesn't, nor + does the Intel C compiler in any version I've + tried. Also, \Cplusplus{} has fixed the problems + in C's const-qualifier rules, so the same code + compiled as \Cplusplus{} shouldn't give any + warnings (and doesn't on the compilers I've tried). + }%%% +{} You may also get compiler warnings about unused variables in +this same file; again don't worry, the code is ok. + +See the \verb|README| file in the source code directory +\verb|LocalInterp/src/UniformCartesian/| for further details on the +implementation. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Examples} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{A Simple Example} + +Here's a simple example in C, of interpolating a \verb|CCTK_REAL| and a +\verb|CCTK_COMPLEX| $10 \times 20$ 2-D array, at 5 interpolation points, +using cubic interpolation. + +{\bf Note:} +Since C allows arrays to be initialized only if the initializer +values are compile-time constants, we have to declare the +\verb|interp_coords[]|, \verb|input_arrays[]|, and \verb|output_arrays[]| +arrays as non-\verb|const|, and set their values with ordinary +(run-time) assignment statements. In \Cplusplus, there's no +restriction on initializer values, so we could declare the arrays +\verb|const| and initialize them as part of their declarations. + +\begin{verbatim} +#define N_DIMS 2 +#define N_INTERP_POINTS 5 +#define N_INPUT_ARRAYS 2 +#define N_OUTPUT_ARRAYS 2 + +/* (x,y) 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) coordinates of interpolation points */ +const CCTK_REAL interp_x[N_INTERP_POINTS]; +const CCTK_REAL interp_y[N_INTERP_POINTS]; +const void* interp_coords[N_DIMS]; /* see note above */ + +/* input arrays */ +/* ... note Cactus uses Fortran storage ordering, i.e. X is contiguous */ +#define NX 10 +#define NY 20 +const CCTK_REAL input_real [NY][NX]; +const CCTK_COMPLEX input_complex[NY][NX]; +const CCTK_INT input_array_dims[N_DIMS] = { NX, NY }; +const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS] + = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX }; +const void* input_arrays[N_INPUT_ARRAYS]; /* see note above */ + +/* 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* output_arrays[N_OUTPUT_ARRAYS]; /* see note above */ + +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!"); + +/* initialize the rest of the parameter arrays */ +interp_coords[0] = (const void *) interp_x; +interp_coords[1] = (const void *) interp_y; +input_arrays [0] = (const void *) input_real; +input_arrays [1] = (const void *) input_complex; +output_arrays[0] = ( void *) output_real; +output_arrays[1] = ( void *) output_complex; + +/* do the actual interpolation, and check for error returns */ +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 Interpolating Derivatives} +\label{AEIThorns/AEILocalInterp/sect-example-derivatives} + +In this example we interpolate the 3-metric and its 1st~derivatives +at a set of interpolation points on a 2-sphere. + +{\bf Note:} +Since we're not using \Cplusplus, we again declare the pointer arrays +non-\verb|const|, and ``initialize'' them with ordinary (run-time) +assignment statements. However, in this example, for greater clarity +we place these assignment statements right after the declarations. +Since C allows declarations only at the start of a \verb|{ }| block, +not in the middle of a block, we nest the rest of the program in extra +blocks (with the \verb|{ }| indented 2 spaces to distinguish them from +``normal'' \verb|{ }| pairs) to allow for further declarations. The +reader will have to decide for herself whether this style is more or +less ugly than separating the declarations and initializations of the +pointer arrays. + +\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* interp_coords[N_DIMS]; /* see note above */ +interp_coords[0] = (const void *) interp_x; +interp_coords[1] = (const void *) interp_y; +interp_coords[2] = (const void *) interp_z; + +/* dimensions of the data grid */ +#define NX 30 +#define NY 40 +#define NZ 50 + + { +/* input arrays */ +/* n.b. we use Fortran storage order: X is contiguous, Z least so */ +#define N_INPUT_ARRAYS 6 +const CCTK_REAL gxx_3D[NZ][NY][NX], gxy_3D[NZ][NY][NX], gxz_3D[NZ][NY][NX], + gyy_3D[NZ][NY][NX], gyz_3D[NZ][NY][NX], + gzz_3D[NZ][NY][NX]; + +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* input_arrays[N_INPUT_ARRAYS]; /* see note above */ +input_arrays[0] = (const void *) gxx_3D; +input_arrays[1] = (const void *) gxy_3D; +input_arrays[2] = (const void *) gxz_3D; +input_arrays[3] = (const void *) gyy_3D; +input_arrays[4] = (const void *) gyz_3D; +input_arrays[5] = (const void *) gzz_3D; + + { +/* output arrays */ +/* (for best efficiency we group all operations on a given input together) */ +#define N_OUTPUT_ARRAYS 24 +CCTK_REAL + gxx[N_INTERP_POINTS], + dx_gxx[N_INTERP_POINTS], dy_gxx[N_INTERP_POINTS], dz_gxx[N_INTERP_POINTS], + gxy[N_INTERP_POINTS], + dx_gxy[N_INTERP_POINTS], dy_gxy[N_INTERP_POINTS], dz_gxy[N_INTERP_POINTS], + gxz[N_INTERP_POINTS], + dx_gxz[N_INTERP_POINTS], dy_gxz[N_INTERP_POINTS], dz_gxz[N_INTERP_POINTS], + gyy[N_INTERP_POINTS], + dx_gyy[N_INTERP_POINTS], dy_gyy[N_INTERP_POINTS], dz_gyy[N_INTERP_POINTS], + gyz[N_INTERP_POINTS], + dx_gyz[N_INTERP_POINTS], dy_gyz[N_INTERP_POINTS], dz_gyz[N_INTERP_POINTS], + gzz[N_INTERP_POINTS], + dx_gzz[N_INTERP_POINTS], dy_gzz[N_INTERP_POINTS], dz_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 }; +void* output_arrays[N_OUTPUT_ARRAYS]; /* see note above */ +output_arrays[ 0] = (void *) gxx; +output_arrays[ 1] = (void *) dx_gxx; +output_arrays[ 2] = (void *) dy_gxx; +output_arrays[ 3] = (void *) dz_gxx; +output_arrays[ 4] = (void *) gxy; +output_arrays[ 5] = (void *) dx_gxy; +output_arrays[ 6] = (void *) dy_gxy; +output_arrays[ 7] = (void *) dz_gxy; +output_arrays[ 8] = (void *) gxz; +output_arrays[ 9] = (void *) dx_gxz; +output_arrays[10] = (void *) dy_gxz; +output_arrays[11] = (void *) dz_gxz; +output_arrays[12] = (void *) gyy; +output_arrays[13] = (void *) dx_gyy; +output_arrays[14] = (void *) dy_gyy; +output_arrays[15] = (void *) dz_gyy; +output_arrays[16] = (void *) gyz; +output_arrays[17] = (void *) dx_gyz; +output_arrays[18] = (void *) dy_gyz; +output_arrays[19] = (void *) dz_gyz; +output_arrays[20] = (void *) gzz; +output_arrays[21] = (void *) dx_gzz; +output_arrays[22] = (void *) dy_gzz; +output_arrays[23] = (void *) dz_gzz; + + { +/* integer codes to specify the derivatives */ +const CCTK_INT operand_indices[N_OUTPUT_ARRAYS] + = { 0, 0, 0, 0, + 1, 1, 1, 1, + 2, 2, 2, 2, + 3, 3, 3, 3, + 4, 4, 4, 4, + 5, 5, 5, 5 }; +#define DERIV(x) x +const CCTK_INT operation_codes[N_OUTPUT_ARRAYS] + = { DERIV(0), DERIV(1), DERIV(2), DERIV(3), + DERIV(0), DERIV(1), DERIV(2), DERIV(3), + DERIV(0), DERIV(1), DERIV(2), DERIV(3), + DERIV(0), DERIV(1), DERIV(2), DERIV(3), + DERIV(0), DERIV(1), DERIV(2), DERIV(3), + DERIV(0), DERIV(1), DERIV(2), DERIV(3) }; + +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} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\section{Acknowledgments} + +Thanks to Thomas Radke for our ongoing collaboration on Cactus +interpolation, +to Ian Hawke and Erik Schnetter for helpful comments on the documentation, +to Tom Goodale and Thomas Radke for many useful design discussions, +to Erik Schnetter for bug reports, +and to all the Cactus crew for a great infrastructure! + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Do not delete next line +% END CACTUS THORNGUIDE + +\end{document} diff --git a/doc/references b/doc/references new file mode 100644 index 0000000..2602c39 --- /dev/null +++ b/doc/references @@ -0,0 +1,15 @@ +This file contains references which are or might in the future relevant +to interpolation. + +Thermodynamic consistency in EOS interpolation +@ARTICLE{2000ApJS..126..501T, + author = {{Timmes}, F.~X. and {Swesty}, F.~D.}, + title = "{The Accuracy, Consistency, and Speed of an Electron-Positron Equation of State Based on Table Interpolation of the Helmholtz Free Energy}", + journal = {Astrophysical Journal Supplement Series}, + year = 2000, + month = feb, + volume = 126, + pages = {501--516}, + url = {http://esoads.eso.org/cgi-bin/nph-bib_query?bibcode=2000ApJS..126..501T&db_key=AST}, + adsnote = {Provided by the NASA Astrophysics Data System} +} |