add a debugging feature: setting the parameter

  AEILocalInterp::log_interp_coords = true
causes this thorn to write a log file on each processor
giving the grid and (some of) the interpolation coordinates


git-svn-id: http://svn.aei.mpg.de/numrel/AEIThorns/AEILocalInterp/trunk@45 0f49ee68-0e4f-0410-9b9c-b2c123ded7ef

7 files changed, 266 insertions, 11 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index 78eb921..15b94e5 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -1392,6 +1392,25 @@ to a positive value turns on various debug printing inside the
 interpolator.  This is intended for debugging the interpolator itself;
 you need to look at the interpolator source code to interpret the output.
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\subsection{Logging}
+
+This thorn has a Boolean parameter \verb|log_interp_coords|.
+If this is set to \verb|true|, this thorn will write a detailed log
+file giving the grid information and interpolation coordinates for
+each interpolator call.  Each processor will write a separate log file,
+with the file name determined via
+\begin{verbatim}
+snprintf(buffer, length,
+         "AEILocalInterp.proc%d.log"
+         CCTK_MyProc(NULL))
+\end{verbatim}
+
+Note that these log files are typically quite large, and writing them         
+may significantly slow the simulation -- you should probably only use
+this option for debugging purposes.
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{Supressing Interpolation}
diff --git a/param.ccl b/param.ccl
index 9f0ce17..a60db7c 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,4 +1,22 @@
 # Parameter definitions for thorn AEILocalInterp
 # $Header$
 
-# there are no parameters for this thorn
+# all parameters are private to this thorn
+private:
+
+#
+# If this logging is done, each processor writes its log data to a file
+# named via
+#	snprintf(buffer, length,
+#		 "AEILocalInterp.proc%d.log"
+#		 CCTK_MyProc())
+# Note that these log files are typically quite large, and writing them
+# may significantly slow the simulation -- you should probably only use
+# this option for debugging purposes.
+#
+Boolean log_interp_coords						\
+  "should we log the grid min(delta)max and the interpolation		\
+   coordinates for each call on the interpolator?"			\
+	STEERABLE=ALWAYS
+{
+} false
diff --git a/src/InterpLocalUniform.c b/src/InterpLocalUniform.c
index 30a5f9a..7cc6ff0 100644
--- a/src/InterpLocalUniform.c
+++ b/src/InterpLocalUniform.c
@@ -70,7 +70,7 @@ CCTK_FILEVERSION(AEITHorns_AEILocalInterp_src_InterpLocalUniform_c)
 /******************************************************************************/
 
 /*
- * *****data structures and functions local to this file *****
+ * ***** data structures and functions local to this file *****
  */
 
 /**************************************/
@@ -1112,6 +1112,63 @@ if (N_dims > MAX_N_DIMS)
 /******************************************************************************/
 
 /*
+ * if logging is requested,
+ *    open a logging file if we haven't already done so
+ */
+  {
+static FILE* log_fp = NULL;
+const bool log_interp_coords
+   = (AEILocalInterp_get_int_param("AEILocalInterp", "log_interp_coords") != 0);
+
+if (log_interp_coords && (log_fp == NULL))
+   then {
+	/* open logging file & write header */
+
+	#define FILE_NAME_BUFFER_SIZE	100
+	char file_name[FILE_NAME_BUFFER_SIZE];
+	Util_snprintf(file_name, FILE_NAME_BUFFER_SIZE,
+		      "AEILocalInterp.proc%d.log",
+		      CCTK_MyProc(NULL));
+
+	log_fp = fopen(file_name, "w");
+	if (log_fp == NULL)
+	   then CCTK_VWarn(ERROR_MSG_SEVERITY_LEVEL,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+"\n"
+"   CCTK_InterpLocalUniform(): can't open coordinate logging file\n"
+"                              \"%s\"!\n"
+			   ,
+			   file_name);				/*NOTREACHED*/
+
+	fprintf(log_fp, "## any missing data (eg if N_dims < 3 or N_interp_points < 3) is printed as 0\n");
+	fprintf(log_fp, "# column  1 = N_dims\n");
+	fprintf(log_fp, "# column  2 = N_interp_points\n");
+	fprintf(log_fp, "# column  3 = grid x_min\n");
+	fprintf(log_fp, "# column  4 = grid delta_x\n");
+	fprintf(log_fp, "# column  5 = grid x_max\n");
+	fprintf(log_fp, "# column  6 = grid y_min\n");
+	fprintf(log_fp, "# column  7 = grid delta_y\n");
+	fprintf(log_fp, "# column  8 = grid y_max\n");
+	fprintf(log_fp, "# column  9 = grid z_min\n");
+	fprintf(log_fp, "# column 10 = grid delta_z\n");
+	fprintf(log_fp, "# column 11 = grid z_max\n");
+	fprintf(log_fp, "# column 12 = interp_x[0]\n");
+	fprintf(log_fp, "# column 13 = interp_y[0]\n");
+	fprintf(log_fp, "# column 14 = interp_z[0]\n");
+	fprintf(log_fp, "# column 15 = interp_x[1]\n");
+	fprintf(log_fp, "# column 16 = interp_y[1]\n");
+	fprintf(log_fp, "# column 17 = interp_z[1]\n");
+	fprintf(log_fp, "# column 18 = interp_x[2]\n");
+	fprintf(log_fp, "# column 19 = interp_y[2]\n");
+	fprintf(log_fp, "# column 20 = interp_z[2]\n");
+	fprintf(log_fp, "# column 21 = interp_x[N_interp_points-1]\n");
+	fprintf(log_fp, "# column 22 = interp_y[N_interp_points-1]\n");
+	fprintf(log_fp, "# column 23 = interp_z[N_interp_points-1]\n");
+	}
+   
+/******************************************************************************/
+
+/*
  * get the parameters from the parameter table, filling in defaults as needed
  */
 
@@ -1567,7 +1624,7 @@ const int return_code
 			 N_output_arrays,
 			    output_array_type_codes, output_arrays,
 			    operand_indices, operation_codes,
-			 debug,
+			 debug, (log_interp_coords ? log_fp : NULL),
 			 &error_info,
 			 &molecule_structure_flags,
 			 p_molecule_min_max_m_info,
@@ -1675,6 +1732,7 @@ return return_code;					/*** NORMAL RETURN ***/
   }
   }
   }
+  }
 }
 
 /******************************************************************************/
diff --git a/src/InterpLocalUniform.h b/src/InterpLocalUniform.h
index d3e23c2..6d84cbb 100644
--- a/src/InterpLocalUniform.h
+++ b/src/InterpLocalUniform.h
@@ -262,3 +262,7 @@ int AEILocalInterp_molecule_posn(fp grid_origin, fp grid_delta,
 
 /* functions in "util.c" */
 int AEILocalInterp_decode_N_parts(int type_code);
+int AEILocalInterp_get_int_param(const char* thorn_or_implementation_name,
+				 const char* parameter_name);
+
+/******************************************************************************/
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;
   }
   }
diff --git a/src/template.h b/src/template.h
index 3290cf2..c98ab4e 100644
--- a/src/template.h
+++ b/src/template.h
@@ -33,7 +33,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,
diff --git a/src/util.c b/src/util.c
index a4bc5eb..ac1a3ee 100644
--- a/src/util.c
+++ b/src/util.c
@@ -78,3 +78,50 @@ case CCTK_VARIABLE_FPOINTER:	return 0;
 default:			return -1;
 	}
 }
+
+/******************************************************************************/
+
+/*@@
+  @routine  AEILocalInterp_get_int_param
+  @date     12 Jan 2007
+  @author   Jonathan Thornburg <jthorn@aei.mpg.de>
+  @desc     This function gets the (scalar) value of a CCTK_INT or Boolean
+            Cactus parameter.
+  @enddesc
+
+  @var      thorn_or_implementation_name
+  @vdesc    The name of the thorn (for private parameters) or implementation
+            (for restricted parameters), eg "AEILocalInterp".
+  @vtype    const char*
+  @endvar
+
+  @var      parameter_name
+  @vdesc    The name of the parameter, e.g. "log_interp_coords"
+  @vtype    const char*
+  @endvar
+
+  @returntype	int
+  @returndesc	This function returns the value of the parameter.
+                If an error occurs, this function calls
+                  CCTK_VWarn(CCTK_WARN_ABORT, ...)
+                (and does not return to the caller).
+  @endreturndesc
+  @@*/
+int AEILocalInterp_get_int_param(const char* const thorn_or_implementation_name,
+				 const char* const parameter_name)
+{
+CCTK_INT data_type;
+const CCTK_INT* const value_ptr
+	= (const CCTK_INT*) CCTK_ParameterGet(parameter_name,
+					      thorn_or_implementation_name,
+                                              NULL);
+
+if (value_ptr == NULL)
+   then CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING,
+"***** AEILocalInterp_decode_N_parts():\n"
+"        can't get value of parameter %s::%s!\n"
+		   ,
+		   thorn_or_implementation_name, parameter_name); /*NOTREACHED*/
+
+return *value_ptr;
+}