Proposal for New Global (Grid Array) Interpolator API ===================================================== Author: Jonathan Thornburg Date: October 2002 Version: $Id: GlobalInterpolation_draft_1.txt,v 1.2 2002/10/22 14:47:18 tradke Exp $ URL: http://www.cactuscode.org/Development/Specs/GlobalInterpolation_draft_1.txt Status: draft for discussion; hopefully to be at least partially implemented in Nov-Dec 2002 Metadata ======== This proposal has several 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. 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. However, 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. 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. 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. 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. 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. 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_INTEPR_POINT_X_RANGE one or more interpolation points is outside (or too close to a boundary of) the grid (cf. the optional parameter-table entry out_of_range_tolerance) *** RFC *** What should the interpolator do about the other points if an out-of-range 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 point). (b) Continue the interpolation with other points, but don't write anything at all to the output arrays for the out-of-range point(s). (c) Continue the interpolation with other points, and write NaNs to the output arrays for the out-of-range 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 points. Some possibilities: (0) Don't report details on the out-of-range 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 point. [LocalInterp's CCTK_InterpLocalUniform() does this, 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) Somehow report on all out-of-range 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. Note that any unrecognized table entries will be passed through to the underlying local interpolator. 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. *** RFC *** What about mesh refinement? If a single GH contains all refinement levels, then maybe there should be optional table arguments to choose which refinment level we want to interpolate from. Or does a GH contain only a single refinment level? 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 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. *** RFC *** Should we define a query API (presumably using the parameter table) to let user code find out whether or not this option is supported? 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_REAL out_of_range_tolerance[2*N_dims]; These are modelled after LocalInterp's CCTK_InterpLocalUniform(), but should presumably only refer to the physical grid boundaires, not to the interprocessor boundaries. See the discussion of the CCTK_ERROR_INTEPR_POINT_X_RANGE error code above for various options on reporting further information about out-of-range points. *** RFC *** Should ghost zones be included or not? Or should the caller be able to select either behavior (eg with another table entry)? const char molecule_family[]; This is exactly the same as in LocalInterp's CCTK_InterpLocalUniform(): if an interpolator supports different shapes of molecules (eg cube vs "sphere"), this parameter may be used to select between them const CCTK_INT operand_indices[N_output_arrays]; const CCTK_INT operation_codes[N_output_arrays]; These are exactly the same as in LocalInterp's CCTK_InterpLocalUniform(): they can be used to have the interpolator take (eg) derivatives. (Jonathan needs this feature for AHFinderDirect.) *** RFC *** There should also be table entries to specify the excision mask. How should we do this? (Maybe delay this until we have a final version of Denis Pollney's new excision-mask scheme?) 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 ---