Proposal for New Global (Grid Array) Interpolator API -- Draft 2 ================================================================ Author: Jonathan Thornburg Date: October 2002 Version: $Header: /cactus/CactusWebSite/Development/Specs/GlobalInterpolation_draft_2.txt,v 1.2 2002/11/29 15:52:44 tradke Exp $ URL: http://www.cactuscode.org/Development/Specs/GlobalInterpolation_draft_2.txt Status: draft for discussion; hopefully to be at least partially implemented in Nov-Dec 2002 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 1 ========================== * 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 interpolation of grid arrays by the CCTK_InterpGV() API. This uses variable-length argument lists, which are rather error-prone (the compiler can't check/warn about extra or missing or wrong-typed arguments). It 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 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 fairly closely modelled after the CCTK_InterpLocalUniform() API, and uses a 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.] Prototype ========= int status = CCTK_InterpGridArrays(const cGH *GH, int N_dims, int operator_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[]); 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. This 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. Details of the operation performed, and what (if any) inputs and/or outputs are specified in the parameter table, depend on the driver thorn and interpolation operator. The first implementation will likely be in PUGHInterp. 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. operator_handle (>= 0) Handle to the interpolation operator as returned by CCTK_InterpHandle(). param_table_handle (>= 0) Handle to a key-value table containing additional parameters for the interpolation. *** 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's new hyperslab routines work similarly, but he has the check in the hyperslab routines (analogous to having it in InterpGridArrays()) rather than down in the table routines as in Erik's suggestion 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 on this processor. N.b. This argument may 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. These coordinates are with respect to the coordinate system defined by coord_system_handle . N.b. These pointers, and the interpolation coordinates themselves, may vary from processor to processor in a multiprocessor Cactus run. N_input_arrays (>= 0) The number of input arrays to be interpolated. Note that if the operand_indices parameter table entry is used to specify a 1-to-many mapping of input arrays to output arrays, only the unique set of input arrays should be given here. Note also that N_input_arrays == 0 is valid, since the main purpose of the call may be (eg) to do some query or setup computation. 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. *** RFC *** We should have some value to specify "skip this interpolation", in the same way as a NULL output pointer. What's a good sentinal value for a Cactus variable index? N_output_arrays (>= 0) The number of output arrays to be returned from the interpolation. Note that this may differ from N_input_arrays , eg if the operand_indices parameter table entry is used to specify a 1-to-many mapping of input arrays to output arrays. Note also that N_output_arrays == 0 is valid, since the main purpose of the call may be (eg) to do some query or setup computation. 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. N.b. These pointers may vary from processor to processor in a multiprocessor Cactus run. N.b. If output_arrays[out] is NULL for some index out , then that interpolation is skipped. (This is useful if the main purpose of the call is (eg) to do some query or setup computation.) 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 (cf. the optional parameter-table entries use_ghost_zones and out_of_range_tolerance) CCTK_ERROR_INTERP_POINT_EXCISED one or more interpolation points is in (or too close to) an excised region of the grid (cf. the optional parameter-table entries excision_gridfn_varindex, excision_gridfn_mask, excision_value, and excision_tolerance) *** 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 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). I (Jonathan) don't like this possibility -- eg it can make it hard to come up with detailed error messages for the Cactus user. (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. Inputs/Outputs in 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. const CCTK_INT refinement_level If present, gives the mesh refinement level from which to interpolate. Note that this must be the same for all the interpolations in a single CCTK_InterpGridArrays() call. *** RFC *** We also need sentinal values (which preferably should be the same for all drivers) for "base grid" and "finest grid where data is available, subject to out_of_range_tolerance ". What should these be? const CCTK_INT local_interpolator_handle; By default the interpolator will choose a suitable local interpolator, but this may be overridden with this table entry. This should be a handle to an interpolator supporting the CCTK_InterpLocalUniform() API (or CCTK_InterpLocalNonUniform() or CCTK_InterpLocalWarped() or whatever, if/when we implement those). 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 its default choice for the local interpolator, as well as what options it requires from the local interpolator. const CCTK_INT local_interpolator_param_table_handle; Handle to a key-value table containing additional parameters for the local interpolation (eg CCTK_InterpLocalUniform() for the current sort of uniform-local-grid Cactus grids). For example, this could be used to specify that CCTK_InterpLocalUniform() should take derivatives, by specifying const CCTK_INT operand_indices[N_output_arrays]; const CCTK_INT operation_codes[N_output_arrays]; *** RFC *** Should we pass any unrecognized entries in param_table_handle table down to the local interpolator? I would say no, we either ignore them altogether, or treat them as an error -- if the user wants to pass parameters to the local interpolator, that's what the local_interpolator_param_table_handle table is for. *** RFC *** CCTK_InterpGridArrays() will typically want to specify some options of its own for the local interpolatyor. Should it be allowed to modify the local_interpolator_param_table_handle table to do this, or should it be required to clone this table and modify the clone? 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 use_ghost_zones This is a Boolean flag (0=false, nonzero=true) specifying whether or not the interpolation may use data from the ghost zones. *** RFC *** Should this be an array like out_of_range_tolerance[] (below) to individually specify this for each face of the current processor's chunk of the grid? If so, should we allow this to differ from one processor to another in a multiprocessor Cactus run? If this parameter is not specified, the default is that yes, the interpolation may use data from the ghost zones. Note that if the interpolation does actually use data from any ghost zone, it is the user's responsibility to ensure that the grid arrays are synchronized before calling the interpolator! const CCTK_REAL out_of_range_tolerance[2*N_dims]; This is just like in CCTK_InterpLocalUniform(): 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, ...]. 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 point is > TOL * coord_delta[axis] outside the grid (or from a ghost zone if use_ghost_zones is false) 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 interpolation 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 (or from a ghost zone if use_ghost_zones is false). Other values of TOL are illegal (reserved for future use). The default value of TOL should be a small positive number, to allow for floating point roundoff errors pushing an interpolation point which was meant to be just on the grid boundary, slightly outside the grid. *** RFC *** Should out_of_range_tolerance[] also depend on refinement_level , eg by making it a [2*N_dims*N_refinement_levels] array with the appropriate semantics? 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 typically it must off-center the interpolation molecule. Figuring out just which direction in which to move 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 help in this: If specified, it should be the Cactus variable index of an "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. const CCTK_REAL excision_tolerance; This is similar to out_of_range_tolerance, but specifies the handling for excision regions: If excision_tolerance == -1.0, then it's an error if an interpolation point is close enough to an excised point to require off-centering the interpolation molecule (or more generally, choosing a centering different from that which would have been chosen if there were no excised points) *** RFC *** What other values of excision_tolerance should be supported? In particular, should we provide some way of specifying how far into the excised region an interpolation point is allowed to be? Just how should we do this, given that the excised region is only known at a "Lego" grid-point granularity? 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 operator_handle = CCTK_InterpHandle("generalized polynomial interpolation"); if (operator_handle < 0) CCTK_WARN(-1, "can't get operator handle!"); /* this sets the interpolation order */ 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, operator_handle, param_table_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 ---