Proposal for New Global (Grid Array) Interpolator API -- Draft 4 ================================================================ Author: Jonathan Thornburg Date: 6.Jan.2003 Version: $Header: /cactus/CactusWebSite/Development/Specs/GlobalInterpolation_draft_4.txt,v 1.2 2003/01/22 15:31:59 tradke Exp $ URL: http://www.cactuscode.org/Development/Specs/GlobalInterpolation_draft_4.txt Status: draft for discussion; implementation is currently (Jan 2004) in progress Metadata ======== This proposal has a number of open points marked by "*** RFC ***" (RFC = Request For Comments). Comments on these, as well as more general discussion on this proposal, should probably be on the developers@cactuscode.org mailing list. Major Changes from draft 3 --> 4 ================================ * merge local and global parameter table; only one table handle is now passed to the grid-array interpolation API Major Changes from draft 2 --> 3 ================================ * add more discussion of the composite nature of the global interpolator, and how it relates to and makes use of the (a) local interpolator * drop global interpolator handle, move local interpolator handle and local interpolator parameter table out of global parameter table into the explicit argument list * revise handling of excision and boundary off-centering as per mailing-list discussion between Jonathan and Erik * slightly revise handling of ghost/boundary zones * define sentinal value -1 for CCTK variable indices Major Changes from draft 1 --> 2 ================================ * Clarify why "no-op" values like N_interp_points == 0, N_input_points == 0, output_arrays[out] == NULL, etc, may be useful * add excision support, based on Denis Pollney's new CactusEinstein/SpaceMask thorn (there are still some RFCs in how we do this) * add separate return code for trying to interpolate in (or too far into) an excised region * clarify documentation of out_of_range_tolerance * split parameter table into separate parameter tables for global interpolation (this API) and to be passed to local interpolator (CCTK_InterpLocalUniform() API for the current flavors of Cactus grids) * remove some optional parameters which would now go in the local interpolator's parameter table * add optional use_ghost_zones parameter to say whether or not the interpolation is allowed to use data from ghost zones * added some parameters to specify refinement levels for Carpet and other mesh-refinement drivers -- Erik, Ian, Scott, and other mesh-refinement gurus, please comment if this is good/bad for you! Motivation ========== Currently Cactus provides "global" (multiprocessor) interpolation of grid arrays by the CCTK_InterpGV() API. This API uses variable-length argument lists, which are rather error-prone (the compiler can't check/warn about extra or missing or wrong-typed arguments). This API also doesn't provide some features various users want (eg Jonathan Thornburg wants derivatives, Erik Schnetter wants hyperslabs), and it lacks extensibility. This spec is for a new global grid-array interpolation API, CCTK_InterpGridArrays(), which whould provide a superset of the CCTK_InterpGV() functionality, and we hope would eventually replace CCTK_InterpGV(). The new API is modelled after, and builds on, the CCTK_InterpLocalUniform() API for local interpolation, and uses Cactus key-value table for extensionality in a similar way to that API. [The CCTK_InterpLocalUniform() API is documented in section E of the Cactus Users' Guide. The key-value table inputs/outputs for this API are currently documented in the thorn guide for thorn LocalInterp.] The semantics of CCTK_InterpGridArrays() are mostly independent of which Cactus driver is being used, but the implementation will depend on, and make use of, driver-specific internals. Thus the Cactus flesh will supply CCTK_InterpGridArrays() via a registration API in same manner as CCTK_InterpGV(), with driver-specific interpolation thorns registering implementations. The first implentation of CCTK_InterpGridArrays() will likely be in PUGHInterp, but we hope is that other drivers (in particular Carpet) will follow soon. C Prototype =========== int status = CCTK_InterpGridArrays(const cGH *GH, int N_dims, int local_interp_handle, int param_table_handle, int coord_system_handle, int N_interp_points, int interp_coords_type_code, const void* const interp_coords[], int N_input_arrays, const CCTK_INT input_array_variable_indices[], int N_output_arrays, const CCTK_INT output_array_type_codes[], void* const output_arrays[]); Fortran Prototype ================= call CCTK_InterpGridArrays(status, . GH, . N_dims, . local_interp_handle, . param_table_handle, . coord_system_handle, . N_interp_points, . interp_coords_type_code, . interp_coords, . N_input_arrays, . input_array_variable_indices, . N_output_arrays, . output_array_type_codes, . output_arrays) integer status CCTK_POINTER_TO_CONST GH integer N_dims integer param_table_handle integer local_interp_handle integer coord_system_handle integer N_interp_points integer interp_coords_type_code CCTK_POINTER_TO_CONST(*) interp_coords integer N_input_arrays CCTK_INT(*) input_array_variable_indices integer N_output_arrays CCTK_INT(*) output_array_type_codes CCTK_POINTER(*) output_arrays This function interpolates a list of Cactus grid arrays (in a multiprocessor Cactus run these are generally distributed over processors). The grid topology and coordinates are implicitly specified via a Cactus coordinate system. Like CCTK_InterpGV(), CCTK_InterpGridArrays() is a collective operation: in a multiprocessor Cactus run the caller *must* call this function synchronously (in parallel) on *each* processor, passing identical arguments except (optionally) for N_interp_points # may differ from processor to processor interp_coords[] # coordinates may differ from processor to processor output_arrays[] # pointers may differ from processor to processor The set of interpolation points may differ from processor to processor in a multiprocessor Cactus run; it may even be empty on some (or all!) processors. The interpolation points may be anywhere in the global Cactus grid. In a multiprocessor Cactus run they may vary from processor to processor; each processor will get whatever interpolated data it asks for, with the interpolator taking care of whatever interprocessor communication may be necessary. The computation is optimized for the case of interpolating a number of arrays at a time; in this case all the interprocessor communication can be done together, and the same interpolation coefficients can be used for all the arrays. Global vs Local Interpolation ============================= CCTK_InterpGridArrays() specifies a global (in general multiprocessor) interpolation. This is a composite operation, combining interprocessor communication with local interpolation. On each processor: (1) Send the coordinates of each interpolation point to the processor which "owns" that part of the grid. (2) If the Cactus driver uses mesh refinement, choose the appropriate grid from which to interpolate. (3) Do a local interpolation using CCTK_InterpLocalUniform(). [This assumes we're interpolating from a uniform grid, which are the only type of grid Cactus supports at present. If/when Cactus adds support for more general grids, then this local interpolation would obviously use the corresponding API for the actual type of grid from which we're interpolating.] The user should be able to specify any compatible local interpolation operator for this step. (4) Send the results of the local interpolation back to the processor(s) which requested those interpolation points. Return Value ============ 0 ==> success (error codes are described below) Arguments ========= GH (!= NULL) Pointer to a valid Cactus grid hierarchy. N_dims (>= 1) Number of dimensions in which to interpolate. This must be <= the dimensionality of the coordinate system defined by coord_system_handle . The default case is that it's ==; see the discussion of the interpolation_hyperslab_handle parameter-table entry for the < case. param_table_handle (>= 0) Handle to a key-value table containing zero or more additional parameters for the interpolation operation. The table can be modified by the local and/or global interpolation routine(s). For example, a table might be used to specify that the local interpolator should take derivatives, by specifying const CCTK_INT operand_indices[N_output_arrays]; const CCTK_INT operation_codes[N_output_arrays]; Also, the global interpolator will typically need to specify some options of its own for the local interpolator. These will override any entries with the same keys in the param_table_handle table. The discussion of individual table entries below says if these are modified in this manner. Finally, the param_table_handle table can be used to pass back arbitrary information by the local and/or global interpolation routine(s) by adding/modifying appropriate key/value pairs. *** RFC *** How should the user indicate that there are no parmameters? - my preferred solution is for the user to pass the handle of an empty table, eg one newly created with Util_CreateTable(). - in a recent discussion on the developers mailing list, Erik Schnetter suggested we have a special "/dev/null" handle which would always mean an empty table (any data written to it would be silently discarded) - Thomas Radke suggests a similar idea, but having the check here in InterpGridArrays() rather than down in the table routines as in Erik's suggestion local_interp_handle; Handle to the local interpolation operator as returned by CCTK_InterpHandle(). (See section B8.3 "Interpolation Operators" in the Cactus Users' Guide for a discussion of interpolation operator handles and registration.) It is the caller's responsibility to ensure that the specified interpolator supports any optional parameter-table entries that CCTK_InterpGridArrays() passes to it. Each thorn providing a CCTK_InterpGridArrays() interpolator should document what options it requires from the local interpolator. coord_system_handle (>= 0) Cactus coordinate system handle (as returned by CCTK_CoordSystemHandle()) defining the mapping between integer grid subscripts and floating-point coordinates. N_interp_points (>= 0) The number of interpolation points requested by this processor. This may (and typically will) vary from processor to processor in a multiprocessor Cactus run. interp_coords_type_code One of the CCTK_VARIABLE_* type codes, giving the data type of the interpolation-point coordinate arrays pointed to by interp_coords[] . interp_coords (!= NULL) (Pointer to) an array of N_dims pointers to 1-D arrays giving the coordinates of the interpolation points requested by this processor. These coordinates are with respect to the coordinate system defined by coord_system_handle . These pointers, and the interpolation coordinates themselves, may (and typically will) vary from processor to processor in a multiprocessor Cactus run. N_input_arrays (>= 0) The number of input arrays to be interpolated. If N_input_arrays == 0 then no interpolation is done; such a call may be useful for setup, interpolator querying, etc. If the operand_indices parameter table entry is used to specify a nontrivial (eg 1-to-many) mapping of input arrays to output arrays, only the unique set of input arrays should be given here. input_array_variable_indices (!= NULL) (Pointer to) an array of N_input_arrays Cactus variable indices (as returned by CCTK_VarIndex() ) specifying the input grid arrays for the interpolation. If input_array_variable_indices[in] == -1 for some index or indices in , then that interpolation is skipped. (This may be useful if the main purpose of the call is (eg) to do some query or setup computation.) N_output_arrays (>= 0) The number of output arrays to be returned from the interpolation. If N_input_arrays == 0 then no interpolation is done; such a call may be useful for setup, interpolator querying, etc. Note that N_output_arrays may differ from N_input_arrays , eg if the operand_indices parameter table entry is used to specify a nontrivial (eg many-to-1) mapping of input arrays to output arrays, If such a mapping is specified, only the unique set of output arrays should be given here. output_array_type_codes (!= NULL) (Pointer to) an array of N_output_arrays CCTK_VARIABLE_* type codes giving the data types of the 1-D output arrays pointed to by output_arrays[] . output_arrays (!= NULL) (Pointer to) an array of N_output_arrays pointers to the (caller-supplied) 1-D output arrays for the interpolation. If output_arrays[out] is NULL for some index or indices out , then that interpolation is skipped. (This may be useful if the main purpose of the call is (eg) to do some query or setup computation.) These pointers may (and typically will) vary from processor to processor in a multiprocessor Cactus run. However, any given pointer must be either NULL on all processors, or non-NULL on all processors. Errors ====== UTIL_ERROR_NULL_POINTER one or more inputs which should be non-NULL pointers, were NULL UTIL_ERROR_BAD_INPUT one or more of the inputs is invalid UTIL_ERROR_NO_MEMORY unable to allocate memory UTIL_ERROR_BAD_HANDLE one or more of the Cactus handles manipulated by this API is invalid CCTK_ERROR_INTERP_POINT_X_RANGE one or more interpolation points is outside (or too close to a boundary of) the grid CCTK_ERROR_INTERP_POINT_EXCISED one or more interpolation points is in (or too close to) an excised region of the grid *** RFC *** What should the interpolator do about the other points if an out-of-range or excised point is detected? Some possibilities: (a) Abort the interpolation (not writing any more values to the output arrays past the first-detected out-of-range or excised point, leaving unchanged whatever garbage might have been in the output arrays) (b) Continue the interpolation with other points, but don't write anything at all to the output arrays for the out-of-range or excised point(s). (c) Continue the interpolation with other points, and write NaNs to the output arrays for the out-of-range or excised point(s). [[if there are Cactus platforms which don't support NaNs, what should we do there? Maybe have a magic "sentinal" value, eg -9.999e99 or suchlike?]] At present LocalInterp's CCTK_InterpLocalUniform() does (a), and also defines the parameter-table outputs CCTK_INT out_of_range_point; # which point is it? CCTK_INT out_of_range_axis; # which coordinate axis? CCTK_INT out_of_range_end; # whidh +/- end of the grid? to report further information about the out-of-range point. For (b) or (c) this "report further information" mechanism would need to be generalized somehow, because there might be multiple out-of-range or excised points. Some possibilities: (0) Don't report details on the out-of-range or excised point(s). (1) Just report details for the first-detected out-of-range or excised point. [LocalInterp's CCTK_InterpLocalUniform() does this for out-of-range points (it doesn't support excision yet), and it's specifically documented that different interpolators may differ on which point they detect first (and hence report upon) if there are multiple out-of-range points.] (2) Report the total number of out-of-range and/or excised points, but only give details for the first-detected one. (3) Report details of all out-of-range or excised points (maybe make this report itself optional, eg conditional on the caller setting a parameter-table entry to say she wants the info). The tricky part here is how to pass back potentially-large arrays whose size the caller wouldn't know in advance. Having the caller pre-allocate arrays of size N_interp_points would be very wasteful in the no-error case. Moreover, N_interp_points might be large, so we wouldn't want to put the arrays themselves directly in the parameter table. About the best solution I can see would be for the interpolator malloc() the arrays and return pointers to them in a parameter-table entry. The caller would be responsible for freeing the arrays and deleting the parameter-table entry when done with them. I (Jonathan) think the combination (c) and (2) is probably the cleanest, maybe with an option to also do (3). (Tht way the user would always get back the details of at least one out-of-range point, even if the malloc() failed in (3).) One more question: what should the error code be if there are both out-of-range and excised points? Inputs/Outputs in global parameter table ======================================== Unless noted otherwise, these must be the same on all processors. const CCTK_INT order; Most interpolators will have this as a mandatory parameter, to give the interpolation order. const CCTK_INT input_array_time_levels[N_input_arrays]; If present, gives the Cactus time levels of the input arrays. If this table entry isn't present, all the inputs are taken from the current (0) time level. const CCTK_INT min_refinement_level, max_refinement_level; If present, these specify a range of the refinement levels from which the interpolator should interpolate. 0 is the base grid, positive values mean refined grids, and -1 means the finest grid available. For each interpolation point, the interpolator will interpolate from the finest grid for which data is available without violating the boundary or excision restrictions specified by the other parameters below. *** RFC *** This means that different points on the same processor could be interpolated from different grids. AMR gurus: is this going to be a problem for efficient implementation? const CCTK_INT interpolation_hyperslab_handle; If this is specified, then interpolation is done only with the subspace defined by this hyperslab. For example, this may be used to do 2-D or 1-D interpolation in 3-D Cactus grid functions (Erik Schnetter wants this for some of his boundary conditions). In this case the N_dims argument, and all others which are dimensioned by N_dims , refer to the dimensionality of the interpolation (= that of the hyperslab), *not* to that of the full grid array. Any given interpolator may or may not support this option; if the interpolator doesn't support it then specifying it gives a UTIL_ERROR_BAD_INPUT error return. Note that if this option is supported, then using it will probably be considerably more efficient than explicitly getting the hyperslab and doing an N_dims-D interpolation in the hyperslab data. (The idea is for the interpolator to use the non-contiguous-input-array features of CCTK_InterpLocalUniform() to essentially do the hyperslabbing on each processor "for free" as part of the interpolation.) const CCTK_INT N_boundary_points_to_omit[2*N_dims] The array elements are matched up with the axes and minimum/maximum ends of the overall Cactus grid (the union over all processors) in the order [x_min, x_max, y_min, y_max, z_min, z_max, ...]. For each such end of the grid, the interpolator will *not* use data from the specified number of points closest to the end. 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 data from any boundary points (either physical or symmetry boundary). * Setting this parameter for some end To the width of a symmetry zone on that end means to not use any points from that symmetry zone. Note that this parameter refers only to *physical* grid boundarys and/or symmetry boundaries, not to interprocessor boundaries in multiprocessor cactus runs. The interpolator is always allowed to use values in interprocessor ghost zones. (I.e. it is always the user's responsibility to ensure that variables are synchronized between processors before interpolating them.) In a multiprocessor Cactus run, the global interpolator will modify this parameter in the parameter table on each processor (to properly describe the boundary handling for that processor's chunk of the grid) before calling the local interpolator. [N.b. this assumes that the local interpolator supports this parameter with the same semantics. LocalInterp's CCTK_InterpLocalUniform() doesn't do this yet, but I (Jonathan) will modify it soon to do so.] 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] const CCTK_REAL excision_extrapolation_tolerance[2*N_dims] These table entries control the interpolator's behavior near grid boundaries and excised regions. We assume that far from any grid boundaries or excised regions, the local interpolator has some default centering for its interpolation molecules. For example, given grid points at integer coordinates, an interpolator using 3-point centered molecules might use molecules containing the points [-1, 0,1] for any interpolation point in [-0.5,0.5). An interpolator using 4-point centered molecules might use molecules containing the points [-1, 0, 1, 2] for any interpolatino point in [0,1). These table entries control the local interpolator's willingness to off-center (change away from the default centering) its interpolation molecules near boundaries and excised points. The boundary_*_tolerance array elements are matched up with the axes and minimum/maximum ends of the overall Cactus grid (the union over all processors) in the order [x_min, x_max, y_min, y_max, z_min, z_max, ...]. The excision_*_tolerance array elements are matched up in the same order, except that "minimum/maximum" are now interpreted as "decreasing/increasing when moving into the excised region". *** RFC *** Should these parameters also depend on the refinement level in a mesh-refinement code, eg by making them [2*N_dims*N_refinement_levels] arrays with the appropriate semantics? 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 excised points; to make this a region we take the lego-block bounding box of the non-excised grid points. *** RFC *** What to do if the set of excised points isn't convex? * The "default-centering region" as that region in interpolation-point space where an interpolation can use the default molecule centering, i.e. where the default-centering molecules require data only from the data-valid region. For example, suppose we have the following 2-D quadrant grid with valid points * and excised points x, and suppose we're using the same molecule-centering strategy as the previous example. Then the lines show the boundary of the data-valid region (marked V) and the default-centering regions for centered 3-, 4-, and 5-point interpolation molecules. (If x increases to the right and y up, then these would be governed by the excision tolerances for x_min and y_min, because x and y are *decreasing* when moving into the excision region.) * * * * * * * * * * * * * * * * 5 ---------------+ 4 -*---*---*---* | * * * * 3 -----------+ | +---+ V -*---*---* | *---* | * * * | +---+ | | x x *---* | * | * * * | | | | x x x * | * | * * * | | | | V 3 4 5 The behavior of the interpolator is then as follows: If off_centering_tolerance == 0 and extrapolation_tolerance == 0 then no off-centering is allowed: the interpolator reports an error (CCTK_ERROR_INTERP_POINT_X_RANGE or CCTK_ERROR_INTERP_POINT_EXCISED as appropriate) if any interpolation point is outside the default- -centering region. If off_centering_tolerance > 0 and extrapolation_tolerance == 0 then the interpolator will report an error if any interpolation point is > off_centering_tolerance * coord_delta[axis] outside the default-centering region in any axis. If an interpolation point is outside the default-centering region by <= this distance in each axis, then the interpolator will off-center the interpolation molecules as necessary. If off_centering_tolerance == 0 and extrapolation_tolerance > 0 then the interpolator will report an error if any interpolation point is > extrapolation_tolerance * coord_delta[axis] outside the valid-data region in any axis. If an interpolation point is outside the valid-data region by <= this distance in each axis, then the interpolator will off-center the interpolation molecules as necessary. If off_centering_tolerance > 0 and extrapolation_tolerance > 0 then the interpolator aborts with a UTIL_ERROR_BAD_INPUT error. If any of these 4 table entries aren't specified, the defaults are off_centering_tolerance = 0 and extrapolation_tolerance = 1.0e-10. In a multiprocessor Cactus run, the global interpolator will modify boundary_off_centering_tolerance boundary_extrapolation_tolerance in the parameter table on each processor (to properly describe the boundary handling for that processor's chunk of the grid) before calling the local interpolator. [N.b. this assumes that the local interpolator supports this parameter with the same semantics. LocalInterp's CCTK_InterpLocalUniform() doesn't do this yet, but I (Jonathan) will modify it soon to do so.] const CCTK_INT excision_gridfn_varindex; const CCTK_INT8 excision_gridfn_mask; const CCTK_INT8 excision_value; Either all of these should be specified, or none of them. If they're specified, they describe an excision mask, used to flag certain grid points as "excised" and thus unavailable for interpolation. The semantics are that excision_gridfn_varindex is the Cactus variable index of the excision mask grid array. A point is taken to be "excised" if an only if SpaceMask_CheckStateBits(ptr, ijk, excision_gridfn_mask, excision_gridfn_mask) is nonzero, where SpaceMask_CheckStateBits() is the macro defined in CactusEinstein/SpaceMask, and ptr is the CCTK_VarDataPtr() of the excision_gridfn_varindex grid array. Note that we will also need to add excision support to the CCTK_InterpLocalUniform() API and implementation. const CCTK_INT excision_table_gridfn_varindex; If the interpolator discovers that a nominally-centered interpolation molecule would require data from an excised region, then (if it doesn't treat this as an error) it must off-center the interpolation molecule. Figuring out just which direction in which to move the molecule in order to have the molecule be entirely in the non-excised part of the grid may be somewhat complicated and/or expensive. excision_table_gridfn_varindex can be used to make this more efficient: If specified, it should be the Cactus variable index of a "space table field" grid array (of integral type), with semantics as documented in CactusEinstein/SpaceMask: for each excised point, the value of this grid array should be the handle of a table giving information about the excised region. If the key const CCTK_REAL centroid[N_dims] is present in this table, it is taken to give the centroid of the exised region; the interpolator will then move the molecule directly away from this centroid to return to the non-excised part of the grid. (Some interpolators may require that excision_table_gridfn_varindex be present and contain this key if any of the other excision table entries are present.) Examples ======== Here's a simple example to do quartic 3-D interpolation of a real and a complex grid array, at 1000 interpolation points (we assume some code somewhere has stored values in the arrays marked with /* input */ comments). (If this is a multiprocessor Cactus run, each processor may have it's own set of 1000 interpolation points.) The example code is legal C++, but not actually legal C, (but only) because automatic arrays are being initialized, and because declarations and assignment statements are interleaved for code clarity and const-correctness. It should be obvious how to transform (obfuscate!) the code into C with explicit assignments and reordering of declarations to precede assignments. --- begin example --- #include "cctk.h" #include "util_Table.h" #define N_DIMS 3 #define N_INTERP_POINTS 1000 #define N_INPUT_ARRAYS 2 #define N_OUTPUT_ARRAYS 2 const cGH *GH; /* input */ /* interpolation points */ const CCTK_REAL interp_x[N_INTERP_POINTS], /* input */ interp_y[N_INTERP_POINTS], /* input */ interp_z[N_INTERP_POINTS]; /* input */ const CCTK_INT interp_coord_type_codes[N_DIMS] = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL }; const void *const interp_coords[N_DIMS] = { (const void *) interp_x, (const void *) interp_y, (const void *) interp_z }; /* input and output arrays */ const CCTK_INT input_array_variable_indices[N_INPUT_ARRAYS] = { CCTK_VarIndex("my_thorn::real_array"), /* no error checking */ CCTK_VarIndex("my_thorn::complex_array") }; /* here :( */ const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS] = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_COMPLEX }; CCTK_REAL output_for_real_array [N_INTERP_POINTS]; /* output stored here * / CCTK_COMPLEX output_for_complex_array[N_INTERP_POINTS]; /* output stored here * / void *const output_arrays[N_OUTPUT_ARRAYS] = { (void *) output_for_real_array, (void *) output_for_complex_array }; const int local_interp_handle = CCTK_InterpHandle("Lagrange polynomial interpolation"); if (local_interp_handle < 0) CCTK_WARN(-1, "can't get operator handle!"); const int param_table_handle = Util_TableCreateFromString("order=4"); const int coord_system_handle = CCTK_CoordSystemHandle("cart3d"); if (coord_system_handle < 0) CCTK_WARN(-1, "can't get coordinate-system handle!"); const int status = CCTK_InterpGridArrays(GH, N_DIMS, param_table_handle, local_interp_handle, coord_system_handle, N_INTERP_POINTS, interp_coord_type_codes, interp_coords, N_INPUT_ARRAYS, input_array_variable_indices, N_OUTPUT_ARRAYS, output_array_type_codes, output_arrays); --- end example ---