diff options
Diffstat (limited to 'src/InterpGridArrays.c')
-rw-r--r-- | src/InterpGridArrays.c | 213 |
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); } |