path: root/src/driver/io.cc

// io.cc -- I/O routines for this thorn
// $Header$
//
// <<<prototypes for functions local to this file>>>
// input_gridfn - read an angular grid function from a data file
// output_gridfn - write an angular grid function to a data file
// output_Jacobians - write a Jacobian matrix or matrices to a data file
// decode_horizon_file_format - decode the horizon_file_format parameter
//
// setup_BH_diagnostics_output_file - create/initialize BH-diagnostics file
// output_BH_diagnostics_fn - append BH diagnostics to file
//
/// io_file_name
//

#include <stdio.h>
#include <assert.h>
#include <math.h>

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

#include "stl_vector.hh"

#include "config.h"
#include "stdc.h"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;

#include "../patch/coords.hh"
#include "../patch/grid.hh"
#include "../patch/fd_grid.hh"
#include "../patch/patch.hh"
#include "../patch/patch_edge.hh"
#include "../patch/patch_interp.hh"
#include "../patch/ghost_zone.hh"
#include "../patch/patch_system.hh"

#include "../elliptic/Jacobian.hh"

#include "../gr/gfns.hh"
#include "../gr/gr.hh"

#include "driver.hh"

//******************************************************************************

//
// ***** prototypes for functions local to this file ***** 
//

namespace {
const char* io_file_name(const struct IO_info& IO_info,
			 const char base_file_name[],
			 int hn, int AHF_iteration = 0);
	  }

//******************************************************************************

//
// This function inputs a gridfn from a data file.
//
// We assume that this gridfn is ghosted, but the ghost zones are *not*
// present in the data file.
//
void input_gridfn(patch_system& ps, int unknown_gfn,
		  const struct IO_info& IO_info, const char base_file_name[],
		  int hn, bool print_msg_flag, int AHF_iteration /* = 0 */)
{
const char* file_name = io_file_name(IO_info, base_file_name,
				     hn, AHF_iteration);
if (print_msg_flag)
   then {
	if ( (unknown_gfn == gfns::gfn__h)
	     && (IO_info.time_iteration == 0) && (AHF_iteration == 0) )
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "   reading initial guess from \"%s\"", file_name);
	   else CCTK_VInfo(CCTK_THORNSTRING,
			   "   reading \"%s\"", file_name);
	}

switch	(IO_info.horizon_file_format)
	{
case horizon_file_format__ASCII_gnuplot:
	ps.read_ghosted_gridfn(unknown_gfn,
			       file_name,
			       false);		// no ghost zones in data file
	break;

case horizon_file_format__HDF5:
	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"input_gridfn(): HDF5 data files not implemented yet!");	/*NOTREACHED*/

default:
	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   input_gridfn(): unknown IO_info.horizon_file_format=(int)%d!\n"
"                   (this should never happen!)"
		   ,
		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
	}
}

//******************************************************************************

//
// This function outputs a gridfn from a data file.
//
// If the gridfn is h, then we also write out the xyz positions of the
// horizon surface points.
//
// FIXME: if the gridfn is not h, we assume that it's nominal-grid.
//
void output_gridfn(patch_system& ps, int unknown_gfn,
		   const struct IO_info& IO_info, const char base_file_name[],
		   int hn, bool print_msg_flag, int AHF_iteration /* = 0 */)
{
const char* file_name = io_file_name(IO_info, base_file_name,
				     hn, AHF_iteration);

switch	(IO_info.horizon_file_format)
	{
case horizon_file_format__ASCII_gnuplot:
	switch	(unknown_gfn)
		{
	case gfns::gfn__h:
		if (print_msg_flag)
		   then CCTK_VInfo(CCTK_THORNSTRING,
				   "   writing h to \"%s\"", file_name);
		ps.print_ghosted_gridfn_with_xyz
		      (unknown_gfn,
		       true, gfns::gfn__h,
		       file_name,
		       IO_info.output_ghost_zones_for_h); // should we include
							  // ghost zones?
		break;
	case gfns::gfn__Theta:
		if (print_msg_flag)
		   then CCTK_VInfo(CCTK_THORNSTRING,
				   "   writing Theta to \"%s\"", file_name);
		ps.print_gridfn(unknown_gfn, file_name);
		break;
	case gfns::gfn__Delta_h:
		if (print_msg_flag)
		   then CCTK_VInfo(CCTK_THORNSTRING,
				   "   writing Delta_h to \"%s\"", file_name);
		ps.print_gridfn(unknown_gfn, file_name);
		break;
	default:
		if (print_msg_flag)
		   then CCTK_VInfo(CCTK_THORNSTRING,
				   "   writing gfn=%d to \"%s\"",
				   unknown_gfn, file_name);
		ps.print_gridfn(unknown_gfn, file_name);
		break;
		}
	break;

case horizon_file_format__HDF5:
	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"output_gridfn(): HDF5 data files not implemented yet!");	/*NOTREACHED*/

default:
	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   output_gridfn(): unknown IO_info.horizon_file_format=(int)%d!\n"
"                    (this should never happen!)"
		   ,
		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
	}
}

//******************************************************************************

//
// This function prints one or two Jacobian matrices (and their difference
// in the latter case) to a named output file.
//
// Bugs:
// - The file format is hardwired to ASCII.
//
// Arguments:
// A null Jacobian pointer means to skip that Jacobian.
//
void print_Jacobians(const patch_system& ps,
		     const Jacobian* Jac_NP, const Jacobian* Jac_SD_FDdr,
		     const struct IO_info& IO_info, const char base_file_name[],
		     int hn, bool print_msg_flag, int AHF_iteration /* = 0 */)
{
if (Jac_NP == NULL)
   then {
	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "print_Jacobinas(): Jac_NP == NULL is not (yet) supported!");
	return;						// *** ERROR RETURN ***
	}

const char* file_name = io_file_name(IO_info, base_file_name,
				     hn, AHF_iteration);
if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "   writing %s to \"%s\"",
		   ((Jac_SD_FDdr == NULL) ? "NP Jacobian"
					  : "NP and SD_FDdr Jacobians"),
		   file_name);

FILE *fileptr = fopen(file_name, "w");
if (fileptr == NULL)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "print_Jacobians(): can't open output file \"%s\"!",
		   file_name);					/*NOTREACHED*/

fprintf(fileptr, "# column 1 = x patch name\n");
fprintf(fileptr, "# column 2 = x patch number\n");
fprintf(fileptr, "# column 3 = x irho\n");
fprintf(fileptr, "# column 4 = x isigma\n");
fprintf(fileptr, "# column 5 = x II\n");
fprintf(fileptr, "# column 6 = y patch name\n");
fprintf(fileptr, "# column 7 = y patch number\n");
fprintf(fileptr, "# column 8 = y irho\n");
fprintf(fileptr, "# column 9 = y isigma\n");
fprintf(fileptr, "# column 10 = y JJ\n");
fprintf(fileptr, "# column 11 = Jac_NP(II,JJ)\n");
if (Jac_SD_FDdr != NULL)
   then {
	fprintf(fileptr, "# column 12 = Jac_SD_FDdr(II,JJ)\n");
	fprintf(fileptr, "# column 13 = abs error in Jac_SD_FDdr\n");
	fprintf(fileptr, "# column 14 = rel error in Jac_SD_FDdr\n");
	}

    for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
    {
    const patch& xp = ps.ith_patch(xpn);

    for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
    {
    for (int x_isigma = xp.min_isigma() ;
	 x_isigma <= xp.max_isigma() ;
	 ++x_isigma)
    {
    const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);

	for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
	{
	const patch& yp = ps.ith_patch(ypn);

	for (int y_irho = yp.min_irho() ;
	     y_irho <= yp.max_irho() ;
	     ++y_irho)
	{
	for (int y_isigma = yp.min_isigma() ;
	     y_isigma <= yp.max_isigma() ;
	     ++y_isigma)
	{
	const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);

	if (! Jac_NP->is_explicitly_stored(II,JJ))
	   then continue;			// skip sparse points

	const fp NP = (*Jac_NP)(II,JJ);
	const fp SD_FDdr = (Jac_SD_FDdr == NULL) ? 0.0 : (*Jac_SD_FDdr)(II,JJ);

	if ((NP == 0.0) && (SD_FDdr == 0.0))
	   then continue;		// skip zero values (if == )

	fprintf(fileptr,
		"%s %d %d %d %d\t%s %d %d %d %d\t%#.10g",
		xp.name(), xpn, x_irho, x_isigma, II,
		yp.name(), ypn, y_irho, y_isigma, JJ,
		double(NP));

	if (Jac_SD_FDdr != NULL)
	   then {
		const fp abs_NP      = jtutil::abs(NP     );
		const fp abs_SD_FDdr = jtutil::abs(SD_FDdr);
		const fp max_abs = jtutil::max(abs_SD_FDdr, abs_NP);
		const fp SD_FDdr_abs_error = SD_FDdr - NP;
		const fp SD_FDdr_rel_error = SD_FDdr_abs_error / max_abs;

		fprintf(fileptr,
			"\t%#.10g\t%e\t%e",
			double(SD_FDdr),
			double(SD_FDdr_abs_error), double(SD_FDdr_rel_error));
		}

	fprintf(fileptr, "\n");
	}
	}
	}
    }
    }
    }

fclose(fileptr);
}

//******************************************************************************

//
// This function decodes the  horizon_file_format  parameter (string) into an
// internal enum for future use.
//
enum horizon_file_format
  decode_horizon_file_format(const char horizon_file_format_string[])
{
if	(STRING_EQUAL(horizon_file_format_string, "ASCII (gnuplot)"))
   then return horizon_file_format__ASCII_gnuplot;
else if (STRING_EQUAL(horizon_file_format_string, "HDF5"))
   then return horizon_file_format__HDF5;
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"decode_horizon_file_format(): unknown horizon_file_format_string=\"%s\"!",
		   horizon_file_format_string);			/*NOTREACHED*/
}

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function creates a BH-diagnostics output file, writes a suitable
// header comment, and returns a stdio file pointer which can be used to
// append data to the file.
//
FILE* setup_BH_diagnostics_output_file(const struct IO_info& IO_info,
				       int hn, int N_horizons)
{
const int N_file_name_buffer = 200;
char file_name_buffer[N_file_name_buffer];

snprintf(file_name_buffer, N_file_name_buffer,
	 "%s.ah%d.%s",
	 IO_info.BH_diagnostics_base_file_name,
	 hn, IO_info.BH_diagnostics_file_name_extension);
FILE *fileptr = fopen(file_name_buffer, "w");
if (fileptr == NULL)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   setup_BH_diagnostics_output_file():\n"
"        can't open BH-diagnostics output file\n"
"        \"%s\"!"
		   ,
		   file_name_buffer);				/*NOTREACHED*/

fprintf(fileptr, "# apparent horizon %d/%d\n", hn, N_horizons);
fprintf(fileptr, "#\n");
fprintf(fileptr, "# column 1 = cctk_iteration\n");
fprintf(fileptr, "# column 2 = cctk_time\n");
fprintf(fileptr, "# column 3 = centroid_x\n");
fprintf(fileptr, "# column 4 = centroid_y\n");
fprintf(fileptr, "# column 5 = centroid_z\n");
fprintf(fileptr, "# column 6 = mean radius\n");
fprintf(fileptr, "# column 7 = xy-plane circumference\n");
fprintf(fileptr, "# column 8 = xz-plane circumference\n");
fprintf(fileptr, "# column 9 = yz-plane circumference\n");
fprintf(fileptr, "# column 10 = area\n");
fprintf(fileptr, "# column 11 = m_irreducible\n");
fflush(fileptr);

return fileptr;
}

//******************************************************************************

//
// This function appends a BH-diagnostics line to a data file.  It
// attempts to ensure that the newly-written line is flushed to disk,
// so the output file can be examined while the Cactus run is still in
// progress.
//
// (The "_fn" in the function name is to distinguish it from the Cactus
// parameter  output_BH_diagnostics .)
//
// Arguments:
// BH_diagnostics = The BH diagnostics to be written
// fileptr = The stdio file pointer to append to
//
void output_BH_diagnostics_fn(const struct BH_diagnostics& BH_diagnostics,
			      const struct IO_info& IO_info,
			      int hn, FILE* fileptr)
{
assert(fileptr != NULL);

fprintf(fileptr,
//  cctk_iteration        mean radius                     area    m_irreducible
//  ==  cctk_time         ======  {xy,xz,yz}-plane        ======  ======
//  ==  ====  centroid_[xyz]      circumferences          ======  ======
//  ==  ====  ==========  ======  ======================  ======  ======
   "%d\t%.3f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n",
	IO_info.time_iteration, double(IO_info.time),
	BH_diagnostics.centroid_x, BH_diagnostics.centroid_y,
				   BH_diagnostics.centroid_z,
	BH_diagnostics.mean_radius,
	BH_diagnostics.circumference_xy, BH_diagnostics.circumference_xz,
					 BH_diagnostics.circumference_yz,
	BH_diagnostics.area,
	BH_diagnostics.m_irreducible);

fflush(fileptr);
}

//******************************************************************************
//******************************************************************************
//******************************************************************************

//
// This function encapsulates our file-naming conventions for "horizon"
// output files (those used for h, H, and other angular grid functions).
//
// Arguments:
// base_file_name[] = from the parameter file
// hn = the horizon number
// AHF_iteration = the apparent horizon finder's internal iteration
//		   number (>= 1) if this is an intermediate iterate,
//		   or the default (0) if this is a final computed
//		   horizon position
//
// Results:
// This function returns (a pointer to) the file name.  The returned
// result points into a private static buffer; the usual caveats apply.
//
namespace {
const char* io_file_name(const struct IO_info& IO_info,
			 const char base_file_name[],
			 int hn, int AHF_iteration /* = 0 */)
{
const int N_file_name_buffer = 200;
static char file_name_buffer[N_file_name_buffer];

const char* file_name_extension;
switch	(IO_info.horizon_file_format)
	{
case horizon_file_format__ASCII_gnuplot:
	file_name_extension = IO_info.ASCII_gnuplot_file_name_extension;
	break;
case horizon_file_format__HDF5:
	file_name_extension = IO_info.HDF5_file_name_extension;
	break;
default:
	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   io_file_name(): unknown IO_info.horizon_file_format=(int)%d!\n"
"                   (this should never happen!)"
		   ,
		   int(IO_info.horizon_file_format));		/*NOTREACHED*/
	}

if (AHF_iteration == 0)
   then snprintf(file_name_buffer, N_file_name_buffer,
		 "%s.t%d.ah%d.%s",
		 base_file_name,
		 IO_info.time_iteration, hn,
		 file_name_extension);
   else snprintf(file_name_buffer, N_file_name_buffer,
		 "%s.t%d.ah%d.it%d.%s",
		 base_file_name,
		 IO_info.time_iteration, hn, AHF_iteration,
		 file_name_extension);

return file_name_buffer;
}
	  }