diff options
| -rw-r--r-- | README | 3 | ||||
| -rw-r--r-- | doc/TODO | 5 | ||||
| -rw-r--r-- | doc/documentation.tex | 933 | ||||
| -rw-r--r-- | doc/future-ideas.tex | 61 | ||||
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c | 293 | ||||
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h | 39 | ||||
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/README | 29 | ||||
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/molecule_posn.c | 177 | ||||
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/template.c | 200 | ||||
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/template.h | 8 | ||||
| -rw-r--r-- | src/GeneralizedPolynomial-Uniform/test_molecule_posn.c | 447 |
11 files changed, 1327 insertions, 868 deletions
@@ -32,6 +32,9 @@ This thorn actually packages together two separate local interpolators: backwards-compatability synonym for "Lagrange polynomial interpolation".) This code lives in the subdirectory src/GeneralizedPolynomial-Uniform/ . +See the README files in the individual interpolators' directories +for more information. + We plan to eventually phase out the CCTK_InterpLocal() API. We may also phase out its interpolator, or we may move this to the newer CCTK_InterpLocalArrays() API. @@ -1,7 +1,6 @@ +$Header$ + Things to do on this thorn: -- test suite for GeneralizedPolynomial-Uniform/Hermite interpolator -- documentation for same (especially note very poor accuracy when - off-centered) - make casts from datatype to CCTK_REAL explicit in GeneralizedPolynomial-Uniform fetch routines - clean up the horrible inefficiency of the diff --git a/doc/documentation.tex b/doc/documentation.tex index 40adc63..fa59fec 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -93,6 +93,11 @@ % 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\thorn#1{{\bf #1}} + \def\Cplusplus{\hbox{C\hspace{-.05em}\raisebox{.4ex}{\tiny\bf ++}}} \def\ie{i.e.\hbox{}} \def\eg{e.g.\hbox{}} @@ -100,9 +105,9 @@ % Add an abstract for this thorn's documentation \begin{abstract} This thorn provides processor-local interpolation of 1-D, 2-D, and 3-D -data arrays. It provides interpolators for both the older -\verb|CCTK_InterpLocal()| API, and the newer and more generic -\verb|CCTK_InterpLocalUniform()| APIs. +data arrays. It supports both the older \verb|CCTK_InterpLocal()| API, +and the newer and more generic \verb|CCTK_InterpLocalUniform()| API. +This thorn guide describes how to use the various interpolators. \end{abstract} % The following sections are suggestive only. @@ -117,26 +122,27 @@ interpolators in the future): \begin{description} \item[Uniform Cartesian] This is the local interpolator which used to live - in the {\bf PUGHInterp} thorn. It was written in C by + in the \thorn{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.) + 2, and~3. \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. + 1~dimension (orders~1-6) and~2 and 3~dimensions (orders~1-4), + and Hermite polynomial interpolation in + 1~dimension (orders~2-4) and~2 and 3~dimensions (orders~2 and~3). + This interpolator supports a number of options, described below. \end{description} -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Terminology} +\subsection{Basic Terminology} +\label{localinterp/sect-basic-terminology} Within Cactus, each interpolator has a character-string name; this is mapped to a Cactus {\bf interpolator handle} by @@ -144,7 +150,7 @@ this is mapped to a Cactus {\bf interpolator handle} by 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 {\bf PUGHInterp}). +thorns such as \thorn{PUGHInterp}). Terminology for interpolation seems to differ a bit from one author to another. Here we refer to the {\bf point-centered} interpolation @@ -164,12 +170,12 @@ algorithm: \begin{verbatim} for each interpolation point { -choose a finite difference molecule position - somewhere near the interpolation point - for each input array +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 + the input data at the molecule points output = polynomial(interpolation point) } } @@ -179,21 +185,37 @@ 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 {\bf order} of the interpolation to be the order of -the fitted polynomial. That is, in our terminology, linear interpolation +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. +interpolator thus has $O\big((\Delta x)^{n+1})\big)$ interpolation errors +for generic smooth input data. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -However, because the interpolating polynomial generally changes if +\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. There are interpolation -algorithms (\eg{} Hermite polynomial and spline interpolation) which -give better smoothness, but none of the present interpolators implement -them. :( +and rises to a local maximum in each grid cell. This is the case, +for example, for Lagrange polynomial interpolation. + +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. See +section~\ref{localinterp/sect-Hermite} +for details. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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 @@ -228,8 +250,13 @@ z = interp_coords[2][pt] \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. +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 {\bf molecule position}, a nominal reference @@ -247,12 +274,12 @@ can be written 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~\ref{localinterp/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 -section~\ref{sect-generic-options}, we use \verb|const| qualifiers +section~\ref{localinterp/sect-generic-options}, 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 @@ -261,126 +288,540 @@ to return outputs back to the caller. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\section{Generic Interpolator Options} -\label{sect-generic-options} +\section{The Uniform Cartesian Interpolator} + +The uniform Cartesian interpolator is the one which used to live +in the \thorn{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{Generalized Polynomial Interpolation} +\label{localinterp/sect-generic-options} -The newer \verb|CCTK_InterpLocalUniform()| and -\verb|CCTK_InterpLocalNonUniform()| APIs specify a {\bf parameter table} -(a Cactus key-value table, manipulated by the \verb|Util_Table*()| APIs) -as a generic mechanism for passing optional or mandatory input/output -to/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. +The generalized polynomial interpolator was written in C +(plus Maple to generate the coefficients) by Jonathan Thornburg in +winter 2001-2002. It provides: +\begin{itemize} +\item Lagrange polynomial interpolation of orders 1--6 for 1-D arrays, + and 1--4 for 2- and 3-D arrays, registered under the operator + name \verb|"Lagrange polynomial interpolation"|.%%% +\footnote{%%% + For backwards compatability, the operator name + \hbox{\tt "generalized polynomial interpolation"} + is also allowed. + }%%% +\item Hermite polynomial interpolation of orders 2-4 for 1-D arrays, + and 2--3 for 2- and 3-D arrays, registered under the operator + name \verb|"Hermite polynomial interpolation"|.%%% +\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. -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. +This interpolator supports a number of options specified by a +{\bf 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{Interpolation Order} -The uniform Cartesian interpolator encodes the order in the operator -name, but other interpolators should use a (mandatory) parameter-table -parameter +The only mandatory parameter for this interpolator is the interpolation +order: \begin{verbatim} const CCTK_INT order; \end{verbatim} +As noted in section~\ref{localinterp/sect-basic-terminology}, +in our terminology linear interpolation is order~1, quadratic is +order~2, cubic is order~3, etc. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\subsection{Handling of Out-of-Range Interpolation Points} +\subsection{Molecule Size and Centering} + +If no grid boundaries or excised points are nearby, the interpolator +centers the molecules around the interpolation point as much as possible. +Figure~\ref{localinterp/fig-molecule-size+centering} gives the molecule size +and details of the centering policy for each interpolation operator +and order. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\begin{figure} + +\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{localinterp/fig-molecule-size+centering} +\end{figure} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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 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} -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.) +More generally, this option may be useful for controlling the extent +to which the interpolator uses input data from ghost and/or symmetry +zones. -Some interpolators may allow this behavior to be changed by the -optional parameter +%%%%%%%%%%%%%%%%%%%% + +\subsubsection{Off-Centering Molecules near Grid Boundaries + and/or Excision Regions} + +The parameter-table entries \begin{verbatim} -const CCTK_REAL out_of_range_tolerance[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 semantics of this are as follows: +control the local interpolator's willingness to off-center +(change away from the default centering) its interpolation molecules +near boundaries and excised points. -The array elements are matched up with the axes and minimum/maximum -ends of the grid in the order -$[x_{\min}, x_{\max}, y_{\min}, y_{\max}, z_{\min}, z_{\max}, \dots]$. -An array value \verb|TOL| is interpreted as follows: +To describe how these parameters work, it's useful to define two +regions in the data coordinate space: \begin{description} -\item[\rm If ${\tt TOL} \ge 0.0$,] - then an interpolation point is considered to be ``out of range'' - if and only if the interpolation point is - $> \verb|TOL| \times \verb|coord_delta[axis]|$ - outside the grid in this coordinate direction. -\item[\rm If ${\tt TOL} = -1.0$,]%%% -\footnote{%%% - Note that this is an {\em exact\/} floating-point - equality test! Although such tests are normally - very dangerous, this one is ok, because -1.0 can - be exactly represented in any reasonable floating-point - format. - }%%% -{} 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. +\item[The ``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 ``default-centering region''] + is that region in interpolation-point space where the + interpolator can use the default molecule centering + (described in figure~\ref{localinterp/fig-molecule-size+centering}), + \ie{} where the default-centering molecules require data + only from the data-valid region. \end{description} -Other values of \verb|TOL| are illegal (reserved for future use). +Figure~\ref{localinterp/fig-valid-data+default-centering} shows an +example of these regions. -To provide the ``fuzz'' noted above, \verb|out_of_range_tolerance[]| -should default to having all elements set to a small positive value, -say $10^4 \epsilon$, where $\epsilon$ is the ``machine epsilon''.%%% -\footnote{%%% - {\bf 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. - }%%% -{} However, at present all interpolators actually set the default -value to $10^4 \epsilon_d$, where $\epsilon_d$ is the machine epsilon -for C \verb|double| values, even though \verb|CCTK_REAL| may actually -be a different data type. (This is a bug, and may get fixed some -day\dots) - -If any interpolation points are out of range (as determined by the -\verb|out_of_range_tolerance[]| critera described above), the -interpolator 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{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{localinterp/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 ($\leq$) \verb|*_off_centering_tolerance[ibndry]| + grid spacings of the default-centering region, {\bf and} +\item within ($\leq$) \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. + +In practice the default-centering region is always a (normally proper) +subset of the valid-data region, so there are three cases of interest: +\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 ($\leq$) + \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%%% +\addtocounter{footnote}{-1}%%% +{} and + $\text{\tt *\!\_extrapolation\_tolerance[ibndry]} > 0.0$:]\mbox{}\\ + The interpolator allows interpolation points to be up to ($\leq$) + \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. +\end{description} + +(Note that setting $\verb|*_off_centering_tolerance[ibndry]| = 0.0$ and\\ +$\verb|*_extrapolation_tolerance[ibndry]| > 0.0$ is almost always a +mistake: 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 if you make such a setting.) + +If any of these table entries aren't specified, the defaults are +$\verb|*_off_centering_tolerance[ibndry]| = +\infty$%%% +\footnotemark%%% +\footnotetext{%%% + In practice any positive value larger than the molecule + radius is fine; $999.0$ is a conventional value here.%%% + }%%% +{} and $\verb|*_extrapolation_tolerance[ibndry]| = 10^{-10}$ for all +\verb|ibndry|, \ie{} interpolation points may be anywhere within +the valid-data region or up to $10^{-10}$ of a grid spacing outside it. +(The latter criterion avoids error returns if floating-point roundoff +moves an interpolation point which was supposed to be just on the +boundary of the valid-data region, slightly outside it.) + +%%%%%%%%%%%%%%%%%%%% + +\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 out of range? */ -/* (value is pt for out-of-range point) */ -CCTK_INT out_of_range_pt; +/* which interpolation point is it? */ +/* (value gives 0-origin pt for the offending point) */ +CCTK_INT error_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; +/* 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; -/* on which end of the grid is the point out of range? */ +/* 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 out_of_range_end; +CCTK_INT error_direction; \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. - }%%% +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{Molecule Family} -\label{sect-generic-options/molecule-family} +\label{localinterp/sect-molecule-family} An interpolator may support different families/shapes of interpolation molecules. Hypercube-shaped molecules are the simplest and most common @@ -389,29 +830,21 @@ or some generalization of this in higher numbers of dimensions. The parameter \begin{verbatim} -const char molecule_family[]; /* set with Util_TableSetString() */ +/* set with Util_TableSetString() */ +const char molecule_family[]; \end{verbatim} -may be used to query or change the strategy for choosing the molecules. +may be used to set or query the molecule family. -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 in the parameter table, this sets the +molecule family/shape based on that string. -If this key is present with a value which is a non-empty string, this -sets the molecule family/shape based on that string. +If this key isn't present in the parameter table, then the +interpolator sets it to the molecule family being used. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Non-Contiguous Input Arrays} -\label{sect-generic-options/non-contiguous-inputs} +\label{localinterp/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 @@ -419,14 +852,8 @@ 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.) -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: +The following (optional) parameter-table entries specify non-contiguous +input arrays: \begin{verbatim} const CCTK_INT input_array_offsets[N_input_arrays]; @@ -436,7 +863,7 @@ 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 +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] @@ -461,25 +888,21 @@ 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} -\label{sect-generic-options/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 +\label{localinterp/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, +running integrals, etc.) + +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]; @@ -495,7 +918,10 @@ 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. +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{localinterp/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 @@ -526,18 +952,19 @@ for the operation codes, by using the macro #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}. +There's an example of this in +section~\ref{localinterp/sect-example-derivatives}. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Jacobian and Domain of Dependence} -\label{sect-generic-options/Jacobian} +\label{localinterp/sect-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} + \label{localinterp/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 @@ -567,7 +994,7 @@ questions: were padded in this way. }%%% \item If this interpolator supports computing derivatives - as described in section~\ref{sect-generic-options/derivatives}, + as described in section~\ref{localinterp/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 @@ -586,19 +1013,17 @@ Jacobian and/or domain of dependence. %%%%%%%%%%%%%%%%%%%% \subsubsection{Determining the Jacobian's structure} -\label{sect-generic-options/Jacobian/structure} +\label{localinterp/sect-Jacobian/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 report its Jacobian structure -using the following parameter-table entries: +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{sect-generic-options/molecule-family}. +section~\ref{localinterp/sect-molecule-family}. If the interpolation molecule size and/or shape vary with the -interpolation coordinates, the interpolator should set the parameter +interpolation coordinates, the interpolator sets the parameter \begin{verbatim} CCTK_INT MSS_is_fn_of_interp_coords; \end{verbatim} @@ -607,9 +1032,9 @@ are independent of the interpolation coordinates) it should set this parameter to 0. If the interpolator supports computing derivatives as described -in section~\ref{sect-generic-options/derivatives}, and if the -interpolation molecule's size and/or shape varies with -\verb|operation_codes[]|, the interpolator should set the +in section~\ref{localinterp/sect-derivatives}, {\em and\/} +if the interpolation molecule's size and/or shape varies with +\verb|operation_codes[]|, the interpolator should sets the parameter \begin{verbatim} CCTK_INT MSS_is_fn_of_which_operation; @@ -636,7 +1061,7 @@ are independent of the input array values; this is a necessary, but not sufficient, condition for the interpolation to be linear), it should set this parameter to 0. -If the actual floating-point values of the Jacobian~$(\ref{eqn-Jacobian})$ +If the actual floating-point values of the Jacobian~$(\ref{localinterp/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 should set the parameter @@ -649,7 +1074,7 @@ this parameter to 0. %%%%%%%%%%%%%%%%%%%% \subsubsection{Fixed-Size Hypercube-Shaped Molecules} -\label{sect-generic-options/Jacobian/fixed-sized-hypercube} +\label{localinterp/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 @@ -658,7 +1083,7 @@ 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{sect-generic-options/Jacobian/structure} +section~\ref{localinterp/sect-Jacobian/structure} returns \begin{verbatim} MSS_is_fn_of_interp_coords = 0 @@ -690,7 +1115,7 @@ 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{eqn-Jacobian})$ itself: +Jacobian~$(\ref{localinterp/eqn-Jacobian})$ itself: \begin{verbatim} CCTK_REAL* const Jacobian_pointer[N_output_arrays]; const CCTK_INT Jacobian_offset [N_output_arrays]; /* optional */ @@ -710,7 +1135,7 @@ Then for each \verb|out| where \verb|Jacobian_pointer[out] != NULL|,%%% supressing the queries in this manner for any remaining output arrays. }%%% -{} the interpolator would store the Jacobian~$(\ref{eqn-Jacobian})$ in +{} the interpolator would store the Jacobian~$(\ref{localinterp/eqn-Jacobian})$ in \begin{verbatim} Jacobian_pointer[out][offset + pt*Jacobian_interp_point_stride @@ -763,153 +1188,8 @@ interp_molecule_positions[0][13] = 9 /* 1.00 [0.8, 0.9, 1.0] */ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -\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 {\bf 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} -\label{sect-generalized-polynomial-interp} - -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|"Lagrange polynomial interpolation"|.%%% -\footnote{%%% - For backwards compatability, the operator name - \hbox{\tt "generalized polynomial interpolation"} - is also allowed. - }%%% -{} (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 features described in -section~\ref{sect-generic-options}: -\begin{itemize} -\item interpolation order (this is a mandatory parameter) -\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}) -\item querying the interpolation's Jacobian and/or domain of dependence -\end{itemize} -It also supports the additional feature: -\begin{itemize} -\item Savitzky-Golay smoothing of the input data - as part of the interpolation process - (described in - section~\ref{sect-generalized-polynomial-interp/smoothing}) -\end{itemize} - -The interpolation order is a mandatory parameter which must be -present in the parameter table when the interpolator is called; -all the other parameter-table oparameters are optional. - -This interpolator does not support any caching features; at present -(May 2002) we have no plans to add these. - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - \subsection{Smoothing} -\label{sect-generalized-polynomial-interp/smoothing} +\label{localinterp/sect-smoothing} The way the generalized polynomial interpolator is implemented it's easy to also do Savitzky-Golay smoothing.%%% @@ -1080,7 +1360,7 @@ if (CCTK_InterpLocalUniform(N_DIMS, %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{An Example of Interpolating Derivatives} -\label{sect-example-derivatives} +\label{localinterp/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 @@ -1276,10 +1556,11 @@ if (CCTK_InterpLocalUniform(N_DIMS, \section{Acknowledgments} -Thanks to Thomas Radke for contributing the former {\bf PUGHInterp} +Thanks to Thomas Radke for contributing the former \thorn{PUGHInterp} \verb|CCTK_InterpLocal()| interpolator, -to Ian Hawke for helpful comments on the documentation, +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! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/doc/future-ideas.tex b/doc/future-ideas.tex new file mode 100644 index 0000000..9495c15 --- /dev/null +++ b/doc/future-ideas.tex @@ -0,0 +1,61 @@ +$Header$ + +% This file contains documentation for features which aren't implemented +% yet, and which are far enough in the future that I (Jonathan) decided +% to remove them from documentation.tex for the time being. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\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} diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c index 56f386a..1f63a03 100644 --- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c +++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.c @@ -38,7 +38,7 @@ #include "Hermite/all_prototypes.h" /* the rcs ID and its dummy function to use it */ -static const char *rcsid = "$Header$"; +static const char* rcsid = "$Header$"; CCTK_FILEVERSION(CactusBase_LocalInterp_src_GeneralizedPolynomialUniform_InterpLocalUniform_c) /******************************************************************************/ @@ -231,8 +231,20 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite, and calls the appropriate subfunction to do the actual interpolation, then finally stores some results back in the parameter table. + + FIXME: + This function is huge, and should be split up into + reasonable-sized subfunctions. @enddesc + @hdate 28.Jan.2003 + @hauthor Jonathan Thornburg <jthorn@aei.mpg.de> + @hdesc Support parameter-table entries + @var N_boundary_points_to_omit, + @var boundary_off_centering_tolerance, and + @var boundary_extrapolation_tolerance. + @endhdesc + ***** misc arguments ***** @var interp_operator @@ -349,26 +361,43 @@ return LocalInterp_InterpLocalUniform(interp_operator_Hermite, @vtype CCTK_INT order @endvar - @var out_of_range_tolerance - @vdesc Specifies how out-of-range interpolation points should - be handled. The array elements are matched up with - the axes and minimum/maximum ends of the grid in the + @var N_boundary_points_to_omit + @vdesc If this key is present, then this array specifies how many + input grid points to omit (not use for the interpolation) + on each grid boundary. The array elements are matched up + with the axes and minimum/maximum ends of the grid in the order [x_min, x_max, y_min, y_max, z_min, z_max, ...]. - An array value TOL is interpreted as follows: - If TOL >= 0.0, - then an interpolation point is considered to be - "out of range" if and only if the interpolation - point is > TOL * coord_delta[axis] - outside the grid in this coordinate direction. - If TOL == -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 - in this direction. - Other values of TOL are illegal. - @vtype const CCTK_REAL out_of_range_tolerance[2*N_dims] + If this key isn't present, it defaults to all zeros, i.e. + to use all the input grid points in the interpolation. + @vtype const CCTK_INT N_boundary_points_to_omit[2*N_dims] + @endvar + + @var boundary_off_centering_tolerance + @vdesc If this key is present, then the interpolator allows + interpolation points to be up to (<=) this many grid + spacings outside the default-centering region on each + grid boundary, off-centering the interpolation molecules + as necessary. + The array elements are matched up with the axes and + minimum/maximum ends of the grid in the order + [x_min, x_max, y_min, y_max, z_min, z_max, ...]. + If this key isn't present, it defaults to all infinities, + i.e. to place no restriction on the interpolation point. + @vtype const CCTK_REAL boundary_off_centering_tolerance[2*N_dims] + @endvar + + @var boundary_extrapolation_tolerance + @vdesc If this key is present, then the interpolator allows + interpolation points to be up to (<=) this many grid + spacings outside the valid-data region on each grid + boundary, doing extrapolation instead of interpolation + as necessary. + The array elements are matched up with the axes and + minimum/maximum ends of the grid in the order + [x_min, x_max, y_min, y_max, z_min, z_max, ...]. + If this key isn't present, it defaults to all 1.0e-10, + i.e. to allow up to 1e-10 grid spacings of extrapolation. + @vtype const CCTK_REAL boundary_extrapolation_tolerance[2*N_dims] @endvar @var input_array_offsets @@ -625,7 +654,7 @@ static * would falsely detect as an out-of-memory condition. * * As a workaround, we pad all our malloc request sizes, i.e. - * CCTK_INT *const p = malloc((N+1) * sizeof(CCTK_INT)); + * CCTK_INT* const p = (CCTK_INT*) malloc((N+1) * sizeof(CCTK_INT)); * if (p == NULL) * then return UTIL_ERROR_NO_MEMORY * This is a kludge, but so are the other plausible solutions. :( :( @@ -657,7 +686,7 @@ if ( (N_dims <= 0) then { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, -"InterpLocalUniform(): invalid arguments\n" +"CCTK_InterpLocalUniform(): invalid arguments\n" " (N_dims <= 0, param_table_handle < 0, N_input_arrays < 0,\n" " N_output_arrays < 0, and/or NULL pointers-that-shouldn't-be-NULL)!"); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ @@ -684,14 +713,14 @@ if (N_dims > MAX_N_DIMS) CCTK_INT order; status = Util_TableGetInt(param_table_handle, &order, "order"); if (status == 1) - then { } /* ok ==> no-op here */ + then { } /* value set from table ==> no-op here */ else { /* n.b. order is a mandatory parameter (no default)!*/ CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad or missing table entry\n" -" for \"order\" parameter!\n" +" CCTK_InterpLocalUniform(): bad or missing table entry for\n" +" \"order\" parameter!\n" " (this is a mandatory parameter)\n" " error status=%d", status); @@ -708,59 +737,122 @@ if ((order < 1) || (order > MAX_ORDER)) MAX_ORDER); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } - + +/******************************************************************************/ /* * out-of-range interpolation-point handling */ { -CCTK_REAL out_of_range_tolerance[2*MAX_N_DIMS]; -const int N_tolerances = 2*N_dims; -status = Util_TableGetRealArray(param_table_handle, - N_tolerances, out_of_range_tolerance, - "out_of_range_tolerance"); +CCTK_INT N_boundary_points_to_omit[MAX_N_BOUNDARIES]; +CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES]; +CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES]; +const int N_boundaries = 2*N_dims; + +status = Util_TableGetIntArray(param_table_handle, + N_boundaries, N_boundary_points_to_omit, + "N_boundary_points_to_omit"); if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) then { - /* default */ - int tol_index; - for (tol_index = 0 ; tol_index < N_tolerances ; ++tol_index) + /* set the default value */ + int ibndry; + for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) { - out_of_range_tolerance[tol_index] - = OUT_OF_RANGE_TOLERANCE_DEFAULT; + N_boundary_points_to_omit[ibndry] = 0; } } -else if (status == N_tolerances) +else if (status == N_boundaries) + then { } /* value set from table ==> no-op here */ +else { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"N_boundary_points_to_omit\" parameter!\n" +" error status=%d", + status); + return status; /*** ERROR RETURN ***/ + } + +/**************************************/ + +status = Util_TableGetRealArray(param_table_handle, + N_boundaries, boundary_off_centering_tolerance, + "boundary_off_centering_tolerance"); +if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) then { - /* check that all values are valid */ - int tol_index; - for (tol_index = 0 ; tol_index < N_tolerances ; ++tol_index) + /* set the default value */ + int ibndry; + for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) { - if (! ( (out_of_range_tolerance[tol_index] >= 0.0) - || (out_of_range_tolerance[tol_index] == -1.0) ) ) - then { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, + boundary_off_centering_tolerance[ibndry] + = BOUNDARY_OFF_CENTER_TOL_DEFAULT; + } + } +else if (status == N_boundaries) + then { } /* value set from table ==> no-op here */ +else { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform():\n" -" invalid out_of_range_tolerance[tol_index=%d] = %g!\n" -" (valid values are -1.0 or >= 0.0)", - tol_index, - out_of_range_tolerance[tol_index]); - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - } +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"boundary_off_centering_tolerance\" parameter!\n" +" error status=%d", + status); + return status; /*** ERROR RETURN ***/ + } + +/**************************************/ + +status = Util_TableGetRealArray(param_table_handle, + N_boundaries, boundary_extrapolation_tolerance, + "boundary_extrapolation_tolerance"); +if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) + then { + /* set the default value */ + int ibndry; + for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) + { + boundary_extrapolation_tolerance[ibndry] + = BOUNDARY_EXTRAP_TOL_DEFAULT; } } +else if (status == N_boundaries) + then { } /* value set from table ==> no-op here */ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"out_of_range_tolerance\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"boundary_extrapolation_tolerance\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } +/* check for almost-always-a-mistake settings, warn if found */ + { +int ibndry; + for (ibndry = 0 ; ibndry < N_boundaries ; ++ibndry) + { + if ( (boundary_off_centering_tolerance[ibndry] == 0.0) + && (boundary_extrapolation_tolerance[ibndry] > 0.0) ) + then CCTK_VWarn(WARNING_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): warning: setting\n" +" boundary_off_centering_tolerance[%d] == 0.0\n" +" and\n" +" boundary_extrapolation_tolerance[%d] > 0.0\n" +" is almost certainly a mistake\n" +" (the boundary_extrapolation_tolerance[%d]\n" +" setting will have no effect because of the\n" +" boundary_off_centering_tolerance[%d] setting)" + , + ibndry, ibndry, ibndry, ibndry); + } + } + /******************************************************************************/ /* @@ -768,23 +860,23 @@ else { */ { #define MOLECULE_FAMILY_BUFSIZ (MAX_MOLECULE_FAMILY_STRLEN+1) -char molecule_family_string[MOLECULE_FAMILY_BUFSIZ]; +char molecule_family_strbuf[MOLECULE_FAMILY_BUFSIZ]; +char* molecule_family_str = molecule_family_strbuf; enum molecule_family molecule_family; status = Util_TableGetString(param_table_handle, - MOLECULE_FAMILY_BUFSIZ, molecule_family_string, + MOLECULE_FAMILY_BUFSIZ, molecule_family_strbuf, "molecule_family"); if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) then { - /* default */ + /* set the default value for code here*/ /* (need enum for main code, string for error messages) */ - strcpy(molecule_family_string, "cube"); - molecule_family = molecule_family_cube; - } -else if (status == 0) - then { - /* this is a query of our default molecule family */ + molecule_family = molecule_family_cube; + molecule_family_str = "cube"; + + /* set this key in the parameter table */ + /* to give the (default) molecule family we're going to use */ status = Util_TableSetString(param_table_handle, - "cube", + molecule_family_str, "molecule_family"); if (status < 0) then { @@ -801,13 +893,13 @@ else if (status == 0) else if (status > 0) then { /* decode molecule family string */ - if (strcmp(molecule_family_string, "cube") == 0) + if (strcmp(molecule_family_strbuf, "cube") == 0) then molecule_family = molecule_family_cube; else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, -"InterpLocalUniform(): unknown molecule_family=\"%s\"!", - molecule_family_string); +"CCTK_InterpLocalUniform(): unknown molecule_family=\"%s\"!", + molecule_family_strbuf); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } } @@ -815,13 +907,15 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"molecule_family\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"molecule_family\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } +/**************************************/ + /* * smoothing */ @@ -829,14 +923,15 @@ else { CCTK_INT smoothing; status = Util_TableGetInt(param_table_handle, &smoothing, "smoothing"); if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) - then smoothing = 0; /* default */ + then smoothing = 0; /* set default value */ else if (status == 1) then { } /* value set from table ==> no-op here */ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry for \"smoothing\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"smoothing\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ @@ -850,8 +945,8 @@ if ((smoothing < 0) || (smoothing > MAX_SMOOTHING)) * input array offsets */ { -CCTK_INT *const input_array_offsets - = malloc(N_input_arrays1 * sizeof(CCTK_INT)); +CCTK_INT* const input_array_offsets + = (CCTK_INT*) malloc(N_input_arrays1 * sizeof(CCTK_INT)); if (input_array_offsets == NULL) then return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ status = Util_TableGetIntArray(param_table_handle, @@ -869,13 +964,15 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"input_array_offsets\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"input_array_offsets\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } +/**************************************/ + /* * input array strides */ @@ -915,13 +1012,15 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"input_array_offsets\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"input_array_offsets\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } +/**************************************/ + /* * input array min/max subscripts */ @@ -942,8 +1041,8 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"input_array_min_subscripts\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"input_array_min_subscripts\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ @@ -981,8 +1080,8 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"input_array_max_subscripts\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"input_array_max_subscripts\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ @@ -994,8 +1093,8 @@ else { * operand indices */ { -CCTK_INT *const operand_indices - = malloc(N_output_arrays1 * sizeof(CCTK_INT)); +CCTK_INT* const operand_indices + = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (operand_indices == NULL) then { free(input_array_offsets); @@ -1059,19 +1158,21 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"operand_indices\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"operand_indices\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ } +/**************************************/ + /* * operation_codes */ { -CCTK_INT *const operation_codes - = malloc(N_output_arrays1 * sizeof(CCTK_INT)); +CCTK_INT* const operation_codes + = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (operation_codes == NULL) then { free(operand_indices); @@ -1110,8 +1211,8 @@ else { CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" -" CCTK_InterpLocalUniform(): bad table entry\n" -" for \"operation_codes\" parameter!\n" +" CCTK_InterpLocalUniform(): bad table entry for\n" +" \"operation_codes\" parameter!\n" " error status=%d", status); return status; /*** ERROR RETURN ***/ @@ -1223,7 +1324,7 @@ if (status) then { /* yes, we're doing Jacobian queries */ Jacobian_info.Jacobian_pointer - = (CCTK_REAL **) malloc(N_output_arrays1 * sizeof(CCTK_REAL *)); + = (CCTK_REAL**) malloc(N_output_arrays1 * sizeof(CCTK_REAL *)); if (Jacobian_info.Jacobian_pointer == NULL) then { free(operation_codes); @@ -1232,8 +1333,8 @@ if (status) return UTIL_ERROR_NO_MEMORY; /*** ERROR RETURN ***/ } { - CCTK_POINTER *Jacobian_pointer_temp - = (CCTK_POINTER *) malloc(N_output_arrays1 * sizeof(CCTK_POINTER)); + CCTK_POINTER* Jacobian_pointer_temp + = (CCTK_POINTER*) malloc(N_output_arrays1 * sizeof(CCTK_POINTER)); if (Jacobian_pointer_temp == NULL) then { free(operation_codes); @@ -1276,7 +1377,7 @@ if (status) } Jacobian_info.Jacobian_offset - = (CCTK_INT *) malloc(N_output_arrays1 * sizeof(CCTK_INT)); + = (CCTK_INT*) malloc(N_output_arrays1 * sizeof(CCTK_INT)); if (Jacobian_info.Jacobian_offset == NULL) then { free(operation_codes); @@ -1550,7 +1651,7 @@ if (p_interp_fn == NULL_INTERP_FN_PTR) " interp_operator=\"%s\", N_dims=%d\n" " molecule_family=\"%s\", order=%d, smoothing=%d", interp_operator_string, N_dims, - molecule_family_string, order, smoothing); + molecule_family_str, order, smoothing); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } @@ -1564,7 +1665,9 @@ const int return_code = (*p_interp_fn)(coord_origin, coord_delta, N_interp_points, interp_coords_type_code, interp_coords, - out_of_range_tolerance, + N_boundary_points_to_omit, + boundary_off_centering_tolerance, + boundary_extrapolation_tolerance, N_input_arrays, input_array_offsets, input_array_strides, input_array_min_subscripts, diff --git a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h index d2990df..851378a 100644 --- a/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h +++ b/src/GeneralizedPolynomial-Uniform/InterpLocalUniform.h @@ -70,11 +70,17 @@ typedef int bool; */ #define MAX_N_DIMS 3 +/* a "boundary" is the combination of a dimension and a min/max "side" */ +#define MAX_N_BOUNDARIES (2*MAX_N_DIMS) + /* * if molecule_family_string is a C string specifying a molecule family * (i.e. value associated with the "molecule_family" key in the parameter * table), we must have * strlen(molecule_family_string) <= MAX_MOLECULE_FAMILY_STRLEN + * n.b. exceeding this won't cause a buffer overflow, it will "just" + * cause the string to be truncated (and probably not recognized + * by the interpolator) */ #define MAX_MOLECULE_FAMILY_STRLEN 20 @@ -103,12 +109,15 @@ enum molecule_family * other compile-time settings */ -/* default for each element of out_of_range_tolerance[] */ -/* FIXME: this is based on C "double", which may not match CCTK_REAL */ -#define OUT_OF_RANGE_TOLERANCE_DEFAULT (10000.0*DBL_EPSILON) +/* default for each element of boundary_off_centering_tolerance[] */ +#define BOUNDARY_OFF_CENTER_TOL_DEFAULT 999.0 + +/* default for each element of boundary_extrapolation_tolerance[] */ +#define BOUNDARY_EXTRAP_TOL_DEFAULT 1.0e-10 -/* CCTK_VWarn() severity level for error messages */ -#define ERROR_MSG_SEVERITY_LEVEL 1 +/* CCTK_VWarn() severity level for error/warning messages */ +#define ERROR_MSG_SEVERITY_LEVEL 0 +#define WARNING_MSG_SEVERITY_LEVEL 1 /******************************************************************************/ @@ -122,6 +131,7 @@ enum molecule_family struct error_flags { int error_pt; + int error_ibndry; int error_axis; int error_end; }; @@ -153,20 +163,16 @@ struct Jacobian_info /* * error codes for LocalInterp_molecule_posn() - * ... these are defined in terms of INT_MIN from <limits.h> */ /* x < minimum allowable x in grid */ -#define MOLECULE_POSN_ERROR_X_LT_MIN (INT_MIN + 0) +#define MOLECULE_POSN_ERROR_X_LT_MIN (-1) /* x > maximum allowable x in grid */ -#define MOLECULE_POSN_ERROR_X_GT_MAX (INT_MIN + 1) +#define MOLECULE_POSN_ERROR_X_GT_MAX (-2) /* grid is smaller than molecule */ -#define MOLECULE_POSN_ERROR_GRID_TINY (INT_MIN + 2) - -/* is a given integer an error code? */ -#define IS_MOLECULE_POSN_ERROR_CODE(x) (x <= MOLECULE_POSN_ERROR_GRID_TINY) +#define MOLECULE_POSN_ERROR_GRID_TINY (-3) /******************************************************************************/ @@ -216,11 +222,12 @@ int LocalInterp_ILU_Hermite(int N_dims, int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta, int grid_i_min, int grid_i_max, int molecule_size, - fp out_of_range_tolerance_min, - fp out_of_range_tolerance_max, + fp boundary_off_centering_tolerance_min, + fp boundary_off_centering_tolerance_max, + fp boundary_extrapolation_tolerance_min, + fp boundary_extrapolation_tolerance_max, fp x, - fp *x_rel, - int *molecule_m_min, int *molecule_m_max); + int* i_center, fp* x_rel); /* functions in util.c */ void LocalInterp_zero_int_array(int N, CCTK_INT array[]); diff --git a/src/GeneralizedPolynomial-Uniform/README b/src/GeneralizedPolynomial-Uniform/README index 0963854..b022489 100644 --- a/src/GeneralizedPolynomial-Uniform/README +++ b/src/GeneralizedPolynomial-Uniform/README @@ -7,7 +7,8 @@ under the name The source code files are as follows: -* startup.c registers the interpolation operator +* startup.c registers the interpolation operators +* InterpLocalUniform.h is an overall header file for the whole interpolator * InterpLocalUniform.c is the top-level driver: it gets the various options from the parameter table, then decodes (N_dims, molecule_family, order, smoothing) @@ -16,6 +17,7 @@ The source code files are as follows: * [123]d.cube.order*.smooth*.c define the individual interpolation subfunctions. Each of them just #defines a whole bunch of macros, then #includes template.c (which has the actual code). +* template.h defines a prototype for the function defined by template.c * template.c is the actual interpolation code. It is written in terms of a large number of macros, which should be #defined before #including template.c. There's a long block comment @@ -26,6 +28,10 @@ The source code files are as follows: Maple-generated code fragments in the [123]d.coeffs/directories; template.c uses various macros to tell it which fragments to #include. +* molecule_posn.c contains the LocalInterp_molecule_posn() function + to compute where in the grid each (an) interpolation molecule should + be centered for each (a given) interpolation point. +* util.c contains some low-level utility routines * [123]d.maple are the top-level Maple code files; they call various functions in interpolate.maple and util.maple to do the actual work. @@ -33,7 +39,19 @@ The source code files are as follows: computing an interpolant and manipulating/printing it in various ways * util.maple contains low-level utility routines +* makefile is a makefile for... +* test_molecule_posn.c is a standalone test driver for the code + in molecule_posn.c . By default it runs a large set of (around 100) + test cases stored in a table in the code. +The subdirectories are as follows: +* [123]d.coeffs/ are empty directory trees left-over from previous + versions of this interpolator (CVS doesn't have any easy way to + remove directories) +* Lagrange/ contains the code for the Lagrange polynomial interpolator. +* Hermite/ contains the code for the Hermite polynomial interpolator. +* common/ contains low-level code common to both the Lagrange and + the Hermite interpolators. To add a new combination of (N_dims, molecule_family, order, smoothing), @@ -60,17 +78,10 @@ to this interpolator, you need to -Other makefile targets: -test_molecule_posn - This makes a standalone test driver for the - LocalInterp_molecule_posn() - function defined in molecule_posn.c - - Copyright ========= -This interpolator is copyright (C) 2002-2003 +This interpolator is copyright (C) 2001-2003 by Jonathan Thornburg <jthorn@aei.mpg.de>. This interpolator is free software; you can redistribute it and/or modify diff --git a/src/GeneralizedPolynomial-Uniform/molecule_posn.c b/src/GeneralizedPolynomial-Uniform/molecule_posn.c index e007396..6d42ef3 100644 --- a/src/GeneralizedPolynomial-Uniform/molecule_posn.c +++ b/src/GeneralizedPolynomial-Uniform/molecule_posn.c @@ -8,7 +8,6 @@ @@*/ #include <math.h> -#include <limits.h> #include <stdlib.h> #include <stdio.h> @@ -73,6 +72,18 @@ static const char *rcsid = "$Header$"; center, i.e. the above diagram has columns labelled with [i+m]. @enddesc + @hdate 27.Jan.2003 + @hauthor Jonathan Thornburg <jthorn@aei.mpg.de> + @hdesc Complete rewrite: now supports + @var boundary_off_centering_tolerance_min, + @var boundary_off_centering_tolerance_max, + @var boundary_extrapolation_tolerance_min, + @var boundary_extrapolation_tolerance_max, + also change to returning status code, + also drop returning @var min_m and @var max_m + since they were never used. + @endhdesc + @var grid_origin @vdesc The floating-point coordinate x of the grid point i=0. @vtype fp grid_origin @@ -98,50 +109,53 @@ static const char *rcsid = "$Header$"; @vtype int molecule_size @endvar - @var out_of_range_tolerance_min, out_of_range_tolerance_max - @vdesc Specifies how out-of-range interpolation points should - be handled for the {minimum,maximum} ends of the grid - respectively. - If out_of_range_tolerance >= 0.0, - then an interpolation point is considered to be - "out of range" if and only if the interpolation - point is > out_of_range_tolerance * grid_delta - outside the grid in any coordinate. - If out_of_range_tolerance == -1.0, - then an interpolation point is considered to be - "out of range" if and only if a centered molecule - would require data from outside the grid. - Other values of out_of_range_tolerance are illegal. - @vtype fp out_of_range_tolerance_min, out_of_range_tolerance_max + @var boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max + @vdesc Specifies how many grid spacings the interpolation point + may be beyond the default-centering region before being + declared "out of range" on the {minimum,maximum} boundaries + of the grid respectively. + @vtype fp boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max + @endvar + + @var boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max + @vdesc Specifies how many grid spacings the interpolation point + may be beyond the grid boundary before being declared + "out of range" on the {minimum,maximum} boundaries + of the grid respectively. + @vtype fp boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max @endvar @var x - @vdesc The floating-point coordinate x at which we would like to center - the molecule. + @vdesc The floating-point coordinate of the interpolation point. @vtype fp x @endvar - @var x_rel - @vdesc A pointer to a floating-point value where this function should - store the x coordinate relative to the molecule center, measured - in units of the grid spacing; or NULL to skip storing this. - @vtype fp *x_rel + @var i_center + @vdesc A pointer to an value where this function should + store the integer coordinate of the molecule center, + or NULL to skip storing this + @vtype int *i_center @vio pointer to out @endvar - @var molecule_m_min, molecule_m_max - @vdesc A pointer to an int where this function should store the - {minimum,maximum} molecule coordinate m of the molecule; + @var x_rel + @vdesc A pointer to where this function should store the + interpolation point's x coordinate relative to the + molecule center, measured in units of the grid spacing; or NULL to skip storing this. - @vtype int *molecule_m_min, *molecule_m_max + @vtype fp *x_rel @vio pointer to out @endvar @returntype int @returndesc - This function returns the integer coordinate i_center at which - the molecule should be centered, or one of the error codes defined - in "InterpLocalUniform.h": + This function returns 0 if the interpolation point is "in range", + or one of the (negative) error codes defined in "InterpLocalUniform.h" + if the interpolation point is "out of range": MOLECULE_POSN_ERROR_X_LT_MIN if x is out-of-range on the min end of the grid (i.e. x < the minimum allowable x) @@ -155,73 +169,70 @@ static const char *rcsid = "$Header$"; int LocalInterp_molecule_posn(fp grid_origin, fp grid_delta, int grid_i_min, int grid_i_max, int molecule_size, - fp out_of_range_tolerance_min, - fp out_of_range_tolerance_max, + fp boundary_off_centering_tolerance_min, + fp boundary_off_centering_tolerance_max, + fp boundary_extrapolation_tolerance_min, + fp boundary_extrapolation_tolerance_max, fp x, - fp *x_rel, - int *molecule_m_min, int *molecule_m_max) + int* i_center, fp* x_rel) { -/* min/max molecule coordinates */ -const int m_max = (molecule_size >> 1); /* a positive number */ -const int m_min = m_max - molecule_size + 1; /* a negative number */ +/* molecule radia (inherently positive) in +/- directions */ +const int mr_plus = (molecule_size >> 1); +const int mr_minus = molecule_size - mr_plus - 1; -/* integer coordinate i of point x, as a floating-point number */ -const fp fp_i = (x - grid_origin) / grid_delta; +/* range of x_rel for which this molecule is centered, cf. diagram in header comment */ +const fp centered_min_x_rel = IS_EVEN(molecule_size) ? 0.0 : -0.5; +const fp centered_max_x_rel = IS_EVEN(molecule_size) ? 1.0 : +0.5; -/* is point x out-of-range? */ -if (out_of_range_tolerance_min >= 0.0) - then { - const fp fp_effective_grid_i_min - = ((fp) grid_i_min) - out_of_range_tolerance_min; - if (fp_i < fp_effective_grid_i_min) - then return INT_MIN; /*** ERROR RETURN ***/ - } -if (out_of_range_tolerance_max >= 0.0) - then { - const fp fp_effective_grid_i_max - = ((fp) grid_i_max) + out_of_range_tolerance_max; - if (fp_i > fp_effective_grid_i_max) - then return INT_MAX; /*** ERROR RETURN ***/ - } +/* range of i where we *could* center the molecule, as floating-point numbers */ +const fp fp_centered_min_possible_i = grid_i_min + mr_minus + centered_min_x_rel; +const fp fp_centered_max_possible_i = grid_i_max - mr_plus + centered_max_x_rel; -/* integer coordinate i where we will center the molecule */ -/* (see diagram in header comment for explanation) */ -/* ... as a floating-point number */ - { -const fp fp_i_center = IS_EVEN(molecule_size) ? floor(fp_i) : JT_ROUND(fp_i); -/* ... as an integer */ -int i_center = (int) fp_i_center; +/* integer coordinate i of interpolation point, as a floating-point number */ +const fp fp_i = (x - grid_origin) / grid_delta; -/* min/max integer coordinates in molecule assuming it's fully centered */ -const int centered_i_min = i_center + m_min; -const int centered_i_max = i_center + m_max; +/* is the molecule larger than the grid? */ +if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max)) + then return MOLECULE_POSN_ERROR_GRID_TINY; /*** ERROR RETURN ***/ -/* check if off-centered molecules are forbidden? */ -if ((out_of_range_tolerance_min == -1.0) && (centered_i_min < grid_i_min)) +/* is interpolation point x beyond the extrapolation tolerance? */ +if (fp_i < (fp)grid_i_min - boundary_extrapolation_tolerance_min) then return MOLECULE_POSN_ERROR_X_LT_MIN; /*** ERROR RETURN ***/ -if ((out_of_range_tolerance_max == -1.0) && (centered_i_max > grid_i_max)) +if (fp_i > (fp)grid_i_max + boundary_extrapolation_tolerance_max) then return MOLECULE_POSN_ERROR_X_GT_MAX; /*** ERROR RETURN ***/ -/* off-center as needed if we're close to the edge of the grid */ +/* is interpolation point x beyond the off-centering tolerance? */ +if (fp_i < fp_centered_min_possible_i - boundary_off_centering_tolerance_min) + then return MOLECULE_POSN_ERROR_X_LT_MIN; /*** ERROR RETURN ***/ +if (fp_i > fp_centered_max_possible_i + boundary_off_centering_tolerance_max) + then return MOLECULE_POSN_ERROR_X_GT_MAX; /*** ERROR RETURN ***/ + +/* now choose the actual molecule position/centering: */ +/* first set up a centered molecule */ { -int shift = 0; -if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max)) - then return MOLECULE_POSN_ERROR_GRID_TINY; /*** ERROR RETURN ***/ -if (centered_i_min < grid_i_min) - then shift = grid_i_min - centered_i_min; /* a positive number */ -if (centered_i_max > grid_i_max) - then shift = grid_i_max - centered_i_max; /* a negative number */ -i_center += shift; +fp fp_i_center = IS_EVEN(molecule_size) /* ... as a floating-point number */ + ? floor(fp_i) + : JT_ROUND(fp_i); +int int_i_center = (int) fp_i_center; /* ... as an integer */ + +/* then clamp the molecule at the grid boundaries */ +if (int_i_center - mr_minus < grid_i_min) + then { + int_i_center = grid_i_min + mr_minus; + fp_i_center = (fp) int_i_center; + } +if (int_i_center + mr_plus > grid_i_max) + then { + int_i_center = grid_i_max - mr_plus; + fp_i_center = (fp) int_i_center; + } /* store the results */ +if (i_center != NULL) + then *i_center = int_i_center; if (x_rel != NULL) - then *x_rel = fp_i - ((fp) i_center); -if (molecule_m_min != NULL) - then *molecule_m_min = m_min; -if (molecule_m_max != NULL) - then *molecule_m_max = m_max; + then *x_rel = fp_i - fp_i_center; -return i_center; - } +return 0; /*** NORMAL RETURN ***/ } } diff --git a/src/GeneralizedPolynomial-Uniform/template.c b/src/GeneralizedPolynomial-Uniform/template.c index 0bb409b..a567119 100644 --- a/src/GeneralizedPolynomial-Uniform/template.c +++ b/src/GeneralizedPolynomial-Uniform/template.c @@ -20,6 +20,13 @@ 3-D molecules. @endhdesc + @hdate 28.Jan.2003 + @hauthor Jonathan Thornburg <jthorn@aei.mpg.de> + @hdesc Support parameter-table entries + @var boundary_off_centering_tolerance and + @var boundary_extrapolation_tolerance. + @endhdesc + @desc The following header files must be #included before #including this file: @@ -301,6 +308,9 @@ only briefly; see the InterpLocalUniform() documentation and/or the thorn guide for this thorn for further details. + If you change the arguments for this function, note that you + must also change the prototype in "template.h". + @var error_flags @vdesc If we encounter an out-of-range point and this pointer is non-NULL, we store a description of the out-of-range @@ -346,13 +356,15 @@ @@*/ int FUNCTION_NAME(/***** coordinate system *****/ - const CCTK_REAL coord_system_origin[], - const CCTK_REAL grid_spacing[], + const CCTK_REAL coord_origin[], + const CCTK_REAL coord_delta[], /***** interpolation points *****/ int N_interp_points, int interp_coords_type_code, const void* const interp_coords[], - const CCTK_REAL out_of_range_tolerance[], + const CCTK_INT N_boundary_points_to_omit[], + const CCTK_REAL boundary_off_centering_tolerance[], + const CCTK_REAL boundary_extrapolation_tolerance[], /***** input arrays *****/ int N_input_arrays, const CCTK_INT input_array_offsets[], @@ -691,8 +703,8 @@ int out; */ { #if N_DIMS >= 1 - const fp origin_x = coord_system_origin[X_AXIS]; - const fp delta_x = grid_spacing[X_AXIS]; + const fp origin_x = coord_origin[X_AXIS]; + const fp delta_x = coord_delta[X_AXIS]; #if defined(HAVE_OP_DX) \ || defined(HAVE_OP_DXY) || defined(HAVE_OP_DXZ) \ || defined(HAVE_OP_DXX) @@ -700,8 +712,8 @@ int out; #endif #endif #if N_DIMS >= 2 - const fp origin_y = coord_system_origin[Y_AXIS]; - const fp delta_y = grid_spacing[Y_AXIS]; + const fp origin_y = coord_origin[Y_AXIS]; + const fp delta_y = coord_delta[Y_AXIS]; #if defined(HAVE_OP_DY) \ || defined(HAVE_OP_DXY) || defined(HAVE_OP_DYZ) \ || defined(HAVE_OP_DYY) @@ -709,8 +721,8 @@ int out; #endif #endif #if N_DIMS >= 3 - const fp origin_z = coord_system_origin[Z_AXIS]; - const fp delta_z = grid_spacing[Z_AXIS]; + const fp origin_z = coord_origin[Z_AXIS]; + const fp delta_z = coord_delta[Z_AXIS]; #if defined(HAVE_OP_DZ) \ || defined(HAVE_OP_DXZ) || defined(HAVE_OP_DYZ) \ || defined(HAVE_OP_DZZ) @@ -764,6 +776,7 @@ int out; /* * interpolate at each point */ +int return_code = 0; int pt; for (pt = 0 ; pt < N_interp_points ; ++pt) { @@ -890,124 +903,84 @@ int pt; * * n.b. we need the final answers in variables with the magic * names (x,y,z) (machine-generated code uses these names), - * but we use temp variables as intermediates for (likely) - * better performance: the temp variables have their addresses - * taken and so may not be register-allocated, whereas the - * final (x,y,z) are "clean" and thus more likely to be - * register-allocated + * but we use temp variables as intermediates for these and + * for center_[ijk] for (likely) better performance: + * the temp variables have their addresses taken and so + * may not be register-allocated, whereas the final variables + * are "clean" and thus more likely to be register-allocated */ { - #if (N_DIMS >= 1) - fp x_temp; - const int center_i - = LocalInterp_molecule_posn(origin_x, delta_x, - input_array_min_subscripts[X_AXIS], - input_array_max_subscripts[X_AXIS], - MOLECULE_SIZE, - out_of_range_tolerance[X_AXIS_MIN], - out_of_range_tolerance[X_AXIS_MAX], - interp_coords_fp[X_AXIS], - &x_temp, - (int *) NULL, (int *) NULL); - const fp x = x_temp; - #endif - #if (N_DIMS >= 2) - fp y_temp; - const int center_j - = LocalInterp_molecule_posn(origin_y, delta_y, - input_array_min_subscripts[Y_AXIS], - input_array_max_subscripts[Y_AXIS], - MOLECULE_SIZE, - out_of_range_tolerance[Y_AXIS_MIN], - out_of_range_tolerance[Y_AXIS_MAX], - interp_coords_fp[Y_AXIS], - &y_temp, - (int *) NULL, (int *) NULL); - const fp y = y_temp; - #endif - #if (N_DIMS >= 3) - fp z_temp; - const int center_k - = LocalInterp_molecule_posn(origin_z, delta_z, - input_array_min_subscripts[Z_AXIS], - input_array_max_subscripts[Z_AXIS], - MOLECULE_SIZE, - out_of_range_tolerance[Z_AXIS_MIN], - out_of_range_tolerance[Z_AXIS_MAX], - interp_coords_fp[Z_AXIS], - &z_temp, - (int *) NULL, (int *) NULL); - const fp z = z_temp; - #endif - #if (N_DIMS > 3) - #error "N_DIMS may not be > 3!" - #endif - - /* is the interpolation point out-of-range? */ - #if (N_DIMS >= 1) - if (IS_MOLECULE_POSN_ERROR_CODE(center_i)) - then { - if (center_i == MOLECULE_POSN_ERROR_GRID_TINY) - then return CCTK_ERROR_INTERP_GRID_TOO_TINY; - if (error_flags != NULL) - then { - error_flags->error_pt = pt; - error_flags->error_axis = X_AXIS; - error_flags->error_end - = (center_i == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1; - } - return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ - } - #endif - #if (N_DIMS >= 2) - if (IS_MOLECULE_POSN_ERROR_CODE(center_j)) - then { - if (center_j == MOLECULE_POSN_ERROR_GRID_TINY) - then return CCTK_ERROR_INTERP_GRID_TOO_TINY; - if (error_flags != NULL) - then { - error_flags->error_pt = pt; - error_flags->error_axis = Y_AXIS; - error_flags->error_end - = (center_j == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1; - } - return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ - } - #endif - #if (N_DIMS >= 3) - if (IS_MOLECULE_POSN_ERROR_CODE(center_k)) - then { - if (center_k == MOLECULE_POSN_ERROR_GRID_TINY) - then return CCTK_ERROR_INTERP_GRID_TOO_TINY; - if (error_flags != NULL) - then { - error_flags->error_pt = pt; - error_flags->error_axis = Z_AXIS; - error_flags->error_end - = (center_k == MOLECULE_POSN_ERROR_X_GT_MAX) ? +1 : -1; - } - return CCTK_ERROR_INTERP_POINT_X_RANGE; /*** ERROR RETURN ***/ - } - #endif + int center_ijk_temp[MAX_N_DIMS]; + fp xyz_temp[MAX_N_DIMS]; + int axis; + for (axis = 0 ; axis < N_DIMS ; ++axis) + { + const int ibndry_min = 2*axis; + const int ibndry_max = 2*axis + 1; + const int status = LocalInterp_molecule_posn + (coord_origin[axis], coord_delta[axis], + input_array_min_subscripts[axis] + + N_boundary_points_to_omit[ibndry_min], + input_array_max_subscripts[axis] + - N_boundary_points_to_omit[ibndry_max], + MOLECULE_SIZE, + boundary_off_centering_tolerance[ibndry_min], + boundary_off_centering_tolerance[ibndry_max], + boundary_extrapolation_tolerance[ibndry_min], + boundary_extrapolation_tolerance[ibndry_max], + interp_coords_fp[axis], + ¢er_ijk_temp[axis], &xyz_temp[axis]); + if (status < 0) + then { + if (status == MOLECULE_POSN_ERROR_GRID_TINY) + then return CCTK_ERROR_INTERP_GRID_TOO_TINY; + /*** ERROR RETURN ***/ + if (error_flags != NULL) + then { + error_flags->error_pt = pt; + error_flags->error_axis = axis; + if (status == MOLECULE_POSN_ERROR_X_LT_MIN) + then { + error_flags->error_ibndry = ibndry_min; + error_flags->error_end = -1; + } + else { + error_flags->error_ibndry = ibndry_max; + error_flags->error_end = +1; + } + } + return CCTK_ERROR_INTERP_POINT_OUTSIDE; /*** ERROR RETURN ***/ + } + } + { + #if (N_DIMS >= 1) + const int center_i = center_ijk_temp[X_AXIS]; + const fp x = xyz_temp [X_AXIS]; + #endif + #if (N_DIMS >= 2) + const int center_j = center_ijk_temp[Y_AXIS]; + const fp y = xyz_temp [Y_AXIS]; + #endif + #if (N_DIMS >= 3) + const int center_k = center_ijk_temp[Z_AXIS]; + const fp z = xyz_temp [Z_AXIS]; + #endif /* * compute 1-d position of molecule center in input arrays */ { - #if (N_DIMS == 1) + #if (N_DIMS == 1) const int molecule_center_posn = stride_i*center_i; - #endif - #if (N_DIMS == 2) + #elif (N_DIMS == 2) const int molecule_center_posn = stride_i*center_i + stride_j*center_j; - #endif - #if (N_DIMS == 3) + #elif (N_DIMS == 3) const int molecule_center_posn = stride_i*center_i + stride_j*center_j + stride_k*center_k; - #endif - #if (N_DIMS > 3) + #else #error "N_DIMS may not be > 3!" #endif @@ -1566,6 +1539,7 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, } } } + } /* end of for (pt = ...) loop */ } diff --git a/src/GeneralizedPolynomial-Uniform/template.h b/src/GeneralizedPolynomial-Uniform/template.h index 3bc2a13..b5eb5b4 100644 --- a/src/GeneralizedPolynomial-Uniform/template.h +++ b/src/GeneralizedPolynomial-Uniform/template.h @@ -8,13 +8,15 @@ @version $Header$ @@*/ int FUNCTION_NAME(/***** coordinate system *****/ - const CCTK_REAL coord_system_origin[], - const CCTK_REAL grid_spacing[], + const CCTK_REAL coord_origin[], + const CCTK_REAL coord_delta[], /***** interpolation points *****/ int N_interp_points, int interp_coords_type_code, const void* const interp_coords[], - const CCTK_REAL out_of_range_tolerance[], + const CCTK_INT N_boundary_points_to_omit[], + const CCTK_REAL boundary_off_centering_tolerance[], + const CCTK_REAL boundary_extrapolation_tolerance[], /***** input arrays *****/ int N_input_arrays, const CCTK_INT input_array_offsets[], diff --git a/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c b/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c index bd3a1f9..4f6a5f5 100644 --- a/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c +++ b/src/GeneralizedPolynomial-Uniform/test_molecule_posn.c @@ -7,11 +7,27 @@ * Usage: * test_molecule_posn # run a preset set of tests * test_molecule_posn molecule_size \ - * out_of_range_tolerance_min \ - * out_of_range_tolerance_max \ + * boundary_off_centering_tolerance_min \ + * boundary_off_centering_tolerance_max \ + * boundary_extrapolation_tolerance_min \ + * boundary_extrapolation_tolerance_max \ + * boundary_off_centering_tolerance_max \ * x # do a single test as specified */ +const char* help_msg = +"usage:\n" +"# run a preset series of tests:\n" +" test_molecule_posn\n" +"# run a single test as specified:\n" +" test_molecule_posn molecule_size \\\n" +" boundary_off_centering_tolerance_min \\\n" +" boundary_off_centering_tolerance_max \\\n" +" boundary_extrapolation_tolerance_min \\\n" +" boundary_extrapolation_tolerance_max \\\n" +" x\n" +; + #include <math.h> #include <string.h> #include <limits.h> @@ -39,8 +55,10 @@ const int grid_i_max = 105; /* x_max = 13.6 */ */ static void run_interactive_test(int molecule_size, - fp out_of_range_tolerance_min, - fp out_of_range_tolerance_max, + fp boundary_off_centering_tolerance_min, + fp boundary_off_centering_tolerance_max, + fp boundary_extrapolation_tolerance_min, + fp boundary_extrapolation_tolerance_max, fp x); static int run_batch_tests(void); @@ -56,152 +74,135 @@ struct args_results { /* args */ int molecule_size; - fp out_of_range_tolerance_min; - fp out_of_range_tolerance_max; + fp boundary_off_centering_tolerance_min; + fp boundary_off_centering_tolerance_max; + fp boundary_extrapolation_tolerance_min; + fp boundary_extrapolation_tolerance_max; fp x; /* results */ + int status; int i_center; fp x_rel; - int m_min, m_max; }; /* test data */ +#define OK 0 +#define X_LT_MIN MOLECULE_POSN_ERROR_X_LT_MIN +#define X_GT_MAX MOLECULE_POSN_ERROR_X_GT_MAX static const struct args_results test_data[] = { - /* molecule size 2 */ - { 2, 1.0, 1.0e-12, 7.19, INT_MIN, -1.1, 0,+1 }, - { 2, 1.0, 1.0e-12, 7.21, 42, -0.9, 0,+1 }, - { 2, 1.0, 1.0e-12, 7.24, 42, -0.6, 0,+1 }, - { 2, 1.0, 1.0e-12, 7.26, 42, -0.4, 0,+1 }, - { 2, 1.0, 1.0e-12, 7.29, 42, -0.1, 0,+1 }, - { 2, -1.0, 1.0e-12, 7.29, INT_MIN, -0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 7.3, 42, 0.0, 0,+1 },/* grid_x_min */ - { 2, 1.0e-12, 1.0e-12, 7.31, 42, +0.1, 0,+1 }, - { 2, -1.0 , 1.0e-12, 7.31, 42, +0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 7.35, 42, +0.5, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 7.39, 42, +0.9, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 7.41, 43, 0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 9.81, 67, +0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 9.85, 67, +0.5, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 9.89, 67, +0.9, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.45, 103, +0.5, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.51, 104, +0.1, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.55, 104, +0.5, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.59, 104, +0.9, 0,+1 }, - { 2, 1.0e-12, -1.0, 13.59, 104, +0.9, 0,+1 }, - { 2, 1.0e-12, 1.0e-12, 13.6, 104, +1.0, 0,+1 },/* grid_x_max */ - { 2, 1.0e-12, 1.0, 13.61, 104, +1.1, 0,+1 }, - { 2, 1.0e-12, -1.0, 13.61, INT_MAX, +1.1, 0,+1 }, - { 2, 1.0e-12, 1.0, 13.65, 104, +1.5, 0,+1 }, - { 2, 1.0e-12, 1.0, 13.69, 104, +1.9, 0,+1 }, - { 2, 1.0e-12, 1.0, 13.71, INT_MAX, +2.1, 0,+1 }, - - /* molecule size 3 */ - { 3, 1.0, 1.0e-12, 7.19, INT_MIN, -2.1, -1,+1 }, - { 3, 1.0, 1.0e-12, 7.21, 43, -1.9, -1,+1 }, - { 3, 1.0, 1.0e-12, 7.24, 43, -1.6, -1,+1 }, - { 3, 1.0, 1.0e-12, 7.26, 43, -1.4, -1,+1 }, - { 3, 1.0, 1.0e-12, 7.29, 43, -1.1, -1,+1 }, - { 3, -1.0, 1.0e-12, 7.29, INT_MIN, -1.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.3, 43, -1.0, -1,+1 },/* grid_x_min */ - { 3, 1.0e-12, 1.0e-12, 7.31, 43, -0.9, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.34, 43, -0.6, -1,+1 }, - { 3, -1.0, 1.0e-12, 7.34, INT_MIN, -0.6, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.36, 43, -0.4, -1,+1 }, - { 3, -1.0, 1.0e-12, 7.36, 43, -0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.39, 43, -0.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.4, 43, 0.0, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.44, 43, +0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 7.46, 44, -0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.8, 67, 0.0, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.81, 67, +0.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.84, 67, +0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.86, 68, -0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.89, 68, -0.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 9.9, 68, 0.0, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.44, 103, +0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.46, 104, -0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.5, 104, 0.0, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.51, 104, +0.1, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.54, 104, +0.4, -1,+1 }, - { 3, 1.0e-12, -1.0, 13.54, 104, +0.4, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.56, 104, +0.6, -1,+1 }, - { 3, 1.0e-12, -1.0, 13.56, INT_MAX, +0.6, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.59, 104, +0.9, -1,+1 }, - { 3, 1.0e-12, 1.0e-12, 13.6, 104, +1.0, -1,+1 },/* grid_x_max */ - { 3, 1.0e-12, 1.0, 13.61, 104, +1.1, -1,+1 }, - { 3, 1.0e-12, 1.0, 13.65, 104, +1.5, -1,+1 }, - { 3, 1.0e-12, 1.0, 13.69, 104, +1.9, -1,+1 }, - { 3, 1.0e-12, 1.0, 13.71, INT_MAX, +2.1, -1,+1 }, - - /* molecule size 4 */ - { 4, 0.2, 1.0e-12, 7.27, INT_MIN, -1.3, -1,+2 }, - { 4, 0.2, 1.0e-12, 7.29, 43, -1.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.3, 43, -1.0, -1,+2 },/* grid_x_min */ - { 4, 1.0e-12, 1.0e-12, 7.33, 43, -0.7, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.39, 43, -0.1, -1,+2 }, - { 4, -1.0, 1.0e-12, 7.39, INT_MIN, -0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.4, 43, 0.0, -1,+2 }, - { 4, -1.0, 1.0e-12, 7.41, 43, +0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.42, 43, +0.2, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.48, 43, +0.8, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 7.51, 44, +0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 9.81, 67, +0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 9.85, 67, +0.5, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 9.89, 67, +0.9, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.39, 102, +0.9, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.41, 103, +0.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.48, 103, +0.8, -1,+2 }, - { 4, 1.0e-12, -1.0, 13.48, 103, +0.8, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.5, 103, +1.0, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.51, 103, +1.1, -1,+2 }, - { 4, 1.0e-12, -1.0, 13.51, INT_MAX, +1.1, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.55, 103, +1.5, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.59, 103, +1.9, -1,+2 }, - { 4, 1.0e-12, 1.0e-12, 13.6, 103, +2.0, -1,+2 },/* grid_x_max */ - { 4, 1.0e-12, 2.0, 13.79, 103, +3.9, -1,+2 }, - { 4, 1.0e-12, 2.0, 13.81, INT_MAX, +4.1, -1,+2 }, - - /* molecule size 5 */ - { 5, 3.0, 1.0e-12, 6.99, INT_MIN, -5.1, -2,+2 }, - { 5, 3.0, 1.0e-12, 7.01, 44, -4.9, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.3, 44, -2.0, -2,+2 },/* grid_x_min */ - { 5, 1.0e-12, 1.0e-12, 7.4, 44, -1.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.44, 44, -0.6, -2,+2 }, - { 5, -1.0, 1.0e-12, 7.44, INT_MIN, -0.6, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.46, 44, -0.4, -2,+2 }, - { 5, -1.0, 1.0e-12, 7.46, 44, -0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.49, 44, -0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.5, 44, 0.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.54, 44, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.56, 45, -0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.6, 45, 0.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 7.64, 45, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.8, 67, 0.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.81, 67, +0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.84, 67, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.86, 68, -0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.89, 68, -0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 9.9, 68, 0.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.34, 102, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.36, 103, -0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.39, 103, -0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.41, 103, +0.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.44, 103, +0.4, -2,+2 }, - { 5, 1.0e-12, -1.0, 13.44, 103, +0.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.46, 103, +0.6, -2,+2 }, - { 5, 1.0e-12, -1.0, 13.46, INT_MAX, +0.6, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.48, 103, +0.8, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.5, 103, +1.0, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.51, 103, +1.1, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.54, 103, +1.4, -2,+2 }, - { 5, 1.0e-12, -1.0, 13.54, INT_MAX, +1.4, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.56, 103, +1.6, -2,+2 }, - { 5, 1.0e-12, -1.0, 13.56, INT_MAX, +1.6, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.59, 103, +1.9, -2,+2 }, - { 5, 1.0e-12, 1.0e-12, 13.6, 103, +2.0, -2,+2 },/* grid_x_max */ - { 5, 1.0e-12, 1.5, 13.74, 103, +3.4, -2,+2 }, - { 5, 1.0e-12, 1.5, 13.76, INT_MAX, +3.6, -2,+2 }, + /*** molecule size 2 ***/ + /* status */ + /* off_centering extrapolation i_center */ + /* min max min max x x_rel */ + { 2, 999.0, 999.0, 1.0, 999.0, 7.19, X_LT_MIN, 0, 0.0 }, + { 2, 999.0, 999.0, 1.0, 999.0, 7.21, OK, 42, -0.9 }, + { 2, 999.0, 999.0, 1.0, 999.0, 7.24, OK, 42, -0.6 }, + { 2, 999.0, 999.0, 1.0, 999.0, 7.26, OK, 42, -0.4 }, + { 2, 0.0, 0.0, 1.0, 999.0, 7.26, X_LT_MIN, 0, 0.0 }, + { 2, 1.0e-6, 999.0, 0.2, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 2, 0.2, 999.0, 1.0e-6, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 2, 0.2, 999.0, 0.2, 999.0, 7.29, OK, 42, -0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.31, OK, 42, +0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.34, OK, 42, +0.4 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.36, OK, 42, +0.6 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.39, OK, 42, +0.9 }, + { 2, 0.0, 0.0, 1.0e-6, 999.0, 7.41, OK, 43, +0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.81, OK, 67, +0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.84, OK, 67, +0.4 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.85, OK, 67, +0.5 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.86, OK, 67, +0.6 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 9.89, OK, 67, +0.9 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 13.45, OK, 103, +0.5 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 13.51, OK, 104, +0.1 }, + { 2, 0.0, 0.0, 1.0e-6, 1.0e-6, 13.59, OK, 104, +0.9 }, + { 2, 0.0, 0.2, 1.0e-6, 1.0e-6, 13.61, X_GT_MAX, 0, 0.0 }, + { 2, 0.0, 1.0e-6, 1.0e-6, 1.0e-6, 13.61, X_GT_MAX, 0, 0.0 }, + { 2, 0.0, 0.2, 0.0, 0.2, 13.61, OK, 104, +1.1 }, + + /*** molecule size 3 ***/ + /* status */ + /* off_centering extrapolation i_center */ + /* min max min max x x_rel */ + { 3, 0.5, 999.0, 999.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 3, 999.0, 999.0, 0.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 3, 0.7, 999.0, 0.2, 999.0, 7.29, OK, 43, -1.1 }, + { 3, 0.3, 999.0, 0.0, 999.0, 7.31, X_LT_MIN, 0, 0.0 }, + { 3, 0.5, 999.0, 0.0, 999.0, 7.31, OK, 43, -0.9 }, + { 3, 0.0, 0.0, 0.0, 0.0, 7.36, OK, 43, -0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 7.44, OK, 43, +0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 7.46, OK, 44, -0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 9.81, OK, 67, +0.1 }, + { 3, 0.0, 0.0, 0.0, 0.0, 9.84, OK, 67, +0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 9.86, OK, 68, -0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 9.89, OK, 68, -0.1 }, + { 3, 0.0, 0.0, 0.0, 0.0, 13.44, OK, 103, +0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 13.46, OK, 104, -0.4 }, + { 3, 0.0, 0.0, 0.0, 0.0, 13.54, OK, 104, +0.4 }, + { 3, 0.0, 0.2, 0.0, 0.0, 13.56, OK, 104, +0.6 }, + { 3, 999.0, 0.0, 999.0, 999.0, 13.56, X_GT_MAX, 0, 0.0 }, + { 3, 0.0, 0.5, 0.0, 0.0, 13.59, OK, 104, +0.9 }, + { 3, 999.0, 0.3, 999.0, 999.0, 13.59, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 0.5, 999.0, 999.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 999.0, 999.0, 0.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 0.7, 999.0, 0.2, 13.61, OK, 104, +1.1 }, + { 3, 999.0, 999.0, 999.0, 0.2, 13.63, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 0.7, 999.0, 999.0, 13.63, X_GT_MAX, 0, 0.0 }, + { 3, 999.0, 0.9, 999.0, 0.4, 13.63, OK, 104, +1.3 }, + + /*** molecule size 4 ***/ + /* status */ + /* off_centering extrapolation i_center */ + /* min max min max x x_rel */ + { 4, 1.0, 999.0, 999.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 4, 999.0, 999.0, 0.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 4, 1.2, 999.0, 0.2, 999.0, 7.29, OK, 43, -1.1 }, + { 4, 0.5, 999.0, 0.0, 999.0, 7.31, X_LT_MIN, 0, 0.0 }, + { 4, 1.0, 999.0, 0.0, 999.0, 7.31, OK, 43, -0.9 }, + { 4, 0.0, 999.0, 0.0, 999.0, 7.39, X_LT_MIN, 0, 0.0 }, + { 4, 0.2, 999.0, 0.0, 999.0, 7.39, OK, 43, -0.1 }, + { 4, 0.0, 0.0, 0.0, 999.0, 7.41, OK, 43, +0.1 }, + { 4, 0.0, 0.0, 0.0, 999.0, 7.49, OK, 43, +0.9 }, + { 4, 0.0, 0.0, 0.0, 999.0, 7.51, OK, 44, +0.1 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.81, OK, 67, +0.1 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.84, OK, 67, +0.4 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.85, OK, 67, +0.5 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.86, OK, 67, +0.6 }, + { 4, 0.0, 0.0, 0.0, 0.0, 9.89, OK, 67, +0.9 }, + { 4, 0.0, 0.0, 0.0, 0.0, 13.44, OK, 103, +0.4 }, + { 4, 0.0, 0.0, 0.0, 0.0, 13.49, OK, 103, +0.9 }, + { 4, 0.0, 0.2, 0.0, 0.0, 13.51, OK, 103, +1.1 }, + { 4, 999.0, 0.0, 999.0, 999.0, 13.51, X_GT_MAX, 0, 0.0 }, + { 4, 999.0, 0.8, 999.0, 0.0, 13.59, X_GT_MAX, 0, 0.0 }, + { 4, 0.0, 1.0, 0.0, 0.0, 13.59, OK, 103, +1.9 }, + { 4, 999.0, 999.0, 999.0, 0.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 4, 0.0, 1.0, 999.0, 999.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 4, 0.0, 1.2, 0.0, 0.2, 13.61, OK, 103, +2.1 }, + + /*** molecule size 5 ***/ + /* status */ + /* off_centering extrapolation i_center */ + /* min max min max x x_rel */ + { 5, 1.5, 999.0, 999.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 5, 999.0, 999.0, 0.0, 999.0, 7.29, X_LT_MIN, 0, 0.0 }, + { 5, 1.7, 999.0, 0.2, 999.0, 7.29, OK, 44, -2.1 }, + { 5, 0.9, 999.0, 0.0, 999.0, 7.35, X_LT_MIN, 0, 0.0 }, + { 5, 1.1, 999.0, 0.0, 999.0, 7.35, OK, 44, -1.5 }, + { 5, 0.0, 999.0, 999.0, 999.0, 7.44, X_LT_MIN, 0, 0.0 }, + { 5, 0.2, 999.0, 0.0, 999.0, 7.44, OK, 44, -0.6 }, + { 5, 0.0, 0.0, 0.0, 0.0, 7.46, OK, 44, -0.4 }, + { 5, 0.0, 0.0, 0.0, 0.0, 9.81, OK, 67, +0.1 }, + { 5, 0.0, 0.0, 0.0, 0.0, 9.84, OK, 67, +0.4 }, + { 5, 0.0, 0.0, 0.0, 0.0, 9.86, OK, 68, -0.4 }, + { 5, 0.0, 0.0, 0.0, 0.0, 9.89, OK, 68, -0.1 }, + { 5, 0.0, 0.0, 0.0, 0.0, 13.44, OK, 103, +0.4 }, + { 5, 999.0, 0.0, 999.0, 999.0, 13.46, X_GT_MAX, 0, 0.0 }, + { 5, 0.0, 0.2, 0.0, 0.0, 13.46, OK, 103, +0.6 }, + { 5, 0.0, 1.1, 0.0, 0.0, 13.55, OK, 103, +1.5 }, + { 5, 999.0, 0.9, 999.0, 999.0, 13.55, X_GT_MAX, 0, 0.0 }, + { 5, 999.0, 999.0, 999.0, 0.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 5, 999.0, 1.5, 999.0, 999.0, 13.61, X_GT_MAX, 0, 0.0 }, + { 5, 999.0, 1.7, 999.0, 0.2, 13.61, OK, 103, +2.1 }, }; #define N_TESTS ((int) (sizeof(test_data)/sizeof(test_data[0]))) @@ -212,7 +213,10 @@ int main(int argc, const char *const argv[]) { bool N_fail; int molecule_size; -double out_of_range_tolerance_min, out_of_range_tolerance_max; +fp boundary_off_centering_tolerance_min; +fp boundary_off_centering_tolerance_max; +fp boundary_extrapolation_tolerance_min; +fp boundary_extrapolation_tolerance_max; double x; switch (argc) @@ -222,44 +226,34 @@ case 1: N_fail = run_batch_tests(); if (N_fail == 0) then { - printf("*** all tests successful ***\n"); + printf("*** all %d tests successful ***\n", N_TESTS); return 0; } else { - printf("*** %d test(s) failed ***\n", N_fail); + printf("*** %d/%d test(s) failed ***\n", N_fail, N_TESTS); return 1; } -case 5: +case 7: if ( (sscanf(argv[1], "%d", &molecule_size) == 1) - && (sscanf(argv[2], "%lf", &out_of_range_tolerance_min) == 1) - && (sscanf(argv[3], "%lf", &out_of_range_tolerance_max) == 1) - && (sscanf(argv[4], "%lf", &x) == 1) ) + && (sscanf(argv[2], "%lf", &boundary_off_centering_tolerance_min) == 1) + && (sscanf(argv[3], "%lf", &boundary_off_centering_tolerance_max) == 1) + && (sscanf(argv[4], "%lf", &boundary_extrapolation_tolerance_min) == 1) + && (sscanf(argv[5], "%lf", &boundary_extrapolation_tolerance_max) == 1) + && (sscanf(argv[6], "%lf", &x) == 1) ) then { run_interactive_test(molecule_size, - out_of_range_tolerance_min, - out_of_range_tolerance_max, + boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max, + boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max, x); - return 0; + return 0; /*** NORMAL EXIT ***/ } - /* fall through */ + /* error ==> fall through */ default: - fprintf(stderr, - "usage:\n" - "# run a single test as specified:\n" - " %s molecule_size \\\n" - " %*s out_of_range_tolerance_min \\\n" - " %*s out_of_range_tolerance_max \\\n" - " %*s x\n" - "# run a preset set of tests:\n" - " %s\n" - , - argv[0], - (int) strlen(argv[0]), "", - (int) strlen(argv[0]), "", - (int) strlen(argv[0]), "", - argv[0]); + fprintf(stderr, help_msg); return 1; } } @@ -271,33 +265,39 @@ default: */ static void run_interactive_test(int molecule_size, - fp out_of_range_tolerance_min, - fp out_of_range_tolerance_max, + fp boundary_off_centering_tolerance_min, + fp boundary_off_centering_tolerance_max, + fp boundary_extrapolation_tolerance_min, + fp boundary_extrapolation_tolerance_max, fp x) { +int i_center; fp x_rel; -int m_min, m_max; printf("testing with molecule_size=%d\n", molecule_size); -printf(" out_of_range_tolerance_[min,max]=[%g,%g]\n", - (double) out_of_range_tolerance_min, - (double) out_of_range_tolerance_max); +printf(" boundary_off_centering_tolerance_[min,max]=[%g,%g]\n", + (double) boundary_off_centering_tolerance_min, + (double) boundary_off_centering_tolerance_max); +printf(" boundary_extrapolation_tolerance_[min,max]=[%g,%g]\n", + (double) boundary_extrapolation_tolerance_min, + (double) boundary_extrapolation_tolerance_max); printf(" x=%g\n", (double) x); { -const int i_center = LocalInterp_molecule_posn(grid_x0, grid_dx, - grid_i_min, grid_i_max, - molecule_size, - out_of_range_tolerance_min, - out_of_range_tolerance_max, - x, - &x_rel, - &m_min, &m_max); - -if ((i_center == INT_MIN) || (i_center == INT_MAX)) - then printf("i_center=%d\n", i_center); - else printf("i_center=%d x_rel=%g m_[min,max]=[%d,%d] i_[min,max]=[%d,%d]\n", - i_center, x_rel, m_min, m_max, i_center+m_min, i_center+m_max); +const int status + = LocalInterp_molecule_posn(grid_x0, grid_dx, + grid_i_min, grid_i_max, + molecule_size, + boundary_off_centering_tolerance_min, + boundary_off_centering_tolerance_max, + boundary_extrapolation_tolerance_min, + boundary_extrapolation_tolerance_max, + x, + &i_center, &x_rel); + +if (status < 0) + then printf("status=%d\n", status); + else printf("i_center=%d x_rel=%g\n", i_center, x_rel); } } @@ -319,41 +319,48 @@ int failure_count = 0; for (i = 0 ; i < N_TESTS ; ++i) { const struct args_results *p = & test_data[i]; - fp x_rel = 0.0; - int m_min = 0, m_max = 0; - int i_center = LocalInterp_molecule_posn(grid_x0, grid_dx, - grid_i_min, grid_i_max, - p->molecule_size, - p->out_of_range_tolerance_min, - p->out_of_range_tolerance_max, - p->x, - &x_rel, - &m_min, &m_max); - bool ok = ((i_center == INT_MIN) || (i_center == INT_MAX)) - ? (i_center == p->i_center) - : ( (i_center == p->i_center) - && fuzzy_EQ(x_rel, p->x_rel) - && (m_min == p->m_min) - && (m_max == p->m_max) ); - - int msglen = - printf("size=%d tol=[%g,%g] x=%g ==> ", - p->molecule_size, - (double) p->out_of_range_tolerance_min, - (double) p->out_of_range_tolerance_max, - (double) p->x); - printf("i_center=%d x_rel=%g m=[%d,%d]\n", - i_center, (double) x_rel, m_min, m_max); - - if (! ok) - then { + int i_center = 999; + fp x_rel = 999.999; + const int status = LocalInterp_molecule_posn + (grid_x0, grid_dx, + grid_i_min, grid_i_max, + p->molecule_size, + p->boundary_off_centering_tolerance_min, + p->boundary_off_centering_tolerance_max, + p->boundary_extrapolation_tolerance_min, + p->boundary_extrapolation_tolerance_max, + p->x, + &i_center, &x_rel); + bool ok = (status == 0) + ? ( (status == p->status) && (i_center == p->i_center) + && fuzzy_EQ(x_rel, p->x_rel) ) + : (status == p->status); + + printf("size=%d off_centering_tol=[%g,%g] extrapolation_tol=[%g,%g] x=%g", + p->molecule_size, + (double) p->boundary_off_centering_tolerance_min, + (double) p->boundary_off_centering_tolerance_max, + (double) p->boundary_extrapolation_tolerance_min, + (double) p->boundary_extrapolation_tolerance_max, + (double) p->x); + printf(": status=%d", status); + if (status == 0) + then printf(" i_center=%d x_rel=%g: ", + i_center, (double) x_rel); + + if (ok) + then printf(": ok\n"); + else { ++failure_count; - { - int error_msglen = printf("***** FAIL: "); - printf("%-*s", msglen-error_msglen, "expected"); - printf("i_center=%d x_rel=%g m_[min,max]=[%d,%d]\n", - p->i_center, (double) p->x_rel, p->m_min, p->m_max); - } + printf(": FAIL\n"); + if (p->status == 0) + then printf("*** expected status=%d i_center=%d x_rel=%g\n", + p->status, p->i_center, (double) p->x_rel); + else printf("*** expected status=%d\n", p->status); + if (status == 0) + then printf("*** got status=%d i_center=%d x_rel=%g\n", + status, i_center, (double) x_rel); + else printf("*** got status=%d\n", status); } } |