xx Thorn Guide for the AEILocalInterp Thorn

Thorn Guide for the AEILocalInterp Thorn

Jonathan Thornburg, incorporating code from Thomas Radke

Date

Abstract

This thorn provides the CCTK_InterpLocalUniform() API for processor-local interpolation of N-dimensional data arrays. The data arrays (in general there may be many of them) must be defined on a uniformly spaced grid.

This thorn provides several variants of Lagrange and Hermite interpolation. At present the Lagrange interpolation operators support ordersx1 throughx6 for 1-D interpolation, andx1 throughx4 for 2-D and 3-D interpolation. The Hermite interpolation operators support ordersx2, 3, andx4 for 1-D interpolation, andx2 andx3 for 2-D and 3-D interpolation.

This thorn supports a number of interpolation options, including parameter-table entries to control the handling of grid boundaries and out-of-range points, non-contiguous input arrays, taking derivatives as part of the interpolation, and querying the Jacobian of the interpolation operation.

1 Introduction

This interpolator provides:

Lagrange polynomial interpolation of orders 1–6 for 1-D interpolation and 1–4 for 2-D and 3-D interpolation. This interpolator actually provides two slightly different flavors of Lagrange interpolation, tensor product and maximum degree, registered under the operator names

xxx"Lagrangexpolynomialxinterpolationx(tensorxproduct)"
xxx"Lagrangexpolynomialxinterpolationx(maximumxdegree)"

Sectionx2.5 of this thorn guide explains the two variants in detail, but for most purposes the tensor-product variant is preferable. For convenience this is also registered under the operator name

xxx"Lagrangexpolynomialxinterpolation"

and this is probably the best interpolation operator for general-purpose

use.1 1For backwards compatability this interpolator is also registered under still another operator name, "generalized polynomial interpolation", but this is deprecated and will probably be removed at some point.

Hermite polynomial interpolation of orders 2-4 for 1-D interpolation, and 2–3 for 2-D and 3-D interpolation, registered under the operator name

     xxx"Hermitexpolynomialxinterpolation"

The code allows arbitrarily-shaped interpolation molecules, but at the moment only hypercube-shaped molecules are implemented. It would be easy to add additional orders and/or dimensions if desired.

This interpolator supports a number of options specified by a parameter table (a Cactus key-value table, manipulated by the 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.

1.1 History

This interpolator was written by Jonathan Thornburg in winter 2001–2002. Between then and Julyx2003 it lived in the LocalInterp thorn in the CactusBase arrangement, but in July 2003 it was moved to the (new) AEILocalInterp thorn in the AEIThorns arrangement so it could stay GPL (Cactus policies forbid GPL code in the CactusBase arrangement).

1.2 Basic Terminology

Within Cactus, each interpolator has a character-string name; this is mapped to a Cactus interpolator handle by CCTK_InterpHandle(). For any given interpolator handle, there may be a separate interpolator defined for each of the interpolation APIs (both the processor-local ones provided by this thorn, and the global ones provided by driver-specific interpolator thorns such as PUGHInterp).

Terminology for interpolation seems to differ a bit from one author to another. Here we refer to the point-centered interpolation of a set of input arrays (defining data on a uniformly or nonuniformly spaced grid of data points) to a set of interpolation points (specified by a corresponding set of interpolation coordinates), with the results being stored in a corresponding set of output arrays. Alternatively, we may refer to cell-centered interpolation, using a grid of data cells and a set of interpolation cells. (This latter terminology is common in hydrodynamics interpolators.)

At present all the interpolators do polynomial interpolation, and allow the interpolation of multiple input arrays (to the same set of interpolation points) in a single interpolator call, using the basic algorithm:

forxeachxinterpolationxpoint
{
choosexanxinterpolationxmoleculexpositionxsomewherexnearxthexinterpolationxpoint
xxxxxxxxforxeachxoutputxarray
xxxxxxxx{
xxxxxxxxcomputexanxinterpolatingxpolynomial
xxxxxxxxxxwhichxapproximatesxthexinputxdataxatxthexmoleculexpoints
xxxxxxxxoutputx=xpolynomial(interpolationxpoint)
xxxxxxxx}
}

In the future, we may add other interpolators where the choice of molecule is data-dependent (and thus may vary from one input array to the next), and/or where the entire input grid is used in each interpolation.

We define the order of the interpolation to be the order of the fitted polynomial. That is, in our terminology linear interpolation is orderx1, quadratic is orderx2, cubic is orderx3, etc. An orderx interpolator thus has interpolation errors for generic smooth input data. Sectionx2.5 explains how the interpolating polynomial is defined for multi-dimensional interpolation.

1.3 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, i.e. its 1stxderivative is generically discontinuous, with 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.2 2For a further discussion of the non-smoothness of interpolation errors, see appendixxF of J. Thornburg, Physical Review D 59(10), 104007. This is the case, for example, for tensor-product Lagrange polynomial interpolation.

[For maximum-degree Lagrange polynomial interpolation, the interpolation error may not be zero at the grid points, and the interpolant itself may not be continuous there. Sectionx2.5 explains this in detail.] This thorn also provides Hermite polynomial interpolation, which guarantees a smooth () interpolant and (for smooth input data) interpolation error. Unfortunately, this comes at the cost of a larger molecule size than the same-order Lagrange interplator, and a much larger interpolation error if the interpolation molecules are off-centered. (By default, the Hermite interpolator doesn’t allow off-centered molecules, though this can be

changed via the boundary_off_centering_tolerance[] and excision_off_centering_tolerance[]

parameter-table entries.)

For most purposes Lagrange polynomial interpolation is better; only use Hermite polynomial interpolation if you need a smooth interpolant.

1.4 More Terminology

As described in the Function Reference section of the Cactus User’s Guide, interp_coords, input_arrays, and output_arrays are actually all pointers to arrays of voidx* pointers, since we support a number of different Cactus data types. Internally, the interpolator casts these voidx* pointers to CCTK_REALx* or whatever the correct Cactus data types are. But when describing how the interpolator accesses the various arrays, for simplicity we usually gloss over this casting, i.e. we pretend that interp_coords, input_arrays, and output_arrays are pointers to arrays of CCTK_REALx* pointers. (This may become clearer once you read the next example.)

We use pt, in, and out as generic 0-origin integer subscripts into the arrays of interpolation points, input arrays, and output arrays respectively. We use (i,j,k) as a generic N_dims-vector of integer subscripts into the input array input_arrays[in]. (Thus (i,j,k) space refers to the grid of data points.) We usually only write array subscripting expressions for the 3-D case; the restrictions/generalizations to other dimensions should be obvious.

For example, for 3-D interpolation, the (x,y,z) coordinates of a typical interpolation point are given by

xx=xinterp_coords[0][pt]
yx=xinterp_coords[1][pt]
zx=xinterp_coords[2][pt]

(Notice here that as described above, we’ve implicitly taken interp_coords to have the Cxtype constxCCTK_REAL*xinterp_coords[], glossing over the casting from its actual Cxtype of

constxvoid*xinterp_coords[].)

We take axis to be a 0-origin integer specifying a coordinate axis (dimension), i.e. 0xforx, 1xforx, 2xforx, x. We take 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.

When describing Jacobians and domains of dependence, it’s useful to introduce the notion of molecule position, a nominal reference point for the interpolation molecule in (i,j,k) coordinate space. (For example, the molecule position might just be the (i,j,k) coordinates of the molecule’s center.) We also introduce molecule coordinates (mi,mj,jk), which are just (i,j,k) coordinates relative to the molecule position. We use m as a generic molecule coordinate. Thus (in notation which should be obvious) a generic molecule operation can be written Note that there is no requirement that the output be semantically located at the position posn! (This may become clearer once you look at the example in sectionx2.9.2.) However, by convention we do assume that is always a valid m coordinate; this simplifies pointer arithmetic.

When describing various entries in the parameter table in, we use const qualifiers on table entries to indicate that the interpolator treats them as const variables/arrays, i.e. it promises not to change them. In contrast, table entries which are not shown as const are considered read-write by the interpolator; it typically uses them to return outputs back to the caller.

2 Interpolator Parameters

The CCTK_InterpLocalUniform() interpolation API (described in detail in the “Function Reference” appendix of the Cactus Users’ Guide) includes a parameter table (a Cactus key-value table, manipulated by the Util_Table*() APIs) whose table handle is passed to the interpolator. This can be used to specify a number of options to the interpolator; it can also be used by the interpolator to return additional information back to the caller.

2.1 Interpolation Order The only mandatory parameter for this interpolator is the interpolation order:

constxCCTK_INTxorder;

As noted in sectionx1.2, in our terminology linear interpolation is orderx1, quadratic is orderx2, cubic is orderx3, etc.

2.2 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. Tablex1 gives the molecule size and details of the centering policy for each interpolation operator and order.

01                                                                                      11     10Molecule xxxxx  00  xxxxx

0101
01

01Interpolation Operator  11Order10 Size  xxxxx  00 Molecule Centering Policy xxxxx
01"Lagrangexpolynomialxinterpolation"                                                    11 1 10 2-point xxxxx  00PICT [figure: \left[\right.; \left)\right.] xxxxx
01"Lagrangexpolynomialxinterpolation"                                                    11 2 10 3-point xxxxx  00PICT [figure: \left[\right.; \left)\right.] xxxxx
01"Lagrangexpolynomialxinterpolation"                                                    11 3 10 4-point xxxxx  00PICT [figure: \left[\right.; \left)\right.] xxxxx
01"Lagrangexpolynomialxinterpolation"                                                    11 4 10 5-point xxxxx  00PICT [figure: \left[\right.; \left)\right.] xxxxx
01"Lagrangexpolynomialxinterpolation"                                                    11 5 10 6-point xxxxx  00PICT [figure: \left[\right.; \left)\right.] xxxxx
01"Lagrangexpolynomialxinterpolation"                                                    11 6 10 7-point xxxxx  00PICT [figure: \left[\right.; \left)\right.] xxxxx
01"Hermitexpolynomialxinterpolation"                                                     11 2 10 4-point xxxxx  00PICT [figure: \left[\right.; \left)\right.] xxxxx
01"Hermitexpolynomialxinterpolation"                                                     11 3 10 6-point xxxxx  00PICT [figure: \left[\right.; \left)\right.] xxxxx
01"Hermitexpolynomialxinterpolation"                                                     11 4 10 6-point xxxxx  00PICT [figure: \left[\right.; \left)\right.] xxxxx

01

Tablex1: This table gives the molecule size and centering for each interpolation operator and order. xshows the input grid points, and xshows 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 for all interpolation points in the range .

2.3 Handling of Grid Boundaries

If an interpolation point is too near a grid boundary and/or an excised region the interpolator can either off-center the interpolation molecule, or refuse to interpolate that point (returning a

CCTK_INTERP_ERROR_POINT_OUTSIDE 3 3For backwards compatability the error code CCTKxINTERPxERRORxPOINTxXxRANGE is also defined; it’s a synonym for CCTKxINTERPxERRORxPOINTxOUTSIDE.

or CCTK_INTERP_ERROR_POINT_EXCISED error code, and possibly also printing a warning message). This behavior is controlled by the (optional) parameter-table entries

constxCCTK_INTxN_boundary_points_to_omit[2*N_dims]
constxCCTK_REALxboundary_off_centering_tolerance[2*N_dims]
constxCCTK_REALxboundary_extrapolation_tolerance[2*N_dims]
constxCCTK_REALxexcision_off_centering_tolerance[2*N_dims]x#xnotximplementedxyet
constxCCTK_REALxexcision_extrapolation_tolerance[2*N_dims]x#xnotximplementedxyet

The elements of these arrays are matched up with the grid boundaries in the order . (For 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 ibndry as a generic 0-origin index into these arrays.

2.3.1 Omitting Boundary Points

For any given grid boundary ibndry, the interpolator doesn’t use input data from the 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:

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.

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.

More generally, this option may be useful for controlling the extent to which the interpolator uses input data from ghost and/or symmetry zones.

2.3.2 Off-Centering Molecules near Grid Boundaries and/or Excision Regions

The (optional) parameter-table entries

constxCCTK_REALxboundary_off_centering_tolerance[2*N_dims]
constxCCTK_REALxboundary_extrapolation_tolerance[2*N_dims]
constxCCTK_REALxexcision_off_centering_tolerance[2*N_dims]x#xnotximplementedxyet
constxCCTK_REALxexcision_extrapolation_tolerance[2*N_dims]x#xnotximplementedxyet

control the local interpolator’s willingness to off-center (change away from the default centering) its interpolation molecules near boundaries and excised points.

To describe how these parameters work, it’s useful to define two regions in the data coordinate space:

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.

The default-centering region

is that region in interpolation-point space where the interpolator can use the default molecule centering (described in tablex1), i.e. where the default-centering molecules require data only from the data-valid region.

Figurex1 shows an example of these regions.

PICT[figure]

Figurex1: This figure shows an example of the valid-data and default-centering regions for a 2-D grid. xmarks 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 boundaryx*xtolerance[ibndry] and excisionx*xtolerance[ibndry] arrays apply to each line

segment of the region boundary:

 10    01 (ibndry = 0)
10    01 (ibndry = 1)
10    01 (ibndry = 2)
10    01 (ibndry = 3)

Near a boundary or excised regions, the interpolator will off-center its molecules (if and) only if the interpolation point is both

within () *_off_centering_tolerance[ibndry] grid spacings of the default-centering region, and

within () *_extrapolation_tolerance[ibndry] grid spacings of the valid-data region.

where we use “*” to denote boundary or excision as appropriate.

There are four cases for these parameters:

and :

No off-centering is allowed: the interpolator reports an error (CCTK_ERROR_INTERP_POINT_OUTSIDE or CCTK_ERROR_INTERP_POINT_EXCISED return code as appropriate) if any interpolation point is outside the default-centering region.

and :

The interpolator allows interpolation points to be up to () *_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.

4 and :

4In practice is a conventional value here. The actual value doesn’t matter so long as it’s larger than any possible molecule radius. The interpolator allows interpolation points to be up to () *_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.

and :

In practice the default-centering region is always a (normally proper) subset of the valid-data region, so this case is nonsensical: the positive value for *_extrapolation_tolerance[ibndry] has no effect because of the setting. The interpolator gives a warning for this case.

boundary_off_centering_tolerance[ibndry]x=x999.0xxx#xeffectivelyxinfinity
boundary_extrapolation_tolerance[ibndry]x=x1.0e-10

If any of these table entries aren’t specified, the defaults are In other words, the interpolation points may be anywhere within the valid-data region or up to grid spacing outside it. (The “fudge factor” helps to avoid spurious errors in case floating-point roundoff moves an interpolation point which was supposed to be just on the boundary of the valid-data region, slightly outside it.)

[In the near future these defaults may be changed: For Lagrange polynomial interpolation the new defaults would be while for Hermite polynomial interpolation they would be

boundary_off_centering_tolerance[ibndry]x=x1.0e-10
boundary_extrapolation_tolerance[ibndry]x=x0.0

This would leave Lagrange interpolation unchanged, while for Hermite interpolation the defaults would forbid any significant off-centering of the interpolation molecules.]

2.3.3 Suppressing Warning Messages about Off-Centering

By default, AEILocalInterp prints a Cactus levelx1 warning message for each interpolation point which causes

it to return a CCTK_INTERP_ERROR_POINT_OUTSIDE or CCTK_INTERP_ERROR_POINT_EXCISED error code. If the

key suppress_warnings is present in the parameter table (it may have any datatype and value), then these messages are not printed.

2.4 Per-Point Status Reporting

By default, an interpolator just returns a single result, either 0 for success or some (negative) error code. If there are multiple interpolation points causing errors, an interpolator should return 0 for success, or the error code for the first point in error (the one with the smallest pt).

However, sometimes you would like to know the status of the interpolation for each point. The following (optional) parameter-table entries may be used to request this, and to query if the interpolator supports this:

/*
x*xToxqueryxwhetherxthexinterpolatorxsupportsxper-pointxstatus,xsetxthis
x*xtoxaxNULLxpointer.xxToxactuallyxrequestxper-pointxstatus,xsetxthisxto
x*xaxnon-NULLxpointerxpointingxtoxaxbufferxofxxN_interp_pointsxxCCTK_INTs
x*xintoxwhichxthexinterpolatorxshouldxstorexthexper-pointxstatus.xxThexstatus
x*xforxpointxxptxxisxdefinedxtoxbexthexresultxwhichxCCTK_InterpLocalUnifom()
x*xwouldxreturnxifxthatxwerexthexonlyxpointxbeingxinterpolated.
x*/
CCTK_POINTERxper_point_status;

If the interpolator supports returning per-point status, and if the key per_point_status is present in the parameter table, then the interpolator will set

CCTK_INTxerror_point_status;

in the parameter table. If , the interpolator will set error_point_status tox0. Otherwise, the interpolator will set error_point_status to the negative of the number of points for which the per-point status is non-zero.

2.5 Multi-Dimensional Interpolation

In multiple dimensions, there are two plausible definitions for the generic interpolating polynomial of a given order (degree) . For convenience I’ll describe these for the 2-D case, but the generalization to any number of dimensions should be obvious:

The generic “tensor product” polynomial of degree is

The generic “maximum degree” polynomial of degree is

Because it has coefficients, the tensor-product polynomial can (and does) pass through all the input data points in a square molecule of size . This implies that the interpolation error vanishes at the input grid points, and that the overall interpolating function is continuous (up to floating-point roundoff errors). However, does have the slightly peculiar property of having terms up to despite being of “degree” . (For example, the “linear” interpolant for would have terms, even though those are formally quadratic in the independent variables and .)

In contrast, the maximum-degree polynomial has only coefficients, so for generic input data (i.e. input data which isn’t actually sampled from a polynomial of the formx) can’t pass through all the input data points in a square molecule of size . The interpolator actually does a least-squares fit of the polymomial to the input data in the molecule, so in general the molecule won’t pass through any of the data points! Moreover, each time the interpolation point crosses a grid square (for oddx) or the center lines of a grid square (for evenx), the set of points used in the molecule changes (cf.xsectionx2.2), so the interpolant generally has a jump discontinuity! For these reasons, the tensor-product choicex is generally preferable.

2.6 Molecule Family

An interpolator may support different families/shapes of interpolation molecules. Hypercube-shaped molecules are the simplest and most common case, but one could also imagine (say) octagon-shaped molecules in 2-D, or some generalization of this in higher numbers of dimensions.

The (optional) parameter

/*xsetxwithxUtil_TableSetString()x*/
constxcharxmolecule_family[];

may be used to set or query the molecule family.

If this key is present in the parameter table, the interpolator sets the molecule family/shape based on the value specified. If this key isn’t present in the parameter table, then the interpolator sets it to the molecule family being used.

At present only hypercubed-shaped molecules are implemented; these are specified by molecule_family set to

"cube".

2.7 Non-Contiguous Input Arrays

Sometimes the input “arrays” used by the interpolator might not be contiguous in memory. For example, we might want to do 2-D interpolation within a plane of a 3-D grid array, but the plane might or might not be contiguous in memory. (Another example would be that the input arrays might be members of a compact group.)

The following (optional) parameter-table entries specify non-contiguous input arrays:

constxCCTK_INTxinput_array_offsets[N_input_arrays];
/*xthexnextx3xtablexentriesxarexsharedxbyxallxinputxarraysx*/
constxCCTK_INTxinput_array_stridesxxxxxxx[N_dims];
constxCCTK_INTxinput_array_min_subscripts[N_dims];
constxCCTK_INTxinput_array_max_subscripts[N_dims];

In general, the interpolator accesses the input using the generic subscripting expression

input_array[in][offsetx+xi*stride_ix+xj*stride_jx+xk*stride_k]

where

offsetx=xinput_array_offsets[in]
(stride_i,stride_j,stride_k)x=xinput_array_strides[]
 and where (i,j,k) run from input_array_min_subscripts[] to input_array_max_subscripts[] inclusive

(n.b.xthis is an inclusive range, i.e. ).

The defaults are that each input array is contiguous in memory, i.e. input_array_offsets[] = 0, stride

determined from input_array_dims[] in the usual Fortran manner, input_array_min_subscripts[] = all 0s,
and input_array_max_subscripts[] = input_array_dims[]-1. If the stride and max subscript are both

specified explicitly, then the explicit input_array_dims[] argument to CCTK_InterpLocalUniform() is ignored.

2.8 Derivatives

If we view the input data as being samples of a smooth function, then instead of estimating values of that function at the interpolation points, the interpolator can instead or additionally estimate values of various derivatives of that function at the interpolation points. (We don’t currently implement it, but one could also imagine interpolating more general molecule operations such as Laplacians.)

The following (optional) parameter-table entries are used to specify this:

constxCCTK_INTxoperand_indices[N_output_arrays];
constxCCTK_INTxoperation_codes[N_output_arrays];

The semantics here are that

output_array[out]x=xop(input_array[operand_indices[out]])

where op is an operator specified by the operation_codes[out] value as described below.

Note that operand_indices[] doesn’t directly name the inputs, but rather gives their (0-origin) subscripts in the list of inputs. This allows for a more efficient implementation in the (common) case where some of the input arrays have many different operations applied to them. (It’s most efficient to group all operations on a given input array together in the operand_indices and operation_codes arrays, as in the example in sectionx5.2.)

Negative operation_codes[out] values are reserved for future use. An operation_codes[out] value which is is taken as a decimal integer encoding a coordinate partial derivative: each decimal digit means to take the coordinate partial derivative along that (1-origin) axis; the order of the digits in a number is ignored. Tablex2 summarizes these resulting derivative codes.

 11operation_codes[out]11What it means

11  0  11interpolate the input array itself (no derivative)
11  1  11interpolate of the input array
11  2  11interpolate of the input array
11  3  11interpolate of the input array
11  12 or 21  11interpolate of the input array
11  13 or 31  11interpolate of the input array
11  23 or 32  11interpolate of the input array
11  11  11interpolate of the input array
11  22  11interpolate of the input array
11  33  11interpolate of the input array

11

Tablex2: This table gives the codes in operationxcodes[out] for each possible 1st or 2ndxderivative in 3-D; for 1-D or 2-D codes referring to nonexistent coordinates are invalid.

At present we do not have any #defines for the operation codes in the Cactus header files. However, you can avoid most of the software-engineering problems of having “magic numbers” for the operation codes, by using the macro

#definexDERIV(op)xxxop

to mark all such operation_codes[] values in your code. There’s an example of this in sectionx5.2.

2.9 Jacobian and Domain of Dependence

The Jacobian of the interpolator is defined as We may want to know the Jacobian itself, and/or “just” the set of (i,j,k) for which this is nonzero (i.e. the Jacobian’s sparsity structure, or equivalently the domain of dependence of the output result, or equivalently the interpolation molecule size and shape) for a given out, in, and

pt.5 5For something like a spline interpolator the Jacobian is generally nonzero for all (i,j,k), but exponentially small for most of them; in this case the Jacobian-query API would probably specify a cutoff to return an approximate Jacobian with reasonable sparsity.

The complexity of doing this depends (strongly!) on the structure of the Jacobian, and in particular on the answers to the following questions:

What molecule family is being used?

Does the interpolation molecule size and/or shape depend on where the interpolation points are in

the grid?6 6We can always make the interpolation molecules be the “same size and shape” by padding them with zero entries to bring them all up to the size of the largest molecule used anywhere in the grid, but this would be very inefficient if many molecules were padded in this way.

If this interpolator supports computing derivatives as described in sectionx2.8, does the interpolation molecule size and/or shape depend on operation_codes[]?

Does the interpolation molecule size and/or shape depend on the actual floating-point values being interpolated? (Examples of this might include ENO (essentially nonoscillatory) and/or TVD (total-variation diminishing) interpolators for hydrodynamics calculations.)

Do the actual floating-point values of the Jacobian depend on the actual floating-point values being interpolated? Equivalently, is the interpolation nonlinear?

Because the different cases differ so much in the amount of information required to describe the Jacobian, it’s hard to define a single API which covers all cases without burdening the simpler cases with excessive complexity. Instead, the interpolator defines a Jacobian-structure–query API to determine which case applies, together with (conceptually) several different APIs for the different cases. (At the moment only a single case is implemented.)

2.9.1 Determining the Jacobian’s structure

The following parameter-table entries may be used to query which of the different Jacobian-structure cases applies:

The parameter molecule_family may may be used to query what molecule family is being used. This is described in detail in sectionx2.6.

If the interpolation molecule size and/or shape vary with the interpolation coordinates, the interpolator sets the parameter

CCTK_INTxMSS_is_fn_of_interp_coords;

to 1. Otherwise (i.e. if the interpolation molecule size and shape are independent of the interpolation coordinates) it should set this parameter to 0.

If the interpolator supports computing derivatives as described in sectionx2.8, and if the interpolation molecule’s size and/or shape varies with operation_codes[], the interpolator sets the parameter

CCTK_INTxMSS_is_fn_of_which_operation;

to 1. Otherwise (i.e. if the interpolator doesn’t support computing derivatives, or if the interpolator does support computing derivatives but the interpolation molecule size and shape are independent of the operation_code[] values), the interpolator sets this parameter to 0. Note that this query tests whether the molecule size and/or shape depend on operation_codes[] in general, independent of whether there are in fact any distinct values (or even any values at all) passed in operation_codes[] in this particular interpolator call. In other words, this is a query about the basic design of the interpolator, not about this particular call.

If the interpolation molecule’s size and/or shape varies with the actual floating-point values of the input arrays, the interpolator sets the parameter

CCTK_INTxMSS_is_fn_of_input_array_values;

to 1. Otherwise (i.e. if the interpolation molecule size and shape are independent of the input array values; this is a necessary, but not sufficient, condition for the interpolation to be linear), the interpolator sets this parameter to 0.

If the actual floating-point values of the Jacobianx (for a given out, in, and pt) depend on the actual floating-point values of the input arrays (i.e. if the interpolation is nonlinear), the interpolator sets the parameter

CCTK_INTxJacobian_is_fn_of_input_array_values;

to 1. Otherwise (i.e. if the interpolation is linear) the interpolator sets this parameter to 0. The current implementation always sets

MSS_is_fn_of_interp_coordsx=x0
MSS_is_fn_of_which_operationx=x0
MSS_is_fn_of_input_array_valuesx=x0
Jacobian_is_fn_of_input_array_valuesx=x0

2.9.2 Fixed-Size Hypercube-Shaped Molecules

The simplest case (and the only one for which we have defined an API at present) is when the molecules are hypercube-shaped and of (typically small) fixed size, independent of the interpolation coordinates and the actual floating-point values in the input arrays (though presumably depending on the interpolation order and on operation_code). In other words, this case applies if (and only if) the Jacobian structure information described in sectionx2.9.1 returns

MSS_is_fn_of_interp_coordsx=x0
MSS_is_fn_of_which_operationx=x0
MSS_is_fn_of_input_array_valuesx=x0
Jacobian_is_fn_of_input_array_valuesx=x0

(These are precisely the values set by the current implementation.) In the rest of this section we describe the query API for this case.

The following parameters may be used to query the molecule size:

CCTK_INTxconstxmolecule_min_m[N_dims];
CCTK_INTxconstxmolecule_max_m[N_dims];

The semantics of these are that if both of these keys are present (the values and datatypes don’t matter), then the interpolator will (re)set the values to give the (inclusive) minimum and maximum mxmolecule coordinates. (Note that either both of these keys should be present, or neither should be present. This simplifies the logic in the interpolator slightly.)

The following parameter may be used to query the molecule positions:

CCTK_INT*xconstxmolecule_positions[N_dims];

The semantics of this is that the caller should set molecule_positions[] to an array of N_dims pointers to (caller-supplied) arrays of N_interp_points CCTK_INTs 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 Jacobianx itself:

CCTK_REAL*xconstxJacobian_pointer[N_output_arrays];
constxCCTK_INTxxxJacobian_offsetx[N_output_arrays];xxx#xoptional,xdefault=allx0
/*xthexnextx3xtablexentriesxarexsharedxbyxallxJacobiansx*/
constxCCTK_INTxJacobian_interp_point_stride;
constxCCTK_INTxJacobian_m_strides[N_dims];
constxCCTK_INTxJacobian_part_stride;xxxxxxxxxxxxxxxxxx#xoptional,xdefault=1
 If Jacobian_pointer is present in the table, then Jacobian_interp_point_stride

and Jacobian_m_strides must also be present. For each out where

Jacobian_pointer[out]x!=xNULL,7 7That is, setting Jacobianxpointer[out] = NULL supresses the Jacobian query. If (as is often the case) the Jacobian is independent of out, then it’s a good idea to avoid unnecessary computations by supressing the queries in this manner for any remaining output arrays. the interpolator would then store the Jacobianx in

Jacobian_pointer[out][offset
xxxxxxxxxxxxxxxxxxxxxx+xpt*Jacobian_interp_point_stride
xxxxxxxxxxxxxxxxxxxxxx+xmi*stride_ix+xmj*stride_jx+xmk*stride_k
xxxxxxxxxxxxxxxxxxxxxx+xpart*Jacobian_part_stride]

where

offsetx=xJacobian_offset[out]
(stride_i,stride_j,stride_k)x=xJacobian_m_strides[]

and where part is 0 for real values and the real parts of complex values, and 1 for the imaginary parts of complex values.

By appropriately setting the various stride parameters, this allows a fairly wide variety of possible storage layouts for the Jacobian.

An example may help to clarify this: Suppose we have a 1-D grid with 11xgrid points, with integer subscriptsx0 throughx10 inclusive, and interpolation coordinates given by coord_originx=x0.0 and coord_deltax=x0.1. Suppose further that we’re doing Lagrange polynomial interpolation, with orderx=x2 and hence (by tablex1) using 3-point molecules. Finally, suppose that we query the Jacobian molecule positions for the N_interp_points=14 interpolation points 0.0, 0.04, 0.06, 0.10, 0.14, 0.16, 0.20, 0.80, 0.84, 0.86, 0.90, 0.94, 0.96, 1.00. Then the queries might return

interp_molecule_min_mx=x-1
interp_molecule_max_mx=x+1
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx/*xinterp_xxxxxxxmoleculexxxxx*/
interp_molecule_positions[0][x0]x=xx1xxxx/*xxx0.00xxxxx[0.0,x0.1,x0.2]x*/
interp_molecule_positions[0][x1]x=xx1xxxx/*xxx0.04xxxxx[0.0,x0.1,x0.2]x*/
interp_molecule_positions[0][x2]x=xx1xxxx/*xxx0.06xxxxx[0.0,x0.1,x0.2]x*/
interp_molecule_positions[0][x3]x=xx1xxxx/*xxx0.10xxxxx[0.0,x0.1,x0.2]x*/
interp_molecule_positions[0][x4]x=xx1xxxx/*xxx0.14xxxxx[0.0,x0.1,x0.2]x*/
interp_molecule_positions[0][x5]x=xx2xxxx/*xxx0.16xxxxx[0.1,x0.2,x0.3]x*/
interp_molecule_positions[0][x6]x=xx2xxxx/*xxx0.20xxxxx[0.1,x0.2,x0.3]x*/
interp_molecule_positions[0][x7]x=xx8xxxx/*xxx0.80xxxxx[0.7,x0.8,x0.9]x*/
interp_molecule_positions[0][x8]x=xx8xxxx/*xxx0.84xxxxx[0.7,x0.8,x0.9]x*/
interp_molecule_positions[0][x9]x=xx9xxxx/*xxx0.86xxxxx[0.8,x0.9,x1.0]x*/
interp_molecule_positions[0][10]x=xx9xxxx/*xxx0.90xxxxx[0.8,x0.9,x1.0]x*/
interp_molecule_positions[0][11]x=xx9xxxx/*xxx0.94xxxxx[0.8,x0.9,x1.0]x*/
interp_molecule_positions[0][12]x=xx9xxxx/*xxx0.96xxxxx[0.8,x0.9,x1.0]x*/
interp_molecule_positions[0][13]x=xx9xxxx/*xxx1.00xxxxx[0.8,x0.9,x1.0]x*/

2.10 Smoothing

The way the generalized polynomial interpolator is implemented it’s easy to also do Savitzky-Golay

smoothing.8 ^{8See sectionx14.8 of Numerical Recipes (2nd edition) for a general discussion of Savitzky-Golay smoothing.} This is best described by way of an example: Suppose we’re doing 1-D cubic interpolation, using (by tablex1) 4-point molecules. In other words, at each interpolation point we use a cubic interpolation polynomial fitted to 4xsurrounding data points. For Savitzky-Golay smoothing, we would instead least-squares fit a cubic polynomial to some larger number of surrounding data points. This combines interpolation with smoothing, so there’s less amplification of noise in the input data in the interpolation outputs.

The optional input

constxCCTK_INTxsmoothing;

specifies how much (how many points) to enlarge the interpolation molecule for this. The default is 0 (no smoothing). 1 would mean to enlarge the molecule by 1 point (e.g. to use a 5-point molecule instead of the usual 4-point one for cubic interpolation). 2 would mean to enlarge by 2 points, (e.g. to use a 6-point molecule for cubic interpolation). Etc etc.

Note that in xdimension, the maximum-degree Lagrange interpolation already uses more data points than the number of free parameters in the interpolating polynomials, i.e. it already does some Savitzky-Golay smoothing. For example, in 2-D a generic cubic polynomial has 10 free parameters , which we least-squares fit to the 16 data points in the molecule.

Savitzky-Golay smoothing is basically free apart from the increase in the molecule size, e.g. a smoothing=2 cubic interpolation has exactly the same cost as any other 6-point–molecule interpolation.

The current implementation has all the framework for smoothing, but only the smoothing=0 case is implemented at the moment.

2.11 Debugging

Setting the optional parameter

constxCCTK_INTxdebug;

to a positive value turns on various debug printing inside the interpolator. This is intended for debugging the interpolator itself; you need to look at the interpolator source code to interpret the output.

2.12 Logging

This thorn has a Boolean parameter log_interp_coords. If this is set to true, this thorn will write a detailed log file giving the grid information and interpolation coordinates for each interpolator call. Each processor will write a separate log file, with the file name determined via

snprintf(buffer,xlength,
xxxxxxxxx"AEILocalInterp.proc%d.log"
xxxxxxxxxCCTK_MyProc(NULL))

Each file begins with a header comment describing the file format in detail. After the header comment, there is one line in the file for each call on CCTK_InterpLocalUniform() which is handled by this thorn (with

log_interp_coords set to true).

Note that these log files are typically quite large, and writing them may significantly slow the simulation – you should probably only use this option for debugging purposes.

3 Supressing Interpolation

Sometimes you may want to call the interpolator, but not actually do any interpolation for some or all of the arrays:

you may be doing a query

you may be just checking values in the parameter table

you may want to interpolate varying subsets of some set of arrays

you may have no interpolation points, yet still want to do the interpolator call for other reasons

In such cases there are several possible ways you can supress some or all of the interpolations:

Set :

This supresses all interpolation. In this case interp_coords may also be NULL if desired.

Set and :

This suppresses all interpolation. However, note that some parameter-table entries like operand_indices and operation_codes are specified as being arrays whose size depends on N_output_arrays, and it’s an error for these table entries to be present with the wrong size (in this case, with any nonzero size).

Set :

This supresses the interpolation for this particular output array. This is probably the cleanest way of selectively turning individual arrays on and off. If and the input arrays aren’t needed for any queries you’re doing, then it’s useful to also set the corresponding , to supress the interpolator fetching data (which in this case will never be used) from the input arrays.

Note that the combination and is an error: conceptually you’re asking for an interpolation (the output array pointer is non-NULL) but not supplying any input data (the input array pointer is

NULL).

4 Implementation

This interpolator’s basic design is to use separate specialized code for each combination of

xxx(N_dims,xmolecule_family,xorder,xsmoothing)

i.e. in practice for each distinct choice of interpolation molecule. Maple is used to generate all the interpolation coefficients. The C preprocessor is then used to generate all the specialized code from a single master “template” file. The template code uses #ifdefs to handle lower dimensions with no extra overhead, e.g. 1-D/2-D interpolation is basically just as efficient as in a purpose-built 1-D/2-D interpolator, with no extra overhead imposed by the interpolator also supporting higher-dimensional interpolation.

The Maple code which generates the interpolation coefficients is quite general-purpose, and can handle an arbitrary dimensionality and molecule size/shape. Generating new coefficients can be rather time-consuming, though, e.g. the current coefficients for 3-D for ordersx1-4 take about 8xcpu minutes to generate using Maplex7 on a 1.7xGHzxP4.

Note that when compiling the code in the directory src/GeneralizedPolynomial-Uniform of this thorn, you may get compiler warnings about casts discarding const qualifiers from pointers in (the #include-ed file) template.c. Don’t worry – the code is actually

ok.9 9The warnings are really a design bug in C’s const-qualifier rules when pointers point to other pointers. See questionx11.10 in the (excellent!) comp.lang.c online FAQ (http://www.eskimo.com/xscs/C-faq/top.html) for further discussion. gccx2.95.3 gives the spurious warnings, but gccx3.2.1 doesn’t, nor does the Intel C compiler in any version I’ve tried. Also, C^{++} has fixed the problems in C’s const-qualifier rules, so the same code compiled as C^{++} shouldn’t give any warnings (and doesn’t on the compilers I’ve tried). You may also get compiler warnings about unused variables in this same file; again don’t worry, the code is ok.

See the README file in the source code directory AEILocalInterp/src/ and in its subdirectories for further details on the implementation.

5 Examples

5.1 A Simple Example Here’s a simple example in C, of interpolating a CCTK_REAL and a CCTK_COMPLEX 2-D array, at 5 interpolation points, using cubic interpolation.

Note: Since C allows arrays to be initialized only if the initializer values are compile-time constants, we have to

declare the interp_coords[], input_arrays[], and output_arrays[] arrays as non-const, and set their

values with ordinary (run-time) assignment statements. In C^{++}, there’s no restriction on initializer values, so we could declare the arrays const and initialize them as part of their declarations.

#definexN_DIMSxxx2
#definexN_INTERP_POINTSxxx5
#definexN_INPUT_ARRAYSxxxx2
#definexN_OUTPUT_ARRAYSxxx2
/*x(x,y)xcoordinatesxofxdataxgridxpointsx*/
#definexX_ORIGINxxx...
#definexX_DELTAxxxx...
#definexY_ORIGINxxx...
#definexY_DELTAxxxx...
constxCCTK_REALxorigin[N_DIMS]x=x{xX_ORIGIN,xY_ORIGINx};
constxCCTK_REALxdeltax[N_DIMS]x=x{xX_DELTA,xxY_DELTAxx};
/*x(x,y)xcoordinatesxofxinterpolationxpointsx*/
constxCCTK_REALxinterp_x[N_INTERP_POINTS];
constxCCTK_REALxinterp_y[N_INTERP_POINTS];
constxvoid*xinterp_coords[N_DIMS];xxxxxxxxxxxxxx/*xseexnotexabovex*/
/*xinputxarraysx*/
/*x...xnotexCactusxusesxFortranxstoragexordering,xi.e.xXxisxcontiguousx*/
#definexNXxxx10
#definexNYxxx20
constxCCTK_REALxxxxinput_realxxx[NY][NX];
constxCCTK_COMPLEXxinput_complex[NY][NX];
constxCCTK_INTxinput_array_dims[N_DIMS]x=x{xNX,xNYx};
constxCCTK_INTxinput_array_type_codes[N_INPUT_ARRAYS]
xxxxxxxx=x{xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_COMPLEXx};
constxvoid*xinput_arrays[N_INPUT_ARRAYS];xxxxxxx/*xseexnotexabovex*/
/*xoutputxarraysx*/
CCTK_REALxxxxoutput_realxxx[N_INTERP_POINTS];
CCTK_COMPLEXxoutput_complex[N_INTERP_POINTS];
constxCCTK_INTxoutput_array_type_codes[N_OUTPUT_ARRAYS]
xxxxxxxx=x{xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_COMPLEXx};
void*xoutput_arrays[N_OUTPUT_ARRAYS];xxxxxxxxxxx/*xseexnotexabovex*/
intxoperator_handle,xparam_table_handle;
operator_handlex=xCCTK_InterpHandle("myxinterpolationxoperator");
ifx(operator_handlex<x0)
xxxxxxxxCCTK_WARN(-1,x"can'txgetxinterpolationxhandle!");
param_table_handlex=xUtil_TableCreateFromString("order=3");
ifx(param_table_handlex<x0)
xxxxxxxxCCTK_WARN(-1,x"can'txcreatexparameterxtable!");
/*xinitializexthexrestxofxthexparameterxarraysx*/

interp_coords[0]x=x(constxvoidx*)xinterp_x;
interp_coords[1]x=x(constxvoidx*)xinterp_y;
input_arraysx[0]x=x(constxvoidx*)xinput_real;
input_arraysx[1]x=x(constxvoidx*)xinput_complex;
output_arrays[0]x=x(xxxxxxvoidx*)xoutput_real;
output_arrays[1]x=x(xxxxxxvoidx*)xoutput_complex;
/*xdoxthexactualxinterpolation,xandxcheckxforxerrorxreturnsx*/
ifx(CCTK_InterpLocalUniform(N_DIMS,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxoperator_handle,xparam_table_handle,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxorigin,xdelta,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxN_INTERP_POINTS,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxCCTK_VARIABLE_REAL,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxinterp_coords,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxN_INPUT_ARRAYS,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxinput_array_dims,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxinput_array_type_codes,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxinput_arrays,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxN_OUTPUT_ARRAYS,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxoutput_array_type_codes,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxoutput_arrays)x<x0)
xxxxxxxxCCTK_WARN(-1,x"errorxreturnxfromxinterpolator!");

5.2 An Example of Interpolating Derivatives

In this example we interpolate the 3-metric and its 1stxderivatives at a set of interpolation points on a 2-sphere.

Note: Since we’re not using C^{++}, we again declare the pointer arrays non-const, and “initialize” them with ordinary (run-time) assignment statements. However, in this example, for greater clarity we place these assignment statements right after the declarations. Since C allows declarations only at the start of a {x} block, not in the middle of a block, we nest the rest of the program in extra blocks (with the {x} indented 2 spaces to distinguish them from “normal” {x} pairs) to allow for further declarations. The reader will have to decide for herself whether this style is more or less ugly than separating the declarations and initializations of the pointer arrays.

#definexN_DIMSxxxxxxxx3
/*xinterpolationxpointsx*/
#definexN_INTERP_POINTSx1000
constxCCTK_REALxinterp_x[N_INTERP_POINTS],
xxxxxxxxxxxxxxxxinterp_y[N_INTERP_POINTS],
xxxxxxxxxxxxxxxxinterp_z[N_INTERP_POINTS];
constxvoid*xinterp_coords[N_DIMS];xxxxxxxxxxxxxx/*xseexnotexabovex*/
interp_coords[0]x=x(constxvoidx*)xinterp_x;
interp_coords[1]x=x(constxvoidx*)xinterp_y;
interp_coords[2]x=x(constxvoidx*)xinterp_z;
/*xdimensionsxofxthexdataxgridx*/
#definexNXxxx30
#definexNYxxx40
#definexNZxxx50
xx{
/*xinputxarraysx*/
/*xn.b.xwexusexFortranxstoragexorder:xXxisxcontiguous,xZxleastxsox*/
#definexN_INPUT_ARRAYSxx6
constxCCTK_REALxgxx_3D[NZ][NY][NX],xgxy_3D[NZ][NY][NX],xgxz_3D[NZ][NY][NX],
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxgyy_3D[NZ][NY][NX],xgyz_3D[NZ][NY][NX],
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxgzz_3D[NZ][NY][NX];
constxCCTK_INTxinput_array_dims[N_DIMS]x=x{NX,xNY,xNZ};
constxCCTK_INTxinput_array_type_codes[N_INPUT_ARRAYS]
xxxxxxxx=x{xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxCCTK_VARIABLE_REALx};
constxvoid*xinput_arrays[N_INPUT_ARRAYS];xxxxxxx/*xseexnotexabovex*/
input_arrays[0]x=x(constxvoidx*)xgxx_3D;
input_arrays[1]x=x(constxvoidx*)xgxy_3D;
input_arrays[2]x=x(constxvoidx*)xgxz_3D;
input_arrays[3]x=x(constxvoidx*)xgyy_3D;
input_arrays[4]x=x(constxvoidx*)xgyz_3D;
input_arrays[5]x=x(constxvoidx*)xgzz_3D;
xx{
/*xoutputxarraysx*/
/*x(forxbestxefficiencyxwexgroupxallxoperationsxonxaxgivenxinputxtogether)x*/
#definexN_OUTPUT_ARRAYSx24
CCTK_REAL
xxxgxx[N_INTERP_POINTS],
xxxxxxdx_gxx[N_INTERP_POINTS],xdy_gxx[N_INTERP_POINTS],xdz_gxx[N_INTERP_POINTS],

xxxgxy[N_INTERP_POINTS],
xxxxxxdx_gxy[N_INTERP_POINTS],xdy_gxy[N_INTERP_POINTS],xdz_gxy[N_INTERP_POINTS],
xxxgxz[N_INTERP_POINTS],
xxxxxxdx_gxz[N_INTERP_POINTS],xdy_gxz[N_INTERP_POINTS],xdz_gxz[N_INTERP_POINTS],
xxxgyy[N_INTERP_POINTS],
xxxxxxdx_gyy[N_INTERP_POINTS],xdy_gyy[N_INTERP_POINTS],xdz_gyy[N_INTERP_POINTS],
xxxgyz[N_INTERP_POINTS],
xxxxxxdx_gyz[N_INTERP_POINTS],xdy_gyz[N_INTERP_POINTS],xdz_gyz[N_INTERP_POINTS],
xxxgzz[N_INTERP_POINTS],
xxxxxxdx_gzz[N_INTERP_POINTS],xdy_gzz[N_INTERP_POINTS],xdz_gzz[N_INTERP_POINTS];
constxCCTK_INTxoutput_array_type_codes[N_OUTPUT_ARRAYS]
xxx=x{xCCTK_VARIABLE_REAL,
xxxxxxxxxxCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,
xxxxxxxCCTK_VARIABLE_REAL,
xxxxxxxxxxCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,
xxxxxxxCCTK_VARIABLE_REAL,
xxxxxxxxxxCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,
xxxxxxxCCTK_VARIABLE_REAL,
xxxxxxxxxxCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,
xxxxxxxCCTK_VARIABLE_REAL,
xxxxxxxxxxCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,
xxxxxxxCCTK_VARIABLE_REAL,
xxxxxxxxxxCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REAL,xCCTK_VARIABLE_REALx};
void*xoutput_arrays[N_OUTPUT_ARRAYS];xxxxxxxxxxx/*xseexnotexabovex*/
output_arrays[x0]x=x(voidx*)xxxxgxx;
output_arrays[x1]x=x(voidx*)xdx_gxx;
output_arrays[x2]x=x(voidx*)xdy_gxx;
output_arrays[x3]x=x(voidx*)xdz_gxx;
output_arrays[x4]x=x(voidx*)xxxxgxy;
output_arrays[x5]x=x(voidx*)xdx_gxy;
output_arrays[x6]x=x(voidx*)xdy_gxy;
output_arrays[x7]x=x(voidx*)xdz_gxy;
output_arrays[x8]x=x(voidx*)xxxxgxz;
output_arrays[x9]x=x(voidx*)xdx_gxz;
output_arrays[10]x=x(voidx*)xdy_gxz;
output_arrays[11]x=x(voidx*)xdz_gxz;
output_arrays[12]x=x(voidx*)xxxxgyy;
output_arrays[13]x=x(voidx*)xdx_gyy;
output_arrays[14]x=x(voidx*)xdy_gyy;
output_arrays[15]x=x(voidx*)xdz_gyy;
output_arrays[16]x=x(voidx*)xxxxgyz;
output_arrays[17]x=x(voidx*)xdx_gyz;
output_arrays[18]x=x(voidx*)xdy_gyz;
output_arrays[19]x=x(voidx*)xdz_gyz;
output_arrays[20]x=x(voidx*)xxxxgzz;
output_arrays[21]x=x(voidx*)xdx_gzz;
output_arrays[22]x=x(voidx*)xdy_gzz;
output_arrays[23]x=x(voidx*)xdz_gzz;
xx{
/*xintegerxcodesxtoxspecifyxthexderivativesx*/
constxCCTK_INTxoperand_indices[N_OUTPUT_ARRAYS]

xxx=x{x0,x0,x0,x0,
xxxxxxx1,x1,x1,x1,
xxxxxxx2,x2,x2,x2,
xxxxxxx3,x3,x3,x3,
xxxxxxx4,x4,x4,x4,
xxxxxxx5,x5,x5,x5x};
#definexDERIV(x)xxxx
constxCCTK_INTxoperation_codes[N_OUTPUT_ARRAYS]
xxx=x{xDERIV(0),xDERIV(1),xDERIV(2),xDERIV(3),
xxxxxxxDERIV(0),xDERIV(1),xDERIV(2),xDERIV(3),
xxxxxxxDERIV(0),xDERIV(1),xDERIV(2),xDERIV(3),
xxxxxxxDERIV(0),xDERIV(1),xDERIV(2),xDERIV(3),
xxxxxxxDERIV(0),xDERIV(1),xDERIV(2),xDERIV(3),
xxxxxxxDERIV(0),xDERIV(1),xDERIV(2),xDERIV(3)x};
intxoperator_handle,xparam_table_handle;
operator_handlex=xCCTK_InterpHandle("myxinterpolationxoperator");
ifx(operator_handlex<x0)
xxxxxxxxCCTK_WARN(-1,x"can'txgetxinterpolationxhandle!");
param_table_handlex=xUtil_TableCreate(UTIL_TABLE_FLAGS_DEFAULT);
ifx(param_table_handlex<x0)
xxxxxxxxCCTK_WARN(-1,x"can'txcreatexparameterxtable!");
ifx(Util_TableSetInt(param_table_handle,x3,x"order")x<x0)
xxxxxxxxCCTK_WARN(-1,x"can'txsetxorderxinxparameterxtable!");
ifx(Util_TableSetIntArray(param_table_handle,
xxxxxxxxxxxxxxxxxxxxxxxxxxN_OUTPUT_ARRAYS,xoperand_indices,
xxxxxxxxxxxxxxxxxxxxxxxxxx"operand_indices")x<x0)
xxxxxxxxCCTK_WARN(-1,x"can'txsetxoperand_indicesxarrayxinxparameterxtable!");
ifx(Util_TableSetIntArray(param_table_handle,
xxxxxxxxxxxxxxxxxxxxxxxxxxN_OUTPUT_ARRAYS,xoperation_codes,
xxxxxxxxxxxxxxxxxxxxxxxxxx"operation_codes")x<x0)
xxxxxxxxCCTK_WARN(-1,x"can'txsetxoperation_codesxarrayxinxparameterxtable!");
ifx(CCTK_InterpLocalUniform(N_DIMS,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxoperator_handle,xparam_table_handle,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxorigin,xdelta,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxN_INTERP_POINTS,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxCCTK_VARIABLE_REAL,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxinterp_coords,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxN_INPUT_ARRAYS,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxinput_array_dims,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxinput_array_type_codes,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxinput_arrays,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxN_OUTPUT_ARRAYS,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxoutput_array_type_codes,
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxoutput_arrays)x<x0)
xxxxxxxxCCTK_WARN(-1,x"errorxreturnxfromxinterpolator!");
xx}
xx}

6 Acknowledgments

Thanks to Thomas Radke for our ongoing collaboration on Cactus interpolation, to Ian Hawke and Erik Schnetter for helpful comments on the documentation, to Tom Goodale and Thomas Radke for many useful design discussions, to Erik Schnetter for bug reports, and to all the Cactus crew for a great infrastructure!

7 Parameters

11logxinterpxcoords 10Scope: privatex01  BOOLEAN

11Description: should we log the grid min(delta)max and the interpolation coordinates for each call on the interpolator?

11             10  x01  Default: false

11

8 Interfaces

General Implements: aeilocalinterp 9 Schedule This section lists all the variables which are assigned storage by thorn Numerical/AEILocalInterp. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage NONE

Scheduled Functions

CCTKxSTARTUP x aeilocalinterpxuxstartup x register cctkxinterplocaluniform() interpolation operators x

11x11After: 11driverxstartup
11 x 11Language:11c
11 x 11Type: 11function