From 618cb394a3e9734b17ac2c428179ffa2e91aa28b Mon Sep 17 00:00:00 2001 From: jthorn Date: Thu, 24 Jul 2003 13:55:35 +0000 Subject: * add support for returning per-point interpolator status * drop support for returning status info on which side of the grid an outside-the-grid point is outside -- this turned out to be awkward to implement in combinatino with per-point status, and I believe it was only used by my apparent horizon finder to produce slightly more informative error messages in a few cases * convert some commented-out debugging code to #ifdef on various DEBUG symbols git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/AEILocalInterp/trunk@9 0f49ee68-0e4f-0410-9b9c-b2c123ded7ef --- src/InterpLocalUniform.c | 191 ++++++++++++++++++++++++++++++---------------- src/InterpLocalUniform.h | 15 ++-- src/template.c | 192 ++++++++++++++++++++++++++++++++--------------- src/template.h | 2 +- 4 files changed, 267 insertions(+), 133 deletions(-) diff --git a/src/InterpLocalUniform.c b/src/InterpLocalUniform.c index d4fafa8..76882d4 100644 --- a/src/InterpLocalUniform.c +++ b/src/InterpLocalUniform.c @@ -10,9 +10,10 @@ ** InterpLocalUniform - main driver routine ** ** check_boundary_tolerances - check boundary tolerances for consistency +** get_error_point_info - get per-point error reporting info from par table ** get_and_decode_molecule_family - get & decode molecule_family from par table ** get_molecule_positions - get molecule_positions from parameter table -** get_Jacobian_query - get Jacobian-query info from parameter table +** get_Jacobian_info - get Jacobian-query info from parameter table ** set_error_info - set error information in parameter table ** set_molecule_structure - set molecule structure info in parameter table ** set_molecule_min_max_m - set molecule size info in parameter table @@ -62,6 +63,9 @@ #include "Lagrange-maximum-degree/all_prototypes.h" #include "Hermite/all_prototypes.h" +#undef AEILOCALINTERP_DEBUG /* define this for verbose debugging output */ +#undef AEILOCALINTERP_DEBUG2 /* define this for even more verbose debugging output */ + /* the rcs ID and its dummy function to use it */ static const char* rcsid = "$Header$"; CCTK_FILEVERSION(AEITHorns_AEILocalInterp_src_InterpLocalUniform_c) @@ -124,6 +128,9 @@ static (int N_boundaries, const CCTK_REAL boundary_off_centering_tolerance[MAX_N_BOUNDARIES], const CCTK_REAL boundary_extrapolation_tolerance[MAX_N_BOUNDARIES]); +static + int get_error_point_info(int param_table_handle, + struct error_info *p_error_info); static int get_and_decode_molecule_family (int param_table_handle, @@ -139,7 +146,7 @@ static struct Jacobian_info* p_Jacobian_info); static int set_error_info(int param_table_handle, - struct error_flags* p_error_flags); + struct error_info* p_error_info); static int set_molecule_structure (int param_table_handle, @@ -151,12 +158,11 @@ static const struct molecule_min_max_m_info* p_molecule_min_max_m_info); static - int get_and_check_INT - (int param_table_handle, const char name[], - bool mandatory_flag, int default_value, - bool check_range_flag, int min_value, int max_value, - const char max_value_string[], - CCTK_INT* p_value); + int get_and_check_INT(int param_table_handle, const char name[], + bool mandatory_flag, int default_value, + bool check_range_flag, int min_value, int max_value, + const char max_value_string[], + CCTK_INT* p_value); static int get_INT_array(int param_table_handle, const char name[], bool default_flag, int default_value, @@ -1150,15 +1156,25 @@ check_boundary_tolerances(N_boundaries, boundary_off_centering_tolerance, boundary_extrapolation_tolerance); +/**************************************/ + + { +struct error_info error_info; +status = get_error_point_info(param_table_handle, + &error_info); +if (status != 0) + then return status; /*** ERROR RETURN ***/ + /**************************************/ { #define MOLECULE_FAMILY_BUFSIZ (MAX_MOLECULE_FAMILY_STRLEN+1) char molecule_family_string[MOLECULE_FAMILY_BUFSIZ]; enum molecule_family molecule_family; -status = get_and_decode_molecule_family(param_table_handle, - MOLECULE_FAMILY_BUFSIZ, molecule_family_string, - &molecule_family); +status = get_and_decode_molecule_family + (param_table_handle, + MOLECULE_FAMILY_BUFSIZ, molecule_family_string, + &molecule_family); if (status != 0) then return status; /*** ERROR RETURN ***/ @@ -1524,8 +1540,11 @@ if (p_interp_fn == NULL_INTERP_FN_PTR) } /* call the subfunction to actually do the interpolation */ +#ifdef AEILOCALINTERP_DEBUG +printf("AEILocalInterp::InterpLocalUniform.c: calling interpolator fn (N_interp_points=%d)\n", N_interp_points); +fflush(stdout); +#endif { -struct error_flags error_flags; struct molecule_structure_flags molecule_structure_flags; const int return_code = (*p_interp_fn)(coord_origin, coord_delta, @@ -1542,33 +1561,34 @@ const int return_code N_output_arrays, output_array_type_codes, output_arrays, operand_indices, operation_codes, - &error_flags, + &error_info, &molecule_structure_flags, p_molecule_min_max_m_info, p_molecule_positions, p_Jacobian_info); +#ifdef AEILOCALINTERP_DEBUG +printf("AEILocalInterp::InterpLocalUniform.c: back from interpolator fn with return_code=%d\n", return_code); +fflush(stdout); +#endif /******************************************************************************/ /* - * is an interpolation point out of range? + * set any further error status in the parameter table */ -if (return_code == CCTK_ERROR_INTERP_POINT_OUTSIDE) +status = set_error_info(param_table_handle, + &error_info); +if (status != 0) then { - status = set_error_info(param_table_handle, - &error_flags); - if (status != 0) + free(input_array_offsets); + free(operand_indices); + free(operation_codes); + if (p_Jacobian_info != NULL) then { - free(input_array_offsets); - free(operand_indices); - free(operation_codes); - if (p_Jacobian_info != NULL) - then { - free(Jacobian_info.Jacobian_offset); - free(Jacobian_info.Jacobian_pointer); - } - return status; /*** ERROR RETURN ***/ + free(Jacobian_info.Jacobian_offset); + free(Jacobian_info.Jacobian_pointer); } + return status; /*** ERROR RETURN ***/ } /******************************************************************************/ @@ -1625,6 +1645,10 @@ if (p_Jacobian_info != NULL) free(Jacobian_info.Jacobian_pointer); } +#ifdef AEILOCALINTERP_DEBUG +printf("AEILocalInterp::InterpLocalUniform.c: returning %d\n", return_code); +fflush(stdout); +#endif return return_code; /*** NORMAL RETURN ***/ } } @@ -1641,6 +1665,7 @@ return return_code; /*** NORMAL RETURN ***/ } } } + } } /******************************************************************************/ @@ -1686,6 +1711,56 @@ int ibndry; /******************************************************************************/ +/* + * This function tries to get + * CCTK_POINTER per_point_status + * from the parameter table, and sets found_per_point_status to + * indicate whether or not that key was found. + * + * If per_point_status is found, this function casts it to a CCTK_INT* + * and stores it in p_error_info->per_point_status. If not (if there's + * no such key in the table), this function sets + * p_error_info->per_point_status to a NULL pointer. + * + * Results: + * This function returns 0 for ok, or the (nonzero) return code for error. + */ +static + int get_error_point_info(int param_table_handle, + struct error_info *p_error_info) +{ +CCTK_POINTER per_point_status; +int status; + +status = Util_TableGetPointer(param_table_handle, + &per_point_status, + "per_point_status"); + +if (status == 1) + then { + p_error_info->found_per_point_status = true; + p_error_info->per_point_status = (CCTK_INT*) per_point_status; + } +else if (status == UTIL_ERROR_TABLE_NO_SUCH_KEY) + then { + p_error_info->found_per_point_status = false; + p_error_info->per_point_status = NULL; + } +else { + CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, +"\n" +" CCTK_InterpLocalUniform(): bad \"per_point_status\" table entry!\n" +" error status=%d", + status); + return status; /*** ERROR RETURN ***/ + } + +return 0; /*** NORMAL RETURN ***/ +} + +/******************************************************************************/ + /* * This function gets * const char molecule_family[] @@ -1931,46 +2006,33 @@ return 0; /*** NORMAL RETURN ***/ /******************************************************************************/ /* - * This function sets the error-info entries - * CCTK_INT error_pt - * CCTK_INT error_ibndry - * CCTK_INT error_axis - * CCTK_INT error_direction - * in the parameter table. + * If p_error_info->found_per_point_status is true, this function + * sets the paramater-table entry + * CCTK_INT error_point_status + * to the negative of p_error_info->error_count. * * Results: * This function returns 0 for ok, or the (nonzero) return code for error. */ static int set_error_info(int param_table_handle, - struct error_flags* p_error_flags) + struct error_info* p_error_info) { -const int status1 = Util_TableSetInt(param_table_handle, - p_error_flags->error_pt, - "error_pt"); -const int status2 = Util_TableSetInt(param_table_handle, - p_error_flags->error_ibndry, - "error_ibndry"); -const int status3 = Util_TableSetInt(param_table_handle, - p_error_flags->error_axis, - "error_axis"); -const int status4 = Util_TableSetInt(param_table_handle, - p_error_flags->error_direction, - "error_direction"); -if ((status1 < 0) || (status2 < 0) || (status3 < 0) || (status4 < 0)) +if (p_error_info->found_per_point_status) then { - CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, - __LINE__, __FILE__, CCTK_THORNSTRING, + const int status = Util_TableSetInt(param_table_handle, + -p_error_info->error_count, + "error_point_status"); + if (status < 0) + then CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL, + __LINE__, __FILE__, CCTK_THORNSTRING, "\n" " CCTK_InterpLocalUniform():\n" -" error setting error-info table entry/entries!\n" -" error status1=%d status2=%d status3=%d status4=%d" - , - status1, status2, status3, status4); - return (status1 < 0) ? status1 /*** ERROR RETURN ***/ - : (status2 < 0) ? status2 /*** ERROR RETURN ***/ - : (status3 < 0) ? status3 /*** ERROR RETURN ***/ - : status4; /*** ERROR RETURN ***/ +" error setting \"error_point_status\" in parameter table!" +" error status=%d" + , + status); + return status; /*** ERROR RETURN ***/ } return 0; /*** NORMAL RETURN ***/ @@ -2097,12 +2159,11 @@ return 0; /*** NORMAL RETURN ***/ * This function returns 0 for ok, or the (nonzero) return code for error. */ static - int get_and_check_INT - (int param_table_handle, const char name[], - bool mandatory_flag, int default_value, - bool check_range_flag, int min_value, int max_value, - const char max_value_string[], - CCTK_INT* p_value) + int get_and_check_INT(int param_table_handle, const char name[], + bool mandatory_flag, int default_value, + bool check_range_flag, int min_value, int max_value, + const char max_value_string[], + CCTK_INT* p_value) { CCTK_INT value; @@ -2281,8 +2342,8 @@ return 0; /*** NORMAL RETURN ***/ */ static int get_REAL_array(int param_table_handle, const char name[], - CCTK_REAL default_value, - int N, CCTK_REAL buffer[]) + CCTK_REAL default_value, + int N, CCTK_REAL buffer[]) { const int status = Util_TableGetRealArray(param_table_handle, diff --git a/src/InterpLocalUniform.h b/src/InterpLocalUniform.h index 28b807a..3515ff0 100644 --- a/src/InterpLocalUniform.h +++ b/src/InterpLocalUniform.h @@ -139,12 +139,17 @@ enum molecule_family /* number of real "parts" in a complex number */ #define COMPLEX_N_PARTS 2 -struct error_flags +struct error_info { - int error_pt; - int error_ibndry; - int error_axis; - int error_direction; + /* did we find error_point_status in the parameter table? */ + bool found_per_point_status; + + /* NULL pointer to skip per-point status, or */ + /* --> array of size N_interp_points to be set to per-point status */ + CCTK_INT* per_point_status; + + /* count of the number of points in error */ + CCTK_INT error_count; }; struct molecule_structure_flags 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], - ¢er_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], + ¢er_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; } } } diff --git a/src/template.h b/src/template.h index b5eb5b4..d1f9d92 100644 --- a/src/template.h +++ b/src/template.h @@ -33,7 +33,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[], -- cgit v1.2.3