From 77a3b86574aa1eb4a09e5093939dd1027902d151 Mon Sep 17 00:00:00 2001 From: tradke Date: Fri, 1 Aug 2003 17:18:09 +0000 Subject: Removed duplicate code for the old CCTK_InterpLocal() API (which is now in CactusBase/LocalInterp). PUGHInterp's implementation of the old CCTK_InterpGV() API is now just a wrapper for CCTK_InterpGridArrays(). git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGH/PUGHInterp/trunk@71 1c20744c-e24a-42ec-9533-f5004cb800e5 --- src/Interpolate.c | 615 ----------------------------- src/Operator.c | 1120 ---------------------------------------------------- src/Startup.c | 93 ++++- src/make.code.defn | 2 +- src/pughInterpGH.h | 27 -- 5 files changed, 91 insertions(+), 1766 deletions(-) delete mode 100644 src/Interpolate.c delete mode 100644 src/Operator.c diff --git a/src/Interpolate.c b/src/Interpolate.c deleted file mode 100644 index b46cb40..0000000 --- a/src/Interpolate.c +++ /dev/null @@ -1,615 +0,0 @@ - /*@@ - @file Interpolate.c - @date Wed 17 Jan 2001 - @author Thomas Radke - @desc - Interpolation of arrays to arbitrary points - - This interpolator is based on the Cactus 3.x Fortran version - written by Paul Walker. It also contains some nice optimization - features from Erik Schnetter. Jonathan Thornburg added some - additional comments in October 2001. - @enddesc - - @history - @date Wed 17 Jan 2001 - @author Thomas Radke - @hdesc Translation from Fortran to C - @date Thu 18 Oct 2001 - @author Jonathan Thornburg - @hdesc Add lots of comments, PUGHINTERP_VERBOSE_DEBUG debugging code - @endhistory - @version $Id$ - @@*/ - -#include -#include -#include -#include -#include - -#include "cctk.h" -#include "cctk_Parameters.h" -#include "pughInterpGH.h" - -/* the rcs ID and its dummy function to use it */ -static const char *rcsid = "$Header$"; -CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Interpolate_c) - -#ifdef PUGHINTERP_VERBOSE_DEBUG - /* if this is >= 0, we print verbose debugging information at this point */ - int PUGHInterp_verbose_debug_n = -1; -#endif - -/* the highest order of interpolation we support so far */ -#define MAXORDER 3 - -/* the highest dimension for variables we can deal with (so far) */ -#define MAXDIM 3 - -/******************************************************************************/ - -/*@@ - @routine INTERPOLATE (macro) - @date 18 Oct 2001 - @author code by ???, these comments by Jonathan Thornburg - @desc - This macro does the interpolation of in_array[] to compute - a single value out_array[n] (actually out_array[n]subpart ; - see the comments below for details). The data to be interpolated - must be real numbers of some type, i.e. if the arrays are - complex this macro must be called separately for the real - and imaginary parts. - @enddesc - - @var cctk_type - @vdesc C type of input and output array elements (might be complex) - @endvar - - @var cctk_subtype - @vdesc C type of actual numbers being interpolated (must be real) - @endvar - - @var subpart - @vdesc string to be suffixed to input/output array element to get - to get real number, i.e. empty string if cctk_type is real, - .Re or .Im as appropriate if cctk_type is complex - @endvar - - @var in_array - @vdesc A pointer to array to be interpolated (strictly speaking, to - the array's [0][0]...[0] element); this is typically passed - as a void * pointer so we typecast it as necessary. - @endvar - - @var out_array - @vdesc A 1-dimensional array where the interpolation result should be - stored; this is typically passed as a void * pointer so we - typecast it as necessary. - @endvar - - @var order - @vdesc The order of the interpolation (1=linear, 2=quadratic, 3=cubic, ...) - @endvar - - @var point - @vdesc [MAXDIM] array of integers giving the integer grid coordinates - of the closest grid point to the interpolation point; the - interpolation molecule/stencil is centered at this point. - @endvar - - @var dims - @vdesc [MAXDIM] array of integers giving the dimensions of in_array . - @endvar - - @var n - @vdesc Position in out_array where we should store the interpolation - result. - @endvar - - @var coeff - @vdesc [MAXDIM][MAX_ORDER+1] array of (floating-point) interpolation - coefficients; detailed semantics are that coeff[axis][m] is the - coefficient of y[m] when the 1-dimensional Lagrange interpolation - polynomial passing through the order+1 points - {(0,y[0]), (1,y[1]), ..., (order,y[order])} - is evaluated at the position x=offset[axis]. - @endvar - @@*/ -/* - * The basic idea here is that conceptually we first interpolate the - * (say) 3D gridfn in the x direction at each y and z grid point, - * then interpolate that 2D plane of values in the y direction at - * each z grid point, and finally interpolate that 1D line of values - * in the z direction. The implementation actually interleaves the - * different directions' interpolations so that only 3 scalar temporaries - * are needed. - */ -#define INTERPOLATE(cctk_type, cctk_subtype, subpart, in_array, out_array, \ - order, point, dims, n, coeff) \ - { \ - int ii, jj, kk; \ - const cctk_type *fi; \ - /* for some reason (probably loop unrolling & pipelining) the compiler \ - produces faster code if arrays are used instead of scalars */ \ - cctk_subtype interp_result, fj[1], fk[1]; \ - \ - \ - interp_result = 0; \ - \ - /* NOTE-MAXDIM: support >3D arrays by adding more loops */ \ - for (kk = 0; kk <= order; kk++) \ - /*if (fabs(coeff[2][kk]) > 1e-6)*/ \ - { \ - fk[0] = 0; \ - for (jj = 0; jj <= order; jj++) \ - /*if (fabs(coeff[1][jj]) > 1e-6)*/ \ - { \ - /* NOTE-MAXDIM: for >3D arrays adapt the index calculation here */ \ - fi = (const cctk_type *) in_array + \ - point[0] + dims[0]*(point[1]+jj + dims[1]*(point[2]+kk)); \ - \ - fj[0] = 0; \ - for (ii = 0; ii <= order; ii++) \ - /*if (fabs(coeff[0][ii]) > 1e-6)*/ \ - { \ - fj[0] += fi[ii]subpart * coeff[0][ii]; \ - } \ - /* at this point we have just computed */ \ - /* fj[0] = in_array[*][jj][kk] interpolated to x=offset[0] */ \ - \ - fk[0] += fj[0] * coeff[1][jj]; \ - } \ - /* at this point we have just computed */ \ - /* fk[0] = fj[0][*][kk] interpolated to y=offset[1] */ \ - /* = in_array[*][*][kk] interpolated to */ \ - /* x=offset[0], y=offset[1] */ \ - \ - interp_result += fk[0] * coeff[2][kk]; \ - } \ - /* at this point we have just computed */ \ - /* interp_result = fk[0][*] interpolated to z=offset[2] */ \ - /* = in_array[*][*][*] interpolated to */ \ - /* x=offset[0], y=offset[1], z=offset[2] */ \ - \ - /* assign the result */ \ - ((cctk_type *) out_array)[n]subpart = interp_result; \ - } /* end of macro */ - -/* this is needed by some preprocessors to pass into INTERPOLATE - as a dummy macro */ -#define NOTHING - -/******************************************************************************/ - -/*@@ - @routine PUGHInterp_Interpolate - @date Wed 17 Jan 2001 - @author Thomas Radke - @desc - This routine interpolates a set of input arrays - to a set of output arrays (one-to-one) at arbitrary points - which are given by their coordinates and the underlying - regular, uniform grid. - - Current limitations of this implementation are: - - arrays up to three (MAXDIM) dimensions only can be handled - - interpolation orders up to three (MAXORDER) only are supported - - coordinates must be given as CCTK_REAL types - - input and output array types must be the same - (no type casting of interpolation results supported) - - Despite of these limitations, the code was programmed almost - generically in that it can easily be extended to support - higher-dimensional arrays or more interpolation orders. - Places where the code would need to be changed to do this, - are marked with NOTE-MAXDIM and/or NOTE-MAXORDER comments - as appropriate. - @enddesc - - @var num_points - @vdesc number of points to interpolate at - @vtype int - @vio in - @endvar - @var num_dims - @vdesc dimensionality of the input arrays - @vtype int - @vio in - @endvar - @var num_arrays - @vdesc number of input/output arrays - @vtype int - @vio in - @endvar - @var dims - @vdesc dimensions of the input arrays - @vtype int[ num_dims ] - @vio in - @endvar - @var coord - @vdesc coordinates to interpolate at - @vtype CCTK_REAL coord[ num_dims * num_points ] - @vio in - @endvar - @var origin - @vdesc origin of the underlying grid - @vtype CCTK_REAL origin[ num_dims ] - @vio in - @endvar - @var delta - @vdesc deltas of the underlying grid - @vtype CCTK_REAL delta[ num_dims ] - @vio in - @endvar - @var in_types - @vdesc CCTK variable types of input arrays - @vtype int in_types[ num_arrays ] - @vio in - @endvar - @var in_arrays - @vdesc list of input arrays - @vtype void *in_arrays[ num_arrays ] - @vio in - @endvar - @var out_types - @vdesc CCTK variable types of output arrays - @vtype int out_types[ num_arrays ] - @vio in - @endvar - @var out_arrays - @vdesc list of output arrays - @vtype void *out_arrays[ num_arrays ] - @vio out - @endvar - - @returntype int - @returndesc - 0 - successful interpolation - -1 - in case of any errors - @endreturndesc -@@*/ -int PUGHInterp_Interpolate (int order, - int num_points, - int num_dims, - int num_arrays, - const int dims[], - const CCTK_REAL coord[], - const CCTK_REAL origin[], - const CCTK_REAL delta[], - const int in_types[], - const void *const in_arrays[], - const int out_types[], - void *const out_arrays[]) -{ - int retval; - int i, a, n, shift; -#if 0 - int out_of_bounds; -#endif - int max_dims[MAXDIM], point[MAXDIM]; - CCTK_REAL delta_inv[MAXDIM]; - CCTK_REAL below[MAXDIM]; - CCTK_REAL offset[MAXDIM]; - CCTK_REAL coeff[MAXDIM][MAXORDER + 1]; - - - /* verify parameters and check against our restrictions */ - retval = -1; - if (num_dims < 1) - { - CCTK_WARN (1, "Number of dimensions must be positive"); - } - else if (num_dims > MAXDIM) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Interpolation of %d-dimensional arrays not implemented", - num_dims); - } - else if (order < 1) - { - CCTK_WARN (1, "Inperpolation order must be positive"); - } - else if (order > MAXORDER) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Interpolation order %d not implemented", order); - } - else if (num_points < 0) - { - CCTK_WARN (1, "Negative number of points given"); - } - else - { - retval = 0; - } - - /* immediately return in case of errors */ - if (retval) - { - return (retval); - } - - /* also immediately return if there's nothing to do */ - if (num_points == 0) - { - return (retval); - } - - /* avoid divisions by delta later on */ - for (i = 0; i < num_dims; i++) - { - delta_inv[i] = 1.0 / delta[i]; - } - - /* duplicate the dims[] vector into one with MAXDIM-1 elements - (with the remaining elements zeroed out) - so that we can use nested loops over MAXDIM dimensions later on */ - memset (max_dims, 0, sizeof (max_dims)); - memcpy (max_dims, dims, (num_dims - 1) * sizeof (int)); - - /* zero out the coefficients and set the elements with index 'order' to one - so that we can use nested loops over MAXDIM dimensions later on */ - memset (coeff, 0, sizeof (coeff)); - for (i = num_dims; i < MAXDIM; i++) - { - coeff[i][0] = 1; - } - - /* zero out the iterator */ - memset (point, 0, sizeof (point)); - - /* loop over all points to interpolate at */ - for (n = 0; n < num_points; n++) - { -#if 0 - /* reset the out-of-bounds flag */ - out_of_bounds = 0; -#endif - - /* loop over all dimensions */ - for (i = 0; i < num_dims; i++) - { - /* closest grid point for stencil/molecule */ - point[i] = floor ((coord[num_dims*n + i] - origin[i]) * delta_inv[i] - - 0.5 * (order - 1)); - -#if 0 - /* test bounds */ - out_of_bounds |= point[i] < 0 || point[i]+order >= dims[i]; -#endif - - /* check that the stencil/molecule isn't bigger than the grid */ - if (order+1 > dims[i]) /* stencil/molecule size = order+1 */ - { - return -1; - } - - /* if beyond lower bound shift the grid point to the right */ - shift = point[i]; - if (shift < 0) - { - point[i] -= shift; - } - - /* if beyond upper bound shift the grid point to the left */ - shift = point[i] + order - (dims[i] - 1); - if (shift > 0) - { - point[i] -= shift; - } - - assert (point[i] >= 0 && point[i]+order < dims[i]); - - /* physical coordinate of that grid point */ - below[i] = origin[i] + point[i] * delta[i]; - - /* offset from that grid point, in fractions of grid points */ - offset[i] = (coord[num_dims*n + i] - below[i]) * delta_inv[i]; - -#if 0 - /* This test ensures that it is really an interpolation, and not - an extrapolation */ - assert (offset[i]>=0.5*(order-1) && offset[i]<=0.5*(order-1)+1.0); -#endif - } - -#ifdef PUGHINTERP_VERBOSE_DEBUG -if (n == PUGHInterp_verbose_debug_n) - { - int ii; - printf("out_of_bounds = %d\n", out_of_bounds); - for (ii = 0 ; ii < num_dims ; ++ii) - { - printf("offset[%d] = %g\n", ii, (double) offset[ii]); - } - } -#endif /* PUGHINTERP_VERBOSE_DEBUG */ - -#if 0 - /* check bounds */ - if (out_of_bounds) - { - char *msg; - - - /* put all information into a single message string for output */ - msg = (char *) malloc (100 + num_dims*(10 + 4*20)); - sprintf (msg, "Interpolation stencil/molecule out of bounds at grid point [%d", - point[0]); - for (i = 1; i < num_dims; i++) - { - sprintf (msg, "%s, %d", msg, point[i]); - } - sprintf (msg, "%s]\nrange would be min/max [%f / %f", msg, - (double) below[0], (double) (below[0] + offset[0])); - for (i = 1; i < num_dims; i++) - { - sprintf (msg, "%s, %f / %f", msg, - (double) below[i], (double) (below[i] + offset[i])); - } - sprintf (msg, "%s]\ngrid is min/max [%f / %f", msg, - (double) origin[0], (double) (origin[0] + (dims[0]-1)*delta[0])); - for (i = 1; i < num_dims; i++) - { - sprintf (msg, "%s, %f / %f", msg, - (double)origin[i], (double)(origin[i] + (dims[i]-1)*delta[i])); - } - sprintf (msg, "%s]", msg); - CCTK_WARN (1, msg); - free (msg); - - continue; - } -#endif - - /* - * *** compute the interpolation coefficients according to the order *** - * - * (Thanks to Erik for formulating these so nicely.) - * - * These formulas are "just" the coefficients of the classical - * Lagrange interpolation polynomials along each dimension. - * For example, in 1 dimension the unique quadratic passing - * through the 3 points {(x0,y0), (x1,y1), (x2,y2)} is: - * ( x-x1)( x-x2) ( x-x0)( x-x2) ( x-x0)( x-x1) - * -------------- y0 + -------------- y1 + -------------- y2 - * (x0-x1)(x0-x2) (x1-x0)(x1-x2) (x2-x0)(x2-x1) - * (It's easy to see this: each of the terms is yi if x=xi, or - * zero if x=any other xj.) To get the formulas below, just negate - * each (x-x) factor, and substitute the values xi=i. - */ - /* - * NOTE-MAXORDER: support higher interpolation orders by adding the - * appropriate coefficients in another else branch - */ - switch (order) - { - case 1: - /* first order (linear) 1D interpolation */ - for (i = 0; i < num_dims; i++) - { - coeff[i][0] = 1 - offset[i]; - coeff[i][1] = offset[i]; - } - break; - case 2: - /* second order (quadratic) 1D interpolation */ - for (i = 0; i < num_dims; i++) - { - coeff[i][0] = (1-offset[i]) * (2-offset[i]) / ( 2 * 1 ); - coeff[i][1] = ( -offset[i]) * (2-offset[i]) / ( 1 * (-1)); - coeff[i][2] = ( -offset[i]) * (1-offset[i]) / ((-1) * (-2)); - } - break; - case 3: - /* third order (cubic) 1D interpolation */ - for (i = 0; i < num_dims; i++) - { - coeff[i][0] = (1-offset[i]) * (2-offset[i]) * (3-offset[i]) / - ( 3 * 2 * 1 ); - coeff[i][1] = ( -offset[i]) * (2-offset[i]) * (3-offset[i]) / - ( 2 * 1 * (-1)); - coeff[i][2] = ( -offset[i]) * (1-offset[i]) * (3-offset[i]) / - ( 1 * (-1) * (-2)); - coeff[i][3] = ( -offset[i]) * (1-offset[i]) * (2-offset[i]) / - ((-1) * (-2) * (-3)); - } - break; - default: - CCTK_WARN (0, "Implementation error"); - } - -#ifdef PUGHINTERP_VERBOSE_DEBUG -if (n == PUGHInterp_verbose_debug_n) - { - int ii,mm; - for (ii = 0 ; ii < num_dims ; ++ii) - { - for (mm = 0 ; mm <= order ; ++mm) - { - printf("coeff[%d][%d] = %g\n", - ii, mm, (double) coeff[ii][mm]); - } - } - } -#endif /* PUGHINTERP_VERBOSE_DEBUG */ - - /* now loop over all arrays to interpolate at the current point */ - for (a = 0; a < num_arrays; a++) - { - /* sorry, for now input and output arrays must be of same type */ - if (in_types[a] != out_types[a]) - { - CCTK_WARN (1, "Type casting of interpolation results not implemented"); - continue; - } - - /* now do the interpolation according to the array type - we support all kinds of CCTK_REAL* and CCTK_COMPLEX* types here */ - if (in_types[a] == CCTK_VARIABLE_REAL) - { - INTERPOLATE (CCTK_REAL, CCTK_REAL, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX) - { - INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX, CCTK_REAL, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#ifdef CCTK_REAL4 - else if (in_types[a] == CCTK_VARIABLE_REAL4) - { - INTERPOLATE (CCTK_REAL4, CCTK_REAL4, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX8) - { - INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX8, CCTK_REAL4, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#endif -#ifdef CCTK_REAL8 - else if (in_types[a] == CCTK_VARIABLE_REAL8) - { - INTERPOLATE (CCTK_REAL8, CCTK_REAL8, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX16) - { - INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX16, CCTK_REAL8, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#endif -#ifdef CCTK_REAL16 - else if (in_types[a] == CCTK_VARIABLE_REAL16) - { - INTERPOLATE (CCTK_REAL16, CCTK_REAL16, NOTHING, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } - else if (in_types[a] == CCTK_VARIABLE_COMPLEX32) - { - INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Re, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - INTERPOLATE (CCTK_COMPLEX32, CCTK_REAL16, .Im, in_arrays[a], - out_arrays[a], order, point, max_dims, n, coeff); - } -#endif - else - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Unsupported variable type %d", (int) in_types[a]); - } - } /* end of loop over all arrays */ - - } /* end of loop over all points to interpolate at */ - - /* we're done */ - return (retval); -} diff --git a/src/Operator.c b/src/Operator.c deleted file mode 100644 index 45c6dad..0000000 --- a/src/Operator.c +++ /dev/null @@ -1,1120 +0,0 @@ - /*@@ - @file Operator.c - @date Tue Apr 15 18:22:45 1997 - @author Paul Walker - @desc - Definition of interpolation operators for regular uniform grids. - @enddesc - - @history - @date Sun Jul 04 1999 - @author Thomas Radke - @hdesc conversion to Cactus 4.0 (copied from pughGetPoints.c) - @date Wed 31 Jan 2001 - @author Thomas Radke - @hdesc translation of fortran interpolators into C - @endhistory - @version $Id$ - @@*/ - -#include -#include -#include /* floor(3) */ - -#include "cctk.h" -#include "CactusPUGH/PUGH/src/include/pugh.h" -#include "pughInterpGH.h" - -/* the rcs ID and its dummy function to use it */ -static const char *rcsid = "$Header$"; -CCTK_FILEVERSION(CactusPUGH_PUGHInterp_Operator_c) - -/* uncomment this to get some debugging output */ -/* #define PUGHINTERP_DEBUG 1 */ - -/* macro do sort interpolation results from a single communication buffer - into their appropriate output arrays */ -#define SORT_TYPED_ARRAY(cctk_type) \ - { \ - int _i; \ - cctk_type *_src, *_dst; \ - \ - \ - _src = (cctk_type *) this->buf; \ - _dst = (cctk_type *) output_arrays[array]; \ - for (_i = 0; _i < myGH->N_points_from[proc]; _i++) \ - { \ - _dst[myGH->indices[_i + offset]] = *_src++; \ - } \ - this->buf = (char *) _src; \ - } - - -#ifdef CCTK_MPI -/* internal structure describing a handle for a single CCTK data type */ -typedef struct -{ - int vtypesize; /* variable type's size in bytes */ - MPI_Datatype mpitype; /* corresponding MPI datatype */ - int N_arrays; /* number of in/out arrays */ - void *sendbuf; /* communication send buffer for this type */ - void *recvbuf; /* communication receive buffer for this type */ - char *buf; /* work pointer for sendbuf */ -} type_desc_t; -#endif - - -/* prototypes of routines defined in this source file */ -static int CheckArguments (const cGH *GH, - int N_dims, - int N_points, - int N_input_arrays, - int N_output_arrays, - const int interp_coord_array_types[]); -static int CheckOutOfBounds (const cGH *GH, const char *coord_system_name, - int order, int N_dims, int N_points, - const int *dims, const CCTK_REAL *const *coords); -#ifdef CCTK_MPI -static int GetLocalCoords (const cGH *GH, - int N_points, - const char *coord_system_name, - const pGExtras *extras, - const CCTK_REAL *const coords[], - int *N_local_points, - CCTK_REAL **local_coords); -#endif - - -/*@@ - @routine PUGHInterp_InterpGV - @date Sun Jul 04 1999 - @author Thomas Radke - @desc - The interpolation operator registered with the CCTK - under the name "regular uniform cartesian". - - Interpolates a list of CCTK variables (domain-decomposed - grid functions or arrays) to a list of output arrays - (one-to-one) at a given number of interpolation points - (indicated by their coordinates). The points are located - on a coordinate system which is assumed to be a uniform - cartesian. - @enddesc - - @var GH - @vdesc Pointer to CCTK grid hierarchy - @vtype cGH * - @vio in - @endvar - @var order - @vdesc interpolation order - @vtype int - @vio in - @endvar - @var coord_system_name - @vdesc name of coordinate system to use for interpolation - @vtype const char * - @vio in - @endvar - @var N_points - @vdesc number of points to be interpolated on this processor - @vtype int - @vio in - @endvar - @var N_input_arrays - @vdesc number of input arrays (given by their indices) - to interpolate from - @vtype int - @vio in - @endvar - @var N_output_arrays - @vdesc number of output arrays to interpolate to - @vtype int - @vio in - @endvar - @var interp_coord_arrays - @vdesc coordinates of points to interpolate at - @vtype void *const [N_dims] - @vio in - @endvar - @var interp_coord_array_types - @vdesc CCTK data type of coordinate arrays - @vtype int [N_dims] - @vio in - @endvar - @var input_arrays - @vdesc list of input arrays to interpolate on - @vtype void *[N_input_arrays] - @vio in - @endvar - @var input_array_types - @vdesc CCTK data types of input arrays - @vtype int [N_input_arrays] - @vio in - @endvar - @var output_arrays - @vdesc list of output arrays to interpolate to - @vtype void *const [N_output_arrays] - @vio out - @endvar - @var output_array_types - @vdesc CCTK data types of output arrays - @vtype int [N_output_arrays] - @vio in - @endvar - - @returntype int - @returndesc - 0 - successful interpolation - -1 - in case of any errors - @endreturndesc -@@*/ -int PUGHInterp_InterpGV (cGH *GH, - int order, - const char *coord_system_name, - int N_points, - int N_input_arrays, - int N_output_arrays, - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const int input_array_indices[], - void *const output_arrays[], - const int output_array_types[]) -{ - int i, nprocs, N_dims, point, array, retval; - CCTK_REAL *interp_local_coords; - CCTK_REAL *origin, *delta; - const CCTK_REAL *const *data; - const void **input_arrays; - const pGH *pughGH; - const pGExtras *extras; - cGroupDynamicData group_data; -#ifdef CCTK_MPI - int offset, proc, type, maxtype, N_local_points; - void **local_output_arrays; - pughInterpGH *myGH; - type_desc_t *this, *type_desc; -#endif - - - /* get dimensionality of the coordinate system */ - N_dims = CCTK_CoordSystemDim (coord_system_name); - if (N_dims <= 0) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Cannot get dimensions of coordinate system '%s'", - coord_system_name); - return (-1); - } - - /* check other arguments */ - retval = CheckArguments (GH, N_dims, N_points, N_input_arrays, - N_output_arrays, interp_coord_array_types); - if (retval <= 0) - { - return (retval); - } - - /* get extension handle for PUGH */ - pughGH = CCTK_GHExtension (GH, "PUGH"); - - /* get the extras pointer of the first coordinate - This is used later on to verify the layout of the input arrays as well - as for mapping points to processors. */ - i = CCTK_CoordIndex (1, NULL, coord_system_name); - extras = ((const pGA *) pughGH->variables[i][0])->extras; - - /* get dimensions, origin, and delta of the processor-local grid */ - /* NOTE: getting the dimensions should be a flesh routine as well - for now we get the dimensions of every coordinate and take the - i'th element - this is inconsistent !! */ - origin = malloc (2 * N_dims * sizeof (CCTK_REAL)); - delta = origin + N_dims; - input_arrays = malloc (N_input_arrays * sizeof (void *)); - for (i = 0; i < N_dims; i++) - { - CCTK_CoordLocalRange (GH, &origin[i], &delta[i], i + 1, NULL, - coord_system_name); - delta[i] = (delta[i] - origin[i]) / extras->lnsize[i]; - } - - /* check that the input arrays dimensions match the coordinate system - (but their dimensionality can be less) */ - for (array = 0; array < N_input_arrays; array++) - { - if (CCTK_GroupDynamicData (GH, - CCTK_GroupIndexFromVarI(input_array_indices[array]), - &group_data) < 0) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Invalid input array index %d", - input_array_indices[array]); - retval = -1; - continue; - } - - if (group_data.dim > N_dims) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Input array variable with index %d has more dimensions " - "than coordinate system '%s'", - input_array_indices[array], coord_system_name); - retval = -1; - continue; - } - - if (memcmp (group_data.lsh, extras->lnsize, group_data.dim * sizeof (int))) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Dimensions of input array variable with index %d " - "doesn't match with coordinate system '%s'", - input_array_indices[array], coord_system_name); - retval = -1; - } - - /* get the data pointer to the input array (use current timelevel) */ - input_arrays[array] = CCTK_VarDataPtrI (GH, 0, input_array_indices[array]); - } - if (retval >= 0) - { - /* check for out-of-bounds points - This check was originally in the local interpolator code but was disabled - there and moved up here instead. Now local arrays can have out-of-bounds - points, grid arrays cannot. */ - retval = CheckOutOfBounds (GH, coord_system_name, order, N_dims, N_points, - extras->nsize, - (const CCTK_REAL *const *) interp_coord_arrays); - } - - if (retval < 0) - { - free (input_arrays); - free (origin); - return (retval); - } - - /* single-processor case is easy: no communication or buffering, just direct - interpolation of interp_coord_arrays from input_arrays into output_arrays */ - nprocs = CCTK_nProcs (GH); - if (nprocs == 1) - { - /* sort the individual interpolation coordinate arrays into a single one */ - interp_local_coords = malloc (N_dims * N_points * sizeof (CCTK_REAL)); - data = (const CCTK_REAL *const *) interp_coord_arrays; - for (point = 0; point < N_points; point++) - { - for (i = 0; i < N_dims; i++) - { - *interp_local_coords++ = data[i][point]; - } - } - interp_local_coords -= N_dims * N_points; - - /* call the interpolator function */ - retval = PUGHInterp_Interpolate (order, - N_points, N_dims, N_output_arrays, - extras->lnsize, interp_local_coords, - origin, delta, - output_array_types, input_arrays, - output_array_types, output_arrays); - - /* free allocated resources */ - free (interp_local_coords); - free (input_arrays); - free (origin); - - return (retval); - } - -#ifdef CCTK_MPI - /*** Here follows the multi-processor case: - All processors locate their points to interpolate at - and exchange the coordinates so that every processor gets - those points which it can process locally. - After interpolation the results have to be send back to the - requesting processors. - For both communications MPI_Alltoallv() is used. - - In order to minimize the total number of MPI_Alltoallv() calls - (which are quite expensive) we collect the interpolation results - for all output arrays of the same CCTK data type into a single - communication buffer. That means, after communication the data - needs to be resorted from the buffer into the output arrays. - ***/ - - /* first of all, set up a structure with information of the - CCTK data types we have to deal with */ - - /* get the maximum value of the output array CCTK data types - NOTE: we assume that CCTK data types are defined as consecutive - positive constants starting from zero */ - for (array = maxtype = 0; array < N_output_arrays; array++) - { - if (maxtype < output_array_types[array]) - { - maxtype = output_array_types[array]; - } - } - - /* now allocate an array of structures for all potential types */ - type_desc = calloc (maxtype + 1, sizeof (type_desc_t)); - - /* count the number of arrays of same type - (the N_arrays element was already initialized to zero by calloc() */ - for (array = 0; array < N_output_arrays; array++) - { - type_desc[output_array_types[array]].N_arrays++; - } - - /* fill in the type description information */ - for (type = retval = 0, this = type_desc; type <= maxtype; type++, this++) - { - if (this->N_arrays > 0) - { - /* get the variable type size in bytes */ - this->vtypesize = CCTK_VarTypeSize (type); - if (this->vtypesize <= 0) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "Invalid variable type %d passed, " - "arrays of such type will be skipped during interpolation", - type); - this->N_arrays = 0; - continue; - } - - /* get the MPI data type to use for communicating such a CCTK data type */ - this->mpitype = PUGH_MPIDataType (pughGH, type); - if (this->mpitype == MPI_DATATYPE_NULL) - { - CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, - "No MPI data type defined for variable type %d, " - "arrays of such type will be skipped during interpolation", - type); - this->N_arrays = 0; - continue; - } - - retval++; - } - } - - /* check that there's at least one array with a valid CCTK data type */ - if (retval <= 0) - { - free (input_arrays); - free (origin); - free (type_desc); - return (-1); - } - - /* map the requested points to interpolate at onto the processors - they belong to and gather the coordinates of all the points - that this processor owns - the number of processor-local points is returned in N_local_points, - their coordinates in interp_local_coords */ - retval = GetLocalCoords (GH, N_points, coord_system_name, extras, - (const CCTK_REAL *const *) interp_coord_arrays, - &N_local_points, &interp_local_coords); - if (retval) - { - free (input_arrays); - free (origin); - free (type_desc); - return (retval); - } - - /* allocate contiguous communication buffers for each CCTK data type - holding the local interpolation results from all input arrays - of that type - If there are no points to process on this processor - set the buffer pointer to an invalid but non-NULL value - otherwise we might get trouble with NULL pointers in MPI_Alltoallv () */ - for (type = 0, this = type_desc; type <= maxtype; type++, this++) - { - if (this->N_arrays > 0 && N_local_points > 0) - { - this->sendbuf = malloc (N_local_points * this->N_arrays *this->vtypesize); - this->buf = this->sendbuf; - } - else - { - /* dereferencing such an address should code crash on most systems */ - this->sendbuf = (void *) this->vtypesize; - } - } - - /* get extension handle for interp */ - myGH = CCTK_GHExtension (GH, "PUGHInterp"); - - /* allocate new output_arrays array for local interpolation results - from this processor */ - local_output_arrays = calloc (N_output_arrays, sizeof (void *)); - - /* now, in a loop over all processors, do the interpolation - and put the results in the communication buffer at the proper offset */ - for (proc = 0; proc < nprocs; proc++) - { - for (type = 0, this = type_desc; type <= maxtype; type++, this++) - { - if (this->N_arrays > 0) - { - for (array = 0; array < N_output_arrays; array++) - { - if (output_array_types[array] == type) - { - local_output_arrays[array] = this->buf; - this->buf += myGH->N_points_to[proc] * this->vtypesize; - } - } - } - } - - /* call the interpolation operator to process all points of all - output arrays for this processor */ - PUGHInterp_Interpolate (order, - myGH->N_points_to[proc], N_dims, N_output_arrays, - extras->lnsize, interp_local_coords, origin, delta, - output_array_types, input_arrays, - output_array_types, local_output_arrays); - - /* have to add offset for this processor to coordinates array */ - interp_local_coords += myGH->N_points_to[proc] * N_dims; - - } /* end of loop over all processors */ - - /* don't need these anymore */ - if (N_local_points > 0) - { - interp_local_coords -= N_local_points * N_dims; - free (interp_local_coords); - } - free (local_output_arrays); - free (input_arrays); - free (origin); - - /* now send the interpolation results back to the processors they were - requested from, also receive my own results that were computed - by other processors - Since all the locally computed results are in a single contiguous buffer - we need to call MPI_Alltoall() only once for each CCTK data type. */ - for (type = 0, this = type_desc; type <= maxtype; type++, this++) - { - /* skip unused types */ - if (this->N_arrays <= 0) - { - continue; - } - - /* set up the communication (this is type-independent) */ - myGH->sendcnt[0] = this->N_arrays * myGH->N_points_to[0]; - myGH->recvcnt[0] = this->N_arrays * myGH->N_points_from[0]; - myGH->senddispl[0] = myGH->recvdispl[0] = 0; - for (proc = 1; proc < pughGH->nprocs; proc++) - { - myGH->sendcnt[proc] = this->N_arrays * myGH->N_points_to[proc]; - myGH->recvcnt[proc] = this->N_arrays * myGH->N_points_from[proc]; - myGH->senddispl[proc] = myGH->senddispl[proc-1] + myGH->sendcnt[proc-1]; - myGH->recvdispl[proc] = myGH->recvdispl[proc-1] + myGH->recvcnt[proc-1]; - } - - /* allocate buffer for receiving my own requested points */ - /* avoid NULL pointers here because MPI_Alltoallv() doesn't like it */ - if (N_points > 0) - { - this->recvbuf = malloc (N_points * this->N_arrays * this->vtypesize); - } - else - { - /* access to such a fake address should crash the code on most systems */ - this->recvbuf = (void *) this->vtypesize; - } - - /* now exchange the data for this CCTK data type */ - CACTUS_MPI_ERROR (MPI_Alltoallv (this->sendbuf, myGH->sendcnt, - myGH->senddispl, this->mpitype, - this->recvbuf, myGH->recvcnt, - myGH->recvdispl, this->mpitype, - pughGH->PUGH_COMM_WORLD)); - - /* now that the data is sent we don't need the buffer anymore */ - if (N_local_points > 0) - { - free (this->sendbuf); - } - - /* no sort neccessary if there are no points */ - if (N_points <= 0) - { - continue; - } - - /* Go back from processor-sorted data to input-ordered data. - The creation of the indices array above makes this not so bad. */ - this->buf = this->recvbuf; - for (proc = offset = 0; proc < nprocs; proc++) - { - for (array = 0; array < N_output_arrays; array++) - { - if (output_array_types[array] != type) - { - continue; - } - - /* now do the sorting according to the CCTK data type */ - if (output_array_types[array] == CCTK_VARIABLE_CHAR) - { - SORT_TYPED_ARRAY (CCTK_BYTE); - } - else if (output_array_types[array] == CCTK_VARIABLE_INT) - { - SORT_TYPED_ARRAY (CCTK_INT); - } - else if (output_array_types[array] == CCTK_VARIABLE_REAL) - { - SORT_TYPED_ARRAY (CCTK_REAL); - } - else if (output_array_types[array] == CCTK_VARIABLE_COMPLEX) - { - SORT_TYPED_ARRAY (CCTK_COMPLEX); - } -#ifdef CCTK_REAL4 - else if (output_array_types[array] == CCTK_VARIABLE_REAL4) - { - SORT_TYPED_ARRAY (CCTK_REAL4); - } - else if (output_array_types[array] == CCTK_VARIABLE_COMPLEX8) - { - SORT_TYPED_ARRAY (CCTK_COMPLEX8); - } -#endif -#ifdef CCTK_REAL8 - else if (output_array_types[array] == CCTK_VARIABLE_REAL8) - { - SORT_TYPED_ARRAY (CCTK_REAL8); - } - else if (output_array_types[array] == CCTK_VARIABLE_COMPLEX16) - { - SORT_TYPED_ARRAY (CCTK_COMPLEX16); - } -#endif -#ifdef CCTK_REAL16 - else if (output_array_types[array] == CCTK_VARIABLE_REAL16) - { - SORT_TYPED_ARRAY (CCTK_REAL16); - } - else if (output_array_types[array] == CCTK_VARIABLE_COMPLEX32) - { - SORT_TYPED_ARRAY (CCTK_COMPLEX32); - } -#endif - else - { - CCTK_WARN (0, "Implementation error"); - } - - } /* end of loop over all output arrays */ - - /* advance the offset into the communication receive buffer */ - offset += myGH->N_points_from[proc]; - - } /* end of loop over all processors */ - - /* this communication receive buffer isn't needed anymore */ - free (this->recvbuf); - - } /* end of loop over all types */ - - /* free remaining resources allocated within this run */ - if (myGH->whichproc) - { - free (myGH->whichproc); - myGH->whichproc = NULL; - } - free (type_desc); -#endif /* MPI */ - - return (0); -} - - -/**************************************************************************/ -/* local routines */ -/**************************************************************************/ - -#ifdef CCTK_MPI - -/*@@ - @routine GetLocalCoords - @date Sun Jul 04 1999 - @author Thomas Radke - @desc - Collect the coordinates of all points to be processed locally - into an array coords[N_dims][N_local_points]. - - This means for the single-processor case to sort - - inCoords1[N_points], inCoords2[npoints], ..., inCoords[npoints] - into coords[N_dims][N_points] - - where N_points == N_local_points. - - In the multiprocessor case all processors map their points' coordinates - to the processor that owns this point and exchange this information - via MPI_Alltoall (). - N_local_points is then the number of all processors' points to be - interpolated locally on this processor. - - This routine returns the number of points to be processed locally and - a pointer to the allocated array of their coordinates. - @enddesc - - @var GH - @vdesc Pointer to CCTK grid hierarchy - @vtype const cGH * - @vio in - @endvar - @var N_points - @vdesc number of points to be interpolated on this processor - @vtype int - @vio in - @endvar - @var isGlobal - @vdesc flag indicating that coordinates are global (and need to be - collected over all processors) - @vtype int - @vio in - @endvar - @var N_dims - @vdesc number of coordinate dimensions for each point - @vtype int - @vio in - @endvar - @var coords - @vdesc coordinates of each point to be interpolated on this processor - @vtype CCTK_REAL array of size N_dims - @vio in - @endvar - @var N_local_points - @vdesc number of points to be processed by this processor - @vtype int * - @vio out - @endvar - @var local_coords - @vdesc coordinates of each point to be processed by this processor - @vtype pointer to CCTK_REAL array of dims[N_dims][N_local_points] - @vio out - @endvar - -@@*/ -static int GetLocalCoords (const cGH *GH, - int N_points, - const char *coord_system_name, - const pGExtras *extras, - const CCTK_REAL *const coords[], - int *N_local_points, - CCTK_REAL **local_coords) -{ - int dim, point; - int tmp, proc, nprocs; - pGH *pughGH; - pughInterpGH *myGH; - CCTK_REAL *range_min, *range_max; - CCTK_REAL *origin, *delta; - CCTK_REAL *proc_coords; -#define FUDGE 0.0 - - - /* get GH extension handles for PUGHInterp and PUGH */ - myGH = CCTK_GHExtension (GH, "PUGHInterp"); - pughGH = CCTK_GHExtension (GH, "PUGH"); - - /* This holds the proccessor for *each* of N_points points */ - if (N_points > 0) - { - myGH->whichproc = malloc (2 * N_points * sizeof (int)); - } - else - { - myGH->whichproc = NULL; - } - /* indices[] is used to make the sorting easier - when receiving the output data */ - myGH->indices = myGH->whichproc + N_points; - - /* initialize whichproc with invalid processor number -1 */ - for (point = 0; point < N_points; point++) - { - myGH->whichproc[point] = -1; - } - - /* initialize N_points_from to 0 for counting it up in the following loop */ - nprocs = CCTK_nProcs (GH); - memset (myGH->N_points_from, 0, nprocs * sizeof (CCTK_INT)); - - /* allocate the ranges for my local coordinates */ - range_min = malloc (4 * extras->dim * sizeof (CCTK_REAL)); - range_max = range_min + 1*extras->dim; - origin = range_min + 2*extras->dim; - delta = range_min + 3*extras->dim; - - /* get the global origin and delta of the coordinate system */ - for (dim = 0; dim < extras->dim; dim++) - { - CCTK_CoordRange (GH, &origin[dim], &delta[dim], dim+1, NULL, coord_system_name); - delta[dim] = (delta[dim] - origin[dim]) / (extras->nsize[dim]-1); - } - - /* locate the points to interpolate at */ - for (proc = 0; proc < nprocs; proc++) - { - for (dim = 0; dim < extras->dim; dim++) - { - /* compute the coordinate ranges */ - /* TODO: use bbox instead -- but the bboxes of other processors - are now known */ - int const has_lower = extras->lb[proc][dim] == 0; - int const has_upper = extras->ub[proc][dim] == GH->cctk_gsh[dim]-1; - range_min[dim] = origin[dim] + (extras->lb[proc][dim] + (!has_lower) * (extras->nghostzones[dim]-0.5) - FUDGE)*delta[dim]; - range_max[dim] = origin[dim] + (extras->ub[proc][dim] - (!has_upper) * (extras->nghostzones[dim]-0.5) + FUDGE)*delta[dim]; - } - - /* and now which point will be processed by what processor */ - for (point = 0; point < N_points; point++) - { - /* skip points which have already been located */ - if (myGH->whichproc[point] >= 0) - { - continue; - } - - /* check whether the point belongs to this processor - (must be within min/max in all dimensions) */ - tmp = 0; - for (dim = 0; dim < extras->dim; dim++) - { - if (coords[dim][point] >= range_min[dim] && - coords[dim][point] <= range_max[dim]) - { - tmp++; - } - } - if (tmp == extras->dim) - { - myGH->whichproc[point] = proc; - myGH->N_points_from[proc]++; - } - } - } - /* don't need this anymore */ - free (range_min); - - /* make sure that all points could be mapped onto a processor */ - for (point = tmp = 0; point < N_points; point++) - { - if (myGH->whichproc[point] < 0) - { - int i; - char *msg = malloc (80 + extras->dim*20); - - - sprintf (msg, "Unable to locate point %d [%f", - point, coords[0][point]); - for (i = 1; i < extras->dim; i++) - { - sprintf (msg, "%s %f", msg, coords[i][point]); - } - sprintf (msg, "%s]", msg); - CCTK_WARN (1, msg); - free (msg); - tmp = 1; /* mark as error */ - } - } - if (tmp) - { - if (myGH->whichproc) - { - free (myGH->whichproc); - myGH->whichproc = NULL; - } - return (-1); - } - - /* Now we want to resolve the N_points_from[]. Currently this is - the form of ( in 2 proc mode ) - P1: Num from P1 NFP2 - P2: NFP1 NFP2 - - and this needs to become - P1: Num to P1 NTP2 - P2: NTP1 NTP1 - - Since NTP1 = NFP2 (and this works in more proc mode too) - this is an all-to-all communication. - */ - CACTUS_MPI_ERROR (MPI_Alltoall (myGH->N_points_from, 1, PUGH_MPI_INT, - myGH->N_points_to, 1, PUGH_MPI_INT, - pughGH->PUGH_COMM_WORLD)); - -#ifdef PUGHINTERP_DEBUG - for (proc = 0; proc < nprocs; proc++) - { - printf ("processor %d <-> %d From: %d To: %d\n", - CCTK_MyProc (GH), proc, myGH->N_points_from[proc], - myGH->N_points_to[proc]); - } -#endif - - /* Great. Now we know how many to expect from each processor, - and how many to send to each processor. So first we have - to send the locations to the processors which hold our data. - This means I send coords[dim][point] to whichproc[point]. - I have N_points_from[proc] to send to each processor. - */ - - /* This is backwards in the broadcast location; the number of points - we are getting is how many everyone else is sending to us, - eg, N_points_to, not how many we get back from everyone else, - eg, N_points_from. The number we are sending, of course, is - all of our locations, eg, N_points */ - *N_local_points = 0; - for (proc = 0; proc < nprocs; proc++) - { - *N_local_points += myGH->N_points_to[proc]; - } - -#ifdef PUGHINTERP_DEBUG - printf ("processor %d gets %d points in total\n", - CCTK_MyProc (GH), *N_local_points); -#endif - - /* allocate the local coordinates array (sorted in processor order) - and the resulting coordinates array that I have to process */ - proc_coords = malloc (extras->dim * N_points * sizeof (CCTK_REAL)); - *local_coords = malloc (extras->dim * *N_local_points * sizeof (CCTK_REAL)); - - /* now sort my own coordinates as tupels of [extras->dim] */ - for (proc = tmp = 0; proc < nprocs; proc++) - { - for (point = 0; point < N_points; point++) - { - if (myGH->whichproc[point] == proc) - { - for (dim = 0; dim < extras->dim; dim++) - { - *proc_coords++ = coords[dim][point]; - } - myGH->indices[tmp++] = point; - } - } - } - proc_coords -= tmp * extras->dim; - - /* So load up the send and recv stuff */ - /* Send extras->dim elements per data point */ - myGH->sendcnt[0] = extras->dim * myGH->N_points_from[0]; - myGH->recvcnt[0] = extras->dim * myGH->N_points_to[0]; - myGH->senddispl[0] = myGH->recvdispl[0] = 0; - for (proc = 1; proc < nprocs; proc++) - { - myGH->sendcnt[proc] = extras->dim * myGH->N_points_from[proc]; - myGH->recvcnt[proc] = extras->dim * myGH->N_points_to[proc]; - myGH->senddispl[proc] = myGH->senddispl[proc-1] + myGH->sendcnt[proc-1]; - myGH->recvdispl[proc] = myGH->recvdispl[proc-1] + myGH->recvcnt[proc-1]; - } - - /* Great, and now exchange the coordinates and collect the ones - that I have to process in *local_coords[] */ - CACTUS_MPI_ERROR (MPI_Alltoallv (proc_coords, myGH->sendcnt, - myGH->senddispl, PUGH_MPI_REAL, - *local_coords, myGH->recvcnt, - myGH->recvdispl, PUGH_MPI_REAL, - pughGH->PUGH_COMM_WORLD)); - - /* don't need this anymore */ - free (proc_coords); - - return (0); -} -#endif /* CCTK_MPI */ - - - /*@@ - @routine CheckArguments - @date Thu 25 Jan 2001 - @author Thomas Radke - @desc - Checks the interpolation arguments passed in via - the flesh's general interpolation calling interface - - This routine also verifies that the parameters meet - the limitations of PUGHInterp's interpolation operators. - @enddesc - - @var GH - @vdesc Pointer to CCTK grid hierarchy - @vtype const cGH * - @vio in - @endvar - @var N_dims - @vdesc dimensionality of the underlying grid - @vtype int - @vio in - @endvar - @var N_points - @vdesc number of points to interpolate at - @vtype int - @vio in - @endvar - @var N_input_arrays - @vdesc number of passed input arrays - @vtype int - @vio in - @endvar - @var N_output_arrays - @vdesc number of passed input arrays - @vtype int - @vio in - @endvar - @var interp_coord_array_types - @vdesc types of passed coordinates to interpolate at - @vtype int [N_dims] - @vio in - @endvar - - @returntype int - @returndesc - +1 for success - 0 for success but nothing to do - -1 for failure (wrong parameters passed or limitations not met) - @endreturndesc -@@*/ -static int CheckArguments (const cGH *GH, - int N_dims, - int N_points, - int N_input_arrays, - int N_output_arrays, - const int interp_coord_array_types[]) -{ - int i; - - - /* check for invalid arguments */ - if (N_dims < 0 || N_points < 0 || N_input_arrays < 0 || N_output_arrays < 0) - { - return (-1); - } - - /* check if there's anything to do at all */ - /* NOTE: N_points can be 0 in a collective call */ - if (N_dims == 0 || (CCTK_nProcs (GH) == 1 && N_points == 0) || - N_input_arrays == 0 || N_output_arrays == 0) - { - return (0); - } - - /* for now we can only deal with coordinates of type CCTK_REAL */ - for (i = 0; i < N_dims; i++) - { - if (interp_coord_array_types[i] != CCTK_VARIABLE_REAL) - { - CCTK_WARN (1, "Interpolation coordinates must be of type CCTK_REAL"); - return (-1); - } - } - - /* PUGHInterp's interpolation operators compute one output array - per input array */ - if (N_input_arrays != N_output_arrays) - { - CCTK_WARN (1, "Number of input arrays must match number of output arrays"); - return (-1); - } - - return (1); -} - - -static int CheckOutOfBounds (const cGH *GH, const char *coord_system_name, - int order, int N_dims, int N_points, - const int *dims, const CCTK_REAL *const *coords) -{ - int i, p, point, out_of_bounds, retval; - CCTK_REAL *origin, *delta, *delta_inv, *below; - char *msg; - - - msg = malloc (100 + N_dims*(10 + 4*30)); - origin = malloc (4 * N_dims * sizeof (CCTK_REAL)); - delta = origin + 1*N_dims; - delta_inv = origin + 2*N_dims; - below = origin + 3*N_dims; - - /* get the global origin and delta of the coordinate system */ - for (i = 0; i < N_dims; i++) - { - CCTK_CoordRange (GH, &origin[i], &delta[i], i+1, NULL, coord_system_name); - delta[i] = (delta[i] - origin[i]) / (dims[i]-1); - - /* avoid expensive divisions by delta later on */ - delta_inv[i] = 1.0 / delta[i]; - } - - retval = 0; - - for (p = 0; p < N_points; p++) - { - /* reset the out-of-bounds flag */ - out_of_bounds = 0; - - /* loop over all dimensions */ - for (i = 0; i < N_dims; i++) - { - /* grid point of the lower-left stencil point */ - point = floor ((coords[i][p] - origin[i]) * delta_inv[i] - - 0.5 * (order - 1)); - - /* test bounds */ - out_of_bounds |= point < 0 || point+order >= dims[i]; - - /* physical coordinate of that grid point */ - below[i] = origin[i] + point * delta[i]; - } - - /* check bounds */ - if (out_of_bounds) - { - /* put all information into a single message string for output */ - sprintf (msg, "Interpolation stencil out of bounds at interpolation " - "coordinate [%f", (double) coords[0][p]); - for (i = 1; i < N_dims; i++) - { - sprintf (msg, "%s, %f", msg, (double) coords[i][p]); - } - sprintf (msg, "%s]\nrange would be min/max [%f / %f", msg, - (double) below[0], (double) (below[0] + order*delta[0])); - for (i = 1; i < N_dims; i++) - { - sprintf (msg, "%s, %f / %f", msg, - (double) below[i], (double) (below[i] + order*delta[i])); - } - sprintf (msg, "%s]\ngrid is min/max [%f / %f", msg, - (double) origin[0], (double) (origin[0] + (dims[0]-1)*delta[0])); - for (i = 1; i < N_dims; i++) - { - sprintf (msg, "%s, %f / %f", msg, - (double)origin[i], (double)(origin[i] + (dims[i]-1)*delta[i])); - } - sprintf (msg, "%s]", msg); - CCTK_WARN (1, msg); - - retval--; - } - } - - /* free allocated resources */ - free (origin); - free (msg); - - return (retval); -} diff --git a/src/Startup.c b/src/Startup.c index f5ff5f2..3cb3226 100644 --- a/src/Startup.c +++ b/src/Startup.c @@ -11,6 +11,7 @@ #include #include "cctk.h" +#include "util_Table.h" #include "cctk_Interp.h" #include "pughInterpGH.h" @@ -61,6 +62,18 @@ static int InterpGV_3rdOrder (cGH *GH, const int in_array_indices[], void *const out_arrays[], const int out_array_types[]); +static int PUGHInterp_InterpGV (cGH *GH, + int order, + const char *local_interpolator, + const char *coord_system_name, + int N_points, + int N_input_arrays, + int N_output_arrays, + const void *const interp_coord_arrays[], + const int interp_coord_array_types[], + const int input_array_indices[], + void *const output_arrays[], + const int output_array_types[]); /*@@ @@ -168,7 +181,8 @@ static int InterpGV_1stOrder (cGH *GH, void *const out_arrays[], const int out_array_types[]) { - return (PUGHInterp_InterpGV (GH, 1, coord_system, num_points, + return (PUGHInterp_InterpGV (GH, 1, "first-order uniform cartesian", + coord_system, num_points, num_in_array_indices, num_out_arrays, interp_coord_arrays, interp_coord_array_types, in_array_indices, out_arrays, out_array_types)); @@ -186,7 +200,8 @@ static int InterpGV_2ndOrder (cGH *GH, void *const out_arrays[], const int out_array_types[]) { - return (PUGHInterp_InterpGV (GH, 2, coord_system, num_points, + return (PUGHInterp_InterpGV (GH, 2, "second-order uniform cartesian", + coord_system, num_points, num_in_array_indices, num_out_arrays, interp_coord_arrays, interp_coord_array_types, in_array_indices, out_arrays, out_array_types)); @@ -204,8 +219,80 @@ static int InterpGV_3rdOrder (cGH *GH, void *const out_arrays[], const int out_array_types[]) { - return (PUGHInterp_InterpGV (GH, 3, coord_system, num_points, + return (PUGHInterp_InterpGV (GH, 3, "third-order uniform cartesian", + coord_system, num_points, num_in_array_indices, num_out_arrays, interp_coord_arrays, interp_coord_array_types, in_array_indices, out_arrays, out_array_types)); } + + +static int PUGHInterp_InterpGV (cGH *GH, + int order, + const char *local_interpolator, + const char *coord_system_name, + int N_points, + int N_input_arrays, + int N_output_arrays, + const void *const interp_coord_arrays[], + const int interp_coord_array_types[], + const int input_array_indices[], + void *const output_arrays[], + const int output_array_types[]) +{ + int i, N_dims, retval; + int local_interp_handle, param_table_handle, coord_system_handle; + char table_string[64]; + CCTK_INT *_input_array_indices, *_output_array_types; + struct + { + const cGH *const_GH; + cGH *non_const_GH; + } _GH; + + + _GH.non_const_GH = GH; + + /* get dimensionality of the coordinate system */ + N_dims = CCTK_CoordSystemDim (coord_system_name); + if (N_dims <= 0) + { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot get dimensions of coordinate system '%s'", + coord_system_name); + return (-1); + } + coord_system_handle = CCTK_CoordSystemHandle (coord_system_name); + local_interp_handle = CCTK_InterpHandle (local_interpolator); + + /* turn native datatype arrays into CCTK_ datatypes */ + _input_array_indices = malloc (2 * N_dims * sizeof (CCTK_INT)); + _output_array_types = _input_array_indices + N_dims; + for (i = 0; i < N_dims; i++) + { + _input_array_indices[i] = input_array_indices[i]; + _output_array_types[i] = output_array_types[i]; + } + + /* set up the parameter table (only supported option is "order") */ + sprintf (table_string, "order = %d", order); + param_table_handle = Util_TableCreateFromString (table_string); + + retval = PUGHInterp_InterpGridArrays (_GH.const_GH, N_dims, + local_interp_handle, param_table_handle, + coord_system_handle, + N_points, + interp_coord_array_types[0], + interp_coord_arrays, + N_input_arrays, + _input_array_indices, + N_output_arrays, + _output_array_types, + output_arrays); + + /* release resources */ + Util_TableDestroy (param_table_handle); + free (_input_array_indices); + + return (retval); +} diff --git a/src/make.code.defn b/src/make.code.defn index eba2a31..f0bfb0c 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -2,4 +2,4 @@ # $Header$ # Source files in this directory -SRCS = Startup.c Operator.c Interpolate.c InterpGridArrays.c +SRCS = Startup.c InterpGridArrays.c diff --git a/src/pughInterpGH.h b/src/pughInterpGH.h index e25b4f6..871cf22 100644 --- a/src/pughInterpGH.h +++ b/src/pughInterpGH.h @@ -47,33 +47,6 @@ int PUGHInterp_InterpGridArrays (const cGH *GH, const CCTK_INT output_array_types[], void *const output_arrays[]); -/* prototypes of interpolation operators to be registered */ -int PUGHInterp_InterpGV (cGH *GH, - int order, - const char *coord_system, - int num_points, - int num_in_array_indices, - int num_out_arrays, - const void *const interp_coord_arrays[], - const int interp_coord_array_types[], - const int in_array_indices[], - void *const out_arrays[], - const int out_array_types[]); - -/* prototypes of routines used internally */ -int PUGHInterp_Interpolate (int order, - int num_points, - int num_dims, - int num_arrays, - const int dims[ /* num_dims */ ], - const CCTK_REAL coord[ /* num_dims*num_points */ ], - const CCTK_REAL origin[ /* num_dims */ ], - const CCTK_REAL delta[ /* num_dims */ ], - const int in_types[ /* num_arrays */ ], - const void *const in_arrays[ /* num_arrays */ ], - const int out_types[ /* num_arrays */ ], - void *const out_arrays[ /* num_arrays */ ]); - #ifdef __cplusplus } // extern "C" #endif -- cgit v1.2.3