aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/InterpGridArrays.c213
1 files changed, 121 insertions, 92 deletions
diff --git a/src/InterpGridArrays.c b/src/InterpGridArrays.c
index d9e9abe..3fc2ab0 100644
--- a/src/InterpGridArrays.c
+++ b/src/InterpGridArrays.c
@@ -185,7 +185,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
void *const output_arrays[])
{
int i, retval;
- CCTK_REAL *origin_local, *delta_local;
+ CCTK_REAL *origin_local, *delta;
CCTK_INT *input_array_dims, *input_array_types, *input_array_time_levels;
const void **input_arrays;
const char *coord_system_name;
@@ -201,8 +201,9 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
char **output_arrays_local;
type_desc_t *this, *type_desc;
CCTK_REAL *range_min, *range_max;
- CCTK_REAL *origin_global, *delta_global;
- CCTK_REAL *interp_coords_proc, *coords;
+ CCTK_REAL *origin_global;
+ CCTK_REAL *interp_coords_proc, *coords, *bbox_local;
+ const CCTK_REAL **bbox_interp_coords;
#endif
@@ -265,7 +266,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
/* allocate some temporary arrays */
origin_local = malloc (2 * N_dims * sizeof (CCTK_REAL));
- delta_local = origin_local + N_dims;
+ delta = origin_local + N_dims;
input_arrays = malloc (N_input_arrays * sizeof (void *));
input_array_dims = malloc ((N_dims + 2*N_input_arrays) * sizeof (CCTK_INT));
input_array_types = input_array_dims + N_dims;
@@ -298,9 +299,9 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
copy the integer dimension array into an CCTK_INT array */
for (i = 0; i < N_dims; i++)
{
- CCTK_CoordLocalRange (GH, &origin_local[i], &delta_local[i], i + 1,
- NULL, coord_system_name);
- delta_local[i] = (delta_local[i] - origin_local[i]) / extras->lnsize[i];
+ CCTK_CoordLocalRange (GH, &origin_local[i], &delta[i], i + 1, NULL,
+ coord_system_name);
+ delta[i] = (delta[i] - origin_local[i]) / extras->lnsize[i];
input_array_dims[i] = extras->lnsize[i];
}
@@ -354,32 +355,12 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
input_array_types[i] = CCTK_VarTypeI (input_array_indices[i]);
}
-#if 0
- /* Check whether the number of ghostzones is sufficient for the local
- interpolator to handle points on the processor's local bounding box.
- This is done as a query, simply by setting N_output_arrays = 0. */
- if (retval == 0)
- {
- for (i = 0; i < N_dims; i++)
- {
- bbox_interp_coords[i] = origin_local;
- retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
- param_table_handle,
- origin_local, delta_local, 2,
- interp_coords_type, bbox_interp_coords,
- N_input_arrays, input_array_dims,
- input_array_types, input_arrays,
- 0, output_array_types,
- output_arrays);
- }
-#endif
-
/* single-processor case directly calls the local interpolator */
if (retval == 0 && pughGH->nprocs == 1)
{
retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
param_table_handle,
- origin_local, delta_local, N_points,
+ origin_local, delta, N_points,
interp_coords_type, interp_coords,
N_input_arrays, input_array_dims,
input_array_types, input_arrays,
@@ -501,22 +482,21 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
{
myGH->whichproc[point] = -1;
}
-
+
/* initialize N_points_from to 0 for counting it up in the following loop */
memset (myGH->N_points_from, 0, pughGH->nprocs * sizeof (CCTK_INT));
/* allocate the ranges for my local coordinates */
- range_min = malloc (4 * N_dims * sizeof (CCTK_REAL));
- range_max = range_min + 1*N_dims;
+ range_min = malloc (5 * N_dims * sizeof (CCTK_REAL));
+ range_max = range_min + 1*N_dims;
origin_global = range_min + 2*N_dims;
- delta_global = range_min + 3*N_dims;
+ bbox_local = range_min + 3*N_dims;
- /* get the global origin and delta of the coordinate system */
+ /* get the global origin of the coordinate system */
for (i = 0; i < N_dims; i++)
{
- CCTK_CoordRange (GH, &origin_global[i], &delta_global[i], i+1, NULL,
+ CCTK_CoordRange (GH, &origin_global[i], &range_max[i], i+1, NULL,
coord_system_name);
- delta_global[i] = (delta_global[i]-origin_global[i]) / (extras->nsize[i]-1);
}
/* locate the points to interpolate at */
@@ -526,15 +506,25 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
{
/* compute the coordinate ranges */
/* TODO: use bbox instead -- but the bboxes of other processors
- are not known */
+ are not known */
int const has_lower = extras->lb[proc][i] == 0;
int const has_upper = extras->ub[proc][i] == GH->cctk_gsh[i]-1;
- range_min[i] = origin_global[i] + (extras->lb[proc][i] + (!has_lower) * (extras->nghostzones[i]-0.5) - FUDGE)*delta_global[i];
- range_max[i] = origin_global[i] + (extras->ub[proc][i] - (!has_upper) * (extras->nghostzones[i]-0.5) + FUDGE)*delta_global[i];
+ range_min[i] = origin_global[i] + (extras->lb[proc][i] - FUDGE +
+ (!has_lower) * (extras->nghostzones[i]-0.5)) * delta[i];
+ range_max[i] = origin_global[i] + (extras->ub[proc][i] + FUDGE -
+ (!has_upper) * (extras->nghostzones[i]-0.5)) * delta[i];
+
#ifdef PUGHINTERP_DEBUG
printf("processor %d: range_min[%d]=%g range_max[%d]=%g\n",
proc, i, (double) range_min[i], i, (double) range_max[i]);
#endif
+
+ /* save this processor's bbox coordinates */
+ if (proc == pughGH->myproc)
+ {
+ bbox_local[2*i + 0] = range_min[i];
+ bbox_local[2*i + 1] = range_max[i];
+ }
}
/* and now which point will be processed by what processor */
@@ -563,29 +553,62 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
}
}
}
- /* don't need this anymore */
- free (range_min);
+
+ /* Check whether the number of ghostzones is sufficient for the local
+ interpolator to handle points on the processor's local bounding box.
+ This is done as a query, simply by passing all arguments as for the
+ real call to the local interpolator, except for the output_arrays[]
+ argument which contains NULL pointers. */
+ bbox_interp_coords = malloc (N_dims * sizeof (void *));
+ for (i = 0; i < N_dims; i++)
+ {
+ bbox_interp_coords[i] = &bbox_local[2 * i];
+ }
+ output_arrays_local = N_output_arrays > 0 ?
+ calloc (N_output_arrays, sizeof (void *)) : NULL;
+
+ retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
+ param_table_handle,
+ origin_local, delta, 2,
+ CCTK_VARIABLE_REAL,
+ (const void *const *) bbox_interp_coords,
+ N_input_arrays, input_array_dims,
+ input_array_types, input_arrays,
+ N_output_arrays, output_array_types,
+ (void *const *) output_arrays_local);
/* make sure that all points could be mapped onto a processor */
- for (point = j = 0; point < N_points; point++)
+ if (retval == 0)
{
- if (myGH->whichproc[point] < 0)
+ for (point = 0; point < N_points; point++)
{
- msg = malloc (80 + N_dims*20);
- sprintf (msg, "Unable to locate point %d [%f",
- point, (double) ((const CCTK_REAL *) interp_coords[0])[point]);
- for (i = 1; i < N_dims; i++)
+ if (myGH->whichproc[point] < 0)
{
- sprintf (msg, "%s %f",
- msg, (double) ((const CCTK_REAL *) interp_coords[i])[point]);
+ msg = malloc (80 + N_dims*20);
+ sprintf (msg, "Unable to locate point %d [%f",
+ point, (double) ((const CCTK_REAL *) interp_coords[0])[point]);
+ for (i = 1; i < N_dims; i++)
+ {
+ sprintf (msg, "%s %f",
+ msg, (double) ((const CCTK_REAL *) interp_coords[i])[point]);
+ }
+ sprintf (msg, "%s]", msg);
+ CCTK_WARN (1, msg);
+ free (msg);
+ retval = UTIL_ERROR_BAD_INPUT;
}
- sprintf (msg, "%s]", msg);
- CCTK_WARN (1, msg);
- free (msg);
- j = 1; /* mark as error */
}
}
- if (j)
+
+ /* don't need this anymore */
+ if (output_arrays_local)
+ {
+ free (output_arrays_local);
+ }
+ free (bbox_interp_coords);
+ free (range_min);
+
+ if (retval)
{
if (myGH->whichproc)
{
@@ -597,7 +620,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
free (input_array_dims);
free (type_desc);
- return (UTIL_ERROR_BAD_INPUT);
+ return (retval);
}
/* Now we want to resolve the N_points_from[]. Currently this is
@@ -631,7 +654,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
This means I send interp_coords[i][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,
@@ -683,7 +706,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
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 *interp_coords_local[] */
coords = NULL;
@@ -738,7 +761,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
the results in the intermediate local output arrays */
retval = CCTK_InterpLocalUniform (N_dims, local_interp_handle,
param_table_handle,
- origin_local, delta_local, N_points_local,
+ origin_local, delta, N_points_local,
interp_coords_type,
(const void **) interp_coords_local,
N_input_arrays, input_array_dims,
@@ -844,7 +867,7 @@ int PUGHInterp_InterpGridArrays (const cGH *GH,
}
/* go back from processor-sorted data to input-ordered data.
- The creation of the indices array above makes this not so bad. */
+ The creation of the indices array above makes this not so bad. */
this->buf = this->recvbuf;
for (proc = offset = 0; proc < pughGH->nprocs; proc++)
{
@@ -967,7 +990,7 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
{
int i, retval;
CCTK_INT *N_boundary_points_to_omit;
- CCTK_REAL *off_centering_tolerance, *extrapolation_tolerance;
+ CCTK_REAL *boundary_off_centering_tolerance,*boundary_extrapolation_tolerance;
retval = 0;
@@ -980,13 +1003,13 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
/* allocate table option arrays and initialize them to their defaults */
N_boundary_points_to_omit = malloc (2 * N_dims * sizeof (CCTK_INT));
- off_centering_tolerance = malloc (2 * 2 * N_dims * sizeof (CCTK_REAL));
- extrapolation_tolerance = off_centering_tolerance + 2 * N_dims;
+ boundary_off_centering_tolerance = malloc (2*2 * N_dims * sizeof (CCTK_REAL));
+ boundary_extrapolation_tolerance = boundary_off_centering_tolerance +2*N_dims;
for (i = 0; i < 2 * N_dims; i++)
{
N_boundary_points_to_omit[i] = 0;
- off_centering_tolerance[i] = 0.0;
- extrapolation_tolerance[i] = 1e-10;
+ boundary_off_centering_tolerance[i] = 999.0;
+ boundary_extrapolation_tolerance[i] = 1e-10;
}
/* read the 'N_boundary_points_to_omit' options array from the table
@@ -996,54 +1019,57 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
"N_boundary_points_to_omit");
if (i < 0 && i != UTIL_ERROR_TABLE_NO_SUCH_KEY)
{
- CCTK_WARN (1, "Array with key 'N_boundary_points_to_omit' must be of "
- "CCTK_INT datatype");
+ CCTK_WARN (1, "Options array with key 'N_boundary_points_to_omit' "
+ "must be of CCTK_INT datatype");
retval = UTIL_ERROR_BAD_INPUT;
}
else if (i >= 0 && i < 2 * N_dims)
{
- CCTK_WARN (1, "Array with key 'N_boundary_points_to_omit' has fewer "
- "elements than interpolation bounding box");
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Options array with key 'N_boundary_points_to_omit' "
+ "must have exactly 2 * N_dims (= %d) elements", N_dims);
retval = UTIL_ERROR_BAD_INPUT;
}
- /* read the 'off_centering_tolerance' options array from the table
+ /* read the 'boundary_off_centering_tolerance' options array from the table
and verify its datatype and element count */
i = Util_TableGetRealArray (param_table_handle, 2 * N_dims,
- off_centering_tolerance,
- "off_centering_tolerance");
+ boundary_off_centering_tolerance,
+ "boundary_off_centering_tolerance");
if (i < 0 && i != UTIL_ERROR_TABLE_NO_SUCH_KEY)
{
- CCTK_WARN (1, "Array with key 'off_centering_tolerance' must be of "
- "CCTK_REAL datatype");
+ CCTK_WARN (1, "Options array with key 'boundary_off_centering_tolerance' "
+ "must be of CCTK_REAL datatype");
retval = UTIL_ERROR_BAD_INPUT;
}
else if (i >= 0 && i < 2 * N_dims)
{
- CCTK_WARN (1, "Array with key 'off_centering_tolerance' has fewer elements "
- "than interpolation bounding box");
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Options array with key 'boundary_off_centering_tolerance' "
+ "must have exactly 2 * N_dims (= %d) elements", N_dims);
retval = UTIL_ERROR_BAD_INPUT;
}
- /* read the 'extrapolation_tolerance' options array from the table
+ /* read the 'boundary_extrapolation_tolerance' options array from the table
and verify its datatype and element count */
i = Util_TableGetRealArray (param_table_handle, 2 * N_dims,
- extrapolation_tolerance,
- "extrapolation_tolerance");
+ boundary_extrapolation_tolerance,
+ "boundary_extrapolation_tolerance");
if (i < 0 && i != UTIL_ERROR_TABLE_NO_SUCH_KEY)
{
- CCTK_WARN (1, "Array with key 'extrapolation_tolerance' must be of "
- "CCTK_REAL datatype");
+ CCTK_WARN (1, "Options array with key 'boundary_extrapolation_tolerance' "
+ "must be of CCTK_REAL datatype");
retval = UTIL_ERROR_BAD_INPUT;
}
- else if (i >= 0 && i < 2 * N_dims)
+ else if (i >= 0 && i != 2 * N_dims)
{
- CCTK_WARN (1, "Array with key 'extrapolation_tolerance' has fewer elements "
- "than interpolation bounding box");
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Options array with key 'boundary_extrapolation_tolerance' "
+ "must have exactly 2 * N_dims (= %d) elements", N_dims);
retval = UTIL_ERROR_BAD_INPUT;
}
- /* read the 'input_array_time_levels' options array from the table
+ /* read the 'input_array_time_levels' options array from the table
and verify its datatype and element count */
i = Util_TableGetIntArray (param_table_handle, N_input_arrays,
input_array_time_levels,
@@ -1055,14 +1081,15 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
}
else if (i < 0)
{
- CCTK_WARN (1, "Array with key 'input_array_time_levels' must be of "
+ CCTK_WARN (1, "Options array with key 'input_array_time_levels' must be of "
"CCTK_INT datatype");
retval = UTIL_ERROR_BAD_INPUT;
}
- else if (i < N_input_arrays)
+ else if (i != N_input_arrays)
{
- CCTK_WARN (1, "Array with key 'input_array_time_levels' has fewer elements "
- "than total number of input arrays");
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Options array with key 'input_array_time_levels' must have "
+ "exactly N_input_arrays (= %d) elements", N_input_arrays);
retval = UTIL_ERROR_BAD_INPUT;
}
@@ -1072,8 +1099,8 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
if (! GH->cctk_bbox[i])
{
N_boundary_points_to_omit[i] = 0;
- off_centering_tolerance[i] = 0.0;
- extrapolation_tolerance[i] = 0.0;
+ boundary_off_centering_tolerance[i] = 0.0;
+ boundary_extrapolation_tolerance[i] = 0.0;
}
}
if (! retval)
@@ -1082,14 +1109,16 @@ static int PrepareParameterTable (const cGH *GH, int param_table_handle,
N_boundary_points_to_omit,
"N_boundary_points_to_omit");
Util_TableSetRealArray (param_table_handle, 2 * N_dims,
- off_centering_tolerance, "off_centering_tolerance");
+ boundary_off_centering_tolerance,
+ "boundary_off_centering_tolerance");
Util_TableSetRealArray (param_table_handle, 2 * N_dims,
- extrapolation_tolerance, "extrapolation_tolerance");
+ boundary_extrapolation_tolerance,
+ "boundary_extrapolation_tolerance");
}
/* free allocated resources */
free (N_boundary_points_to_omit);
- free (off_centering_tolerance);
+ free (boundary_off_centering_tolerance);
return (retval);
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 08:59:15 +0000