From 57f2e0c8839578e57770cac1beca01880e8ae607 Mon Sep 17 00:00:00 2001 From: jthorn Date: Sun, 13 Aug 2006 18:15:09 +0000 Subject: * change control of debugging code from #ifdef to a new key "debug" in the parameter table (default is still no debugging) * add some additional debugging code git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/AEILocalInterp/trunk@34 0f49ee68-0e4f-0410-9b9c-b2c123ded7ef --- doc/documentation.tex | 12 +++ src/InterpLocalUniform.c | 39 +++++---- src/InterpLocalUniform.h | 1 + src/molecule_posn.c | 41 ++++++++++ src/template.c | 206 +++++++++++++++++++++++++---------------------- src/template.h | 4 +- src/test_molecule_posn.c | 4 + 7 files changed, 196 insertions(+), 111 deletions(-) diff --git a/doc/documentation.tex b/doc/documentation.tex index bd32575..e02a644 100644 --- a/doc/documentation.tex +++ b/doc/documentation.tex @@ -1380,6 +1380,18 @@ exactly the same cost as any other 6-point--molecule interpolation. The current implementation has all the framework for smoothing, but only the \verb|smoothing=0| case is implemented at the moment. +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\subsection{Debugging} + +Setting the optional parameter +\begin{verbatim} +const CCTK_INT debug; +\end{verbatim} +to a positive value turns on various debug printing inside the +interpolator. This is likely only for debugging the interpolator; +you need to look at the interpolator source code to interpret the output. + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Supressing Interpolation} diff --git a/src/InterpLocalUniform.c b/src/InterpLocalUniform.c index 9e7653e..96a2cb7 100644 --- a/src/InterpLocalUniform.c +++ b/src/InterpLocalUniform.c @@ -63,9 +63,6 @@ #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) @@ -1118,6 +1115,14 @@ if (N_dims > MAX_N_DIMS) * get the parameters from the parameter table, filling in defaults as needed */ +CCTK_INT debug; +status = get_and_check_INT(param_table_handle, "debug", + false, 0, /* optional parameter defaulting to 0 */ + false, 0, 0, NULL, + &debug); +if (status != 0) + then return status; /*** ERROR RETURN ***/ + { CCTK_INT order; status = get_and_check_INT(param_table_handle, "order", @@ -1540,10 +1545,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 +if (debug > 0) + then { + printf("AEILocalInterp::InterpLocalUniform.c: calling interpolator fn (N_interp_points=%d)\n", N_interp_points); + fflush(stdout); + } { struct molecule_structure_flags molecule_structure_flags; const int return_code @@ -1561,15 +1567,17 @@ const int return_code N_output_arrays, output_array_type_codes, output_arrays, operand_indices, operation_codes, + debug, &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 +if (debug > 0) + then { + printf("AEILocalInterp::InterpLocalUniform.c: back from interpolator fn with return_code=%d\n", return_code); + fflush(stdout); + } /******************************************************************************/ @@ -1645,10 +1653,11 @@ 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 +if (debug > 0) + then { + printf("AEILocalInterp::InterpLocalUniform.c: returning %d\n", return_code); + fflush(stdout); + } return return_code; /*** NORMAL RETURN ***/ } } diff --git a/src/InterpLocalUniform.h b/src/InterpLocalUniform.h index 30368ea..a111f02 100644 --- a/src/InterpLocalUniform.h +++ b/src/InterpLocalUniform.h @@ -251,6 +251,7 @@ int AEILocalInterp_molecule_posn(fp grid_origin, fp grid_delta, fp boundary_extrapolation_tolerance_min, fp boundary_extrapolation_tolerance_max, fp x, + int debug, int* i_center, fp* x_rel); /* functions in "util.c" */ diff --git a/src/molecule_posn.c b/src/molecule_posn.c index 603d3e1..663d917 100644 --- a/src/molecule_posn.c +++ b/src/molecule_posn.c @@ -137,6 +137,11 @@ static const char *rcsid = "$Header$"; @vtype fp x @endvar + @var debug + @vdesc A debugging flag (0 = no debug output, > 0 = print debug output) + @vtype int x + @endvar + @var i_center @vdesc A pointer to an value where this function should store the integer coordinate of the molecule center, @@ -177,8 +182,25 @@ int AEILocalInterp_molecule_posn(fp grid_origin, fp grid_delta, fp boundary_extrapolation_tolerance_min, fp boundary_extrapolation_tolerance_max, fp x, + int debug, int* i_center, fp* x_rel) { +if (debug >= 8) + then { + printf("AEILocalInterp_molecule_posn():\n"); + printf(" grid_origin=%g grid_delta=%g\n", + (double) grid_origin, (double) grid_delta); + printf(" grid_i_[min,max]=[%d,%d] molecule_size=%d\n", + grid_i_min, grid_i_max, molecule_size); + printf(" boundary_off_centering_tolerance_[min,max]=[%g,%g]\n", + (double) boundary_off_centering_tolerance_min, + (double) boundary_off_centering_tolerance_max); + printf(" boundary_extrapolation_tolerance_[min,max]=[%g,%g]\n", + (double) boundary_extrapolation_tolerance_min, + (double) boundary_extrapolation_tolerance_max); + printf(" x=%g\n", (double) x); + } + /* molecule radia (inherently positive) in +/- directions */ const int mr_plus = (molecule_size >> 1); const int mr_minus = molecule_size - mr_plus - 1; @@ -194,6 +216,17 @@ 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; +if (debug > 8) + then { + printf(" mr_{plus,minus}={%d,%d}\n", mr_plus, mr_minus); + printf(" centered_[min,max]_x_rel=[%g,%g]\n", + (double) centered_min_x_rel, (double) centered_max_x_rel); + printf(" fp_centered_[min,max]_possible_i=[%g,%g]\n", + (double) fp_centered_min_possible_i, + (double) fp_centered_max_possible_i); + 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 ***/ @@ -217,6 +250,9 @@ fp fp_i_center = IS_EVEN(molecule_size) /* ... as a floating-point number */ ? floor(fp_i) : JT_ROUND(fp_i); int int_i_center = (int) fp_i_center; /* ... as an 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; @@ -232,6 +268,11 @@ if (int_i_center > grid_i_max - mr_plus) fp_i_center = (fp) int_i_center; } +if (debug > 8) + then printf( + " final result after clamping: fp_i_center=%g int_i_center=%d\n", + (double) fp_i_center, int_i_center); + /* store the results */ if (i_center != NULL) then *i_center = int_i_center; diff --git a/src/template.c b/src/template.c index 48c1a5c..c1ad45f 100644 --- a/src/template.c +++ b/src/template.c @@ -209,9 +209,6 @@ @version $Header$ @@*/ -#undef AEILOCALINTERP_DEBUG /* define this for verbose debugging output */ -#undef AEILOCALINTERP_DEBUG2 /* define this for even more verbose debugging output */ - /******************************************************************************/ /* @@ -353,9 +350,11 @@ int FUNCTION_NAME(/***** coordinate system *****/ int N_output_arrays, const CCTK_INT output_array_type_codes[], void* const output_arrays[], - /***** operation info */ + /***** operation info *****/ const CCTK_INT operand_indices[], const CCTK_INT operation_codes[], + /***** debugging *****/ + int debug, /***** other return results *****/ struct error_info* error_info, struct molecule_structure_flags* molecule_structure_flags, @@ -757,10 +756,13 @@ int out; 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 + +if (debug > 0) + then { + printf("AEILocalInterp:: in interpolator fn: N_DIMS=%d N_interp_points=%d\n", (int) N_DIMS, N_interp_points); + fflush(stdout); + } + for (pt = 0 ; pt < N_interp_points ; ++pt) { /* this struct holds a molecule-sized piece of a single */ @@ -880,24 +882,25 @@ fflush(stdout); } } -#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 - fflush(stdout); -#endif /* AEILOCALINTERP_DEBUG2 */ + if (debug >= 6) + then { + #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 + fflush(stdout); + } /* * compute position of interpolation molecules with respect to @@ -921,20 +924,21 @@ fflush(stdout); { const int ibndry_min = 2*axis; const int ibndry_max = 2*axis + 1; - #ifdef AEILOCALINTERP_DEBUG2 - printf("axis=%d ibndry_{min,max}=%d,%d MOLECULE_SIZE=%d\n", - axis, ibndry_min, ibndry_max, MOLECULE_SIZE); - printf(" coord_origin[%d]=%g coord_delta[%d]=%g\n", - axis, (double)coord_origin[axis], - axis, (double)coord_delta[axis]); - printf(" input_array_[min,max]_subscripts[%d]=[%d,%d]\n", - axis, - (int)input_array_min_subscripts[axis], - (int)input_array_max_subscripts[axis]); - printf(" N_boundary_points_to_omit[ibndry_[min,max]]=[%d,%d]\n", - (int)N_boundary_points_to_omit[ibndry_min], - (int)N_boundary_points_to_omit[ibndry_max]); - #endif + if (debug >= 8) + then { + printf("axis=%d ibndry_{min,max}=%d,%d MOLECULE_SIZE=%d\n", + axis, ibndry_min, ibndry_max, MOLECULE_SIZE); + printf(" coord_origin[%d]=%g coord_delta[%d]=%g\n", + axis, (double)coord_origin[axis], + axis, (double)coord_delta[axis]); + printf(" input_array_[min,max]_subscripts[%d]=[%d,%d]\n", + axis, + (int)input_array_min_subscripts[axis], + (int)input_array_max_subscripts[axis]); + printf(" N_boundary_points_to_omit[ibndry_[min,max]]=[%d,%d]\n", + (int)N_boundary_points_to_omit[ibndry_min], + (int)N_boundary_points_to_omit[ibndry_max]); + } const int mp_status = AEILocalInterp_molecule_posn (coord_origin[axis], coord_delta[axis], input_array_min_subscripts[axis] @@ -947,7 +951,10 @@ fflush(stdout); boundary_extrapolation_tolerance[ibndry_min], boundary_extrapolation_tolerance[ibndry_max], interp_coords_fp[axis], + debug, ¢er_ijk_temp[axis], &xyz_temp[axis]); + if (debug >= 8) + then printf("==> got mp_status=%d\n", mp_status); /* * unfortunately, the status codes from AEILocalInterp_molecule_posn() @@ -1054,10 +1061,11 @@ fflush(stdout); 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 + if (debug >= 6) + then { + printf("AEILocalInterp:: pt=%d has this_point_status=%d ==> skipping it\n", pt, this_point_status); + fflush(stdout); + } /* this point is in error (eg outside grid) */ /* ==> don't try to interpolate at this point! */ @@ -1078,22 +1086,25 @@ fflush(stdout); const fp z = xyz_temp [Z_AXIS]; #endif -#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); - #else - #error "N_DIMS may not be > 3!" - #endif - fflush(stdout); -#endif + if (debug >= 6) + then { + #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); + #else + #error "N_DIMS may not be > 3!" + #endif + fflush(stdout); + } /* * compute 1-d position of molecule center in input arrays @@ -1277,28 +1288,30 @@ case CCTK_VARIABLE_REAL: LOAD_DATA_REAL(input_array_ptr_real, STRIDE_IJK, &data); - #if defined(AEILOCALINTERP_DEBUG2) \ - && (N_DIMS == 2) && (MOLECULE_SIZE == 4) - /* we assume a cubical molecule :( */ - printf("AEILocalInterp:: loaded data is:\n"); - printf(" data_m1_m1 = %g\n", (double) data.data_m1_m1); - printf(" data_0_m1 = %g\n", (double) data.data_0_m1); - printf(" data_p1_m1 = %g\n", (double) data.data_p1_m1); - printf(" data_p2_m1 = %g\n", (double) data.data_p2_m1); - printf(" data_m1_0 = %g\n", (double) data.data_m1_0); - printf(" data_0_0 = %g\n", (double) data.data_0_0); - printf(" data_p1_0 = %g\n", (double) data.data_p1_0); - printf(" data_p2_0 = %g\n", (double) data.data_p2_0); - printf(" data_m1_p1 = %g\n", (double) data.data_m1_p1); - printf(" data_0_p1 = %g\n", (double) data.data_0_p1); - printf(" data_p1_p1 = %g\n", (double) data.data_p1_p1); - printf(" data_p2_p1 = %g\n", (double) data.data_p2_p1); - printf(" data_m1_p2 = %g\n", (double) data.data_m1_p2); - printf(" data_0_p2 = %g\n", (double) data.data_0_p2); - printf(" data_p1_p2 = %g\n", (double) data.data_p1_p2); - printf(" data_p2_p2 = %g\n", (double) data.data_p2_p2); - fflush(stdout); - #endif + #if (N_DIMS == 2) && (MOLECULE_SIZE == 4) + if (debug >= 10) + then { + /* we assume a cubical molecule :( */ + printf("AEILocalInterp:: loaded data is:\n"); + printf(" data_m1_m1 = %g\n", (double) data.data_m1_m1); + printf(" data_0_m1 = %g\n", (double) data.data_0_m1); + printf(" data_p1_m1 = %g\n", (double) data.data_p1_m1); + printf(" data_p2_m1 = %g\n", (double) data.data_p2_m1); + printf(" data_m1_0 = %g\n", (double) data.data_m1_0); + printf(" data_0_0 = %g\n", (double) data.data_0_0); + printf(" data_p1_0 = %g\n", (double) data.data_p1_0); + printf(" data_p2_0 = %g\n", (double) data.data_p2_0); + printf(" data_m1_p1 = %g\n", (double) data.data_m1_p1); + printf(" data_0_p1 = %g\n", (double) data.data_0_p1); + printf(" data_p1_p1 = %g\n", (double) data.data_p1_p1); + printf(" data_p2_p1 = %g\n", (double) data.data_p2_p1); + printf(" data_m1_p2 = %g\n", (double) data.data_m1_p2); + printf(" data_0_p2 = %g\n", (double) data.data_0_p2); + printf(" data_p1_p2 = %g\n", (double) data.data_p1_p2); + printf(" data_p2_p2 = %g\n", (double) data.data_p2_p2); + fflush(stdout); + } + #endif break; } @@ -1510,13 +1523,15 @@ 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 */ + if (debug >= 10) then { - printf("AEILocalInterp:: pt=%d out=%d --> storing result %g\n", pt, out, (double) result); - fflush(stdout); + 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]) { @@ -1524,16 +1539,17 @@ case CCTK_VARIABLE_REAL: { CCTK_REAL *const output_array_ptr_real = (CCTK_REAL *) output_arrays[out]; - #ifdef AEILOCALINTERP_DEBUG2 - if ((pt & (pt-1)) == 0) /* pt is 0 or a power of 2 */ + if (debug >= 10) then { - printf(" result addr is %p\n", - (void *) &output_array_ptr_real[pt]); - printf(" previous value there was %g\n", - output_array_ptr_real[pt]); - fflush(stdout); + if ((pt & (pt-1)) == 0) /* pt is 0 or a power of 2 */ + then { + printf(" result addr is %p\n", + (void *) &output_array_ptr_real[pt]); + printf(" previous value there was %g\n", + output_array_ptr_real[pt]); + fflush(stdout); + } } - #endif output_array_ptr_real[pt] = (CCTK_REAL) result; break; } diff --git a/src/template.h b/src/template.h index d1f9d92..3290cf2 100644 --- a/src/template.h +++ b/src/template.h @@ -29,9 +29,11 @@ int FUNCTION_NAME(/***** coordinate system *****/ int N_output_arrays, const CCTK_INT output_array_type_codes[], void* const output_arrays[], - /***** operation info */ + /***** operation info *****/ const CCTK_INT operand_indices[], const CCTK_INT operation_codes[], + /***** debugging *****/ + int debug, /***** other return results *****/ struct error_info* error_info, struct molecule_structure_flags* molecule_structure_flags, diff --git a/src/test_molecule_posn.c b/src/test_molecule_posn.c index f8ff218..831ad14 100644 --- a/src/test_molecule_posn.c +++ b/src/test_molecule_posn.c @@ -284,6 +284,7 @@ printf(" boundary_extrapolation_tolerance_[min,max]=[%g,%g]\n", printf(" x=%g\n", (double) x); { +const int debug = 0; const int status = AEILocalInterp_molecule_posn(grid_x0, grid_dx, grid_i_min, grid_i_max, @@ -293,6 +294,7 @@ const int status boundary_extrapolation_tolerance_min, boundary_extrapolation_tolerance_max, x, + debug, &i_center, &x_rel); if (status < 0) @@ -321,6 +323,7 @@ int failure_count = 0; const struct args_results *p = & test_data[i]; int i_center = 999; fp x_rel = 999.999; + const int debug = 0; const int status = AEILocalInterp_molecule_posn (grid_x0, grid_dx, grid_i_min, grid_i_max, @@ -330,6 +333,7 @@ int failure_count = 0; p->boundary_extrapolation_tolerance_min, p->boundary_extrapolation_tolerance_max, p->x, + debug, &i_center, &x_rel); bool ok = (status == 0) ? ( (status == p->status) && (i_center == p->i_center) -- cgit v1.2.3