diff options
Diffstat (limited to 'src/template.c')
-rw-r--r-- | src/template.c | 206 |
1 files changed, 111 insertions, 95 deletions
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; } |