aboutsummaryrefslogtreecommitdiff
path: root/src/template.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/template.c')
-rw-r--r--src/template.c206
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,
&center_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;
}