diff options
Diffstat (limited to 'archive/api2')
| -rw-r--r-- | archive/api2 | 460 |
1 files changed, 460 insertions, 0 deletions
diff --git a/archive/api2 b/archive/api2 new file mode 100644 index 0000000..422be0d --- /dev/null +++ b/archive/api2 @@ -0,0 +1,460 @@ +This is version 2 of my proposal for the new interpolator API. The +main changes since version 1 are +* function names have changed (longer, more descriptive) +* output arrays may now be N-dimensional (before they were only 1-D) +* a bunch more optional stuff in the table to specify array strides etc + +Don't be scared by the length, for most uses it's not that complicated! +There are some examples below... + + + +Synopsis +======== + + int status = CCTK_InterpLocalArrays(arguments described below) + int status = CCTK_InterpGridArrays(arguments described below) + +return is 0 for ok, various -ve values for error codes + + + +Function Arguments +================== + +for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays(): + /***** misc arguments *****/ + const cGH *GH; + /* note N_dims is the number of dimensions in the *interpolation*; */ + /* this may be smaller than the number of dimensions of the input arrays */ + /* if the storage indexing is set up appropriately (eg to interpolate */ + /* in 1-D lines or 2-D planes of 3-D grid arrays) */ + int N_dims; + int operator_handle, int coord_system_handle; + int param_table_handle; /* handle to "parameter table", a key-value */ + /* table, see below for table entries */ + + /***** arguments specifying the interpolation points *****/ + int N_interp_points; + /* array of CCTK_VARIABLE_* codes giving the data types */ + /* of the arrays pointed to by interp_coords[] */ + const int interp_coord_type_codes[N_dims]; + /* array[N_dims] of pointers to arrays[N_interp_points] */ + /* giving x,y,z,... coordinates of interpolation points */ + const void *const interp_coords[N_dims]; + + /***** arguments specifying the inputs (the data to be interpolated) *****/ + int N_input_arrays; +for CCTK_InterpLocalArrays(): + /* array of CCTK_VARIABLE_* codes giving data types of input arrays */ + const int input_array_type_codes[N_input_arrays]; + /* array of input array dimensions (common to all input arrays) */ + const int input_array_dims[N_dims]; + /* array of pointers to input arrays */ + const void *const input_arrays[N_input_arrays]; +for CCTK_InterpGridArrays(): + /* array of CCTK variable indices of input arrays */ + const int input_array_variable_indices[N_input_arrays]; + +for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays() again: + /***** arguments specifying the outputs (the interpolation results) *****/ + int N_output_arrays; + /* array of CCTK_VARIABLE_* codes giving data types of output arrays */ + const int output_array_type_codes[N_output_arrays]; + /* array[N_output_arrays] of pointers to output arrays[N_interp_points] */ + void *const output_arrays[N_output_arrays]; + + + +Information Passed in the Parameter Table +========================================= + +The "parameter table" may be used to specify non-default storage indexing +for input or output arrays, and/or various options for the interpolation +itself. Some interpolators may not implement all of these options. + + +Storage Indexing Options +------------------------ + +Sometimes one of the "arrays" used by the interpolator isn't contiguous +in memory. For example, we might want to do 2-D interpolation within a +plane of a 3-D grid array, and/or the grid array might be a member of a +compact group. To support this, we use several optional table entries +(these should be supported by all interpolation operators): + +For the input arrays, we use + + const int input_array_offsets[N_input_arrays]; + /* next 3 table entries are shared by all input arrays */ + const int input_array_strides [N_dims]; + const int input_array_min_subscripts[N_dims]; + const int input_array_max_subscripts[N_dims]; + +Then for input array number a, the generic subscripting expression for +the 3-D case is + data_pointer[offset + i*istride + j*jstride + k*kstride] +where + data_pointer = input_arrays[a] + offset = input_array_offsets[a] + (istride,jstride,kstride) = input_array_stride[] +and where (i,j,k) run from input_array_min_subscripts[] to +input_array_max_subscripts[] inclusive. + +The defaults are offset=0, stride=determined from input_array_dims[] +in the usual Fortran manner, input_array_min_subscripts[] = 0, +input_array_max_subscripts[] = input_array_dims[]-1. If the stride +and max subscript are both specified explicitly, then the +input_array_dims[] function argument is ignored. + +For CCTK_InterpGridArrays() operating on a member of a compact group +the offset and strides are interpreted in units of _grid_points_. This +has the advantage that interpolator calls need not be changed if a grid +array is changed from being simple to/from compact. In terms of actual +memory addressing, then, the internal subscripting expression for this +case would be + group_data_pointer[offset + member_number + i*istride*N_members + + j*jstride*N_members + + k*kstride*N_members] + +By default the interpolation-point coordinates and the output arrays +are all contiguous 1-D arrays. This may be changed with the table +entries + + const int interp_coords_offsets[N_dims]; + const int output_array_offsets[N_output_arrays]; + /* next 4 table entries are shared by all interp coords and output arrays */ + const int interp_point_N_dims; + const int interp_point_strides [interp_point_N_dims]; + const int interp_point_min_subscripts[interp_point_N_dims]; + const int interp_point_max_subscripts[interp_point_N_dims]; + +For example, if we wanted to do 3-D interpolation, interpolating a value +at each non-ghost-zone point of a 2-D grid of points, with the grid point +coordinates stored as 2-D arrays, we would use + N_dims = 3 + interp_point_N_dims = 2 + interp_point_strides[] = set up from the full size of the 2-D grid + interp_point_{min,max}_subscripts[] = specify the non-ghost-zone points + of the 2-D grid + +Excision Options +---------------- + +Some interpolators may specifically support excision, where a mask array +(same dimensionality and indexing as the input arrays) is used to mark +some grid points as valid (ok to use data there) and other grid points +as invalid (the interpolator isn't allowed to use data there). + +If an interpolator supports this, it should use the following optional +parameters: + +for CCTK_InterpLocalArrays(); + const int mask_type_code; /* one of the CCTK_VARIABLE_* codes */ + const void *const mask_array; +for CCTK_InterpGridArrays(): + const int mask_variable_index; + +for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays(): + /* we consider a grid point to be valid if and only if the mask */ + /* has a value in the closed interval [mask_valid_min,mask_valid_max] */ + /* n.b. the caller should beware of possible rounding problems here; */ + /* it may be appropriate to widen the valid interval slightly */ + /* if the endpoints aren't exactly-representable floating-point */ + /* values */ + const mask_type mask_valid_min, mask_valid_max; + +The same type of storage options supported for the input and/or output +arrays, are also supported for the mask; the mask may have its own offset, +but shares any input-array stride and/or min/max subscript specification: + + const int mask_offset; + + +The remaining parameter-table options are specific to the new interpolator +I'm currently implementing for PUGHInterp. This registers (only) a single +operator, "generalized polynomial interpolation". + + +Interpolation Order and Molecule Family +--------------------------------------- + +The mandatory parameter + + const int order; + +sets the order of the interpolating polynomial (1=linear, 2=quadratic, +3=cubic, etc). Thus the simplest call can just use (eg) + Util_TableCreateFromString("order=3") +for cubic interpolation. + +All the remaining parameters in the table are optional; if they're +omitted defaults will be supplied. + + /* this selects one of a family of related operators */ + /* the default (and the only one I'm implementing right now) */ + /* is "cube" to use the usual hypercube-shaped molecules */ + const char *const molecule_family; + +Smoothing +--------- + +The way I'm implementing the interpolation it's easy to also do +Savitzky-Golay type smoothing (= moving least-square fitting, cf +Numerical Recipes 2nd edition section 14.8). This is controlled by +the parameter + + const int smoothing; + +which says 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. + +This type of 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. + +Derivatives +----------- + +This interpolator can optionally (and again at no extra cost) take +partial derivatives as part of the interpolation: + const int operand_indices[N_output_arrays]; + const int opcodes [N_output_arrays]; +The semantics here are that + output array[i] = op(input array[ operand_indices[i] ]) +where op is specified as an integer operation code as described below. + +Note that the array operand_indices[] doesn't directly name the inputs, +but rather gives the indices (0-origin) in the list of inputs. This +allows for a more efficient implementation in the case where some of +the input arrays have many different operations applied to them. + +The operations are coded based on the decimal digits of the integer: +each decimal digit means to take the derivative in that direction; +the order of the digits in a number is ignored. For example: + 0 = no derivative, "just" interpolate + 1 = interpolate d/dx1 (derivative along x1 coordinate) + 2 = interpolate d/dx2 (derivative along x2 coordinate) + 11 = interpolate d^2/dx1^2 (2nd derivative along x1 coordinate) + 22 = interpolate d^2/dx2^2 (2nd derivative along x2 coordinate) + 12 = 21 = interpolate d^2/dx1 dx2 (mixed 2nd partial derivative in x1 and x2) + 122 = 212 = 221 = interpolate d^3/dx1 dx2^2 (mixed 3rd partial derivative) + 222 = interpolate d^3/dx2^3 (3rd derivative along x2 coordinate) + 123 = 132 = 213 = 231 = 312 = 321 + = interpolate d^3/dx1 dx2 dx3 (mixed 3rd partial derivative) + + + +Pointers in Fortran +=================== + +One possible problem area with this API is that it requires creating +arrays of pointers pointing to other arrays. In C this is no problem, +but in Fortran 77 this is difficult. So, I propose adding two new Cactus +functions to make this easier for Fortran users: + + CCTK_POINTER Util_PointerTo(any Fortran variable or array) + CCTK_POINTER Util_NullPointer() + +Util_PointerTo would be #defined to %loc on those compilers which have +that extension to standard Fortran, or would be a Cactus-provided utility +routine for other cases. It's trivial to write the latter case in C so +long as the Fortran compiler actually uses call by reference; I've never +heard of a Fortran compiler doing otherwise for arrays. (And even for +Fortran scalar variables it would be very hard for a compiler to do otherwise +in light of separate compilation and 1-element arrays being allowed to be +passed to/from scalar variables.) + + + +A Simple Example +================ + +Here's a simple example, written in Fortran 77, to do quadratic interpolation +of a real and a complex local array in 3-D: + +c input arrays: + integer ni, nj, nk + parameter (ni=..., nj=..., nk=...) + CCTK_REAL real_gridfn (ni,nj,nk) + CCTK_COMPLEX complex_gridfn(ni,nj,nk) + +c interpolation coordinates + integer N_interp + parameter (N_interp = ...) + CCTK_REAL xcoord(N_interp), y_coord(N_interp), z_coord(N_interp) + +c output arrays: + CCTK_REAL real_at_xyz (N_interp) + CCTK_COMPLEX complex_at_xyz(N_interp) + + integer status, dummy + integer input_array_type_codes(2) + data input_array_type_codes /CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_COMPLEX/ + integer input_array_dims(3) + CCTK_POINTER input_arrays(2) + integer interp_coord_type_codes(3) + data interp_coord_type_codes /CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL/ + CCTK_POINTER interp_coords(3) + integer output_array_type_codes(2) + data output_array_type_codes /CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_COMPLEX/ + CCTK_POINTER output_arrays(2) + + input_array_dims(1) = ni + input_array_dims(2) = nj + input_array_dims(3) = nk + interp_coords(1) = Util_PointerTo(xcoord) + interp_coords(2) = Util_PointerTo(ycoord) + interp_coords(3) = Util_PointerTo(zcoord) + output_arrays(1) = Util_PointerTo(real_at_xyz) + output_arrays(2) = Util_PointerTo(complex_at_xyz) + + call CCTK_InterpLocalArrays + $ (status, ! return code + cctk_GH, + 3, ! number of dimensions + operator_handle, coord_system_handle, + Util_TableCreateFromString("order=2"), + N_interp, + interp_coord_type_codes, interp_coords, + 2, ! number of input arrays + input_array_type_codes, input_array_dims, input_arrays, + 2, ! number of output arrays + output_array_type_codes, output_arrays) + + if (status .lt. 0) then + call CCTK_WARN(status, "Error return from interpolator!") + call CCTK_Exit(dummy, GH, status) + end if + + + +A More Complicated Example +========================== + +Here's a more complicated example, written in C++. (I'm really only using +C++ to get cleaner initialization of the various arrays, this is still +"almost C".) This example is a simplified form of what I will be doing +in my new apparent horizon finder: + +// +// input grid functions (12 of them, all 3-D CCTK_REAL): +// gxx, gxy, gxz, gyy, gyz, gzz, +// Kxx, Kxy, Kxz, Kyy, Kyz, Kzz +// +// interpolation coordinates: +// xcoord, ycoord, zcoord (all CCTK_REAL[N_interp_points]) +// +// we want to interpolate the gij and Kij, and also interpolate all the +// first derivatives of the gij, so the output arrays are +// (30 of them, all CCTK_REAL[N_interp_points]) +// I_gxx, dx_gxx, dy_gxx, dz_gxx, +// I_gxy, dx_gxy, dy_gxy, dz_gxy, +// I_gxz, dx_gxz, dy_gxz, dz_gxz, +// I_gyy, dx_gyy, dy_gyy, dz_gyy, +// I_gyz, dx_gyz, dy_gyz, dz_gyz, +// I_gzz, dx_gzz, dy_gzz, dz_gzz, +// I_Kxx, I_Kxy, I_Kxz, I_Kyy, I_Kyz, I_Kzz +// + +#define VP(x) static_cast<void *>(x) + +const int N_dims = 3; +const int interp_coord_type_codes[N_dims] + = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL }; +const void *const interp_coords[N_dims] + = { VP(xcoord), VP(ycoord), VP(zcoord) }; + +const int N_input_arrays = 12; +const int input_array_types[N_input_arrays] + = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, + CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL }; + +const int input_array_variable_indices[N_input_arrays] + = { CCTK_VarIndex("somethorn::gxx"), + CCTK_VarIndex("somethorn::gxy"), + CCTK_VarIndex("somethorn::gxz"), + CCTK_VarIndex("somethorn::gyy"), + CCTK_VarIndex("somethorn::gyz"), + CCTK_VarIndex("somethorn::gzz"), + CCTK_VarIndex("somethorn::Kxx"), + CCTK_VarIndex("somethorn::Kxy"), + CCTK_VarIndex("somethorn::Kxz"), + CCTK_VarIndex("somethorn::Kyy"), + CCTK_VarIndex("somethorn::Kyz"), + CCTK_VarIndex("somethorn::Kzz") }; + +const int N_output_arrays = 30; +int output_array_type_codes[N_output_arrays]; + for (int oi = 0 ; oi < N_output_arrays ; ++oi) + { + output_array_type_codes[oi] = CCTK_VARIABLE_REAL; + } + +void *const output_arrays[N_output_arrays] + = { + VP(I_gxx), VP(dx_gxx), VP(dy_gxx), VP(dz_gxx), + VP(I_gxy), VP(dx_gxy), VP(dy_gxy), VP(dz_gxy), + VP(I_gxz), VP(dx_gxz), VP(dy_gxz), VP(dz_gxz), + VP(I_gyy), VP(dx_gyy), VP(dy_gyy), VP(dz_gyy), + VP(I_gyz), VP(dx_gyz), VP(dy_gyz), VP(dz_gyz), + VP(I_gzz), VP(dx_gzz), VP(dy_gzz), VP(dz_gzz), + VP(I_Kxx), VP(I_Kxy), VP(I_Kxz), VP(I_Kyy), VP(I_Kyz), VP(I_Kzz) + }; + +const int operand_indices[N_output_arrays]; + = { + 0, 0, 0, 0, // gxx + 1, 1, 1, 1, // gxy + 2, 2, 2, 2, // gxz + 3, 3, 3, 3, // gyy + 4, 4, 4, 4, // gyz + 5, 5, 5, 5, // gzz + 6, 7, 8, 9, 10, 11 // Kxx-Kzz + }; + +const int opcodes[N_output_arrays] + = { + 0, 1, 2, 3, // I, dx, dy, dz + 0, 1, 2, 3, // I, dx, dy, dz + 0, 1, 2, 3, // I, dx, dy, dz + 0, 1, 2, 3, // I, dx, dy, dz + 0, 1, 2, 3, // I, dx, dy, dz + 0, 1, 2, 3, // I, dx, dy, dz + 0, 0, 0, 0, 0, 0 // all I + }; + +int param_table_handle = Util_TableCreate(UTIL_TABLE_DEFAULT); +Util_TableSet1Int(param_table_handle, + 3, "order"); +Util_TableSetInt(param_table_handle, + N_output_arrays, operand_indices, "operand_indices"); +Util_TableSetInt(param_table_handle, + N_output_arrays, opcodes, "opcodes"); + +int status = CCTK_InterpGridArrays(GH, + N_dims, + operator_handle, coord_system_handle, + param_table_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); +if (status < 0) + { + CCTK_WARN(status, "error return from CCTK_InterpGridArrays()!"); + CCTK_Exit(GH, status); /*NOTREACHED*/ + } +Util_TableDestroy(param_table_handle); |