/*@@ @file template.c @date 23.Oct.2001 @author Jonathan Thornburg @desc This file is a template for generalized interpolation for 1-, 2-, or 3-d molecules. It's used by defining various C preprocessor macros, then #including this file. Each such inclusion defines a single interpolation function (plus other supporting functions local to this file), which does generalized interpolation for a single molecule, or more accurately for a single molecule for each possible operation_codes[] value. @enddesc @hdate 29.Aug.2002 @hauthor Jonathan Thornburg @hdesc Restructure code and break it up into smaller functions to reduce excessive compilation cpu/memory with big 3-D molecules. @endhdesc @hdate 28.Jan.2003 @hauthor Jonathan Thornburg @hdesc Support parameter-table entries @var boundary_off_centering_tolerance and @var boundary_extrapolation_tolerance. @endhdesc @desc The following header files must be #included before #including this file: "util_ErrorCodes.h" "cctk.h" "InterpLocalUniform.h" The following preprocessor macros must be defined before #including this file (note the macros which are file names must all include the proper quotes for a #include in their expansion, e.g. #define DATA_VAR_DCL_FILE_NAME \ "coeffs/2d.cube.size3/data-var.dcl.c" @enddesc @var FUNCTION_NAME @vdesc The name of the interpolation function, e.g. #define FUNCTION_NAME AEILocalInterp_U_LagTP_2cube_20 @endvar @var N_DIMS @vdesc The number of dimensions in which to interpolate, e.g. #define N_DIMS 2 The present implementation restricts this to 1, 2, or 3, but this could easily be changed if needed. Note that MAX_N_DIMS (defined in "InterpLocalUniform.c") is a compile-time upper bound for N_DIMS, useful for sizing arrays etc. @endvar @var XYZ @desc A comma-separated list of the (x,y,z) coordinates in the interpolation, e.g. #define XYZ x,y @endvar @var FP_XYZ @desc A comma-separated list of function-prototype declarations of the (x,y,z) coordinates in the interpolation, e.g. #define FP_XYZ fp x, fp y @endvar @var STRIDE_IJK @desc A comma-separated list of the (i,j,k) array strides for the interpolation, e.g. #define STRIDE_IJK stride_i, stride_j @endvar @var JACOBIAN_MIJK_STRIDE @desc A comma-separated list of the (i,j,k) strides for the Jacobian, e.g. #define JACOBIAN_MIJK_STRIDE \ Jacobian_mi_stride, Jacobian_mj_stride @endvar @var MOLECULE_MIN_M @vdesc The minimum m coordinate in the molecule, e.g. #define MOLECULE_MIN_M -1 The present implementation takes this to be the same in each dimension, but this could easily be changed if needed. @endvar @var MOLECULE_MAX_M @vdesc The maximum m coordinate in the molecule, e.g. #define MOLECULE_MAX_M 1 The present implementation takes this to be the same in each dimension, but this could easily be changed if needed. @endvar @var MOLECULE_SIZE @vdesc The diameter of (number of points in) the molecules to be used, e.g. #define MOLECULE_SIZE 3 The present implementation takes this to be the same in each dimension, but this could easily be changed if needed. @endvar @var HAVE_OP_{I,DX,DY,DXX,DXY,DYY,...} @vdesc Each of these symbols should be defined or not defined according as if the corresponding derivative operator is to be supported or not supported by this function, e.g. #define HAVE_OP_I #define HAVE_OP_DX #define HAVE_OP_DY #define HAVE_OP_DXX #define HAVE_OP_DXY #define HAVE_OP_DYY if we support all 1st and 2nd derivatives in 2-D. @endvar @var DATA_STRUCT @vdesc The name of a C struct containing all the data variables, e.g. #define DATA_STRUCT data_struct_2d_cube_size3 @endvar @var COEFF_STRUCT @vdesc The name of a C struct containing all the molecule coefficients, e.g. #define COEFFS_STRUCT coeffs_struct_2d_cube_size3 @endvar @var LOAD_DATA_{REAL,REAL{4,8,16},COMPLEX,COMPLEX{8,16,32}} @vdesc The name of a C function to load data from the input arrays into the DATA_STRUCT structure. Typically these will be chosem from among the functions defined in "common/load-template.[ch]", e.g. #define LOAD_DATA_REAL AEILocalInterp_load_2dcube3_r #define LOAD_DATA_REAL4 AEILocalInterp_load_2dcube3_r4 #define LOAD_DATA_REAL8 AEILocalInterp_load_2dcube3_r8 #define LOAD_DATA_REAL16 AEILocalInterp_load_2dcube3_r16 #define LOAD_DATA_COMPLEX AEILocalInterp_load_2dcube3_c #define LOAD_DATA_COMPLEX8 AEILocalInterp_load_2dcube3_c8 #define LOAD_DATA_COMPLEX16 AEILocalInterp_load_2dcube3_c16 #define LOAD_DATA_COMPLEX32 AEILocalInterp_load_2dcube3_c32 @endvar @var EVALUATE_MOLECULE @vdesc The name of a C function to evaluate the dot product of the DATA_STRUCT and COEFF_STRUCT structures. Typically these will be chosem from among the functions defined in "common/evaluate.[ch]", e.g. #define EVALUATE_MOLECULE LocalInterp_eval_2d_cube6 @endvar @var STORE_COEFFS @vdesc The name of a C function to store a COEDFF_STRUCT structure of the DATA_STRUCT and COEFF_STRUCT structures. Typically this will be chosem from among the functions defined in "common/store.[ch]", e.g. #define STORE_COEFFS LocalInterp_store_2d_cube6 @endvar @var COEFFS_{I,DX,DY,DXX,DXY,DYY,...}_COMPUTE_FILE_NAME @vdesc Each of these macros should be the name of a file (presumably machine-generated) containing a sequence of C assignment statements to compute all the coefficients in a corresponding coeffs_* structure as polynomials in the variables (x,y,z), e.g. (for 2D size-3 molecules, dx operator) fp t36; fp t42; fp t41; fp t40; fp t39; fp t37; t36 = RATIONAL(1.0,3.0)*x; t42 = RATIONAL(1.0,4.0)*y+t36; t41 = t36+RATIONAL(-1.0,4.0)*y; t40 = RATIONAL(1.0,6.0); t39 = RATIONAL(-1.0,6.0); t37 = RATIONAL(-2.0,3.0)*x; coeffs_dx->coeff_m1_m1 = t39+t42; coeffs_dx->coeff_0_m1 = t37; coeffs_dx->coeff_p1_m1 = t40+t41; coeffs_dx->coeff_m1_0 = t36+t39; coeffs_dx->coeff_0_0 = t37; coeffs_dx->coeff_p1_0 = t36+t40; coeffs_dx->coeff_m1_p1 = t39+t41; coeffs_dx->coeff_0_p1 = t37; coeffs_dx->coeff_p1_p1 = t40+t42; As illustrated, the code may use the macro RATIONAL (defined later in this file) to represent rational-number coefficients, and it may also declare temporary variables as needed. Note that this is the only included code which depends on the actual interpolation scheme used; all the rest just depends on the interpolation dimension and molecule family and size. @endvar @version $Header$ @@*/ /******************************************************************************/ /* * ***** prototypes for functions local to this file ***** */ #ifdef HAVE_OP_I static void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *coeffs_I); #endif #ifdef HAVE_OP_DX static void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dx); #endif #ifdef HAVE_OP_DY static void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dy); #endif #ifdef HAVE_OP_DZ static void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dz); #endif #ifdef HAVE_OP_DXX static void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxx); #endif #ifdef HAVE_OP_DXY static void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxy); #endif #ifdef HAVE_OP_DXZ static void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxz); #endif #ifdef HAVE_OP_DYY static void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyy); #endif #ifdef HAVE_OP_DYZ static void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyz); #endif #ifdef HAVE_OP_DZZ static void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dzz); #endif /******************************************************************************/ /*@@ @routine FUNCTION_NAME @date 23.Oct.2001 @author Jonathan Thornburg @desc This function does generalized interpolation of one or more 2d arrays to arbitrary points. For details, see the header comments for InterpLocalUniform() (in "InterpLocalUniform.c" in this same directory). This function's arguments are mostly all a subset of those of InterpLocalUniform() ; the difference is that this function takes all its arguments explicitly, whereas InputLocalArrays() takes some of them indirectly via a key/value parameter table. Here we document only the "new" explicit arguments, and these only briefly; see the InterpLocalUniform() documentation and/or the thorn guide for this thorn for further details. If you change the arguments for this function, note that you must also change the prototype in "template.h". @var log_fp @vdesc This is NULL for no logging, or a valid stdio stream pointer to enable logging of the grid and interpolation coordinates. @vtype FILE* @vio in @var error_info @vdesc If error_info->per_point_status is non-NULL, then we store per-point interpolator status for each point in error_info->per_point_status[pt] . We set error_info->error_count to the number of points in error. We set the remaining elements of error_info to a detailed description of the first point in error (if any). @vtype struct error_info* error_info; @vio out @var molecule_structure_flags @vdesc If this pointer is non-NULL, we store flags describing the interpolation molecule's structure in the pointed-to structure. @vtype struct molecule_structure_flags* molecule_structure_flags; @vio out @var molecule_min_max_m_info @vdesc If this pointer is non-NULL, we store the interpolation molecule's min/max m in the pointed-to structure. @vtype struct molecule_min_max_m_info* molecule_min_max_m_info; @vio out @var molecule_positions @vdesc If this pointer is non-NULL, we store the interpolation molecule's positions in the (caller-supplied) arrays pointed to by this pointer. @vtype CCTK_INT* const molecule_positions[]; @vio out @var Jacobian_info @vdesc If this pointer is non-NULL, we store the Jacobian of the interpolation in the arrays (and in the manner) pointed to by this structure. @vtype struct Jacobian_info* Jacobian_info; @vio out @returntype int @returndesc This function's return results are a subset of those of InterpLocalUniform(): 0 success UTIL_ERROR_BAD_INPUT one of the input arguments is invalid CCTK_ERROR_INTERP_POINT_OUTSIDE interpolation point is out of range @endreturndesc @@*/ int FUNCTION_NAME(/***** coordinate system *****/ const CCTK_REAL coord_origin[], const CCTK_REAL coord_delta[], /***** interpolation points *****/ int N_interp_points, int interp_coords_type_code, const void* const interp_coords[], const CCTK_INT N_boundary_points_to_omit[], const CCTK_REAL boundary_off_centering_tolerance[], const CCTK_REAL boundary_extrapolation_tolerance[], /***** input arrays *****/ 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[], /***** output arrays *****/ int N_output_arrays, const CCTK_INT output_array_type_codes[], void* const output_arrays[], /***** operation info *****/ const CCTK_INT operand_indices[], const CCTK_INT operation_codes[], /***** debugging *****/ int debug, FILE* log_fp, /***** other return results *****/ struct error_info* error_info, struct molecule_structure_flags* molecule_structure_flags, struct molecule_min_max_m_info* molecule_min_max_m_info, CCTK_INT* const molecule_positions[], struct Jacobian_info* Jacobian_info) { /* * ***** Naming conventions ***** * * in, out = 0-origin indices each selecting an input/output array * pt = 0-origin index selecting an interpolation point * part = 0-origin index selecting real/imaginary part of a complex number */ /* * ***** Implementation notes: ***** * * The basic outline of this function is as follows: * * store molecule structure flags * if (querying molecule min/max m) * then store molecule min/max m * compute "which derivatives are wanted" flags * precompute 1/dx factors * set up input array strides * set up Jacobian strides (or dummy values if no Jacobian info) * for (int pt = 0 ; pt < N_interp_point ; ++pt) * { * struct DATA_STRUCT data; * struct COEFFS_STRUCT coeffs_I, coeffs_dx, ..., coeffs_dzz; * for (int axis = 0 ; axis < N_dims ; ++axis) * { * switch (interp_coords_type_codes) * { * case CCTK_VARIABLE_REAL: * load interp coords for this axis from * const CCTK_REAL *const interp_coords[] * break; * ... * } * } * compute position of interpolation molecules * with respect to the grid into the XYZ variables * and store per-point status if requested * if (this point is outside the grid) * then continue; * if (querying molecule positions) * then store this molecule position * compute molecule center 1-d position in input arrays * * /# compute interpolation coefficients at this point #/ * if (want_I) * then compute_coeffs_I(XYZ, &coeffs_I); * if (want_dx) * then compute_coeffs_dx(XYZ, &coeffs_dx); * ... * if (want_dzz) * then compute_coeffs_dzz(XYZ, &coeffs_dx); * * for (int out = 0 ; out < N_output_arrays ; ++out) * { * const int in = operand_indices[out]; * ***decode*** the output array datatype * to determine whether it's real or complex * const int N_output_parts = output array is complex ? 2 : 1; * for (int part = 0 ; part < N_output_parts ; ++part) * { * if ( (input_arrays[in] != NULL) * && ( (input_arrays[in] != value at last load) * || (part != value at last load) ) ) * then { * save input_arrays[in] and part for * "previous value" test above * ***decode*** the input array datatype * to determine whether it's real or complex * const int N_input_parts * = input array is complex ? 2 : 1; * if (N_input_parts != N_output_parts) * then error(...) * switch (input_array_type_codes[in]) * { * case CCTK_VARIABLE_REAL: * LOAD_DATA_REAL(input_array_ptr, * STRIDE_IJK, * &data); * break; * ... * case CCTK_VARIABLE_COMPLEX: * LOAD_DATA_COMPLEX(input_array_ptr, * STRIDE_IJK, part, * &data); * break; * ... * } * } * * if (output_arrays[out] != NULL) * then { * /# interpolate at this point #/ * fp result; * switch (operation_codes[out]) * { * case 0: * result = EVALUATE_MOLECULE(&coeffs_I, * &data); * break; * case 1: * result = inverse_dx * * EVALUATE_MOLECULE(&coeffs_dx, * &data); * break; * ... * case 33: * result = inverse_dz * inverse_dz * * EVALUATE_MOLECULE(&coeffs_dzz, * &data); * break; * } * * /# store result in output array #/ * switch (output_array_type_codes[out]) * { * case CCTK_VARIABLE_REAL: * store result in * CCTK_REAL *const output_arrays[] * break; * ... * case CCTK_VARIABLE_COMPLEX: * store result in * CCTK_COMPLEX *const output_arrays[][part] * break; * ... * } * } * if (querying Jacobian && (Jacobian_pointer[out] != NULL)) * then { * /# store Jacobian coefficients at this point #/ * CCTK_REAL *const Jacobian_ptr * = Jacobian_info->Jacobian_pointer[out] * + Jacobian_info->Jacobian_offset[out]; * switch (operation_code) * { * case 0: * STORE_COEFFS(1.0, &coeffs_I, * Jacobian_ptr, JACOBIAN_MIJK_STRIDE, * pt, part, * break; * case 1: * STORE_COEFFS(inverse_dx, &coeffs_dx, * Jacobian_ptr, JACOBIAN_MIJK_STRIDE, * pt, part, * break; * ... * case 33: * STORE_COEFFS(inverse_dz*inverse_dz, &coeffs_dzz, * Jacobian_ptr, JACOBIAN_MIJK_STRIDE, * pt, part, * break; * } * } * } * } * } * * For complex datatypes we pointer-alias the N-dimensional input array * to a 2-dimensional array where the 1st axis is the 1-D subscript * corresponding to the N input axes, and the 2nd axes has 2 elements * subscripted by [part] for the (real,imaginary) components of the * complex values. * * 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 * interpolation errors. */ /* layout of axes in N_dims-element arrays */ #define X_AXIS 0 #define Y_AXIS 1 #define Z_AXIS 2 #if (N_DIMS > 3) #error "N_DIMS may not be > 3!" #endif /* layout of axes and min/max ends in out_of_range_tolerance[] array */ #define X_AXIS_MIN 0 #define X_AXIS_MAX 1 #define Y_AXIS_MIN 2 #define Y_AXIS_MAX 3 #define Z_AXIS_MIN 4 #define Z_AXIS_MAX 5 #if (N_DIMS > 3) #error "N_DIMS may not be > 3!" #endif /* basic sanity check on molecule size */ #define MOLECULE_M_COUNT (MOLECULE_MAX_M - MOLECULE_MIN_M + 1) #if (MOLECULE_SIZE != MOLECULE_M_COUNT) #error "MOLECULE_SIZE inconsistent with MOLECULE_{MIN,MAX}_M!" #endif /* macros used by machine-generated interpolation coefficient expressions */ /* * FIXME: right now this is used as (eg) RATIONAL(1.0,2.0); * it might be cleaner if it were RATIONAL(1,2) with the * preprocessor ## operator used to glue on the .0 * (I _think_ this is portable, but is it really?) */ #define RATIONAL(num,den) (num/den) /* * store molecule structure flags, molecule min/max m (if requested) */ if (molecule_structure_flags != NULL) then { molecule_structure_flags->MSS_is_fn_of_interp_coords = 0; molecule_structure_flags->MSS_is_fn_of_which_operation = 0; molecule_structure_flags->MSS_is_fn_of_input_array_values = 0; molecule_structure_flags->Jacobian_is_fn_of_input_array_values = 0; } if (molecule_min_max_m_info != NULL) then { int axis; for (axis = 0 ; axis < N_DIMS ; ++axis) { molecule_min_max_m_info->molecule_min_m[axis] = MOLECULE_MIN_M; molecule_min_max_m_info->molecule_max_m[axis] = MOLECULE_MAX_M; } } /* * compute "which derivatives are wanted" flags */ { #ifdef HAVE_OP_I bool want_I = false; #endif #ifdef HAVE_OP_DX bool want_dx = false; #endif #ifdef HAVE_OP_DY bool want_dy = false; #endif #ifdef HAVE_OP_DZ bool want_dz = false; #endif #ifdef HAVE_OP_DXX bool want_dxx = false; #endif #ifdef HAVE_OP_DXY bool want_dxy = false; #endif #ifdef HAVE_OP_DXZ bool want_dxz = false; #endif #ifdef HAVE_OP_DYY bool want_dyy = false; #endif #ifdef HAVE_OP_DYZ bool want_dyz = false; #endif #ifdef HAVE_OP_DZZ bool want_dzz = false; #endif { int out; for (out = 0 ; out < N_output_arrays ; ++out) { switch (operation_codes[out]) { #ifdef HAVE_OP_I case 0: want_I = true; break; #endif #ifdef HAVE_OP_DX case 1: want_dx = true; break; #endif #ifdef HAVE_OP_DY case 2: want_dy = true; break; #endif #ifdef HAVE_OP_DZ case 3: want_dz = true; break; #endif #ifdef HAVE_OP_DXX case 11: want_dxx = true; break; #endif #ifdef HAVE_OP_DXY case 12: case 21: want_dxy = true; break; #endif #ifdef HAVE_OP_DXZ case 13: case 31: want_dxz = true; break; #endif #ifdef HAVE_OP_DYY case 22: want_dyy = true; break; #endif #ifdef HAVE_OP_DYZ case 23: case 32: want_dyz = true; break; #endif #ifdef HAVE_OP_DZZ case 33: want_dzz = true; break; #endif default: CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): generalized interpolation operation_code %d\n" " is not supported!\n" " (0-origin) output #out=%d" , (int) operation_codes[out], out); /*NOTREACHED*/ return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } } } /* * save origin/delta variables, precompute 1/delta factors * ... in theory we could compute only those factors we're going to use, * but it's not worth the trouble, so we just compute them all */ { #if N_DIMS >= 1 #if defined(HAVE_OP_DX) \ || defined(HAVE_OP_DXY) || defined(HAVE_OP_DXZ) \ || defined(HAVE_OP_DXX) const fp inverse_dx = 1.0 / coord_delta[X_AXIS]; #endif #endif #if N_DIMS >= 2 #if defined(HAVE_OP_DY) \ || defined(HAVE_OP_DXY) || defined(HAVE_OP_DYZ) \ || defined(HAVE_OP_DYY) const fp inverse_dy = 1.0 / coord_delta[Y_AXIS]; #endif #endif #if N_DIMS >= 3 #if defined(HAVE_OP_DZ) \ || defined(HAVE_OP_DXZ) || defined(HAVE_OP_DYZ) \ || defined(HAVE_OP_DZZ) const fp inverse_dz = 1.0 / coord_delta[Z_AXIS]; #endif #endif #if (N_DIMS > 3) #error "N_DIMS may not be > 3!" #endif /* * set up input array strides */ #if (N_DIMS >= 1) const int stride_i = input_array_strides[X_AXIS]; #endif #if (N_DIMS >= 2) const int stride_j = input_array_strides[Y_AXIS]; #endif #if (N_DIMS >= 3) const int stride_k = input_array_strides[Z_AXIS]; #endif #if (N_DIMS > 3) #error "N_DIMS may not be > 3!" #endif /* * set up Jacobian m strides, or dummy values if no Jacobian info */ #if N_DIMS >= 1 const int Jacobian_mi_stride = (Jacobian_info == NULL) ? 0 : Jacobian_info->Jacobian_m_strides[X_AXIS]; #endif #if (N_DIMS >= 2) const int Jacobian_mj_stride = (Jacobian_info == NULL) ? 0 : Jacobian_info->Jacobian_m_strides[Y_AXIS]; #endif #if (N_DIMS >= 3) const int Jacobian_mk_stride = (Jacobian_info == NULL) ? 0 : Jacobian_info->Jacobian_m_strides[Z_AXIS]; #endif #if (N_DIMS > 3) #error "N_DIMS may not be > 3!" #endif /* * if requested, write the interpolation grid to the log file * (we'll write (some of) the individual interpolation coordinates * as we interpolate them (below)) */ if (log_fp != NULL) then { fprintf(log_fp, "%d\t%d", (int) N_DIMS, N_interp_points); { fp grid_min_xyz[MAX_N_DIMS], grid_max_xyz[MAX_N_DIMS]; int axis; for (axis = 0 ; axis < N_DIMS ; ++axis) { grid_min_xyz[axis] = coord_origin[axis] + input_array_min_subscripts[axis]*coord_delta[axis]; grid_max_xyz[axis] = coord_origin[axis] + input_array_max_subscripts[axis]*coord_delta[axis]; } #if (N_DIMS == 1) fprintf(log_fp, "\t%g\t%g\t%g\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0", (double) grid_min_xyz[X_AXIS], (double) coord_delta [X_AXIS], (double) grid_max_xyz[X_AXIS]); #elif (N_DIMS == 2) fprintf(log_fp, "\t%g\t%g\t%g\t%g\t%g\t%g\t0.0\t0.0\t0.0", (double) grid_min_xyz[X_AXIS], (double) coord_delta [X_AXIS], (double) grid_max_xyz[X_AXIS], (double) grid_min_xyz[Y_AXIS], (double) coord_delta [Y_AXIS], (double) grid_max_xyz[Y_AXIS]); #elif (N_DIMS == 3) fprintf(log_fp, "\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g", (double) grid_min_xyz[X_AXIS], (double) coord_delta [X_AXIS], (double) grid_max_xyz[X_AXIS], (double) grid_min_xyz[Y_AXIS], (double) coord_delta [Y_AXIS], (double) grid_max_xyz[Y_AXIS], (double) grid_min_xyz[Z_AXIS], (double) coord_delta [Z_AXIS], (double) grid_max_xyz[Z_AXIS]); #else #error "N_DIMS may not be > 3!" #endif } } /* * interpolate at each point */ { int pt; int return_status = 0; error_info->error_count = 0; if (debug > 0) then { printf("AEILocalInterp:: in interpolator fn: N_DIMS=%d N_interp_points=%d\n", (int) N_DIMS, N_interp_points); fflush(stdout); } for (pt = 0 ; pt < N_interp_points ; ++pt) { /* this struct holds a molecule-sized piece of a single */ /* real input array, or of a real/complex component of */ /* a complex input array */ struct DATA_STRUCT data; /* each of these structs holds a molecule-sized set of */ /* interpolation coefficients -- recall we are representing */ /* the interpolant as a linear combination of the data values */ #ifdef HAVE_OP_I struct COEFFS_STRUCT coeffs_I; #endif #ifdef HAVE_OP_DX struct COEFFS_STRUCT coeffs_dx; #endif #ifdef HAVE_OP_DY struct COEFFS_STRUCT coeffs_dy; #endif #ifdef HAVE_OP_DZ struct COEFFS_STRUCT coeffs_dz; #endif #ifdef HAVE_OP_DXX struct COEFFS_STRUCT coeffs_dxx; #endif #ifdef HAVE_OP_DXY struct COEFFS_STRUCT coeffs_dxy; #endif #ifdef HAVE_OP_DXZ struct COEFFS_STRUCT coeffs_dxz; #endif #ifdef HAVE_OP_DYY struct COEFFS_STRUCT coeffs_dyy; #endif #ifdef HAVE_OP_DYZ struct COEFFS_STRUCT coeffs_dyz; #endif #ifdef HAVE_OP_DZZ struct COEFFS_STRUCT coeffs_dzz; #endif /* * load the interpolation point coordinates * from the interp_coords[] arrays * FIXME: Maybe it would be better (faster) to do this * with N_DIMS open-coded calls on a function? * But then we'd have to have a sentinal value * return for the unknown-type-code error case. * Yuk! :( :( */ 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_coords_type_code) { 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[pt]; break; } #ifdef HAVE_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[pt]; break; } #endif #ifdef HAVE_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[pt]; break; } #endif #ifdef HAVE_CCTK_REAL16 case CCTK_VARIABLE_REAL16: { /* FIXME: maybe we should warn (once per cactus run) */ /* that we may be doing arithmetic in lower */ /* precision than the interp coords? */ const CCTK_REAL16 *const interp_coords_ptr_real16 = (const CCTK_REAL16 *) interp_coords_ptr; interp_coords_fp[axis] = interp_coords_ptr_real16[pt]; break; } #endif default: CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): interpolation coordinates datatype %d\n" " is not supported!" , interp_coords_type_code); /*NOTREACHED*/ return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ /* end of switch (interp_coords_type_code) */ } /* end of for (axis = ...) loop to get interp coords */ } } if (debug >= 6) then { #if (N_DIMS == 1) printf("pt=%d interp coord %.16g\n", pt, (double) interp_coords_fp[X_AXIS]); #elif (N_DIMS == 2) printf("pt=%d interp coords %.16g %.16g\n", pt, (double) interp_coords_fp[X_AXIS], (double) interp_coords_fp[Y_AXIS]); #elif (N_DIMS == 3) printf("pt=%d interp coords %.16g %.16g %.16g\n", pt, (double) interp_coords_fp[X_AXIS], (double) interp_coords_fp[Y_AXIS], (double) interp_coords_fp[Z_AXIS]); #else #error "N_DIMS may not be > 3!" #endif fflush(stdout); } /* * if requested, write this point's interpolation coords to the log file */ if ( (log_fp != NULL) && ((pt < 3) || (pt == N_interp_points-1)) ) then { #if (N_DIMS == 1) fprintf(log_fp, "\t%g\t0.0\t0.0", (double) interp_coords_fp[X_AXIS]); #elif (N_DIMS == 2) fprintf(log_fp, "\t%g\t%g\t0.0", (double) interp_coords_fp[X_AXIS], (double) interp_coords_fp[Y_AXIS]); #elif (N_DIMS == 3) fprintf(log_fp, "\t%g\t%g\t%g", (double) interp_coords_fp[X_AXIS], (double) interp_coords_fp[Y_AXIS], (double) interp_coords_fp[Z_AXIS]); #else #error "N_DIMS may not be > 3!" #endif } /* * compute position of interpolation molecules with respect to * the grid, i.e. compute (x,y,z) coordinates of interpolation * point relative to molecule center, in units of the grid spacing * * n.b. we need the final answers in variables with the magic * names (x,y,z) (machine-generated code uses these names), * but we use temp variables as intermediates for these and * for center_[ijk] for (likely) better performance: * the temp variables have their addresses taken and so * may not be register-allocated, whereas the final variables * are "clean" and thus more likely to be register-allocated */ { int this_point_status = 0; int center_ijk_temp[MAX_N_DIMS]; fp xyz_temp[MAX_N_DIMS]; int axis; for (axis = 0 ; axis < N_DIMS ; ++axis) { const int ibndry_min = 2*axis; const int ibndry_max = 2*axis + 1; if (debug >= 8) then { printf("axis=%d ibndry_{min,max}=%d,%d MOLECULE_SIZE=%d\n", axis, ibndry_min, ibndry_max, MOLECULE_SIZE); printf(" coord_origin[%d]=%g coord_delta[%d]=%g\n", axis, (double)coord_origin[axis], axis, (double)coord_delta[axis]); printf(" input_array_[min,max]_subscripts[%d]=[%d,%d]\n", axis, (int)input_array_min_subscripts[axis], (int)input_array_max_subscripts[axis]); printf(" N_boundary_points_to_omit[ibndry_[min,max]]=[%d,%d]\n", (int)N_boundary_points_to_omit[ibndry_min], (int)N_boundary_points_to_omit[ibndry_max]); } { const int mp_status = AEILocalInterp_molecule_posn (coord_origin[axis], coord_delta[axis], input_array_min_subscripts[axis] + N_boundary_points_to_omit[ibndry_min], input_array_max_subscripts[axis] - N_boundary_points_to_omit[ibndry_max], MOLECULE_SIZE, boundary_off_centering_tolerance[ibndry_min], boundary_off_centering_tolerance[ibndry_max], boundary_extrapolation_tolerance[ibndry_min], boundary_extrapolation_tolerance[ibndry_max], interp_coords_fp[axis], debug, ¢er_ijk_temp[axis], &xyz_temp[axis]); if (debug >= 8) then printf("==> got mp_status=%d\n", mp_status); /* * unfortunately, the status codes from AEILocalInterp_molecule_posn() * are different from the ones we need ==> we have to decode them * to get the status for this axis */ { int this_axis_status = -1; switch (mp_status) { case 0: this_axis_status = 0; break; case MOLECULE_POSN_ERROR_GRID_TINY: this_axis_status = CCTK_ERROR_INTERP_GRID_TOO_SMALL; break; case MOLECULE_POSN_ERROR_X_LT_MIN: case MOLECULE_POSN_ERROR_X_GT_MAX: this_axis_status = CCTK_ERROR_INTERP_POINT_OUTSIDE; break; case MOLECULE_POSN_ERROR_NAN: this_axis_status = CCTK_ERROR_INTERP_COORD_NAN; break; case MOLECULE_POSN_ERROR_DX_ZERO: this_axis_status = CCTK_ERROR_INTERP_DELTA_X_ZERO; break; default: CCTK_VWarn(BUG_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): internal error!\n" " unexpected result %d from AEILocalInterp_molecule_posn()\n" " at pt=%d axis=%d!\n" , mp_status, pt, axis); /*NOTREACHED*/ } if (this_axis_status < this_point_status) then this_point_status = this_axis_status; /* end of for (axis = ...) loop to locate interp points in the grid */ } } } if (error_info->per_point_status != NULL) then error_info->per_point_status[pt] = this_point_status; if (this_point_status < 0) then { ++error_info->error_count; if (error_info->print_warning_msg) then { fp grid_min_xyz[MAX_N_DIMS], grid_max_xyz[MAX_N_DIMS]; for (axis = 0 ; axis < N_DIMS ; ++axis) { grid_min_xyz[axis] = coord_origin[axis] + input_array_min_subscripts[axis]*coord_delta[axis]; grid_max_xyz[axis] = coord_origin[axis] + input_array_max_subscripts[axis]*coord_delta[axis]; } CCTK_VWarn(WARNING_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform():\n" " interpolation point is either outside the grid,\n" " or inside but too close to the grid boundary!\n" " (this may be caused by a global interpolation with\n" " driver::ghost_size too small)\n" " 0-origin interpolation point number pt=%d of N_interp_points=%d\n" #if (N_DIMS == 1) " interpolation point x=%g\n" " grid x_min(delta_x)x_max = %g(%g)%g\n" , pt, N_interp_points, (double) interp_coords_fp[X_AXIS], (double) grid_min_xyz[X_AXIS], (double) coord_delta [X_AXIS], (double) grid_max_xyz[X_AXIS]); #elif (N_DIMS == 2) " interpolation point (x,y)=(%g,%g)\n" " grid x_min(delta_x)x_max = %g(%g)%g\n" " grid y_min(delta_y)y_max = %g(%g)%g\n" , pt, N_interp_points, (double) interp_coords_fp[X_AXIS], (double) interp_coords_fp[Y_AXIS], (double) grid_min_xyz[X_AXIS], (double) coord_delta [X_AXIS], (double) grid_max_xyz[X_AXIS], (double) grid_min_xyz[Y_AXIS], (double) coord_delta [Y_AXIS], (double) grid_max_xyz[Y_AXIS]); #elif (N_DIMS == 3) " interpolation point (x,y,z)=(%g,%g,%g)\n" " grid x_min(delta_x)x_max = %g(%g)%g\n" " grid y_min(delta_y)y_max = %g(%g)%g\n" " grid z_min(delta_z)z_max = %g(%g)%g\n" , pt, N_interp_points, (double) interp_coords_fp[X_AXIS], (double) interp_coords_fp[Y_AXIS], (double) interp_coords_fp[Z_AXIS], (double) grid_min_xyz[X_AXIS], (double) coord_delta [X_AXIS], (double) grid_max_xyz[X_AXIS], (double) grid_min_xyz[Y_AXIS], (double) coord_delta [Y_AXIS], (double) grid_max_xyz[Y_AXIS], (double) grid_min_xyz[Z_AXIS], (double) coord_delta [Z_AXIS], (double) grid_max_xyz[Z_AXIS]); #else #error "N_DIMS may not be > 3!" #endif } if (this_point_status < return_status) then return_status = this_point_status; if (debug >= 6) then { printf("AEILocalInterp:: pt=%d has this_point_status=%d ==> skipping it\n", pt, this_point_status); fflush(stdout); } /* this point is in error (eg outside grid) */ /* ==> don't try to interpolate at this point! */ continue; } { #if (N_DIMS >= 1) const int center_i = center_ijk_temp[X_AXIS]; const fp x = xyz_temp [X_AXIS]; #endif #if (N_DIMS >= 2) const int center_j = center_ijk_temp[Y_AXIS]; const fp y = xyz_temp [Y_AXIS]; #endif #if (N_DIMS >= 3) const int center_k = center_ijk_temp[Z_AXIS]; const fp z = xyz_temp [Z_AXIS]; #endif if (debug >= 6) then { #if (N_DIMS == 1) printf("interp center_i = %d\n", center_i); printf("interp grid-relative x = %.16g\n", (double) x); #elif (N_DIMS == 2) printf("interp center_ij = %d %d\n", center_i, center_j); printf("interp grid-relative xy = %.16g %.16g\n", (double) x, (double) y); #elif (N_DIMS == 3) printf("interp center_ijk = %d %d %d\n", center_i, center_j, center_k); printf("interp grid-relative xyz = %.16g %.16g %.16g\n", (double) x, (double) y, (double) z); #else #error "N_DIMS may not be > 3!" #endif fflush(stdout); } /* * compute 1-d position of molecule center in input arrays */ { #if (N_DIMS == 1) const int molecule_center_posn = stride_i*center_i; #elif (N_DIMS == 2) const int molecule_center_posn = stride_i*center_i + stride_j*center_j; #elif (N_DIMS == 3) const int molecule_center_posn = stride_i*center_i + stride_j*center_j + stride_k*center_k; #else #error "N_DIMS may not be > 3!" #endif /* * molecule position queries */ if (molecule_positions != NULL) then { #if (N_DIMS >= 1) molecule_positions[X_AXIS][pt] = center_i; #endif #if (N_DIMS >= 2) molecule_positions[Y_AXIS][pt] = center_j; #endif #if (N_DIMS >= 3) molecule_positions[Z_AXIS][pt] = center_k; #endif } /* * compute the coefficients at this point for whichever * operation_codes[] values (derivatatives) are wanted * (these are polynomials in the variables (x,y,z)) */ #ifdef HAVE_OP_I if (want_I) then compute_coeffs_I(XYZ, &coeffs_I); #endif #ifdef HAVE_OP_DX if (want_dx) then compute_coeffs_dx(XYZ, &coeffs_dx); #endif #ifdef HAVE_OP_DY if (want_dy) then compute_coeffs_dy(XYZ, &coeffs_dy); #endif #ifdef HAVE_OP_DZ if (want_dz) then compute_coeffs_dz(XYZ, &coeffs_dz); #endif #ifdef HAVE_OP_DXX if (want_dxx) then compute_coeffs_dxx(XYZ, &coeffs_dxx); #endif #ifdef HAVE_OP_DXY if (want_dxy) then compute_coeffs_dxy(XYZ, &coeffs_dxy); #endif #ifdef HAVE_OP_DXZ if (want_dxz) then compute_coeffs_dxz(XYZ, &coeffs_dxz); #endif #ifdef HAVE_OP_DYY if (want_dyy) then compute_coeffs_dyy(XYZ, &coeffs_dyy); #endif #ifdef HAVE_OP_DYZ if (want_dyz) then compute_coeffs_dyz(XYZ, &coeffs_dyz); #endif #ifdef HAVE_OP_DZZ if (want_dzz) then compute_coeffs_dzz(XYZ, &coeffs_dzz); #endif /* * compute each output array at this point */ { int out; /* * next 2 initializers must be invalid values to make sure we * execute the ***load*** the first time in the test at the * top of the part loop below */ const void* input_array_ptr__last_load = NULL; int part__last_load = -1; for (out = 0 ; out < N_output_arrays ; ++out) { const int in = operand_indices[out]; const void* const input_array_ptr = input_arrays[in]; /* * ***decode*** the output array datatype * to determine whether it's real or complex, */ const int N_output_parts = AEILocalInterp_decode_N_parts(output_array_type_codes[out]); if (! ((N_output_parts == 1) || (N_output_parts == 2))) then { CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform():\n" " output array doesn't seem to be a real or complex number,\n" " or more precisely, output array has number of \"real parts\"\n" " (1=real, 2=complex) which isn't 1 or 2!\n" " 0-origin output #out=%d\n" " datatype code=%d\n" " (datatype codes are defined by the Cactus flesh,\n" " see src/include/cctk_Constants.h)\n" " ==> N_parts=%d" , out, (int) output_array_type_codes[out], N_output_parts); /*NOTREACHED*/ return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } { int part; for (part = 0 ; part < N_output_parts ; ++part) { if ( (input_array_ptr != NULL) && ((input_array_ptr != input_array_ptr__last_load) || (part != part__last_load)) ) then { /* remember when we did the following load */ input_array_ptr__last_load = input_array_ptr; part__last_load = part; /* * ***decode*** the input array datatype * to determine whether it's real or complex, */ { const int N_input_parts = AEILocalInterp_decode_N_parts(input_array_type_codes[in]); if (N_input_parts != N_output_parts) then { CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform():\n" " data types are incompatible between input and output arrays, or\n" " more precisely, number of \"real parts\" (1=real, 2=complex) differ!\n" " 0-origin input #in =%d datatype code=%d N_parts=%d\n" " 0-origin output #out=%d datatype code=%d N_parts=%d\n" " (datatype codes are defined by the Cactus flesh,\n" " see src/include/cctk_Constants.h)" , in, (int) input_array_type_codes[in], N_input_parts, out, (int) output_array_type_codes[out], N_output_parts); /*NOTREACHED*/ return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ } /* * load the molecule-sized piece of * input_arrays[in][part] at this molecule * position, into the data struct */ { const int input_posn = molecule_center_posn + input_array_offsets[in]; switch (input_array_type_codes[in]) { case CCTK_VARIABLE_REAL: { const CCTK_REAL *const input_array_ptr_real = ((const CCTK_REAL *) input_array_ptr) + input_posn; LOAD_DATA_REAL(input_array_ptr_real, STRIDE_IJK, &data); #if (N_DIMS == 2) && (MOLECULE_SIZE == 4) if (debug >= 10) then { /* we assume a cubical molecule :( */ printf("AEILocalInterp:: loaded data is:\n"); printf(" data_m1_m1 = %g\n", (double) data.data_m1_m1); printf(" data_0_m1 = %g\n", (double) data.data_0_m1); printf(" data_p1_m1 = %g\n", (double) data.data_p1_m1); printf(" data_p2_m1 = %g\n", (double) data.data_p2_m1); printf(" data_m1_0 = %g\n", (double) data.data_m1_0); printf(" data_0_0 = %g\n", (double) data.data_0_0); printf(" data_p1_0 = %g\n", (double) data.data_p1_0); printf(" data_p2_0 = %g\n", (double) data.data_p2_0); printf(" data_m1_p1 = %g\n", (double) data.data_m1_p1); printf(" data_0_p1 = %g\n", (double) data.data_0_p1); printf(" data_p1_p1 = %g\n", (double) data.data_p1_p1); printf(" data_p2_p1 = %g\n", (double) data.data_p2_p1); printf(" data_m1_p2 = %g\n", (double) data.data_m1_p2); printf(" data_0_p2 = %g\n", (double) data.data_0_p2); printf(" data_p1_p2 = %g\n", (double) data.data_p1_p2); printf(" data_p2_p2 = %g\n", (double) data.data_p2_p2); fflush(stdout); } #endif break; } #ifdef HAVE_CCTK_REAL4 case CCTK_VARIABLE_REAL4: { const CCTK_REAL4 *const input_array_ptr_real4 = ((const CCTK_REAL4 *) input_array_ptr) + input_posn; LOAD_DATA_REAL4(input_array_ptr_real4, STRIDE_IJK, &data); break; } #endif #ifdef HAVE_CCTK_REAL8 case CCTK_VARIABLE_REAL8: { const CCTK_REAL8 *const input_array_ptr_real8 = ((const CCTK_REAL8 *) input_array_ptr) + input_posn; LOAD_DATA_REAL8(input_array_ptr_real8, STRIDE_IJK, &data); break; } #endif #ifdef HAVE_CCTK_REAL16 case CCTK_VARIABLE_REAL16: { /* FIXME: maybe we should warn (once per cactus run) that we may be */ /* doing arithmetic in lower precision than the data type? */ const CCTK_REAL16 *const input_array_ptr_real16 = ((const CCTK_REAL16 *) input_array_ptr) + input_posn; LOAD_DATA_REAL16(input_array_ptr_real16, STRIDE_IJK, &data); break; } #endif case CCTK_VARIABLE_COMPLEX: { const CCTK_REAL (*const input_array_ptr_complex)[COMPLEX_N_PARTS] = ((const CCTK_REAL (*)[COMPLEX_N_PARTS]) input_array_ptr) + input_posn; LOAD_DATA_COMPLEX(input_array_ptr_complex, STRIDE_IJK, part, &data); break; } #ifdef HAVE_CCTK_COMPLEX8 case CCTK_VARIABLE_COMPLEX8: { const CCTK_REAL4 (*const input_array_ptr_complex8)[COMPLEX_N_PARTS] = ((const CCTK_REAL4 (*)[COMPLEX_N_PARTS]) input_array_ptr) + input_posn; LOAD_DATA_COMPLEX8(input_array_ptr_complex8, STRIDE_IJK, part, &data); break; } #endif #ifdef HAVE_CCTK_COMPLEX16 case CCTK_VARIABLE_COMPLEX16: { const CCTK_REAL8 (*const input_array_ptr_complex16)[COMPLEX_N_PARTS] = ((const CCTK_REAL8 (*)[COMPLEX_N_PARTS]) input_array_ptr) + input_posn; LOAD_DATA_COMPLEX16(input_array_ptr_complex16, STRIDE_IJK, part, &data); break; } #endif #ifdef HAVE_CCTK_COMPLEX32 case CCTK_VARIABLE_COMPLEX32: { /* FIXME: maybe we should warn (once per cactus run) that we may be */ /* doing arithmetic in lower precision than the data type? */ const CCTK_REAL16 (*const input_array_ptr_complex32)[COMPLEX_N_PARTS] = ((const CCTK_REAL16 (*)[COMPLEX_N_PARTS]) input_array_ptr) + input_posn; LOAD_DATA_COMPLEX32(input_array_ptr_complex32, STRIDE_IJK, part, &data); break; } #endif default: CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): input datatype %d not supported!\n" " (0-origin) input #in=%d" , (int) input_array_type_codes[in], in); /*NOTREACHED*/ return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ /* end of switch (input_array_type_codes[in]) */ } } } /* end of load input array values */ } if (output_arrays[out] != NULL) then { /* * interpolate at this point */ fp result; switch (operation_codes[out]) { #ifdef HAVE_OP_I case 0: result = EVALUATE_MOLECULE(&coeffs_I, &data); break; #endif #ifdef HAVE_OP_DX case 1: result = inverse_dx * EVALUATE_MOLECULE(&coeffs_dx, &data); break; #endif #ifdef HAVE_OP_DY case 2: result = inverse_dy * EVALUATE_MOLECULE(&coeffs_dy, &data); break; #endif #ifdef HAVE_OP_DZ case 3: result = inverse_dz * EVALUATE_MOLECULE(&coeffs_dz, &data); break; #endif #ifdef HAVE_OP_DXX case 11: result = inverse_dx * inverse_dx * EVALUATE_MOLECULE(&coeffs_dxx, &data); break; #endif #ifdef HAVE_OP_DXY case 12: case 21: result = inverse_dx * inverse_dy * EVALUATE_MOLECULE(&coeffs_dxy, &data); break; #endif #ifdef HAVE_OP_DXZ case 13: case 31: result = inverse_dx * inverse_dz * EVALUATE_MOLECULE(&coeffs_dxz, &data); break; #endif #ifdef HAVE_OP_DYY case 22: result = inverse_dy * inverse_dy * EVALUATE_MOLECULE(&coeffs_dyy, &data); break; #endif #ifdef HAVE_OP_DYZ case 23: case 32: result = inverse_dy * inverse_dz * EVALUATE_MOLECULE(&coeffs_dyz, &data); break; #endif #ifdef HAVE_OP_DZZ case 33: result = inverse_dz * inverse_dz * EVALUATE_MOLECULE(&coeffs_dzz, &data); break; #endif default: CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): generalized interpolation operation_code %d\n" " is not supported!\n" " (0-origin) output #out=%d" , (int) operation_codes[out], out); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ /* end of switch (operation_codes[out]) */ } /* * store result in output array */ if (debug >= 10) then { if ((pt & (pt-1)) == 0) /* pt is 0 or a power of 2 */ then { printf("AEILocalInterp:: pt=%d out=%d --> storing result %g\n", pt, out, (double) result); fflush(stdout); } } switch (output_array_type_codes[out]) { case CCTK_VARIABLE_REAL: { CCTK_REAL *const output_array_ptr_real = (CCTK_REAL *) output_arrays[out]; if (debug >= 10) then { if ((pt & (pt-1)) == 0) /* pt is 0 or a power of 2 */ then { printf(" result addr is %p\n", (void *) &output_array_ptr_real[pt]); printf(" previous value there was %g\n", output_array_ptr_real[pt]); fflush(stdout); } } output_array_ptr_real[pt] = (CCTK_REAL) result; break; } #ifdef HAVE_CCTK_REAL4 case CCTK_VARIABLE_REAL4: { CCTK_REAL4 *const output_array_ptr_real4 = (CCTK_REAL4 *) output_arrays[out]; output_array_ptr_real4[pt] = (CCTK_REAL4) result; break; } #endif #ifdef HAVE_CCTK_REAL8 case CCTK_VARIABLE_REAL8: { CCTK_REAL8 *const output_array_ptr_real8 = (CCTK_REAL8 *) output_arrays[out]; output_array_ptr_real8[pt] = (CCTK_REAL8) result; break; } #endif #ifdef HAVE_CCTK_REAL16 case CCTK_VARIABLE_REAL16: { CCTK_REAL16 *const output_array_ptr_real16 = (CCTK_REAL16 *) output_arrays[out]; output_array_ptr_real16[pt] = (CCTK_REAL16) result; break; } #endif case CCTK_VARIABLE_COMPLEX: { CCTK_REAL (*const output_array_ptr_complex)[COMPLEX_N_PARTS] = (CCTK_REAL (*)[COMPLEX_N_PARTS]) output_arrays[out]; output_array_ptr_complex[pt][part] = (CCTK_REAL) result; break; } #ifdef HAVE_CCTK_COMPLEX8 case CCTK_VARIABLE_COMPLEX8: { CCTK_REAL4 (*const output_array_ptr_complex8)[COMPLEX_N_PARTS] = (CCTK_REAL4 (*)[COMPLEX_N_PARTS]) output_arrays[out]; output_array_ptr_complex8[pt][part] = (CCTK_REAL4) result; break; } #endif /* HAVE_CCTK_COMPLEX8 */ #ifdef HAVE_CCTK_COMPLEX16 case CCTK_VARIABLE_COMPLEX16: { CCTK_REAL8 (*const output_array_ptr_complex16)[COMPLEX_N_PARTS] = (CCTK_REAL8 (*)[COMPLEX_N_PARTS]) output_arrays[out]; output_array_ptr_complex16[pt][part] = (CCTK_REAL8) result; break; } #endif /* HAVE_CCTK_COMPLEX16 */ #ifdef HAVE_CCTK_COMPLEX32 case CCTK_VARIABLE_COMPLEX32: { CCTK_REAL16 (*const output_array_ptr_complex32)[COMPLEX_N_PARTS] = (CCTK_REAL16 (*)[COMPLEX_N_PARTS]) output_arrays[out]; output_array_ptr_complex32[pt][part] = (CCTK_REAL16) result; break; } #endif /* HAVE_CCTK_COMPLEX32 */ default: CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): output datatype %d not supported!\n" " (0-origin) output #out=%d" , (int) output_array_type_codes[out], out); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ /* end of switch (output type code) */ } /* end of if (output_arrays[out] != NULL) */ } /* * handle querying the Jacobian */ if ( (Jacobian_info != NULL) && (Jacobian_info->Jacobian_pointer[out] != NULL)) then { CCTK_REAL *const Jacobian_ptr = Jacobian_info->Jacobian_pointer[out] + Jacobian_info->Jacobian_offset[out] + Jacobian_info->Jacobian_interp_point_stride*pt + Jacobian_info->Jacobian_part_stride*part; switch (operation_codes[out]) { #ifdef HAVE_OP_I case 0: STORE_COEFFS(1.0, &coeffs_I, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif #ifdef HAVE_OP_DX case 1: STORE_COEFFS(inverse_dx, &coeffs_dx, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif #ifdef HAVE_OP_DY case 2: STORE_COEFFS(inverse_dy, &coeffs_dy, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif #ifdef HAVE_OP_DZ case 3: STORE_COEFFS(inverse_dz, &coeffs_dz, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif #ifdef HAVE_OP_DXX case 11: STORE_COEFFS(inverse_dx*inverse_dx, &coeffs_dxx, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif #ifdef HAVE_OP_DXY case 12: case 21: STORE_COEFFS(inverse_dx*inverse_dy, &coeffs_dxy, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif #ifdef HAVE_OP_DXZ case 13: case 31: STORE_COEFFS(inverse_dx*inverse_dz, &coeffs_dxz, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif #ifdef HAVE_OP_DYY case 22: STORE_COEFFS(inverse_dy*inverse_dy, &coeffs_dyy, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif #ifdef HAVE_OP_DYZ case 23: case 32: STORE_COEFFS(inverse_dy*inverse_dz, &coeffs_dyz, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif #ifdef HAVE_OP_DZZ case 33: STORE_COEFFS(inverse_dz*inverse_dz, &coeffs_dzz, Jacobian_ptr, JACOBIAN_MIJK_STRIDE); break; #endif default: CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform(): generalized interpolation operation_code %d\n" " is not supported!\n" " (0-origin) output #out=%d" , (int) operation_codes[out], out); return UTIL_ERROR_BAD_INPUT; /*** ERROR RETURN ***/ /* end of switch(operation_codes[out])*/ } /* end of Jacobian-query code */ } /* end of for (part = ...) loop */ } } /* end of for (out = ...) loop */ } } } } } /* end of for (pt = ...) loop */ } /* * if we're writing logging information, finish this */ if (log_fp != NULL) then { switch (N_interp_points) { case 0: /* we printed 0 points, but need 4 */ fprintf(log_fp, "\t0.0\t0.0\t0.0"); fprintf(log_fp, "\t0.0\t0.0\t0.0"); /* fall through */ case 1: /* we printed 2 points (0, N_interp_points-1), but need 4 */ fprintf(log_fp, "\t0.0\t0.0\t0.0"); /* fall through */ case 2: /* we printed 3 points (0, 1, N_interp_points-1), but need 4 */ fprintf(log_fp, "\t0.0\t0.0\t0.0"); /* fall through */ } fprintf(log_fp, "\n"); fflush(log_fp); } return return_status; } } } } /******************************************************************************/ /*@@ @routine compute_coeffs_I @routine compute_coeffs_dx @routine compute_coeffs_dy @routine compute_coeffs_dz @routine compute_coeffs_dxx @routine compute_coeffs_dxy @routine compute_coeffs_dxz @routine compute_coeffs_dyy @routine compute_coeffs_dyz @routine compute_coeffs_dzz @date 29.Aug.2002 @author Jonathan Thornburg @desc Each of these functions computes the corresponding set of interpolation coefficients at the point given by the XYZ variables, using the machine-generated experessions in the files COEFFS_{I,DX,DY,DXX,DXY,DYY,...}_COMPUTE_FILE_NAME These functions *must* be static, because this whole file is #included (and hence compiled) multiple times. @enddesc @@*/ #ifdef HAVE_OP_I static void compute_coeffs_I(FP_XYZ, struct COEFFS_STRUCT *coeffs_I) { #include COEFFS_I_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DX static void compute_coeffs_dx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dx) { #include COEFFS_DX_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DY static void compute_coeffs_dy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dy) { #include COEFFS_DY_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DZ static void compute_coeffs_dz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dz) { #include COEFFS_DZ_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DXX static void compute_coeffs_dxx(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxx) { #include COEFFS_DXX_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DXY static void compute_coeffs_dxy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxy) { #include COEFFS_DXY_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DXZ static void compute_coeffs_dxz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dxz) { #include COEFFS_DXZ_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DYY static void compute_coeffs_dyy(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyy) { #include COEFFS_DYY_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DYZ static void compute_coeffs_dyz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dyz) { #include COEFFS_DYZ_COMPUTE_FILE_NAME } #endif #ifdef HAVE_OP_DZZ static void compute_coeffs_dzz(FP_XYZ, struct COEFFS_STRUCT *coeffs_dzz) { #include COEFFS_DZZ_COMPUTE_FILE_NAME } #endif