diff options
Diffstat (limited to 'archive')
| -rw-r--r-- | archive/2d.cube.order2.smooth0.c | 534 | ||||
| -rw-r--r-- | archive/C_str.maple | 42 | ||||
| -rw-r--r-- | archive/README | 6 | ||||
| -rw-r--r-- | archive/api1 | 323 | ||||
| -rw-r--r-- | archive/api2 | 460 | ||||
| -rw-r--r-- | archive/api2.1 | 488 | ||||
| -rw-r--r-- | archive/api2.2 | 496 | ||||
| -rw-r--r-- | archive/api3 | 491 |
8 files changed, 0 insertions, 2840 deletions
diff --git a/archive/2d.cube.order2.smooth0.c b/archive/2d.cube.order2.smooth0.c deleted file mode 100644 index 34e17ca..0000000 --- a/archive/2d.cube.order2.smooth0.c +++ /dev/null @@ -1,534 +0,0 @@ - /*@@ - @file 2d.cube.order2.smooth0.c - @date 23 Oct 2001 - @author Jonathan Thornburg <jthorn@aei.mpg.de> - @desc - Generalized interpolation for 2d, hypercube-shaped molecules, - order=2, smoothing=0. For details, see the header comments - for "InterpLocalArrays.c" in this directory. - @enddesc - - @version $Id$ - @@*/ - -#include <math.h> -#include <limits.h> -#include <stdlib.h> -#include <string.h> -#include <stdio.h> - -#include "util_ErrorCodes.h" -#include "cctk.h" -#include "InterpLocalArrays.h" - -/* the rcs ID and its dummy function to use it */ -static const char *rcsid = "$Header$"; -CCTK_FILEVERSION( - CactusPUGH_LocalInterp_GeneralizedPolynomial_2d_cube_order2_smooth0_c - ) - -/******************************************************************************/ - -/*@@ - @routine LocalInterp_ILA_2d_cube_ord2_sm0 - @date 23 Oct 2001 - @author Jonathan Thornburg <jthorn@aei.mpg.de> - @desc - This function does generalized interpolation of one or more - 2d arrays to arbitrary points. For details, see the header - comments for InterpLocalArrays() (in "InterpLocalArrays.c" - in this same directory). - - This function's arguments are all a subset of those of - InterpLocalArrays() ; the only difference is that this function - takes all its arguments explicitly, whereas InputLocalArrays() - takes some of them indirectly via a key/value parameter table. - - @returntype int - @returndesc This function's return result is the same as that of - InterpLocalArrays(): - 0 ==> successful interpolation - -1 ==> in case of any errors - @endreturndesc - - @@*/ -int LocalInterp_ILA_2d_cube_ord2_sm0 - (int param_table_handle, - const CCTK_REAL coord_system_origin[], - const CCTK_REAL grid_spacing[], - int N_interp_points, - const CCTK_INT interp_coord_type_codes[], - const void *const interp_coords[], - int N_input_arrays, - const CCTK_INT input_array_offsets[], - const CCTK_INT input_array_strides[], - const CCTK_INT input_array_min_subscripts[], - const CCTK_INT input_array_max_subscripts[], - const CCTK_INT input_array_type_codes[], - const void *const input_arrays[], - int N_output_arrays, - const CCTK_INT output_array_type_codes[], - void *const output_arrays[], - const CCTK_INT operand_indices[], const CCTK_INT opcodes[]) -{ -/* - * Implementation notes: - * - * The basic outline of this function is as follows: - * - * compute "which derivatives are wanted" flags - * for each interpolation point - * { - * declare all the coefficients - * declare all the data-values variables - * ***fetch*** interpolation point coordinates - * compute coefficients for all derivatives which are wanted - * for each output array - * { - * int part; - * for (part = 0 ; part <= 1 ; ++part) - * { - * if (this output array is computed - * using a different input array - * than the previous one || part != 0) - * then ***fetch*** the input array values - * into local variables - * { - * fp result; - * switch (opcode) - * { - * case 0: - * result = evaluate the interpolant - * break; - * case 1: - * result = evaluate the interpolant - * break; - * case ... - * } - * ***store*** result in output array - * bool complex_flag = is datatype complex? - * if (! complex_flag) - * then break; - * } - * } - * } - * } - * - * Here "***fetch***" and "***store***" are all actually switches on - * the various array datatypes. For complex datatypes they offset the - * 1D array position by part to handle real/imaginary components of - * the complex values. - * - * At present we do all floating-point computations in type "fp" - * (typically a typedef for CCTK_REAL), so arrays of higher precision - * than this will incur extra rounding errors. In practice these should - * be negligible compared to the "truncation" interpolation errors. - */ - -/* - * Naming conventions: - * input, output = 0-origin indices each selecting an input/output array - * point = 0-origin index selecting an interpolation point - */ - -/* - * these are compile-time constants here; InterpLocalArrays() decoded - * them and called us (as opposed to another function) based in part - * on these values - */ -#define N_DIMS 2 -#define MOLECULE_SIZE 3 - -/* layout of axes in N_dims-element arrays */ -#define X_AXIS 0 -#define Y_AXIS 1 - -/* input array size, strides, and subscripting computation */ -const int stride_i = input_array_strides[X_AXIS]; -const int stride_j = input_array_strides[Y_AXIS]; -#define SUB2(i,j) (i*stride_i + j*stride_j) - -/* macros used by machine-generated interpolation coefficient expressions */ -#define RATIONAL(num,den) (num/den) - -/* - * compute flags specifying which derivatives are wanted - */ -bool want_I = false; -bool want_dx = false, want_dy = false; - { -int output; - for (output = 0 ; output < N_output_arrays ; ++output) - { - switch (opcodes[output]) - { - case 0: want_I = true; break; - case 1: want_dx = true; break; - case 2: want_dy = true; break; - default: - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Generalized interpolation opcode %d not supported", - opcodes[output]); /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - } - } - } - -/* - * interpolate at each point - */ - { -int point; - for (point = 0 ; point < N_interp_points ; ++point) - { - /* declare all the interpolation coefficients */ - #include "coeffs/2d.cube.order2.smooth0.I.dcl.c" - #include "coeffs/2d.cube.order2.smooth0.dx.dcl.c" - #include "coeffs/2d.cube.order2.smooth0.dy.dcl.c" - - /* declare all the data-values variables */ - #include "coeffs/2d.cube.size3.data-var.dcl.c" - - /* - * ***fetch*** interpolation point coordinates - */ - fp interp_coords_fp[N_DIMS]; - int axis; - for (axis = 0 ; axis < N_DIMS ; ++axis) - { - /* pointer to array of interp coords for this axis */ - const void *const interp_coords_ptr = interp_coords[axis]; - - switch (interp_coord_type_codes[axis]) - { -#ifdef CCTK_REAL - case CCTK_VARIABLE_REAL: - { - const CCTK_REAL *const interp_coords_ptr_real - = (const CCTK_REAL *) interp_coords_ptr; - interp_coords_fp[axis] = interp_coords_ptr_real[point]; - break; - } -#endif -#ifdef CCTK_REAL4 - case CCTK_VARIABLE_REAL4: - { - const CCTK_REAL4 *const interp_coords_ptr_real4 - = (const CCTK_REAL4 *) interp_coords_ptr; - interp_coords_fp[axis] = interp_coords_ptr_real4[point]; - break; - } -#endif -#ifdef CCTK_REAL8 - case CCTK_VARIABLE_REAL8: - { - const CCTK_REAL8 *const interp_coords_ptr_real8 - = (const CCTK_REAL8 *) interp_coords_ptr; - interp_coords_fp[axis] = interp_coords_ptr_real8[point]; - break; - } -#endif - default: - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "interp-coords datatype %d not supported", - interp_coord_type_codes[axis]); - /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - } - } - - /* - * locate interpolation molecules with respect to the grid, - * i.e. compute interp_rel_(x,y) - */ - { - fp interp_rel_x, interp_rel_y; /* (x,y) coordinates of interpolation */ - /* point relative to molecule center */ - /* (in units of the grid spacing) */ - const int center_i - = LocalInterp_molecule_posn(coord_system_origin[X_AXIS], - grid_spacing[X_AXIS], - input_array_min_subscripts[X_AXIS], - input_array_max_subscripts[X_AXIS], - MOLECULE_SIZE, - interp_coords_fp[X_AXIS], - &interp_rel_x, - (int *) NULL, (int *) NULL); - const int center_j - = LocalInterp_molecule_posn(coord_system_origin[Y_AXIS], - grid_spacing[Y_AXIS], - input_array_min_subscripts[Y_AXIS], - input_array_max_subscripts[Y_AXIS], - MOLECULE_SIZE, - interp_coords_fp[Y_AXIS], - &interp_rel_y, - (int *) NULL, (int *) NULL); - const int center_sub = SUB2(center_i, center_j); - - /* - * compute the coefficients for whichever derivatives are wanted - * using machine-generated coefficient expressions - * ... these expressions are polynomials in (x,y) - * ==> we need these names for the relative coordinates - * (copying to fresh local variables will likely also - * give better register allocation, since [xy]_rel - * had their addresses taken and so probably won't be - * register-allocated) - */ - const fp x = interp_rel_x; - const fp y = interp_rel_y; - if (want_I) - then { - #include "coeffs/2d.cube.order2.smooth0.I.coeff.c" - } - if (want_dx) - then { - #include "coeffs/2d.cube.order2.smooth0.dx.coeff.c" - } - if (want_dy) - then { - #include "coeffs/2d.cube.order2.smooth0.dy.coeff.c" - } - - /* - * compute each output array at this point - */ - { - int output; - const void *input_array_ptr = NULL; - for (output = 0 ; output < N_output_arrays ; ++output) - { - const int input = operand_indices[output]; - const int input_offset = input_array_offsets[input]; - - /* - * for each real/imag part of complex data values - * ... for real we'll break out of this loop at the bottom - * after only a single iteration; - * for complex we'll do both iterations - */ - int part; - for (part = 0 ; part <= 1 ; ++part) - { - if ( (input_arrays[input] != input_array_ptr) - || (part != 0) ) - then { - /* - * ***fetch*** the input array values - * into local variables - */ - input_array_ptr = input_arrays[input]; - switch (input_array_type_codes[input]) - { -#ifdef CCTK_REAL -case CCTK_VARIABLE_REAL: - { - const CCTK_REAL *const input_array_ptr_real - = (const CCTK_REAL *) input_array_ptr; - #undef DATA - #define DATA(i,j) \ - input_array_ptr_real[input_offset + center_sub + SUB2(i,j)] - #include "coeffs/2d.cube.size3.data-var.assign.c" - break; - } -#endif -#ifdef CCTK_REAL4 -case CCTK_VARIABLE_REAL4: - { - const CCTK_REAL4 *const input_array_ptr_real4 - = (const CCTK_REAL4 *) input_array_ptr; - #undef DATA - #define DATA(i,j) \ - input_array_ptr_real4[input_offset + center_sub + SUB2(i,j)] - #include "coeffs/2d.cube.size3.data-var.assign.c" - break; - } -#endif -#ifdef CCTK_REAL8 -case CCTK_VARIABLE_REAL8: - { - const CCTK_REAL8 *const input_array_ptr_real8 - = (const CCTK_REAL8 *) input_array_ptr; - #undef DATA - #define DATA(i,j) \ - input_array_ptr_real8[input_offset + center_sub + SUB2(i,j)] - #include "coeffs/2d.cube.size3.data-var.assign.c" - break; - } -#endif -#ifdef CCTK_COMPLEX -case CCTK_VARIABLE_COMPLEX: - { - const CCTK_COMPLEX *const input_array_ptr_complex - = (const CCTK_COMPLEX *) input_array_ptr; - #undef DATA - #define DATA(i,j) \ - ( (const CCTK_REAL *) \ - & input_array_ptr_complex[input_offset + center_sub + SUB2(i,j)] ) [part] - #include "coeffs/2d.cube.size3.data-var.assign.c" - break; - } -#endif -#ifdef CCTK_COMPLEX8 -case CCTK_VARIABLE_COMPLEX8: - { - const CCTK_COMPLEX8 *const input_array_ptr_complex8 - = (const CCTK_COMPLEX8 *) input_array_ptr; - #undef DATA - #define DATA(i,j) \ - ( (const CCTK_REAL4 *) \ - & input_array_ptr_complex8[input_offset + center_sub + SUB2(i,j)] ) [part] - #include "coeffs/2d.cube.size3.data-var.assign.c" - break; - } -#endif -#ifdef CCTK_COMPLEX16 -case CCTK_VARIABLE_COMPLEX16: - { - const CCTK_COMPLEX16 *const input_array_ptr_complex16 - = (const CCTK_COMPLEX16 *) input_array_ptr; - #undef DATA - #define DATA(i,j) \ - ( (const CCTK_REAL8 *) \ - & input_array_ptr_complex16[input_offset + center_sub + SUB2(i,j)] ) [part] - #include "coeffs/2d.cube.size3.data-var.assign.c" - break; - } -#endif -default: - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "input datatype %d not supported", - input_array_type_codes[input]); /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - } - } - - /* - * evaluate the interpolant - * as a linear combination of the data variables - */ - { - fp result; - switch (opcodes[output]) - { - case 0: - result = - #include "coeffs/2d.cube.order2.smooth0.I.eval.c" - break; - case 1: - result = - #include "coeffs/2d.cube.order2.smooth0.dx.eval.c" - break; - case 2: - result = - #include "coeffs/2d.cube.order2.smooth0.dy.eval.c" - break; - default: - CCTK_VWarn(1, __LINE__, __FILE__, - CCTK_THORNSTRING, - "opcode %d not supported", - opcodes[output]); /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; - /*** ERROR RETURN ***/ - } - - /* - * ***store*** the result in the output array - */ - switch (output_array_type_codes[output]) - { -case CCTK_VARIABLE_REAL: - { - CCTK_REAL *const output_array_ptr_real - = (CCTK_REAL *) output_arrays[output]; - output_array_ptr_real[point] = result; - break; - } -case CCTK_VARIABLE_REAL4: - { - CCTK_REAL4 *const output_array_ptr_real4 - = (CCTK_REAL4 *) output_arrays[output]; - output_array_ptr_real4[point] = result; - break; - } -case CCTK_VARIABLE_REAL8: - { - CCTK_REAL8 *const output_array_ptr_real8 - = (CCTK_REAL8 *) output_arrays[output]; - output_array_ptr_real8[point] = result; - break; - } -case CCTK_VARIABLE_COMPLEX: - { - CCTK_COMPLEX *const output_array_ptr_complex - = (CCTK_COMPLEX *) output_arrays[output]; - ((CCTK_REAL *) & output_array_ptr_complex[point]) [part] - = result; - break; - } -case CCTK_VARIABLE_COMPLEX8: - { - CCTK_COMPLEX8 *const output_array_ptr_complex8 - = (CCTK_COMPLEX8 *) output_arrays[output]; - ((CCTK_REAL4 *) & output_array_ptr_complex8[point]) [part] - = result; - break; - } -case CCTK_VARIABLE_COMPLEX16: - { - CCTK_COMPLEX16 *const output_array_ptr_complex16 - = (CCTK_COMPLEX16 *) output_arrays[output]; - ((CCTK_REAL8 *) & output_array_ptr_complex16[point]) [part] - = result; - break; - } -default: - CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, - "output datatype %d not supported", - output_array_type_codes[output]); /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ - } - - /* decode datatype: is it real or complex? */ - { - bool complex_flag; - switch (output_array_type_codes[output]) - { - case CCTK_VARIABLE_REAL: - case CCTK_VARIABLE_REAL4: - case CCTK_VARIABLE_REAL8: - complex_flag = false; - break; - case CCTK_VARIABLE_COMPLEX: - case CCTK_VARIABLE_COMPLEX8: - case CCTK_VARIABLE_COMPLEX16: - complex_flag = true; - break; - default: - CCTK_VWarn(1, __LINE__, __FILE__, - CCTK_THORNSTRING, - "output datatype %d not supported", - output_array_type_codes[output]); - /*NOTREACHED*/ - return UTIL_ERROR_BAD_INPUT; - /*** ERROR RETURN ***/ - } - - /* skip part=1 (imaginary part) for real datatypes */ - if (! complex_flag) - then break; - } - } - /* end of for (part = ...) loop */ - } - /* end of for (output = ...) loop */ - } - } - } - - /* end of for (point = ...) loop */ - } - -return 0; /*** NORMAL RETURN ***/ - } -} diff --git a/archive/C_str.maple b/archive/C_str.maple deleted file mode 100644 index e0335a6..0000000 --- a/archive/C_str.maple +++ /dev/null @@ -1,42 +0,0 @@ -# -# This function is a wrapper around codegen[C]() which returns the -# genenerated code explictly as a Maple string. -# -# Arguments: -# expr = (in) The expression for which code is to be generated. -# ... = (in) Any further arguments are taken as options to be passed -# to codegen[C]() -# -# Results: -# The function returns a maple string of C code. -# -C_str := -proc(expr::algebraic) -local tempname, str, temp; - -# name of temp file -# FIXME: should use process number to ensure uniqueness -tempname := "/tmp/C_str.tmp.c"; - -# truncate temp file to zero length -fopen(tempname, WRITE); -fclose(tempname); - -# generate the code -codegen[C](args, filename=tempname); - -# read the code back in -str := ""; - while true - do - temp := readline(tempname); - if (temp = 0) - then break; - end if; - str := cat(str, temp); - end do; -fclose(tempname); - -# strip off the leading " t0 = " -return op(2,sscanf(str, "%s = %[^;];")); -end proc; diff --git a/archive/README b/archive/README deleted file mode 100644 index 6601932..0000000 --- a/archive/README +++ /dev/null @@ -1,6 +0,0 @@ -This directory contains stuff that's no longer used, but that I -thought might be useful for future reference, so I didn't want to -throw it away. - -Compared to the CVS attic, having a separate directory like this -seems more convenient -- you can 'cat, awk, and grep' more easily. diff --git a/archive/api1 b/archive/api1 deleted file mode 100644 index 5bc5f22..0000000 --- a/archive/api1 +++ /dev/null @@ -1,323 +0,0 @@ -As I've explained in a recent message to the cactususers mailing list, -I'm working on a new "generalized interpolator" for Cactus, which will -be able to interpolate derivatives of grid functions as well as the -grid functions themselves. Corresponding to this, I (in consultation -with Thomas Radke and Tom Goodale) have also come up with a new design -for the Cactus interpolation API. - -The purpose of this message is to solicit comments and/or suggestions -about the new API. - - -Problems with the Current Interpolator API -========================================== - -I see two major problems with the current interpolator API -(CCTK_InterpGV() and CCTK_InterpLocal()): - -The current API uses variable argument lists to allow interpolating -many arrays in a single call. This has turned out to be somewhat -error-prone. (There is *no* compiler checking of variable argument lists!) -To allow better compiler error checking, we propose replacing all the -variable argument lists with explicit arrays of pointers. - -The current interpolator interface encodes all parameters for the -interpolator (eg things like the order) into the character-string -interpolator name. This is a bit awkward, but managable if there are -only a few parameters, all integers. But it's completely unworkable for -a generalized interpolator, which needs a whole bunch of extra parameters -to tell it what derivatives to take. Instead, we propose that a -generalized interpolator get a handle to one of the new key-value -tables, and that any extra parameters (such as "which derivatives to -take" info) be passed in the table. This makes it easy to add new -interpolation operators which take additional parameters (eg a rational -interpolator might have two order parameters, one for numerator and -one for denominator), without having to change the API. - - -The Generalized Interpolator API (version 1) -============================================ - -In more detail, here is my proposal for the new API: -(don't be scared by the length, it's not that complicated!) -(there are some examples below) - - int status = CCTK_GInterpGV(arguments described below) - int status = CCTK_GInterpLocal(arguments described below) - -arguments (both CCTK_GInterpGV() and CCTK_GInterpLocal()) - const cGH *GH; - int N_dims; - int operator_handle, int coord_system_handle; - int param_table_handle; /* handle to key-value table, see below */ - - int N_interp_points; - /* array of CCTK_VARIABLE_* codes giving data types */ - /* of arrays pointed to by interp_coords[] arrays */ - const int interp_coord_types[N_dims]; - /* array of pointers to arrays giving x,y,... coordinates of interp points */ - const void *const interp_coords[N_dims]; - - int N_input_arrays, N_output_arrays; - -for CCTK_GInterpGV(): - /* array of CCTK variable indices of input arrays */ - const int input_array_indices[N_input_arrays]; - -for CCTK_GInterpLocal(): - /* array of CCTK_VARIABLE_* codes giving data types of input arrays */ - const int input_array_types[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 both functions again: - /* array of CCTK_VARIABLE_* codes giving data types of output arrays */ - const int output_array_types[N_output_arrays]; - /* array of pointers to output arrays */ - void *const output_arrays[N_output_arrays]; - -The information passed in the table would depend on the interpolation -operator. For what I am currently implementing, I will register only -a single operator, "generalized polynomial interpolation", and take the -following parameters in the table: - - int order; /* mandatory */ - -Thus the simplest call can just use Util_TableCreateFromString("order=3") -for a 3rd-order interpolation, for example. - -All the remaining parameters in the table are optional; if they're -omitted defaults will be supplied: - - const char *molecule_type; /* this selects one of a family of related */ - /* operators; the default (and the only one */ - /* I'm implementing right now) is to use the */ - /* usual hypercube-shaped molecules */ - - int smoothing; /* the way I'm implementing the interpolation */ - /* it's easy to also do Savitzky-Golay type */ - /* smoothing; this parameter gives how much */ - /* to enlarge the interpolation molecule for */ - /* this; the default is 0 (no smoothing) */ - -Optionally, the caller can specify a mask gridfn -for CCTK_GInterpGV(): - int mask_gridfn_index; -for CCTK_GInterpLocal(); - int mask_type; /* one of the CCTK_VARIABLE_* codes */ - const void *const mask_array; -and the range of values for the mask which correspond to valid points: - CCTK_REAL mask_min, mask_max; /* valid <--> closed interval [min,max] */ - -And last but not least, optionally, the caller can specify taking -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 -(we'll put #defines for these in one of the cctk header files). -Note that the array operand_indices[] doesn't directly name the inputs, -but rather gives the indices 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. - - - -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 a new Cactus -function pointer_to() to make this easier for Fortran users: - - CCTK_POINTER pointer_to(any array in your code) - -This 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, as -almost all Fortran compilers do anyway for arrays. - - - -A Simple Example -================ - -Here's a simple example, written in Fortran 77, interpolating a real -and a complex grid function in 3D: - -c -c input grid functions: -c real_gridfn (CCTK_REAL) -c complex_gridfn (CCTK_COMPLEX) -c -c interpolation coordinates -c xcoord, ycoord, zcoord (CCTK_REAL arrays of size N) -c -c output arrays: -c real_at_xyz (CCTK_REAL array of size N) -c complex_at_xyz (CCTK_COMPLEX array of size N) -c - integer status - CCTK_POINTER input_gridfn_indices(2) - integer interp_coord_types(3) - CCTK_POINTER interp_coords(3) - integer output_array_types(2) - CCTK_POINTER output_arrays(2) - - call CCTK_VarIndex(input_gridfn_indices(1), "mythorn::real_gridfn") - call CCTK_VarIndex(input_gridfn_indices(2), "mythorn::complex_gridfn") - -c could also initialize this with a DATA statement - interp_coord_types(1) = CCTK_VARIABLE_REAL - interp_coord_types(2) = CCTK_VARIABLE_REAL - interp_coord_types(3) = CCTK_VARIABLE_REAL - interp_coords(1) = pointer_to(xcoord) - interp_coords(2) = pointer_to(ycoord) - interp_coords(3) = pointer_to(zcoord) -c could also initialize this with a DATA statement - output_array_types(1) = CCTK_VARIABLE_REAL - output_array_types(2) = CCTK_VARIABLE_COMPLEX - output_arrays(1) = pointer_to(real_at_xyz) - output_arrays(2) = pointer_to(complex_at_xyz) - call CCTK_InterpGV(status, - cctk_GH, - 3, # number of dimensions - operator_handle, coord_system_handle, - Util_TableCreateFromString("order=3"), - N, interp_coord_types, interp_coords, - 2, 2, # number of input/output arrays - input_gridfn_indices, - output_array_types, output_arrays); - if (status .lt. 0) then -c help, an error occured! - 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 3D CCTK_REAL): -// gxx, gxy, gxz, gyy, gyz, gzz, -// Kxx, Kxy, Kxz, Kyy, Kyz, Kzz -// -// interpolation coordinates: -// xcoord, ycoord, zcoord (CCTK_REAL arrays of size N) -// -// output arrays (30 of them, all 3D CCTK_REAL) -// (interpolate the gij and Kij, also interpolate all first derivs of the gij) -// 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, -// - -const int gxx_index = CCTK_VarIndex("somethorn::gxx"); -const int gxy_index = CCTK_VarIndex("somethorn::gxy"); -const int gxz_index = CCTK_VarIndex("somethorn::gxz"); -const int gyy_index = CCTK_VarIndex("somethorn::gyy"); -const int gyz_index = CCTK_VarIndex("somethorn::gyz"); -const int gzz_index = CCTK_VarIndex("somethorn::gzz"); -const int Kxx_index = CCTK_VarIndex("somethorn::Kxx"); -const int Kxy_index = CCTK_VarIndex("somethorn::Kxy"); -const int Kxz_index = CCTK_VarIndex("somethorn::Kxz"); -const int Kyy_index = CCTK_VarIndex("somethorn::Kyy"); -const int Kyz_index = CCTK_VarIndex("somethorn::Kyz"); -const int Kzz_index = CCTK_VarIndex("somethorn::Kzz"); - -const int N_dims = 3; -int interp_coord_types[N_dims] - = { CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL }; -const void *const interp_coords[N_dims] - = { (const void *) &xcoord, - (const void *) &ycoord, - (const void *) &zcoord }; - -const int N_input_arrays = 12; -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 N_output_arrays = 30; -int output_array_types[N_output_arrays]; - for (int oi = 0 ; oi < N_output_arrays ; ++oi) - { - output_array_types[oi] = CCTK_VARIABLE_REAL; - } - -#define VP(x) static_cast<void *>(x) -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), - }; - -int operand_indices[N_output_arrays]; - = { - gxx_index, gxx_index, gxx_index, gxx_index, - gxy_index, gxy_index, gxy_index, gxy_index, - gxz_index, gxz_index, gxz_index, gxz_index, - gyy_index, gyy_index, gyy_index, gyy_index, - gyz_index, gyz_index, gyz_index, gyz_index, - gzz_index, gzz_index, gzz_index, gzz_index, - Kxx_index, Kxy_index, Kxz_index, Kyy_index, Kyz_index, Kzz_index, - }; - -const int op_I = CCTK_INTERP_OPCODE_I; -const int op_dx = CCTK_INTERP_OPCODE_DX; -const int op_dy = CCTK_INTERP_OPCODE_DY; -const int op_dz = CCTK_INTERP_OPCODE_DZ; -int opcodes[N_output_arrays] - = { - op_I, op_dx, op_dy, op_dz, - op_I, op_dx, op_dy, op_dz, - op_I, op_dx, op_dy, op_dz, - op_I, op_dx, op_dy, op_dz, - op_I, op_dx, op_dy, op_dz, - op_I, op_dx, op_dy, op_dz, - op_I, op_I, op_I, op_I, op_I, op_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_GInterpGV(GH, - N_dims, - operator_handle, coord_system_handle, - param_table_handle, - N_interp_points, interp_coord_types, interp_coords, - N_input_arrays, N_output_arrays, - input_array_indices, - output_array_types, output_arrays); -if (status < 0) - { /* something bad happened */ } -Util_TableDestroy(param_table_handle); - - - diff --git a/archive/api2 b/archive/api2 deleted file mode 100644 index 422be0d..0000000 --- a/archive/api2 +++ /dev/null @@ -1,460 +0,0 @@ -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); diff --git a/archive/api2.1 b/archive/api2.1 deleted file mode 100644 index e316aba..0000000 --- a/archive/api2.1 +++ /dev/null @@ -1,488 +0,0 @@ -Author: Jonathan Thornburg <jthorn@aei.mpg.de> -Date: December 2001 - -This is version 2.1 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 -* no #defines in the header file for derivative operation codes -* int arrays changed to CCTK_INT - -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 - -(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 *****/ -for CCTK_InterpGridArrays() only: - const cGH *GH; -for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays(): - /* 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 "parameter table", a key-value */ - /* table, see below for table entries */ - - /***** arguments specifying the coordinate system *****/ -for CCTK_InterpLocalArrays(): - /* the local coordinate system is specified as a generic linear mapping */ - /* from (integer) local input array subscripts --> (global) coordinates: */ - /* coordinate = coord_system_origin[axis] + subscript*grid_spacing[axis] */ - const CCTK_REAL coord_system_origin[N_dims]; /* coords of subscript 0 */ - const CCTK_REAL grid_spacing[N_dims]; -for CCTK_InterpGridArrays(): - int coord_system_handle; /* handle to Cactus coordinate system */ - /* specifying mapping of */ - /* (integer) input array subscripts */ - /* <--> (floating point) coordinates */ - - /***** arguments specifying the interpolation points *****/ - int N_interp_points; - /* (pointer to) array of CCTK_VARIABLE_* codes giving the */ - /* data types of the arrays pointed to by interp_coords[] */ - const CCTK_INT interp_coord_type_codes[N_dims]; - /* (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; -for CCTK_InterpLocalArrays(): - /* 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]; -for CCTK_InterpGridArrays(): - /* array of CCTK variable indices of input arrays */ - const CCTK_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 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. - - -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] - -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, -but shares any input-array stride and/or min/max subscript specification: - - const CCTK_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 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 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) - -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 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_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, 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); diff --git a/archive/api2.2 b/archive/api2.2 deleted file mode 100644 index 35eca3c..0000000 --- a/archive/api2.2 +++ /dev/null @@ -1,496 +0,0 @@ -Author: Jonathan Thornburg <jthorn@aei.mpg.de> -Date: 30 December 2001 - -This is version 2.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, - time levels, etc -* no #defines in the header file for derivative operation codes -* int arrays changed to CCTK_INT - -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 - -(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 *****/ -for CCTK_InterpGridArrays() only: - const cGH *GH; -for both CCTK_InterpLocalArrays() and CCTK_InterpGridArrays(): - /* 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 "parameter table", a key-value */ - /* table, see below for table entries */ - - /***** arguments specifying the coordinate system *****/ -for CCTK_InterpLocalArrays(): - /* the local coordinate system is specified as a generic linear mapping */ - /* from (integer) local input array subscripts --> (global) coordinates: */ - /* coordinate = coord_system_origin[axis] + subscript*grid_spacing[axis] */ - const CCTK_REAL coord_system_origin[N_dims]; /* coords of subscript 0 */ - const CCTK_REAL grid_spacing[N_dims]; -for CCTK_InterpGridArrays(): - int coord_system_handle; /* handle to Cactus coordinate system */ - /* specifying mapping of */ - /* (integer) input array subscripts */ - /* <--> (floating point) coordinates */ - - /***** arguments specifying the interpolation points *****/ - int N_interp_points; - /* (pointer to) array of CCTK_VARIABLE_* codes giving the */ - /* data types of the arrays pointed to by interp_coords[] */ - const CCTK_INT interp_coord_type_codes[N_dims]; - /* (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; -for CCTK_InterpLocalArrays(): - /* 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]; -for CCTK_InterpGridArrays(): - /* array of CCTK variable indices of input arrays */ - const CCTK_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 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 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) - -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 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_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, 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); 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); |