From 8bf581ab72264463f2da5136e5385011d431ffe2 Mon Sep 17 00:00:00 2001 From: jthorn Date: Sun, 13 Aug 2006 19:54:28 +0000 Subject: fix a nasty bug where, if the interpolation coordinate was close to a grid boundary, the actual interpolation was sometimes done as if the interpolation coordinate was displaced by one grid spacing :( :( also add more debugging code (introduced to help track down this bug, but there's no harm in leaving it in for possible future use) also add explicit checks for a couple of error cases: * interpolation coordinate is not a finite number (eg NaN or infinity) * grid origin/delta not a finite number (eg NaN or infinity) * grid spacing is 0 (we need to divide by this!) also add explicit note about the test driver in this directory git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/AEILocalInterp/trunk@35 0f49ee68-0e4f-0410-9b9c-b2c123ded7ef --- src/InterpLocalUniform.h | 5 +++- src/molecule_posn.c | 65 ++++++++++++++++++++++++++++++++++-------------- src/template.c | 2 ++ 3 files changed, 53 insertions(+), 19 deletions(-) diff --git a/src/InterpLocalUniform.h b/src/InterpLocalUniform.h index a111f02..bb412d0 100644 --- a/src/InterpLocalUniform.h +++ b/src/InterpLocalUniform.h @@ -35,7 +35,7 @@ typedef int bool; /* round floating-point value to nearest integer */ /* ... result is expressed as floating point! */ /* ... needs for floor() */ -#define JT_ROUND(x) floor((x) + 0.5) +#define ROUND_TO_INTEGER__RESULT_IS_FP(x) floor((x) + 0.5) /******************************************************************************/ @@ -175,6 +175,9 @@ struct Jacobian_info /* grid is smaller than molecule */ #define MOLECULE_POSN_ERROR_GRID_TINY (-3) +/* someone passed us a NaN or other non-finite floating-point number */ +#define MOLECULE_POSN_ERROR_NAN (-4) + /******************************************************************************/ /* diff --git a/src/molecule_posn.c b/src/molecule_posn.c index 663d917..89debab 100644 --- a/src/molecule_posn.c +++ b/src/molecule_posn.c @@ -78,10 +78,10 @@ static const char *rcsid = "$Header$"; @hdate 27.Jan.2003 @hauthor Jonathan Thornburg @hdesc Complete rewrite: now supports - @var boundary_off_centering_tolerance_min, - @var boundary_off_centering_tolerance_max, - @var boundary_extrapolation_tolerance_min, - @var boundary_extrapolation_tolerance_max, + @var boundary_off_centering_tolerance_min, + @var boundary_off_centering_tolerance_max, + @var boundary_extrapolation_tolerance_min, + @var boundary_extrapolation_tolerance_max, also change to returning status code, also drop returning @var min_m and @var max_m since they were never used. @@ -185,6 +185,17 @@ int AEILocalInterp_molecule_posn(fp grid_origin, fp grid_delta, int debug, int* i_center, fp* x_rel) { +/* + * ***** IMPORTANT ***** + * + * This code is ++tricky. There's a fairly rigorous test driver for + * it in the file test_molecule_posn.c in this directory. This is a + * standalone test driver, which can be compiled via + * gmake -f Makefile.standalone + * in this directory. If you make any changes to this code, please + * verify that all the tests still pass! + */ + if (debug >= 8) then { printf("AEILocalInterp_molecule_posn():\n"); @@ -201,17 +212,38 @@ if (debug >= 8) printf(" x=%g\n", (double) x); } +/* + * basic sanity checks + */ + +/* is the molecule larger than the grid? */ +if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max)) + then return MOLECULE_POSN_ERROR_GRID_TINY; /*** ERROR RETURN ***/ + +/* has someone passed us NaNs for coordinates etc? */ +if (!finite(grid_origin) || !finite(grid_delta) || !finite(x)) + then return MOLECULE_POSN_ERROR_NAN; + +/* ditto if the grid spacing is zero (we'll need to divide by it)! */ +if (grid_delta == 0.0) + then return MOLECULE_POSN_ERROR_NAN; + + /* molecule radia (inherently positive) in +/- directions */ const int mr_plus = (molecule_size >> 1); const int mr_minus = molecule_size - mr_plus - 1; -/* range of x_rel for which this molecule is centered, cf. diagram in header comment */ +/* range of x_rel for which this molecule is centered, */ +/* cf. diagram in header comment */ const fp centered_min_x_rel = IS_EVEN(molecule_size) ? 0.0 : -0.5; const fp centered_max_x_rel = IS_EVEN(molecule_size) ? 1.0 : +0.5; -/* range of i where we *could* center the molecule, as floating-point numbers */ -const fp fp_centered_min_possible_i = grid_i_min + mr_minus + centered_min_x_rel; -const fp fp_centered_max_possible_i = grid_i_max - mr_plus + centered_max_x_rel; +/* range of i where a centered molecule would fit within the grid, */ +/* as floating-point numbers */ +const fp fp_centered_min_possible_i = grid_i_min + mr_minus + + centered_min_x_rel; +const fp fp_centered_max_possible_i = grid_i_max - mr_plus + + centered_max_x_rel; /* integer coordinate i of interpolation point, as a floating-point number */ const fp fp_i = (x - grid_origin) / grid_delta; @@ -227,10 +259,6 @@ if (debug > 8) printf(" fp_i=%g\n", (double) fp_i); } -/* is the molecule larger than the grid? */ -if (molecule_size > HOW_MANY_IN_RANGE(grid_i_min,grid_i_max)) - then return MOLECULE_POSN_ERROR_GRID_TINY; /*** ERROR RETURN ***/ - /* is interpolation point x beyond the extrapolation tolerance? */ if (fp_i < (fp)grid_i_min - boundary_extrapolation_tolerance_min) then return MOLECULE_POSN_ERROR_X_LT_MIN; /*** ERROR RETURN ***/ @@ -246,22 +274,23 @@ if (fp_i > fp_centered_max_possible_i + boundary_off_centering_tolerance_max) /* now choose the actual molecule position/centering: */ /* first set up a centered molecule */ { -fp fp_i_center = IS_EVEN(molecule_size) /* ... as a floating-point number */ +/* compute as a floating point number */ +/* (because that's what the C math library provides), */ +/* but this value is guaranteed to be integral */ +fp fp_i_center = IS_EVEN(molecule_size) ? floor(fp_i) - : JT_ROUND(fp_i); -int int_i_center = (int) fp_i_center; /* ... as an integer */ + : ROUND_TO_INTEGER__RESULT_IS_FP(fp_i); +int int_i_center = (int) fp_i_center; /* convert to integer */ if (debug > 8) then printf(" initial: fp_i_center=%g int_i_center=%d\n", (double) fp_i_center, int_i_center); -/* then clamp the molecule at the grid boundaries */ -if (int_i_center < grid_i_min) int_i_center = grid_i_min; +/* clamp the molecule at the grid boundaries */ if (int_i_center < grid_i_min + mr_minus) then { int_i_center = grid_i_min + mr_minus; fp_i_center = (fp) int_i_center; } -if (int_i_center > grid_i_max) int_i_center = grid_i_max; if (int_i_center > grid_i_max - mr_plus) then { int_i_center = grid_i_max - mr_plus; diff --git a/src/template.c b/src/template.c index c1ad45f..09ed7b3 100644 --- a/src/template.c +++ b/src/template.c @@ -974,6 +974,8 @@ if (debug > 0) case MOLECULE_POSN_ERROR_X_GT_MAX: this_axis_status = CCTK_ERROR_INTERP_POINT_OUTSIDE; break; + case MOLECULE_POSN_ERROR_NAN: + this_axis_status = CCTK_ERROR_INTERP_COORD_NAN; default: CCTK_VWarn(BUG_MSG_SEVERITY_LEVEL, __LINE__, __FILE__, CCTK_THORNSTRING, -- cgit v1.2.3