diff options
Diffstat (limited to 'archive/api3')
| -rw-r--r-- | archive/api3 | 491 |
1 files changed, 0 insertions, 491 deletions
diff --git a/archive/api3 b/archive/api3 deleted file mode 100644 index 9aeedd5..0000000 --- a/archive/api3 +++ /dev/null @@ -1,491 +0,0 @@ -Author: Jonathan Thornburg <jthorn@aei.mpg.de> -Date: 21 January 2002 - -This is version 3.0 of my proposal for the new interpolator API. The -main changes since previous versions are -* This only describes processor-local interpolation -- global (or grid) - interpolation is trickier, and we're deferring it till after the - processor-local code is finished -* The API has split again to separate out {uniform, nonuniform, irregular} - grids. This proposal only addresses the first two; irregular is hard - and needs lots of work all over Cactus before we can support it. -* Since noone seems to care about it, I've dropped support for different - coordinates being of different types (i.e. x is CCTK_REAL4 but y is - CCTK_REAL16); now all the coordinates must be the same type. -* Ditto for the interpolation coordinates. - -Don't be scared by the length, for most uses it's not that complicated! -There are some examples below... - - - -Synopsis -======== - - int status = CCTK_InterpLocalUniform(arguments described below) - int status = CCTK_InterpLocalNonUniform(arguments described below) - -return is 0 for ok, various -ve values for error codes - -(N.b. the flesh APIs to register interpolation operators will also need -their C function prototypes changed to match the changes here.) - - - -Function Arguments -================== - - /***** misc arguments *****/ - /* 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 param_table_handle; /* handle to a key-value table used to pass */ - /* additional parameters to the interpolator */ - - /***** arguments specifying the local coordinate system *****/ -for CCTK_InterpLocalUniform() - /* the local coordinate system is specified as a generic linear mapping */ - /* from (integer) local input array subscripts --> (global) coordinates: */ - /* coordinate = origin[axis] + subscript*delta[axis] */ - const CCTK_REAL origin[N_dims]; /* coords of subscript 0 */ - const CCTK_REAL delta[N_dims]; -for CCTK_InterpLocalNonuniform() - /* the local coordinate system is specified by user-supplied arrays: */ - /* coordinate = array[subscript] (separate array for each axis) */ - /* FIXME: what if subscript=0 is out-of-range??? */ - const void *const coords[N_dims]; /* coords of subscript 0 */ - - /***** arguments specifying the interpolation points *****/ - int N_interp_points; - int interp_coords_type_code; /* CCTK_VARIABLE_* type code for */ - /* interp_coords arrays */ - /* (pointer to) 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; - /* array of input array dimensions (common to all input arrays) */ - const CCTK_INT input_array_dims[N_dims]; - /* array of CCTK_VARIABLE_* codes giving data types of input arrays */ - const CCTK_INT input_array_type_codes[N_input_arrays]; - /* array of pointers to input arrays */ - const void *const input_arrays[N_input_arrays]; - - /***** arguments specifying the outputs (the interpolation results) *****/ - int N_output_arrays; - /* array of CCTK_VARIABLE_* codes giving data types of output arrays */ - const CCTK_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. - - -Array Addressing/Subscripting 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 CCTK_INT input_array_offsets[N_input_arrays]; - /* next 3 table entries are shared by all input arrays */ - const CCTK_INT input_array_strides [N_dims]; - const CCTK_INT input_array_min_subscripts[N_dims]; - const CCTK_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] - -For CCTK_InterpGridArrays(), by default the input (grid) arrays are at -the "current" Cactus time level (level 0). This may be changed with the -table entry - const CCTK_INT input_array_time_levels[N_input_arrays]; - -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 CCTK_INT interp_coords_offsets[N_dims]; - const CCTK_INT output_array_offsets[N_output_arrays]; - /* next 4 table entries are shared by all interp coords and output arrays */ - const CCTK_INT interp_point_N_dims; - const CCTK_INT interp_point_strides [interp_point_N_dims]; - const CCTK_INT interp_point_min_subscripts[interp_point_N_dims]; - const CCTK_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 CCTK_INT mask_type_code; /* one of the CCTK_VARIABLE_* codes */ - const void *const mask_array; /* same dimensions/indexing as input arrays */ -for CCTK_InterpGridArrays(): - const CCTK_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 -and/or time level, but shares any input-array stride and/or min/max -subscript specification: - - const CCTK_INT mask_offset; - const CCTK_INT mask_time_level; - - -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 CCTK_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 CCTK_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 CCTK_INT operand_indices[N_output_arrays]; - const CCTK_INT operation_codes[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) - -After discussion with Tom Goodale, we have decided *not* to put #defines -for the operation codes in any of the interpolator header files -- the -operation codes are specific to this particular interpolation operator, -not common to all operators, so they don't belong in the overall -common-to-all header files. - - - -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 - CCTK_INT input_array_type_codes(2) - data input_array_type_codes /CCTK_VARIABLE_REAL, - $ CCTK_VARIABLE_COMPLEX/ - CCTK_INT input_array_dims(3) - CCTK_POINTER input_arrays(2) - CCTK_INT interp_coord_type_codes(3) - data interp_coord_type_codes /CCTK_VARIABLE_REAL, - $ CCTK_VARIABLE_REAL, - $ CCTK_VARIABLE_REAL/ - CCTK_POINTER interp_coords(3) - CCTK_INT 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 - 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, Util_NullPointer(), 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 CCTK_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 CCTK_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 CCTK_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; -CCTK_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 CCTK_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 CCTK_INT operation_codes[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_TableSetInt(param_table_handle, 3, "order"); -Util_TableSetIntArray(param_table_handle, - N_output_arrays, operand_indices, - "operand_indices"); -Util_TableSetIntArray(param_table_handle, - N_output_arrays, operation_codes, - "operation_codes"); - -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); |