diff options
Diffstat (limited to 'src/Operator.c')
| -rw-r--r-- | src/Operator.c | 168 |
1 files changed, 126 insertions, 42 deletions
diff --git a/src/Operator.c b/src/Operator.c index 8064fd3..0ae230d 100644 --- a/src/Operator.c +++ b/src/Operator.c @@ -19,6 +19,7 @@ #include <stdlib.h> #include <string.h> +#include <math.h> /* floor(3) */ #include "cctk.h" #include "CactusPUGH/PUGH/src/include/pugh.h" @@ -70,6 +71,9 @@ static int CheckArguments (const cGH *GH, 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, @@ -177,26 +181,18 @@ int PUGHInterp_InterpGV (cGH *GH, void *const output_arrays[], const int output_array_types[]) { - const CCTK_REAL **data; + int i, nprocs, N_dims, point, array, retval; CCTK_REAL *interp_local_coords; CCTK_REAL *origin, *delta; - int nprocs; - int dim, N_dims, *dims; - int array; - int point; - int retval; + const CCTK_REAL **data; const void **input_arrays; - pGH *pughGH; - pGExtras *extras; + const pGH *pughGH; + const pGExtras *extras; cGroupDynamicData group_data; #ifdef CCTK_MPI - int N_local_points; - int offset; - int proc; - int type; + int offset, proc, type, maxtype, N_local_points; void **local_output_arrays; pughInterpGH *myGH; - int maxtype; type_desc_t *this, *type_desc; #endif @@ -219,29 +215,29 @@ int PUGHInterp_InterpGV (cGH *GH, 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 !! */ - dims = malloc (2 * N_dims * sizeof (int)); origin = malloc (2 * N_dims * sizeof (CCTK_REAL)); delta = origin + N_dims; input_arrays = malloc (N_input_arrays * sizeof (void *)); - for (dim = 0; dim < N_dims; dim++) + for (i = 0; i < N_dims; i++) { - CCTK_GrouplshVI (GH, N_dims, dims + N_dims, - CCTK_CoordIndex (dim + 1, NULL, coord_system_name)); - dims[dim] = dims[dim + N_dims]; - - CCTK_CoordLocalRange (GH, &origin[dim], &delta[dim], dim + 1, NULL, + CCTK_CoordLocalRange (GH, &origin[i], &delta[i], i + 1, NULL, coord_system_name); - delta[dim] = (delta[dim] - origin[dim]) / dims[dim]; + delta[i] = (delta[i] - origin[i]) / extras->lnsize[i]; } - /* get extension handle for PUGH */ - pughGH = CCTK_GHExtension (GH, "PUGH"); - extras = NULL; - /* check that the input arrays dimensions match the coordinate system (but their dimensionality can be less) */ for (array = 0; array < N_input_arrays; array++) @@ -267,7 +263,7 @@ int PUGHInterp_InterpGV (cGH *GH, continue; } - if (memcmp (group_data.lsh, dims, group_data.dim * sizeof (int))) + 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 " @@ -276,21 +272,24 @@ int PUGHInterp_InterpGV (cGH *GH, retval = -1; } - /* the extras pointer will refer to the first input array with the - highest dimension (need this information later on) */ - if (! extras || extras->dim < group_data.dim) - { - extras = ((pGA *) pughGH->variables[input_array_indices[array]][0])->extras; - } - /* 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 **) interp_coord_arrays); + } + if (retval < 0) { free (input_arrays); free (origin); - free (dims); return (retval); } @@ -304,9 +303,9 @@ int PUGHInterp_InterpGV (cGH *GH, data = (const CCTK_REAL **) interp_coord_arrays; for (point = 0; point < N_points; point++) { - for (dim = 0; dim < N_dims; dim++) + for (i = 0; i < N_dims; i++) { - *interp_local_coords++ = data[dim][point]; + *interp_local_coords++ = data[i][point]; } } interp_local_coords -= N_dims * N_points; @@ -314,7 +313,7 @@ int PUGHInterp_InterpGV (cGH *GH, /* call the interpolator function */ retval = PUGHInterp_Interpolate (order, N_points, N_dims, N_output_arrays, - dims, interp_local_coords, + extras->lnsize, interp_local_coords, origin, delta, output_array_types, input_arrays, output_array_types, output_arrays); @@ -323,7 +322,6 @@ int PUGHInterp_InterpGV (cGH *GH, free (interp_local_coords); free (input_arrays); free (origin); - free (dims); return (retval); } @@ -405,7 +403,6 @@ int PUGHInterp_InterpGV (cGH *GH, if (retval <= 0) { free (input_arrays); - free (dims); free (origin); free (type_desc); return (-1); @@ -423,7 +420,6 @@ int PUGHInterp_InterpGV (cGH *GH, { free (input_arrays); free (origin); - free (dims); free (type_desc); return (retval); } @@ -478,7 +474,7 @@ int PUGHInterp_InterpGV (cGH *GH, output arrays for this processor */ PUGHInterp_Interpolate (order, myGH->N_points_to[proc], N_dims, N_output_arrays, - dims, interp_local_coords, origin, delta, + extras->lnsize, interp_local_coords, origin, delta, output_array_types, input_arrays, output_array_types, local_output_arrays); @@ -496,7 +492,6 @@ int PUGHInterp_InterpGV (cGH *GH, free (local_output_arrays); free (input_arrays); free (origin); - free (dims); /* now send the interpolation results back to the processors they were requested from, also receive my own results that were computed @@ -1038,3 +1033,92 @@ static int CheckArguments (const cGH *GH, 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, *offset; + char *msg; + + + msg = malloc (100 + N_dims*(10 + 4*30)); + origin = malloc (5 * N_dims * sizeof (CCTK_REAL)); + delta = origin + 1*N_dims; + delta_inv = origin + 2*N_dims; + below = origin + 3*N_dims; + offset = origin + 4*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++) + { + /* closest grid point for stencil */ + 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]; + + /* offset from that grid point, in fractions of grid points */ + offset[i] = (coords[i][p] - below[i]) * delta_inv[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 grid point [%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] + offset[0])); + for (i = 1; i < N_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 < 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); +} |