diff options
Diffstat (limited to 'CarpetExtra/CarpetIOSer/src/ioser.cc')
-rw-r--r-- | CarpetExtra/CarpetIOSer/src/ioser.cc | 781 |
1 files changed, 781 insertions, 0 deletions
diff --git a/CarpetExtra/CarpetIOSer/src/ioser.cc b/CarpetExtra/CarpetIOSer/src/ioser.cc new file mode 100644 index 000000000..431b6a295 --- /dev/null +++ b/CarpetExtra/CarpetIOSer/src/ioser.cc @@ -0,0 +1,781 @@ +#include <alloca.h> +#include <assert.h> +#include <limits.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <sys/stat.h> +#include <sys/types.h> + +#include <fstream> +#include <iomanip> +#include <vector> + +#include "cctk.h" +#include "cctk_Parameters.h" + +#include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h" + +#include "Carpet/CarpetLib/src/data.hh" +#include "Carpet/CarpetLib/src/gdata.hh" +#include "Carpet/CarpetLib/src/gf.hh" +#include "Carpet/CarpetLib/src/ggf.hh" +#include "Carpet/CarpetLib/src/vect.hh" + +#include "carpet.hh" + +#include "ioser.hh" + +#include "bbhutil.h" +extern "C" { +int vsxynt(const char *name, double t, double *x, double *y, int n); +} + + +static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/CarpetIOSer/src/ioser.cc,v 1.1.1.1 2002/02/08 14:55:31 shawley Exp $"; + +// That's a hack +namespace Carpet { + void UnsupportedVarType (const int vindex); +} + + +using namespace std; +using namespace Carpet; + + + +static bool CheckForVariable (const cGH* const cgh, + const char* const varlist, const int vindex); +static void SetFlag (int index, const char* optstring, void* arg); + +//void assign_ser_output_element (double& range, const ivect& index); +//void assign_ser_output_element (double& range, const vect<int,dim> index); + +template<int D,int DD> +static void write_ser ( const char *gfname, + const generic_data<D>* const gfdata, + const bbox<int,D>& gfext, + const int vi, + const int time, + const vect<int,D>& org, + const vect<int,DD>& dirs, + const int tl, + const int rl, + const int c, + const int ml, + const CCTK_REAL coord_time, + const vect<CCTK_REAL,D>& coord_lower, + const vect<CCTK_REAL,D>& coord_upper, + const int strip_left, const int strip_right); + + +int CarpetIOSerStartup() +{ + CarpetIOSer<1>::Startup(); + CarpetIOSer<2>::Startup(); + CarpetIOSer<3>::Startup(); + + return 0; +} + + + +// Definition of static members +template<int outdim> int CarpetIOSer<outdim>::GHExtension; +template<int outdim> int CarpetIOSer<outdim>::IOMethod; +template<int outdim> vector<bool> CarpetIOSer<outdim>::do_truncate; +template<int outdim> vector<vector<int> > CarpetIOSer<outdim>::last_output; + + + +template<int outdim> +int CarpetIOSer<outdim>::Startup() +{ + char msg[1000]; + sprintf (msg, "AMR %dD Ser I/O provided by CarpetIOSer", outdim); + CCTK_RegisterBanner (msg); + + char extension_name[1000]; + sprintf (extension_name, "CarpetIOSer_%dD", outdim); + GHExtension = CCTK_RegisterGHExtension(extension_name); + CCTK_RegisterGHExtensionSetupGH (GHExtension, SetupGH); + + char method_name[1000]; + sprintf (method_name, "IOSer_%dD", outdim); + IOMethod = CCTK_RegisterIOMethod (method_name); + CCTK_RegisterIOMethodOutputGH (IOMethod, OutputGH); + CCTK_RegisterIOMethodOutputVarAs (IOMethod, OutputVarAs); + CCTK_RegisterIOMethodTimeToOutput (IOMethod, TimeToOutput); + CCTK_RegisterIOMethodTriggerOutput (IOMethod, TriggerOutput); + + return 0; +} + + + +template<int outdim> +void* CarpetIOSer<outdim> +::SetupGH (tFleshConfig* const fc, const int convLevel, cGH* const cgh) +{ + DECLARE_CCTK_PARAMETERS; + + // Truncate all files if this is not a restart + do_truncate.resize(CCTK_NumVars(), true); + + // No iterations have yet been output + last_output.resize(maxreflevels); + for (int rl=0; rl<maxreflevels; ++rl) { + last_output[rl].resize(CCTK_NumVars(), INT_MIN); + } + + // We register only once, ergo we get only one handle. We store + // that statically, so there is no need to pass anything to Cactus. + return 0; +} + + + +template<int outdim> +int CarpetIOSer<outdim> +::OutputGH (const cGH* const cgh) +{ + for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) { + if (TimeToOutput(cgh, vindex)) { + TriggerOutput(cgh, vindex); + } + } + return 0; +} + + + +template<int outdim> +int CarpetIOSer<outdim> +::OutputVarAs (const cGH* const cgh, + const char* const varname, const char* const alias) +{ + DECLARE_CCTK_PARAMETERS; + + const int n = CCTK_VarIndex(varname); + assert (n>=0 && n<CCTK_NumVars()); + const int group = CCTK_GroupIndexFromVarI (n); + assert (group>=0 && group<(int)Carpet::arrdata.size()); + const int n0 = CCTK_FirstVarIndexI(group); + assert (n0>=0 && n0<CCTK_NumVars()); + const int var = n - n0; + assert (var>=0 && var<CCTK_NumVars()); + const int tl = 0; + + // Check for storage + if (! CCTK_QueryGroupStorageI(cgh, group)) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "Cannot output variable \"%s\" because it has no storage", + varname); + return 0; + } + + // Create the output directory + const char* myoutdir = GetStringParameter("outdir%dD", outdir); + if (CCTK_MyProc(cgh)==0) { + CCTK_CreateDirectory (0755, myoutdir); + } + + // Loop over all direction combinations + vect<int,outdim> dirs; + for (int d=0; d<outdim; ++d) dirs[d] = 0; + + bool done; + do { + + // Output each combination only once + bool ascending = true; + for (int d1=0; d1<outdim; ++d1) { + for (int d2=d1+1; d2<outdim; ++d2) { + ascending = ascending && dirs[d1] < dirs[d2]; + } + } + + // Skip output if the dimensions are not ascending + if (ascending) { + + // If this output is desired at all + bool desired; + switch (outdim) { + case 1: + switch (dirs[0]) { + case 0: + desired = out1D_x; + break; + case 1: + desired = out1D_y; + break; + case 2: + desired = out1D_z; + break; + default: + abort(); + } + break; + default: + abort(); + } + + // Skip output if not desired + if (desired) { + + // Invent a file name + char* const filename + = (char*)alloca(strlen(myoutdir)+strlen(alias)+100); + sprintf (filename, "%s/%s.", myoutdir, alias); + for (int d=0; d<outdim; ++d) { + assert (dirs[d]>=0 && dirs[d]<3); + const char* const coords = "xyz"; + sprintf (filename, "%s%c", filename, coords[dirs[d]]); + } + const char* const suffixes = "lpv"; + sprintf (filename, "%s%c", filename, suffixes[outdim-1]); + + ofstream file; + + if (CCTK_MyProc(cgh)==0) { + // If this is the first time, then write a nice header on + // the root processor + if (do_truncate[n]) { + struct stat fileinfo; + if (! IOUtil_RestartFromRecovery(cgh) + || stat(filename, &fileinfo)!=0) { + file.open (filename, ios::out | ios::trunc); + assert (file.good()); + file << "# " << varname; + for (int d=0; d<outdim; ++d) { + file << " " << "xyz"[dirs[d]]; + } + file << " (" << alias << ")" << endl; + file << "#" << endl; + assert (file.good()); + } + } + if (! file.is_open()) { + file.open (filename, ios::app); + assert (file.good()); + } + file << setprecision(15); + assert (file.good()); + } + + assert (outdim <= CCTK_GroupDimI(group)); + + // Find the output offset + vect<int,dim> offset(0); + switch (outdim) { + case 1: + switch (dirs[0]) { + case 0: + offset[1] = GetGridOffset (cgh, 2, + "out%dD_xline_yi", "out_xline_yi", + "out%dD_xline_y", "out_xline_y", + out_xline_y); + offset[2] = GetGridOffset (cgh, 3, + "out%dD_xline_zi", "out_xline_zi", + "out%dD_xline_z", "out_xline_z", + out_xline_z); + break; + case 1: + offset[0] = GetGridOffset (cgh, 1, + "out%dD_yline_xi", "out_yline_xi", + "out%dD_yline_x", "out_yline_x", + out_yline_x); + offset[2] = GetGridOffset (cgh, 3, + "out%dD_yline_zi", "out_yline_zi", + "out%dD_yline_z", "out_yline_z", + out_yline_z); + break; + case 2: + offset[0] = GetGridOffset (cgh, 1, + "out%dD_zline_xi", "out_zline_xi", + "out%dD_zline_x", "out_zline_x", + out_zline_x); + offset[1] = GetGridOffset (cgh, 2, + "out%dD_zline_yi", "out_zline_yi", + "out%dD_zline_y", "out_zline_y", + out_zline_y); + break; + default: + abort(); + } + break; + default: + abort(); + } + + // Traverse all components on this refinement and multigrid + // level + BEGIN_COMPONENT_LOOP(cgh) { + + const generic_gf<dim>* ff = 0; + + assert (var < (int)arrdata[group].data.size()); + ff = (generic_gf<dim>*)arrdata[group].data[var]; + + const generic_data<dim>* const data + = (*ff) (tl, reflevel, component, mglevel); + const bbox<int,dim> ext = data->extent(); + const vect<int,dim> offset1 = offset * ext.stride(); + + // coordinates + const double coord_time = cgh->cctk_time; + vect<double,dim> global_lower, global_upper; + for (int d=0; d<dim; ++d) { + const int ierr = CCTK_CoordRange + (cgh, &global_lower[d], &global_upper[d], d+1, 0, "cart3d"); + if (ierr<0) { + global_lower[d] = 0; + global_upper[d] = 0; + } + } + const vect<double,dim> global_extent (hh->baseextent.upper() - hh->baseextent.lower() + hh->baseextent.stride() * (dd->lghosts + dd->ughosts)); + vect<double,dim> coord_delta; + for (int d=0; d<dim; ++d) { + if (global_extent[d] != 0) { + coord_delta[d] = (global_upper[d] - global_lower[d]) / global_extent[d]; + } else { + coord_delta[d] = 0; + } + } + // Note: don't permute the "coord_delta" and "data->extent().lower()" + // (you'll pick up the integer operator* then) + const vect<double,dim> coord_lower = global_lower + coord_delta * vect<double,dim>(data->extent().lower()); + const vect<double,dim> coord_upper = global_lower + coord_delta * vect<double,dim>(data->extent().upper()); + + char *gfname = new char[strlen(alias)+strlen(out1D_tag)+1]; + sprintf(gfname,"%s%s",alias,out1D_tag); + write_ser (gfname, data, ext, n, cgh->cctk_iteration, offset1, dirs, + tl, reflevel, component, mglevel, + coord_time, coord_lower, coord_upper, + strip_left, strip_right); + delete [] gfname; + + // Append EOL after every component + if (CCTK_MyProc(cgh)==0) { + if (separate_components) { + assert (file.good()); + file << endl; + } + } + assert (file.good()); + + } END_COMPONENT_LOOP(cgh); + + // Append EOL after every complete set of components + if (CCTK_MyProc(cgh)==0) { + if (separate_grids) { + assert (file.good()); + file << endl; + } + file.close(); + assert (file.good()); + } + + assert (! file.is_open()); + + } // if (desired) + + } // if (ascending) + + // Next direction combination + done = true; + for (int d=0; d<outdim; ++d) { + ++dirs[d]; + if (dirs[d]<CCTK_GroupDimI(group)) { + done = false; + break; + } + dirs[d] = 0; + } + + } while (! done); // all directions + + // Don't truncate again + do_truncate[n] = false; + + return 0; +} + + + +template<int outdim> +int CarpetIOSer<outdim> +::TimeToOutput (const cGH* const cgh, const int vindex) +{ + DECLARE_CCTK_PARAMETERS; + + assert (vindex>=0 && vindex<CCTK_NumVars()); + + const int myoutevery = GetIntParameter("out%dD_every", out_every); + + if (myoutevery < 0) { + // Nothing should be output at all + return 0; + } + + if (cgh->cctk_iteration % myoutevery != 0) { + // Nothing should be output during this iteration + return 0; + } + + if (! CheckForVariable(cgh, GetStringParameter("out%dD_vars", ""), vindex)) { + // This variable should not be output + return 0; + } + + if (last_output[reflevel][vindex] == cgh->cctk_iteration) { + // Has already been output during this iteration + char* varname = CCTK_FullName(vindex); + CCTK_VWarn (5, __LINE__, __FILE__, CCTK_THORNSTRING, + "Skipping output for variable \"%s\", because this variable " + "has already been output during the current iteration -- " + "probably via a trigger during the analysis stage", + varname); + free (varname); + return 0; + } + + assert (last_output[reflevel][vindex] < cgh->cctk_iteration); + + // Should be output during this iteration + return 1; +} + + + +template<int outdim> +int CarpetIOSer<outdim> +::TriggerOutput (const cGH* const cgh, const int vindex) +{ + assert (vindex>=0 && vindex<CCTK_NumVars()); + + char* varname = CCTK_FullName(vindex); + const int retval = OutputVarAs (cgh, varname, CCTK_VarName(vindex)); + free (varname); + + last_output[reflevel][vindex] = cgh->cctk_iteration; + + return retval; +} + + + +template<int outdim> +int CarpetIOSer<outdim> +::GetGridOffset (const cGH* const cgh, const int dir, + const char* const itempl, const char* const iglobal, + const char* const ctempl, const char* const cglobal, + const CCTK_REAL cfallback) +{ + // First choice: explicit coordinate + char cparam[1000]; + sprintf (cparam, ctempl, outdim); + if (CCTK_ParameterQueryTimesSet (cparam, CCTK_THORNSTRING) > 0) { + int ptype; + const CCTK_REAL* const pcoord + = (const CCTK_REAL*)CCTK_ParameterGet (cparam, CCTK_THORNSTRING, &ptype); + assert (pcoord); + const CCTK_REAL coord = *pcoord; + assert (ptype == PARAMETER_REAL); + return CoordToOffset (cgh, dir, coord, 0); + } + + // Second choice: explicit index + char iparam[1000]; + sprintf (iparam, itempl, outdim); + if (CCTK_ParameterQueryTimesSet (iparam, CCTK_THORNSTRING) > 0) { + int ptype; + const int* const pindex + = (const int*)CCTK_ParameterGet (iparam, CCTK_THORNSTRING, &ptype); + assert (pindex); + const int index = *pindex; + assert (ptype == PARAMETER_INT); + return index; + } + + // Third choice: explicit global coordinate + if (CCTK_ParameterQueryTimesSet (cglobal, "IO") > 0) { + int ptype; + const CCTK_REAL* const pcoord + = (const CCTK_REAL*)CCTK_ParameterGet (cglobal, "IO", &ptype); + assert (pcoord); + const CCTK_REAL coord = *pcoord; + assert (ptype == PARAMETER_REAL); + return CoordToOffset (cgh, dir, coord, 0); + } + + // Fourth choice: explicit global index + if (CCTK_ParameterQueryTimesSet (iglobal, "IO") > 0) { + int ptype; + const int* const pindex + = (const int*)CCTK_ParameterGet (iglobal, "IO", &ptype); + assert (pindex); + const int index = *pindex; + assert (ptype == PARAMETER_INT); + return index; + } + + // Fifth choice: default coordinate + return CoordToOffset (cgh, dir, cfallback, 0); +} + + + +template<int outdim> +int CarpetIOSer<outdim> +::CoordToOffset (const cGH* cgh, const int dir, const CCTK_REAL coord, + const int ifallback) +{ + assert (dir>=1 && dir<=dim); + + CCTK_REAL lower, upper; + CCTK_CoordRange (cgh, &lower, &upper, dir, 0, "cart3d"); + + const int npoints = cgh->cctk_gsh[dir-1]; + + const CCTK_REAL rindex = (coord - lower) / (upper - lower) * (npoints-1); + const int levfac = cgh->cctk_levfac[dir-1]; + int cindex = (int)floor(rindex / levfac + 0.5 + 1e-6) * levfac; + + if (cindex<0 || cindex>=npoints) { + CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING, + "The specified coordinate value %g is not within the grid range [%g,%g]", + coord, lower, upper); + cindex = ifallback; + } + + assert (cindex>=0 && cindex<npoints); + + return cindex; +} + + + +template<int outdim> +const char* CarpetIOSer<outdim> +::GetStringParameter (const char* const parametertemplate, + const char* const fallback) +{ + char parametername[1000]; + sprintf (parametername, parametertemplate, outdim); + if (CCTK_ParameterQueryTimesSet (parametername, CCTK_THORNSTRING) > 0) { + int ptype; + const char* const* const ppval + = (const char* const*)CCTK_ParameterGet (parametername, CCTK_THORNSTRING, + &ptype); + assert (ppval); + const char* const pval = *ppval; + assert (ptype == PARAMETER_STRING); + return pval; + } + + return fallback; +} + + + +template<int outdim> +int CarpetIOSer<outdim> +::GetIntParameter (const char* const parametertemplate, const int fallback) +{ + char parametername[1000]; + sprintf (parametername, parametertemplate, outdim); + if (CCTK_ParameterQueryTimesSet (parametername, CCTK_THORNSTRING) > 0) { + int ptype; + const int* const ppval + = (const int*)CCTK_ParameterGet (parametername, CCTK_THORNSTRING, + &ptype); + assert (ppval); + const int pval = *ppval; + assert (ptype == PARAMETER_INT); + return pval; + } + + return fallback; +} + + + +bool CheckForVariable (const cGH* const cgh, + const char* const varlist, const int vindex) +{ + const int numvars = CCTK_NumVars(); + assert (vindex>=0 && vindex<numvars); + + bool* const flags = (bool*)alloca(numvars * sizeof(bool)); + + for (int i=0; i<numvars; ++i) { + flags[i] = false; + } + + CCTK_TraverseString (varlist, SetFlag, flags, CCTK_GROUP_OR_VAR); + + return flags[vindex]; +} + +void SetFlag (int index, const char* optstring, void* arg) +{ + ((bool*)arg)[index] = true; +} + + +// Output +/* +template<class T, int D> +void data<T,D>::assign_ser_output_element (double& range, const vect<int,dim> index) +{ + range = (*this)[index]; +} +*/ + + +// Output +template<int D,int DD> +void write_ser ( const char *gfname, + const generic_data<D>* const gfdata, + const bbox<int,D>& gfext, + const int vi, + const int time, + const vect<int,D>& org, + const vect<int,DD>& dirs, + const int tl, + const int rl, + const int c, + const int ml, + const CCTK_REAL coord_time, + const vect<CCTK_REAL,D>& coord_lower, + const vect<CCTK_REAL,D>& coord_upper, + const int strip_left, const int strip_right) + +{ + assert (DD<=D); + + if (gfdata->proc()==0) { + // output on processor 0 + + int rank; + MPI_Comm_rank (dist::comm, &rank); + if (rank == 0) { + + + assert (D>=1 && D<=3); + assert (DD==1); + const char* const coords = "xyz"; + + const vect<int,DD> lo = gfext.lower()[dirs]; + const vect<int,DD> up = gfext.upper()[dirs]; + const vect<int,DD> str = gfext.stride()[dirs]; + const bbox<int,DD> ext(lo,up,str); + + // Check whether the output origin is contained in the extent of + // the data + vect<int,dim> org1(org); + for (int d=0; d<DD; ++d) org1[dirs[d]] = ext.lower()[d]; + if (gfext.contains(org1)) { + + const int nelems = ext.num_points(); + double *domain = new double[nelems]; + double *range = new double[nelems]; + + int elem = 0; + for (bbox<int,DD>::iterator it=ext.begin(); it!=ext.end(); ++it) { + vect<int,dim> index(org); + for (int d=0; d<DD; ++d) index[dirs[d]] = (*it)[d]; + const int d = dirs[0]; + if (gfext.upper()[d] - gfext.lower()[d] != 0) { + domain[elem] = coord_lower[d] + (index[d] - gfext.lower()[d]) * +(coord_upper[d] - coord_lower[d]) / (gfext.upper()[d] - gfext.lower()[d]); } else { + domain[elem] = coord_lower[d]; + } + +// assign_ser_output_element (range[elem], index); +// range[elem] = (*(data<T,D>*)gfdata)[index]; +// elem++; + +// cout << "CarpetIOSer, before switch" << endl << flush; + + switch (CCTK_VarTypeI(vi)) { +#define TYPECASE(N,T) \ + case N: \ + range[elem] = (*(data<T,D>*)gfdata)[index]; \ + elem++;\ + break; +#include "Carpet/Carpet/src/typecase" +#undef TYPECASE + default: + UnsupportedVarType(vi); + } + // cout << "CarpetIOSer, after switch" << endl << flush; + + } + assert(elem == nelems); + + char *window_name = new char[strlen(gfname)+100]; + sprintf(window_name,"%s_l%d.%c",gfname,rl,"xyz"[dirs[0]]); + vsxynt(window_name, coord_time, domain+strip_left, range+strip_left, + nelems-strip_left-strip_right); + + delete [] window_name; + delete [] domain; + delete [] range; + + } // if extent().contains(org1) + + } + + } else { + // copy to processor 0 and output there + +/* + generic_data* const tmp = make_typed(); + tmp->allocate(extent(), 0); + tmp->copy_from (this, extent()); + tmp->write_ser (gfname, time, org, dirs, tl, rl, c, ml, + coord_time, coord_lower, coord_upper); + delete tmp; +*/ + + generic_data<D>* const tmp = gfdata->make_typed(); + tmp->allocate(gfdata->extent(), 0); + tmp->copy_from (gfdata, gfdata->extent()); + write_ser (gfname, tmp, gfext, vi, time, org, dirs, tl, rl, c, ml, + coord_time, coord_lower, coord_upper, strip_left, strip_right); + delete tmp; + + + + } +} + + + + + +// Explicit instantiation for all output dimensions +template class CarpetIOSer<1>; +template class CarpetIOSer<2>; +template class CarpetIOSer<3>; + +template +void write_ser ( const char *gfname, + const generic_data<3>* const gfdata, + const bbox<int,3>& gfext, + const int vi, + const int time, + const vect<int,3>& org, + const vect<int,1>& dirs, + const int tl, + const int rl, + const int c, + const int ml, + const CCTK_REAL coord_time, + const vect<CCTK_REAL,3>& coord_lower, + const vect<CCTK_REAL,3>& coord_upper, + const int strip_left, const int strip_right); + |