aboutsummaryrefslogtreecommitdiff
path: root/src/Operator.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/Operator.c')
-rw-r--r--src/Operator.c1293
1 files changed, 1293 insertions, 0 deletions
diff --git a/src/Operator.c b/src/Operator.c
new file mode 100644
index 0000000..76a71f5
--- /dev/null
+++ b/src/Operator.c
@@ -0,0 +1,1293 @@
+ /*@@
+ @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 <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "CactusPUGH/PUGH/src/include/pugh.h"
+#include "pughInterpGH.h"
+
+/* the rcs ID and its dummy function to use it */
+static 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, source, dest, points, index_array, offset)\
+ { \
+ int i; \
+ cctk_type *src, *dst; \
+ \
+ \
+ src = (cctk_type *) source; \
+ dst = (cctk_type *) dest; \
+ for (i = 0; i < points; i++) \
+ { \
+ dst[index_array[i + offset]] = *src++; \
+ } \
+ source = (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 mpi_type; /* corresponding MPI datatype */
+ int num_arrays; /* number of in/out arrays */
+ int array; /* iterator for num_arrays */
+ void *send_buffer; /* communication send buffer for this type */
+ void *recv_buffer; /* communication receive buffer for this type */
+ char *buffer; /* work pointer for send_buffer */
+} type_desc_t;
+#endif
+
+
+/* prototypes of routines defined in this source file */
+static int PUGHInterp_CheckArguments (cGH *GH,
+ int num_dims,
+ int num_points,
+ int num_in_arrays,
+ int num_out_arrays,
+ int interp_coord_array_types[]);
+#ifdef CCTK_MPI
+static int GetLocalCoords (cGH *GH,
+ int num_points,
+ const char *coord_system,
+ pGExtras *extras,
+ CCTK_REAL *coords[],
+ int *num_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
+ @vdesc name of coordinate system to use for interpolation
+ @vtype const char *
+ @vio in
+ @endvar
+ @var num_points
+ @vdesc number of points to be interpolated on this processor
+ @vtype int
+ @vio in
+ @endvar
+ @var num_in_array_indices
+ @vdesc number of input arrays (given by their indices)
+ to interpolate from
+ @vtype int
+ @vio in
+ @endvar
+ @var num_out_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 *[num_dims]
+ @vio in
+ @endvar
+ @var interp_coord_array_types
+ @vdesc CCTK data type of coordinate arrays
+ @vtype int [num_dims]
+ @vio in
+ @endvar
+ @var in_arrays
+ @vdesc list of input arrays to interpolate on
+ @vtype void *[num_in_arrays]
+ @vio in
+ @endvar
+ @var in_array_types
+ @vdesc CCTK data types of input arrays
+ @vtype int [num_in_arrays]
+ @vio in
+ @endvar
+ @var out_arrays
+ @vdesc list of output arrays to interpolate to
+ @vtype void *[num_out_arrays]
+ @vio out
+ @endvar
+ @var out_array_types
+ @vdesc CCTK data types of output arrays
+ @vtype int [num_out_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,
+ int num_points,
+ int num_in_array_indices,
+ int num_out_arrays,
+ void *interp_coord_arrays[],
+ int interp_coord_array_types[],
+ int in_array_indices[],
+ void *out_arrays[],
+ int out_array_types[])
+{
+ CCTK_REAL **data, *interp_local_coords;
+ CCTK_REAL *origin, *delta;
+ int nprocs;
+ int timelevel;
+ int dim, num_dims, *dims;
+ int array;
+ int point;
+ int retval;
+ void **in_arrays;
+ pGH *pughGH;
+ pGExtras *extras;
+ cGroupDynamicData group_data;
+#ifdef CCTK_MPI
+ int num_local_points;
+ int offset;
+ int proc;
+ int type;
+ void **local_out_arrays;
+ pughInterpGH *myGH;
+ int max_type;
+ type_desc_t *type_desc;
+#endif
+
+
+ /* get dimensionality of the coordinate system */
+ num_dims = CCTK_CoordSystemDim (coord_system);
+ if (num_dims <= 0)
+ {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Cannot get dimensions of coordinate system '%s'",
+ coord_system);
+ return (-1);
+ }
+
+ /* check other arguments */
+ retval = PUGHInterp_CheckArguments (GH, num_dims, num_points,
+ num_in_array_indices, num_out_arrays,
+ interp_coord_array_types);
+ if (retval <= 0)
+ {
+ return (retval);
+ }
+
+ /* 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 = (int *) malloc (2 * num_dims * sizeof (int));
+ origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL));
+ delta = origin + num_dims;
+ in_arrays = malloc (num_in_array_indices * sizeof (void *));
+ for (dim = 0; dim < num_dims; dim++)
+ {
+ CCTK_GrouplshVI (GH, num_dims, dims + num_dims,
+ CCTK_CoordIndex (dim + 1, NULL, coord_system));
+ dims[dim] = dims[dim + num_dims];
+
+ CCTK_CoordLocalRange (GH, &origin[dim], &delta[dim], dim + 1, NULL,
+ coord_system);
+ delta[dim] = (delta[dim] - origin[dim]) / dims[dim];
+ }
+
+ /* get extension handle for PUGH */
+ pughGH = (pGH *) 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 < num_in_array_indices; array++)
+ {
+ if (CCTK_GroupDynamicData (GH,
+ CCTK_GroupIndexFromVarI(in_array_indices[array]),
+ &group_data) < 0)
+ {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Invalid input array index %d",
+ in_array_indices[array]);
+ retval = -1;
+ continue;
+ }
+
+ if (group_data.dim > num_dims)
+ {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Input array variable with index %d has more dimensions "
+ "than coordinate system '%s'",
+ in_array_indices[array], coord_system);
+ retval = -1;
+ continue;
+ }
+
+ if (memcmp (group_data.lsh, dims, 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'",
+ in_array_indices[array], coord_system);
+ 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[in_array_indices[array]][0])->extras;
+ }
+
+ /* get the data pointer to the input array (use current timelevel) */
+ timelevel = CCTK_NumTimeLevelsFromVarI (in_array_indices[array]) - 1;
+ if (timelevel > 0)
+ {
+ timelevel--;
+ }
+ in_arrays[array] = CCTK_VarDataPtrI(GH, timelevel, in_array_indices[array]);
+ }
+ if (retval < 0)
+ {
+ free (in_arrays);
+ free (origin);
+ free (dims);
+ return (retval);
+ }
+
+ /* single-processor case is easy: no communication or buffering, just direct
+ interpolation of interp_coord_arrays from in_arrays into out_arrays */
+ nprocs = CCTK_nProcs (GH);
+ if (nprocs == 1)
+ {
+ /* sort the individual interpolation coordinate arrays into a single one */
+ interp_local_coords = (CCTK_REAL *) malloc (num_dims * num_points *
+ sizeof (CCTK_REAL));
+ data = (CCTK_REAL **) interp_coord_arrays;
+ for (point = 0; point < num_points; point++)
+ {
+ for (dim = 0; dim < num_dims; dim++)
+ {
+ *interp_local_coords++ = data[dim][point];
+ }
+ }
+ interp_local_coords -= num_dims * num_points;
+
+ /* call the interpolator function */
+ retval = PUGHInterp_Interpolate (order,
+ num_points, num_dims, num_out_arrays,
+ dims, interp_local_coords,
+ origin, delta,
+ out_array_types, in_arrays,
+ out_array_types, out_arrays);
+
+ /* free allocated resources */
+ free (interp_local_coords);
+ free (in_arrays);
+ free (origin);
+ free (dims);
+
+ 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 */
+ max_type = 0;
+ for (array = 0; array < num_out_arrays; array++)
+ {
+ if (max_type < out_array_types[array])
+ {
+ max_type = out_array_types[array];
+ }
+ }
+
+ /* now allocate an array of structures for all potential types */
+ type_desc = (type_desc_t *) calloc (max_type + 1, sizeof (type_desc_t));
+
+ /* count the number of arrays of same type
+ (the num_arrays element was already initialized to zero by calloc() */
+ for (array = 0; array < num_out_arrays; array++)
+ {
+ type_desc[out_array_types[array]].num_arrays++;
+ }
+
+ /* fill in the type description information */
+ retval = 0;
+ for (type = 0; type <= max_type; type++)
+ {
+ if (type_desc[type].num_arrays > 0)
+ {
+ /* get the variable type size in bytes */
+ type_desc[type].vtypesize = CCTK_VarTypeSize (type);
+ if (type_desc[type].vtypesize <= 0)
+ {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Invalid variable type %d passed, "
+ "arrays of such type will be skipped during interpolation",
+ type);
+ type_desc[type].num_arrays = 0;
+ continue;
+ }
+
+ /* get the MPI data type to use for communicating such a CCTK data type */
+ type_desc[type].mpi_type = PUGH_MPIDataType (pughGH, type);
+ if (type_desc[type].mpi_type == 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);
+ type_desc[type].num_arrays = 0;
+ continue;
+ }
+
+ retval++;
+ }
+ }
+
+ /* check that there's at least one array with a valid CCTK data type */
+ if (retval <= 0)
+ {
+ free (in_arrays);
+ free (dims);
+ 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 num_local_points,
+ their coordinates in interp_local_coords */
+ retval = GetLocalCoords (GH, num_points, coord_system, extras,
+ (CCTK_REAL **) interp_coord_arrays,
+ &num_local_points, &interp_local_coords);
+ if (retval)
+ {
+ free (in_arrays);
+ free (origin);
+ free (dims);
+ 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; type <= max_type; type++)
+ {
+ if (type_desc[type].num_arrays > 0 && num_local_points > 0)
+ {
+ type_desc[type].send_buffer = malloc (num_local_points *
+ type_desc[type].num_arrays *
+ type_desc[type].vtypesize);
+ type_desc[type].buffer = (char *) type_desc[type].send_buffer;
+ type_desc[type].array = 0;
+ }
+ else
+ {
+ /* dereferencing such an address should code crash on most systems */
+ type_desc[type].send_buffer = (void *) 1;
+ }
+ }
+
+ /* get extension handle for interp */
+ myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp");
+
+ /* allocate new out_arrays array for local interpolation results
+ from this processor */
+ local_out_arrays = calloc (num_out_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; type <= max_type; type++)
+ {
+ if (type_desc[type].num_arrays > 0)
+ {
+ for (array = 0; array < num_out_arrays; array++)
+ {
+ if (out_array_types[array] == type)
+ {
+ local_out_arrays[array] = type_desc[type].buffer;
+ type_desc[type].buffer += myGH->num_points_to[proc] *
+ type_desc[type].vtypesize;
+ }
+ }
+ }
+ }
+
+ /* call the interpolation operator to process all points of all
+ output arrays for this processor */
+ PUGHInterp_Interpolate (order,
+ myGH->num_points_to[proc], num_dims, num_out_arrays,
+ dims, interp_local_coords, origin, delta,
+ out_array_types, in_arrays,
+ out_array_types, local_out_arrays);
+
+ /* have to add offset for this processor to coordinates array */
+ interp_local_coords += myGH->num_points_to[proc] * num_dims;
+
+ } /* end of loop over all processors */
+
+ /* don't need these anymore */
+ if (num_local_points > 0)
+ {
+ interp_local_coords -= num_local_points * num_dims;
+ free (interp_local_coords);
+ }
+ free (local_out_arrays);
+ free (in_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
+ 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; type <= max_type; type++)
+ {
+ /* skip unused types */
+ if (type_desc[type].num_arrays <= 0)
+ {
+ continue;
+ }
+
+ /* set up the communication (this is type-independent) */
+ myGH->send_count[0] = type_desc[type].num_arrays * myGH->num_points_to[0];
+ myGH->recv_count[0] = type_desc[type].num_arrays * myGH->num_points_from[0];
+ myGH->send_displ[0] = myGH->recv_displ[0] = 0;
+ for (proc = 1; proc < pughGH->nprocs; proc++)
+ {
+ myGH->send_count[proc] = type_desc[type].num_arrays *
+ myGH->num_points_to[proc];
+ myGH->recv_count[proc] = type_desc[type].num_arrays *
+ myGH->num_points_from[proc];
+ myGH->send_displ[proc] = myGH->send_displ[proc-1] +
+ myGH->send_count[proc-1];
+ myGH->recv_displ[proc] = myGH->recv_displ[proc-1] +
+ myGH->recv_count[proc-1];
+ }
+
+ /* allocate buffer for receiving my own requested points */
+ /* avoid NULL pointers here because MPI_Alltoallv() doesn't like it */
+ if (num_points > 0)
+ {
+ type_desc[type].recv_buffer = malloc (num_points *
+ type_desc[type].num_arrays *
+ type_desc[type].vtypesize);
+ }
+ else
+ {
+ /* access to such a fake address should crash the code on most systems */
+ type_desc[type].recv_buffer = (void *) 1;
+ }
+
+ /* now exchange the data for this CCTK data type */
+ CACTUS_MPI_ERROR (MPI_Alltoallv (type_desc[type].send_buffer,
+ myGH->send_count, myGH->send_displ,
+ type_desc[type].mpi_type,
+ type_desc[type].recv_buffer,
+ myGH->recv_count, myGH->recv_displ,
+ type_desc[type].mpi_type,
+ pughGH->PUGH_COMM_WORLD));
+
+ /* now that the data is sent we don't need the buffer anymore */
+ if (num_local_points > 0)
+ {
+ free (type_desc[type].send_buffer);
+ }
+
+ /* no sort neccessary if there are no points */
+ if (num_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. */
+ offset = 0;
+ type_desc[type].buffer = (char *) type_desc[type].recv_buffer;
+ for (proc = 0; proc < nprocs; proc++)
+ {
+ for (array = 0; array < num_out_arrays; array++)
+ {
+ if (out_array_types[array] != type)
+ {
+ continue;
+ }
+
+ /* now do the sorting according to the CCTK data type */
+ if (out_array_types[array] == CCTK_VARIABLE_CHAR)
+ {
+ SORT_TYPED_ARRAY (CCTK_CHAR, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+ else if (out_array_types[array] == CCTK_VARIABLE_INT)
+ {
+ SORT_TYPED_ARRAY (CCTK_INT, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+ else if (out_array_types[array] == CCTK_VARIABLE_REAL)
+ {
+ SORT_TYPED_ARRAY (CCTK_REAL, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+ else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX)
+ {
+ SORT_TYPED_ARRAY (CCTK_COMPLEX, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+#ifdef CCTK_REAL4
+ else if (out_array_types[array] == CCTK_VARIABLE_REAL4)
+ {
+ SORT_TYPED_ARRAY (CCTK_REAL4, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+ else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX8)
+ {
+ SORT_TYPED_ARRAY (CCTK_COMPLEX8, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+#endif
+#ifdef CCTK_REAL8
+ else if (out_array_types[array] == CCTK_VARIABLE_REAL8)
+ {
+ SORT_TYPED_ARRAY (CCTK_REAL8, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+ else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX16)
+ {
+ SORT_TYPED_ARRAY (CCTK_COMPLEX16, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+#endif
+#ifdef CCTK_REAL16
+ else if (out_array_types[array] == CCTK_VARIABLE_REAL16)
+ {
+ SORT_TYPED_ARRAY (CCTK_REAL16, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+ else if (out_array_types[array] == CCTK_VARIABLE_COMPLEX32)
+ {
+ SORT_TYPED_ARRAY (CCTK_COMPLEX32, type_desc[type].buffer,
+ out_arrays[array], myGH->num_points_from[proc],
+ myGH->indices, offset);
+ }
+#endif
+ else
+ {
+ CCTK_WARN (0, "Implementation error");
+ }
+
+ } /* end of loop over all output arrays */
+
+ /* advance the offset into the communication receive buffer */
+ offset += myGH->num_points_from[proc];
+
+ } /* end of loop over all processors */
+
+ /* this communication receive buffer isn't needed anymore */
+ free (type_desc[type].recv_buffer);
+
+ } /* 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);
+}
+
+
+/*@@
+ @routine PUGHInterp_InterpLocal
+ @date Sun Jul 04 1999
+ @author Thomas Radke
+ @desc
+ The interpolation operator registered with the CCTK
+ under the name "uniform cartesian".
+
+ Interpolates a list of non-distributed, processor-local
+ input 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 num_points
+ @vdesc number of points to be interpolated on this processor
+ @vtype int
+ @vio in
+ @endvar
+ @var num_dims
+ @vdesc dimensionality of the underlying grid
+ @vtype int
+ @vio in
+ @endvar
+ @var num_in_arrays
+ @vdesc number of input arrays to interpolate from
+ @vtype int
+ @vio in
+ @endvar
+ @var num_out_arrays
+ @vdesc number of output arrays to interpolate to
+ @vtype int
+ @vio in
+ @endvar
+ @var coord_dims
+ @vdesc dimensions of the underlying grid
+ @vtype int [num_dims]
+ @vio in
+ @endvar
+ @var coord_arrays
+ @vdesc list of grid coordinate arrays
+ @vtype void *[num_dims]
+ @vio in
+ @endvar
+ @var coord_array_types
+ @vdesc CCTK data type of coordinate arrays
+ @vtype int [num_dims]
+ @vio int
+ @endvar
+ @var interp_coord_arrays
+ @vdesc coordinates of points to interpolate at
+ @vtype void *[num_dims]
+ @vio in
+ @endvar
+ @var interp_coord_array_types
+ @vdesc CCTK data type of coordinate arrays to interpolate at
+ @vtype int [num_dims]
+ @vio out
+ @endvar
+ @var in_arrays
+ @vdesc list of input arrays to interpolate on
+ @vtype void *[num_in_arrays]
+ @vio in
+ @endvar
+ @var in_array_types
+ @vdesc CCTK data types of input arrays
+ @vtype int [num_in_arrays]
+ @vio in
+ @endvar
+ @var out_arrays
+ @vdesc list of output arrays to interpolate to
+ @vtype void *[num_out_arrays]
+ @vio out
+ @endvar
+ @var out_array_types
+ @vdesc CCTK data types of output arrays
+ @vtype int [num_out_arrays]
+ @vio in
+ @endvar
+
+ @returntype int
+ @returndesc
+ 0 - successful interpolation
+ -1 - in case of any errors
+ @endreturndesc
+@@*/
+int PUGHInterp_InterpLocal (cGH *GH,
+ int order,
+ int num_points,
+ int num_dims,
+ int num_in_arrays,
+ int num_out_arrays,
+ int coord_dims[],
+ void *coord_arrays[],
+ int coord_array_types[],
+ void *interp_coord_arrays[],
+ int interp_coord_array_types[],
+ void *in_arrays[],
+ int in_array_types[],
+ void *out_arrays[],
+ int out_array_types[])
+{
+ int dim, point, retval;
+ CCTK_REAL *coords, *origin, *delta, **data;
+
+
+ /* check arguments */
+ retval = PUGHInterp_CheckArguments (GH, num_dims, num_points,
+ num_in_arrays, num_out_arrays,
+ interp_coord_array_types);
+ if (retval <= 0)
+ {
+ return (retval);
+ }
+ for (dim = 0; dim < num_dims; dim++)
+ {
+ if (coord_array_types[dim] != CCTK_VARIABLE_REAL)
+ {
+ CCTK_WARN (1, "Coordinates should be of type CCTK_REAL");
+ return (-1);
+ }
+ if (interp_coord_array_types[dim] != CCTK_VARIABLE_REAL)
+ {
+ CCTK_WARN (1, "Interpolation coordinates should be of type CCTK_REAL");
+ return (-1);
+ }
+ }
+
+ /* get the grid spacings - this assumes a cartesian grid */
+ origin = (CCTK_REAL *) malloc (2 * num_dims * sizeof (CCTK_REAL));
+ delta = origin + num_dims;
+ data = (CCTK_REAL **) coord_arrays;
+ for (dim = 0; dim < num_dims; dim++)
+ {
+ origin[dim] = data[dim][0];
+ delta[dim] = data[dim][1] - data[dim][0];
+ }
+
+ /* sort the individual interpolation coordinate arrays into a single one */
+ coords = (CCTK_REAL *) malloc (num_dims * num_points * sizeof (CCTK_REAL));
+ data = (CCTK_REAL **) interp_coord_arrays;
+ for (point = 0; point < num_points; point++)
+ {
+ for (dim = 0; dim < num_dims; dim++)
+ {
+ *coords++ = data[dim][point];
+ }
+ }
+ coords -= num_dims * num_points;
+
+ /* call the interpolator function */
+ retval = PUGHInterp_Interpolate (order,
+ num_points, num_dims, num_out_arrays,
+ coord_dims, coords, origin, delta,
+ in_array_types, in_arrays,
+ out_array_types, out_arrays);
+
+ /* free allocated resources */
+ free (coords);
+ free (origin);
+
+ return (retval);
+}
+
+
+/**************************************************************************/
+/* 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[num_dims][num_local_points].
+ <B>
+ This means for the single-processor case to sort
+
+ inCoords1[num_points], inCoords2[npoints], ..., inCoords<num_dims>[npoints]
+ into coords[num_dims][num_points]
+
+ where num_points == num_local_points.
+ <B>
+ In the multiprocessor case all processors map their points' coordinates
+ to the processor that owns this point and exchange this information
+ via MPI_Alltoall ().
+ num_local_points is then the number of all processors' points to be
+ interpolated locally on this processor.
+ <B>
+ 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 cGH
+ @vio in
+ @endvar
+ @var num_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 num_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 num_dims
+ @vio in
+ @endvar
+ @var num_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[num_dims][num_local_points]
+ @vio out
+ @endvar
+
+@@*/
+static int GetLocalCoords (cGH *GH,
+ int num_points,
+ const char *coord_system,
+ pGExtras *extras,
+ CCTK_REAL *coords[],
+ int *num_local_points,
+ CCTK_REAL **local_coords)
+{
+ DECLARE_CCTK_PARAMETERS
+ 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.00001
+
+
+ /* Do a fliperooni in octant mode */
+ /* This flips the sign around the origin only */
+ /* if the left side (cx0) is bigger than the position to interpolate at */
+ /* BoxInBox needs this for interpolating onto grids */
+ if (CCTK_Equals (domain, "octant"))
+ {
+ for (point = 0; point < num_points; point++)
+ {
+ tmp = 0;
+ for (dim = 0; dim < extras->dim; dim++)
+ {
+ if (coords[dim][point] <= GH->cctk_origin_space[dim])
+ {
+ tmp = 1; /* mark flipping */
+ break;
+ }
+ }
+/*** FIXME *** really change input arrays here ??? */
+ if (tmp)
+ {
+ for (dim = 0; dim < extras->dim; dim++)
+ {
+ coords[dim][point] = fabs (coords[dim][point]);
+ }
+ }
+ }
+ }
+
+ /* get GH extension handles for PUGHInterp and PUGH */
+ myGH = (pughInterpGH *) CCTK_GHExtension (GH, "PUGHInterp");
+ pughGH = (pGH *) CCTK_GHExtension (GH, "PUGH");
+
+ /* This holds the proccessor for *each* of num_points points */
+ if (num_points > 0)
+ {
+ myGH->whichproc = (int *) malloc (2 * num_points * sizeof (int));
+ }
+ else
+ {
+ myGH->whichproc = NULL;
+ }
+ /* indices[] is used to make the sorting easier
+ when receiving the output data */
+ myGH->indices = myGH->whichproc + num_points;
+
+ /* initialize whichproc with invalid processor number -1 */
+ for (point = 0; point < num_points; point++)
+ {
+ myGH->whichproc[point] = -1;
+ }
+
+ /* initialize num_points_from to 0 for counting it up in the following loop */
+ nprocs = CCTK_nProcs (GH);
+ memset (myGH->num_points_from, 0, nprocs * sizeof (CCTK_INT));
+
+ /* allocate the ranges for my local coordinates */
+ range_min = (CCTK_REAL *) 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);
+ delta[dim] = (delta[dim] - origin[dim]) / extras->nsize[dim];
+ }
+
+ /* locate the points to interpolate at */
+ for (proc = 0; proc < nprocs; proc++)
+ {
+ for (dim = 0; dim < extras->dim; dim++)
+ {
+ /* compute the coordinate ranges */
+ range_min[dim] = origin[dim] + (extras->lb[proc][dim] - FUDGE)*delta[dim];
+ range_max[dim] = origin[dim] + (extras->ub[proc][dim] + FUDGE)*delta[dim];
+ }
+
+ /* and now which point will be processed by what processor */
+ for (point = 0; point < num_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->num_points_from[proc]++;
+ }
+ }
+ }
+ /* don't need this anymore */
+ free (range_min);
+
+ /* make sure that all points could be mapped onto a processor */
+ tmp = 0;
+ for (point = 0; point < num_points; point++)
+ {
+ if (myGH->whichproc[point] < 0)
+ {
+ int i;
+ char *msg = (char *) 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 num_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->num_points_from, 1, PUGH_MPI_INT,
+ myGH->num_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->num_points_from[proc],
+ myGH->num_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 num_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, num_points_to, not how many we get back from everyone else,
+ eg, num_points_from. The number we are sending, of course, is
+ all of our locations, eg, num_points */
+ *num_local_points = 0;
+ for (proc = 0; proc < nprocs; proc++)
+ {
+ *num_local_points += myGH->num_points_to[proc];
+ }
+
+#ifdef PUGHINTERP_DEBUG
+ printf ("processor %d gets %d points in total\n",
+ CCTK_MyProc (GH), *num_local_points);
+#endif
+
+ /* allocate the local coordinates array (sorted in processor order)
+ and the resulting coordinates array that I have to process */
+ proc_coords = (CCTK_REAL *) malloc (extras->dim * num_points *
+ sizeof (CCTK_REAL));
+ *local_coords = (CCTK_REAL *) malloc (extras->dim * *num_local_points *
+ sizeof (CCTK_REAL));
+
+ /* now sort my own coordinates as tupels of [extras->dim] */
+ tmp = 0;
+ for (proc = 0; proc < nprocs; proc++)
+ {
+ for (point = 0; point < num_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->send_count[0] = extras->dim * myGH->num_points_from[0];
+ myGH->recv_count[0] = extras->dim * myGH->num_points_to[0];
+ myGH->send_displ[0] = myGH->recv_displ[0] = 0;
+ for (proc = 1; proc < nprocs; proc++)
+ {
+ myGH->send_count[proc] = extras->dim * myGH->num_points_from[proc];
+ myGH->recv_count[proc] = extras->dim * myGH->num_points_to[proc];
+ myGH->send_displ[proc] = myGH->send_displ[proc-1] +
+ myGH->send_count[proc-1];
+ myGH->recv_displ[proc] = myGH->recv_displ[proc-1] +
+ myGH->recv_count[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->send_count,
+ myGH->send_displ, PUGH_MPI_REAL,
+ *local_coords, myGH->recv_count,
+ myGH->recv_displ, PUGH_MPI_REAL,
+ pughGH->PUGH_COMM_WORLD));
+
+ /* don't need this anymore */
+ free (proc_coords);
+
+ return (0);
+}
+#endif /* CCTK_MPI */
+
+
+ /*@@
+ @routine PUGHInterp_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 cGH *
+ @vio in
+ @endvar
+ @var num_dims
+ @vdesc dimensionality of the underlying grid
+ @vtype int
+ @vio in
+ @endvar
+ @var num_points
+ @vdesc number of points to interpolate at
+ @vtype int
+ @vio in
+ @endvar
+ @var num_in_arrays
+ @vdesc number of passed input arrays
+ @vtype int
+ @vio in
+ @endvar
+ @var num_out_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 [num_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 PUGHInterp_CheckArguments (cGH *GH,
+ int num_dims,
+ int num_points,
+ int num_in_arrays,
+ int num_out_arrays,
+ int interp_coord_array_types[])
+{
+ int i;
+
+
+ /* check for invalid arguments */
+ if (num_dims < 0 || num_points < 0 || num_in_arrays < 0 || num_out_arrays < 0)
+ {
+ return (-1);
+ }
+
+ /* check if there's anything to do at all */
+ /* NOTE: num_points can be 0 in a collective call */
+ if (num_dims == 0 || (CCTK_nProcs (GH) == 1 && num_points == 0) ||
+ num_in_arrays == 0 || num_out_arrays == 0)
+ {
+ return (0);
+ }
+
+ /* for now we can only deal with coordinates of type CCTK_REAL */
+ for (i = 0; i < num_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 (num_in_arrays != num_out_arrays)
+ {
+ CCTK_WARN (1, "Number of input arrays must match number of output arrays");
+ return (-1);
+ }
+
+ return (1);
+}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-07 06:58:40 +0000