aboutsummaryrefslogtreecommitdiff
path: root/src/template.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/template.c')
-rw-r--r--src/template.c192
1 files changed, 130 insertions, 62 deletions
diff --git a/src/template.c b/src/template.c
index fe16818..6cce81c 100644
--- a/src/template.c
+++ b/src/template.c
@@ -209,6 +209,9 @@
@version $Header$
@@*/
+#undef AEILOCALINTERP_DEBUG /* define this for verbose debugging output */
+#undef AEILOCALINTERP_DEBUG2 /* define this for even more verbose debugging output */
+
/******************************************************************************/
/*
@@ -279,11 +282,16 @@
If you change the arguments for this function, note that you
must also change the prototype in "template.h".
- @var error_flags
- @vdesc If we encounter an out-of-range point and this pointer
- is non-NULL, we store a description of the out-of-range
- point in the pointed-to structure.
- @vtype struct error_flags* error_flags;
+ @var error_info
+ @vdesc If error_info->per_point_status is non-NULL,
+ then we store per-point interpolator status for
+ each point in error_info->per_point_status[pt] .
+ We set error_info->error_count to the number of
+ points in error.
+ We set the remaining elements of error_info
+ to a detailed description of the first point in error
+ (if any).
+ @vtype struct error_info* error_info;
@vio out
@var molecule_structure_flags
@@ -349,7 +357,7 @@ int FUNCTION_NAME(/***** coordinate system *****/
const CCTK_INT operand_indices[],
const CCTK_INT operation_codes[],
/***** other return results *****/
- struct error_flags* error_flags,
+ struct error_info* error_info,
struct molecule_structure_flags* molecule_structure_flags,
struct molecule_min_max_m_info* molecule_min_max_m_info,
CCTK_INT* const molecule_positions[],
@@ -392,6 +400,9 @@ int FUNCTION_NAME(/***** coordinate system *****/
* }
* compute position of interpolation molecules
* with respect to the grid into the XYZ variables
+ * and store per-point status if requested
+ * if (this point is outside the grid)
+ * then continue;
* if (querying molecule positions)
* then store this molecule position
* compute molecule center 1-d position in input arrays
@@ -744,6 +755,12 @@ int out;
* interpolate at each point
*/
int pt;
+int return_status = 0;
+error_info->error_count = 0;
+#ifdef AEILOCALINTERP_DEBUG
+printf("AEILocalInterp:: in interpolator fn: N_DIMS=%d N_interp_points=%d\n", (int) N_DIMS, N_interp_points);
+fflush(stdout);
+#endif
for (pt = 0 ; pt < N_interp_points ; ++pt)
{
/* this struct holds a molecule-sized piece of a single */
@@ -859,23 +876,28 @@ int pt;
/* end of switch (interp_coords_type_code) */
}
- /* end of for (axis = ...) loop */
+ /* end of for (axis = ...) loop to get interp coords */
}
}
-#ifdef PICKLE
- #if (N_DIMS == 1)
-printf("pt=%d interp coord %.16g\n",
- pt, (double) interp_coords_fp[0]);
- #endif
- #if (N_DIMS == 3)
-printf("pt=%d interp coords %.16g %.16g %.16g\n",
- pt, (double) interp_coords_fp[0],
- (double) interp_coords_fp[1],
- (double) interp_coords_fp[2]);
+#ifdef AEILOCALINTERP_DEBUG2
+ #if (N_DIMS == 1)
+ printf("pt=%d interp coord %.16g\n",
+ pt, (double) interp_coords_fp[0]);
+ #elif (N_DIMS == 2)
+ printf("pt=%d interp coords %.16g %.16g\n",
+ pt, (double) interp_coords_fp[0],
+ (double) interp_coords_fp[1],
+ #elif (N_DIMS == 3)
+ printf("pt=%d interp coords %.16g %.16g %.16g\n",
+ pt, (double) interp_coords_fp[0],
+ (double) interp_coords_fp[1],
+ (double) interp_coords_fp[2]);
+ #else
+ #error "N_DIMS may not be > 3!"
#endif
-#endif
-
+ fflush(stdout);
+#endif /* AEILOCALINTERP_DEBUG2 */
/*
* compute position of interpolation molecules with respect to
@@ -891,6 +913,7 @@ printf("pt=%d interp coords %.16g %.16g %.16g\n",
* are "clean" and thus more likely to be register-allocated
*/
{
+ int this_point_status = 0;
int center_ijk_temp[MAX_N_DIMS];
fp xyz_temp[MAX_N_DIMS];
int axis;
@@ -898,41 +921,75 @@ printf("pt=%d interp coords %.16g %.16g %.16g\n",
{
const int ibndry_min = 2*axis;
const int ibndry_max = 2*axis + 1;
- const int status = AEILocalInterp_molecule_posn
- (coord_origin[axis], coord_delta[axis],
- input_array_min_subscripts[axis]
- + N_boundary_points_to_omit[ibndry_min],
- input_array_max_subscripts[axis]
- - N_boundary_points_to_omit[ibndry_max],
- MOLECULE_SIZE,
- boundary_off_centering_tolerance[ibndry_min],
- boundary_off_centering_tolerance[ibndry_max],
- boundary_extrapolation_tolerance[ibndry_min],
- boundary_extrapolation_tolerance[ibndry_max],
- interp_coords_fp[axis],
- &center_ijk_temp[axis], &xyz_temp[axis]);
- if (status < 0)
- then {
- if (status == MOLECULE_POSN_ERROR_GRID_TINY)
- then return CCTK_ERROR_INTERP_GRID_TOO_TINY;
- /*** ERROR RETURN ***/
- if (error_flags != NULL)
- then {
- error_flags->error_pt = pt;
- error_flags->error_axis = axis;
- if (status == MOLECULE_POSN_ERROR_X_LT_MIN)
- then {
- error_flags->error_ibndry = ibndry_min;
- error_flags->error_direction = -1;
- }
- else {
- error_flags->error_ibndry = ibndry_max;
- error_flags->error_direction = +1;
- }
- }
- return CCTK_ERROR_INTERP_POINT_OUTSIDE; /*** ERROR RETURN ***/
+ const int mp_status = AEILocalInterp_molecule_posn
+ (coord_origin[axis], coord_delta[axis],
+ input_array_min_subscripts[axis]
+ + N_boundary_points_to_omit[ibndry_min],
+ input_array_max_subscripts[axis]
+ - N_boundary_points_to_omit[ibndry_max],
+ MOLECULE_SIZE,
+ boundary_off_centering_tolerance[ibndry_min],
+ boundary_off_centering_tolerance[ibndry_max],
+ boundary_extrapolation_tolerance[ibndry_min],
+ boundary_extrapolation_tolerance[ibndry_max],
+ interp_coords_fp[axis],
+ &center_ijk_temp[axis], &xyz_temp[axis]);
+
+ /*
+ * unfortunately, the status codes from AEILocalInterp_molecule_posn()
+ * are different from the ones we need ==> we have to decode them
+ * to get the status for this axis
+ */
+ int this_axis_status;
+ switch (mp_status)
+ {
+ case 0:
+ this_axis_status = 0;
+ break;
+ case MOLECULE_POSN_ERROR_GRID_TINY:
+ this_axis_status = CCTK_ERROR_INTERP_GRID_TOO_SMALL;
+ break;
+ case MOLECULE_POSN_ERROR_X_LT_MIN:
+ case MOLECULE_POSN_ERROR_X_GT_MAX:
+ this_axis_status = CCTK_ERROR_INTERP_POINT_OUTSIDE;
+ break;
+ default:
+ CCTK_VWarn(BUG_MSG_SEVERITY_LEVEL,
+ __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+" CCTK_InterpLocalUniform(): internal error!\n"
+" unexpected result %d from AEILocalInterp_molecule_posn()\n"
+" at pt=%d axis=%d!\n"
+ ,
+ mp_status, pt, axis); /*NOTREACHED*/
}
+
+ if (this_axis_status < this_point_status)
+ then this_point_status = this_axis_status;
+
+ /* end of for (axis = ...) loop to locate interp points in the grid */
}
+
+ if (error_info->per_point_status != NULL)
+ then error_info->per_point_status[pt] = this_point_status;
+
+ if (this_point_status < 0)
+ then {
+ ++error_info->error_count;
+
+ if (this_point_status < return_status)
+ then return_status = this_point_status;
+
+#ifdef AEILOCALINTERP_DEBUG
+printf("AEILocalInterp:: pt=%d has this_point_status=%d ==> skipping it\n", pt, this_point_status);
+fflush(stdout);
+#endif
+
+ /* this point is in error (eg outside grid) */
+ /* ==> don't try to interpolate at this point! */
+ continue;
+ }
+
{
#if (N_DIMS >= 1)
const int center_i = center_ijk_temp[X_AXIS];
@@ -947,16 +1004,20 @@ printf("pt=%d interp coords %.16g %.16g %.16g\n",
const fp z = xyz_temp [Z_AXIS];
#endif
-#ifdef PICKLE
- #if (N_DIMS == 1)
-printf("interp center_i = %d\n", center_i);
-printf("interp grid-relative x = %.16g\n", (double) x);
- #endif
- #if (N_DIMS == 3)
-printf("interp center_ijk = %d %d %d\n", center_i, center_j, center_k);
-printf("interp grid-relative xyz = %.16g %.16g %.16g\n",
- (double) x, (double) y, (double) z);
+#ifdef AEILOCALINTERP_DEBUG
+ #if (N_DIMS == 1)
+ printf("interp center_i = %d\n", center_i);
+ printf("interp grid-relative x = %.16g\n", (double) x);
+ #elif (N_DIMS == 2)
+ printf("interp center_ij = %d %d\n", center_i, center_j);
+ printf("interp grid-relative xy = %.16g %.16g\n", (double) x, (double) y);
+ #elif (N_DIMS == 3)
+ printf("interp center_ijk = %d %d %d\n", center_i, center_j, center_k);
+ printf("interp grid-relative xyz = %.16g %.16g %.16g\n",
+ (double) x, (double) y, (double) z);
+ #error "N_DIMS may not be > 3!"
#endif
+ fflush(stdout);
#endif
/*
@@ -1352,6 +1413,13 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
/*
* store result in output array
*/
+ #ifdef AEILOCALINTERP_DEBUG
+ if ((pt & (pt-1)) == 0) /* pt is 0 or a power of 2 */
+ then {
+ printf("AEILocalInterp:: pt=%d out=%d --> storing result %g\n", pt, out, (double) result);
+ fflush(stdout);
+ }
+ #endif
switch (output_array_type_codes[out])
{
@@ -1563,7 +1631,7 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
/* end of for (pt = ...) loop */
}
-return 0; /*** NORMAL RETURN ***/
+return return_status;
}
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 18:03:31 +0000