diff options
Diffstat (limited to 'src/template.c')
| -rw-r--r-- | src/template.c | 123 |
1 files changed, 116 insertions, 7 deletions
diff --git a/src/template.c b/src/template.c index fe7a267..5cb439d 100644 --- a/src/template.c +++ b/src/template.c @@ -279,6 +279,12 @@ If you change the arguments for this function, note that you must also change the prototype in "template.h". + @var log_fp + @vdesc This is NULL for no logging, or a valid stdio stream pointer + to enable logging of the grid and interpolation coordinates. + @vtype FILE* + @vio in + @var error_info @vdesc If error_info->per_point_status is non-NULL, then we store per-point interpolator status for @@ -354,7 +360,7 @@ int FUNCTION_NAME(/***** coordinate system *****/ const CCTK_INT operand_indices[], const CCTK_INT operation_codes[], /***** debugging *****/ - int debug, + int debug, FILE* log_fp, /***** other return results *****/ struct error_info* error_info, struct molecule_structure_flags* molecule_structure_flags, @@ -750,6 +756,60 @@ int out; #error "N_DIMS may not be > 3!" #endif + +/* + * if requested, write the interpolation grid to the log file + * (we'll write (some of) the individual interpolation coordinates + * as we interpolate them (below)) + */ +if (log_fp != NULL) + then { + fprintf(log_fp, "%d\t%d", (int) N_DIMS, N_interp_points); + + { + fp grid_min_xyz[MAX_N_DIMS], grid_max_xyz[MAX_N_DIMS]; + int axis; + for (axis = 0 ; axis < N_DIMS ; ++axis) + { + grid_min_xyz[axis] + = coord_origin[axis] + + input_array_min_subscripts[axis]*coord_delta[axis]; + grid_max_xyz[axis] + = coord_origin[axis] + + input_array_max_subscripts[axis]*coord_delta[axis]; + } + + #if (N_DIMS == 1) + fprintf(log_fp, "\t%g\t%g\t%g\t0.0\t0.0\t0.0\t0.0\t0.0\t0.0", + (double) grid_min_xyz[X_AXIS], + (double) coord_delta [X_AXIS], + (double) grid_max_xyz[X_AXIS]); + #elif (N_DIMS == 2) + fprintf(log_fp, "\t%g\t%g\t%g\t%g\t%g\t%g\t0.0\t0.0\t0.0", + (double) grid_min_xyz[X_AXIS], + (double) coord_delta [X_AXIS], + (double) grid_max_xyz[X_AXIS], + (double) grid_min_xyz[Y_AXIS], + (double) coord_delta [Y_AXIS], + (double) grid_max_xyz[Y_AXIS]); + #elif (N_DIMS == 3) + fprintf(log_fp, "\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g", + (double) grid_min_xyz[X_AXIS], + (double) coord_delta [X_AXIS], + (double) grid_max_xyz[X_AXIS], + (double) grid_min_xyz[Y_AXIS], + (double) coord_delta [Y_AXIS], + (double) grid_max_xyz[Y_AXIS], + (double) grid_min_xyz[Z_AXIS], + (double) coord_delta [Z_AXIS], + (double) grid_max_xyz[Z_AXIS]); + #else + #error "N_DIMS may not be > 3!" + #endif + } + } + + /* * interpolate at each point */ @@ -886,16 +946,16 @@ if (debug > 0) then { #if (N_DIMS == 1) printf("pt=%d interp coord %.16g\n", - pt, (double) interp_coords_fp[0]); + pt, (double) interp_coords_fp[X_AXIS]); #elif (N_DIMS == 2) printf("pt=%d interp coords %.16g %.16g\n", - pt, (double) interp_coords_fp[0], - (double) interp_coords_fp[1]); + pt, (double) interp_coords_fp[X_AXIS], + (double) interp_coords_fp[Y_AXIS]); #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]); + pt, (double) interp_coords_fp[X_AXIS], + (double) interp_coords_fp[Y_AXIS], + (double) interp_coords_fp[Z_AXIS]); #else #error "N_DIMS may not be > 3!" #endif @@ -903,6 +963,29 @@ if (debug > 0) } /* + * if requested, write this point's interpolation coords to the log file + */ + if ( (log_fp != NULL) + && ((pt < 3) || (pt == N_interp_points-1)) ) + then { + #if (N_DIMS == 1) + fprintf(log_fp, "\t%g\t0.0\t0.0", + (double) interp_coords_fp[X_AXIS]); + #elif (N_DIMS == 2) + fprintf(log_fp, "\t%g\t%g\t0.0", + (double) interp_coords_fp[X_AXIS], + (double) interp_coords_fp[Y_AXIS]); + #elif (N_DIMS == 3) + fprintf(log_fp, "\t%g\t%g\t%g", + (double) interp_coords_fp[X_AXIS], + (double) interp_coords_fp[Y_AXIS], + (double) interp_coords_fp[Z_AXIS]); + #else + #error "N_DIMS may not be > 3!" + #endif + } + + /* * compute position of interpolation molecules with respect to * the grid, i.e. compute (x,y,z) coordinates of interpolation * point relative to molecule center, in units of the grid spacing @@ -1760,6 +1843,32 @@ CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, /* end of for (pt = ...) loop */ } +/* + * if we're writing logging information, finish this + */ +if (log_fp != NULL) + then { + switch (N_interp_points) + { + case 0: + /* we printed 0 points, but need 4 */ + fprintf(log_fp, "\t0.0\t0.0\t0.0"); + fprintf(log_fp, "\t0.0\t0.0\t0.0"); + /* fall through */ + case 1: + /* we printed 2 points (0, N_interp_points-1), but need 4 */ + fprintf(log_fp, "\t0.0\t0.0\t0.0"); + /* fall through */ + case 2: + /* we printed 3 points (0, 1, N_interp_points-1), but need 4 */ + fprintf(log_fp, "\t0.0\t0.0\t0.0"); + /* fall through */ + } + + fprintf(log_fp, "\n"); + fflush(log_fp); + } + return return_status; } } |