aboutsummaryrefslogtreecommitdiff
path: root/CarpetAttic/CarpetIOFlexIOCheckpoint/src/ioflexio.cc
diff options
context:
space:
mode:
Diffstat (limited to 'CarpetAttic/CarpetIOFlexIOCheckpoint/src/ioflexio.cc')
-rw-r--r--CarpetAttic/CarpetIOFlexIOCheckpoint/src/ioflexio.cc1134
1 files changed, 1134 insertions, 0 deletions
diff --git a/CarpetAttic/CarpetIOFlexIOCheckpoint/src/ioflexio.cc b/CarpetAttic/CarpetIOFlexIOCheckpoint/src/ioflexio.cc
new file mode 100644
index 000000000..1fba5c365
--- /dev/null
+++ b/CarpetAttic/CarpetIOFlexIOCheckpoint/src/ioflexio.cc
@@ -0,0 +1,1134 @@
+#include <assert.h>
+#include <limits.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/stat.h>
+#include <sys/types.h>
+
+#include <algorithm>
+#include <fstream>
+#include <sstream>
+#include <vector>
+
+
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "AMRwriter.hh"
+#include "AmrGridReader.hh"
+#ifdef HDF4
+# include "HDFIO.hh"
+#endif
+#ifdef HDF5
+# include "H5IO.hh"
+#endif
+#include "IEEEIO.hh"
+#include "IO.hh"
+
+// Hack to stop FlexIO data type clash with LAM MPI
+#undef BYTE
+#undef CHAR
+
+#include "CactusBase/IOUtil/src/ioGH.h"
+
+#include "bbox.hh"
+#include "data.hh"
+#include "gdata.hh"
+#include "ggf.hh"
+#include "vect.hh"
+
+#include "carpet.hh"
+
+#include "ioflexio.hh"
+
+
+extern "C" {
+ static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/CarpetAttic/CarpetIOFlexIOCheckpoint/src/ioflexio.cc,v 1.20 2004/01/13 15:46:52 cott Exp $";
+ CCTK_FILEVERSION(Carpet_CarpetIOFlexIO_ioflexio_cc);
+}
+
+
+
+namespace CarpetIOFlexIO {
+
+ using namespace std;
+ using namespace Carpet;
+ using namespace CarpetIOFlexIOUtil;
+ using namespace CarpetCheckpointRestart;
+
+ // Variable definitions
+ // int GHExtension;
+ int IOMethod;
+ vector<bool> do_truncate;
+ vector<vector<int> > last_output;
+
+
+
+ static const char* GetStringParameter (const char* const parametername,
+ const char* const fallback);
+ static int GetIntParameter (const char* const parametername, int fallback);
+ static bool CheckForVariable (const cGH* const cgh,
+ const char* const varlist, const int vindex);
+ static void SetFlag (int index, const char* optstring, void* arg);
+
+
+
+ int CarpetIOFlexIO_Startup ()
+ {
+ CCTK_RegisterBanner ("AMR 3D FlexIO I/O provided by CarpetIOFlexIO");
+
+ // GHExtension = CCTK_RegisterGHExtension("CarpetIOFlexIO");
+ CCTK_RegisterGHExtensionSetupGH (CCTK_RegisterGHExtension("CarpetIOFlexIO"),SetupGH);
+
+ IOMethod = CCTK_RegisterIOMethod ("CarpetIOFlexIO");
+ CCTK_RegisterIOMethodOutputGH (IOMethod, OutputGH);
+ CCTK_RegisterIOMethodOutputVarAs (IOMethod, OutputVarAs);
+ CCTK_RegisterIOMethodTimeToOutput (IOMethod, TimeToOutput);
+ CCTK_RegisterIOMethodTriggerOutput (IOMethod, TriggerOutput);
+
+ /* register the CarpetIOFlexIO recovery routine to thorn IOUtil */
+ if (IOUtil_RegisterRecover ("CarpetIOFlexIO recovery", CarpetIOFlexIO_Recover) < 0)
+ {
+ CCTK_WARN (1, "Failed to register IOFlexIO recovery routine");
+ }
+
+
+ return 0;
+ }
+
+
+
+ void* SetupGH (tFleshConfig* const fc, const int convLevel, cGH* const cgh)
+ {
+ DECLARE_CCTK_PARAMETERS;
+ CarpetIOFlexIOGH* myGH;
+ CCTK_INT i;
+
+ // 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.
+
+ /* allocate a new GH extension structure */
+
+
+ CCTK_INT numvars = CCTK_NumVars ();
+ myGH = (CarpetIOFlexIOGH*) malloc (sizeof (CarpetIOFlexIOGH));
+ myGH->out_last = (int *) malloc (numvars * sizeof (int));
+ myGH->requests = (ioRequest **) calloc (numvars, sizeof (ioRequest *));
+ myGH->cp_filename_list = (char **) calloc (abs (checkpoint_keep), sizeof (char *));
+ myGH->cp_filename_index = 0;
+ myGH->out_vars = strdup ("");
+ myGH->out_every_default = out_every - 1;
+
+ for (i = 0; i < numvars; i++)
+ {
+ myGH->out_last [i] = -1;
+ }
+
+ myGH->open_output_files = NULL;
+
+
+ return (myGH);
+
+ return 0;
+ }
+
+
+ int OutputGH (const cGH* const cgh) {
+ for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) {
+ if (TimeToOutput(cgh, vindex)) {
+ TriggerOutput(cgh, vindex);
+ }
+ }
+ return 0;
+ }
+
+
+
+ int WriteGF (const cGH* const cgh, IObase* writer, AMRwriter* amrwriter, ioRequest* request, const int called_from_checkpoint)
+ {
+
+ DECLARE_CCTK_PARAMETERS;
+
+ const int varindex = request->vindex;
+
+ const int group = CCTK_GroupIndexFromVarI (varindex);
+ const int n0 = CCTK_FirstVarIndexI(group);
+ assert (n0>=0 && n0<CCTK_NumVars());
+ const int var = varindex - n0;
+ assert (var>=0 && var<CCTK_NumVars());
+ int tl = 0;
+ const int grouptype = CCTK_GroupTypeI(group);
+
+ // let's get the correct Carpet time level (which is the (-1) * Cactus timelevel):
+ if (request->timelevel==0)
+ tl = 0;
+ else
+ tl = - request->timelevel;
+
+ assert (! ( (grouptype != CCTK_GF) && reflevel>0));
+
+ if (CCTK_MyProc(cgh)==0) {
+
+ amrwriter->setType (FlexIODataType(CCTK_VarTypeI(varindex)));
+
+ int gpdim = CCTK_GroupDimI(group);
+
+ // need gpdim=1 if scalar (flexio wants this)
+ if(gpdim == 0)
+ if(grouptype == CCTK_SCALAR)
+ gpdim = 1;
+ else
+ CCTK_WARN(0,"Non-scalar variable with dimension 0!!!");
+
+
+ // Set coordinate information
+ CCTK_REAL lower[dim], upper[dim];
+ double origin[dim], delta[dim], timestep;
+ for (int d=0; d<dim; ++d) {
+ const int ierr = CCTK_CoordRange
+ (cgh, &lower[d], &upper[d], d+1, 0, "cart3d");
+ assert (ierr==0);
+ origin[d] = lower[d];
+ delta[d] = cgh->cctk_delta_space[d];
+ }
+ timestep = cgh->cctk_delta_time;
+ amrwriter->setTopLevelParameters
+ (gpdim, origin, delta, timestep, maxreflevels);
+
+ // Set refinement information
+ int interlevel_timerefinement;
+ int interlevel_spacerefinement[dim];
+ int initial_gridplacementrefinement[dim];
+ interlevel_timerefinement = hh->reffact;
+ for (int d=0; d<dim; ++d) {
+ interlevel_spacerefinement[d] = hh->reffact;
+ initial_gridplacementrefinement[d] = 1;
+ }
+ amrwriter->setRefinement
+ (interlevel_timerefinement, interlevel_spacerefinement,
+ initial_gridplacementrefinement);
+
+ // Set level
+ amrwriter->setLevel (reflevel);
+
+ // Set current time
+ amrwriter->setTime (cgh->cctk_iteration);
+ }
+
+
+
+ // Traverse all components on this refinement and multigrid level
+ BEGIN_COMPONENT_LOOP(cgh, grouptype) {
+
+ const ggf<dim>* ff = 0;
+
+ assert (var < (int)arrdata[group].data.size());
+ ff = (ggf<dim>*)arrdata[group].data[var];
+
+ const gdata<dim>* const data = (*ff) (tl, reflevel, component, mglevel);
+
+ // get some more group information
+ cGroupDynamicData gdyndata;
+
+ int ierr = CCTK_GroupDynamicData(cgh,group,&gdyndata);
+ assert(ierr==0);
+
+ cGroup cgdata;
+ ierr = CCTK_GroupData(group,&cgdata);
+ assert(ierr==0);
+
+ /* handle CCTK_DISTRIB_CONSTANT scalar and arrays */
+#if 0
+ if (cgdata.disttype == CCTK_DISTRIB_CONSTANT) {
+ assert(grouptype == CCTK_ARRAY || grouptype == CCTK_SCALAR);
+ if (hh->processors[reflevel][component] == 0) {
+ if (grouptype == CCTK_SCALAR) {
+ CCTK_VInfo (CCTK_THORNSTRING, "dumping SCALAR distrib const");
+ int rank=1;
+ int dim[1]={1};
+ writer -> write(FlexIODataType(CCTK_VarTypeI(varindex)),rank,dim,CCTK_VarDataPtrI(cgh,tl,varindex));
+ DumpCommonAttributes(cgh,writer,request);
+ continue;
+ }
+ else {
+ writer -> write(FlexIODataType(CCTK_VarTypeI(varindex)),cgdata.dim,gdyndata.lsh,CCTK_VarDataPtrI(cgh,tl,varindex));
+ DumpCommonAttributes(cgh,writer,request);
+ continue;
+ }
+ }
+ else {
+ continue;
+ }
+
+ }
+#endif
+
+ // Make temporary copy on processor 0
+ bbox<int,dim> ext = data->extent();
+ vect<int,dim> lo = ext.lower();
+ vect<int,dim> hi = ext.upper();
+ vect<int,dim> str = ext.stride();
+
+ // Ignore ghost zones if desired
+
+ const int out3D_output_outer_boundary_var = (called_from_checkpoint) ? -1 : out3D_output_outer_boundary;
+ const int out3D_max_num_lower_ghosts_var = (called_from_checkpoint) ? -1 : out3D_max_num_lower_ghosts;
+ const int out3D_max_num_upper_ghosts_var = (called_from_checkpoint) ? -1 : out3D_max_num_upper_ghosts;
+
+
+ for (int d=0; d<dim; ++d) {
+ const int max_lower_ghosts = (gdyndata.bbox[2*d ] && out3D_output_outer_boundary_var) ? -1 : out3D_max_num_lower_ghosts_var;
+ const int max_upper_ghosts = (gdyndata.bbox[2*d+1] && out3D_output_outer_boundary_var) ? -1 : out3D_max_num_upper_ghosts_var;
+
+ const int num_lower_ghosts = max_lower_ghosts == -1 ? gdyndata.nghostzones[d] : min(out3D_max_num_lower_ghosts_var, gdyndata.nghostzones[d]);
+ const int num_upper_ghosts = max_upper_ghosts == -1 ? gdyndata.nghostzones[d] : min(out3D_max_num_upper_ghosts_var, gdyndata.nghostzones[d]);
+
+ lo[d] += (gdyndata.nghostzones[d] - num_lower_ghosts) * str[d];
+ hi[d] -= (gdyndata.nghostzones[d] - num_upper_ghosts) * str[d];
+ }
+
+ ext = bbox<int,dim>(lo,hi,str);
+
+ gdata<dim>* const tmp = data->make_typed (varindex);
+ tmp->allocate (ext, 0);
+ //fprintf(stderr,"\n writing1 %d\n",CCTK_MyProc(cgh));
+ if ( !((cgdata.disttype == CCTK_DISTRIB_CONSTANT) && (hh->processors[reflevel][component]!=0))) {
+
+ if (cgdata.disttype == CCTK_DISTRIB_CONSTANT) {
+ assert(grouptype == CCTK_ARRAY || grouptype == CCTK_SCALAR);
+ //fprintf(stderr,"\n scalar %d %d comp: %d\n",CCTK_MyProc(cgh),varindex,component);
+ int origin[dim], dims[dim];
+ for (int d=0; d<dim; ++d) {
+ origin[d] = (ext.lower() / ext.stride())[d];
+ dims[d] = (ext.shape() / ext.stride())[d];
+ }
+ if (CCTK_MyProc(cgh)==0) {
+ amrwriter->write (origin, dims, (void*)data->storage());
+ DumpCommonAttributes(cgh,writer,request);
+ }
+ delete tmp;
+ continue;
+ } else {
+
+ for (comm_state<dim> state; !state.done(); state.step()) {
+ tmp->copy_from (state, data, ext);
+ }
+
+ //fprintf(stderr,"\n writing2 %d component: %d varindex: %d distrib_const: %d\n",CCTK_MyProc(cgh),component,varindex,(cgdata.disttype == CCTK_DISTRIB_CONSTANT));
+ // Write data
+ if (CCTK_MyProc(cgh)==0) {
+ int origin[dim], dims[dim];
+ for (int d=0; d<dim; ++d) {
+ origin[d] = (ext.lower() / ext.stride())[d];
+ dims[d] = (ext.shape() / ext.stride())[d];
+ }
+
+ amrwriter->write (origin, dims, (void*)tmp->storage());
+
+ // dump attributes
+ DumpCommonAttributes(cgh,writer,request);
+
+ }
+ // Delete temporary copy
+
+ delete tmp;
+
+ }
+ }
+ } END_COMPONENT_LOOP;
+
+
+ return 0;
+ }
+
+
+ int 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;
+ }
+
+ const int grouptype = CCTK_GroupTypeI(group);
+ if (grouptype != CCTK_GF && reflevel>0) return 0;
+
+ int first_vindex = CCTK_FirstVarIndexI (group);
+ /* get the default I/O request for this group */
+ ioRequest* request = IOUtil_DefaultIORequest (cgh, first_vindex, 1);
+
+ // Get grid hierarchy extentsion from IOUtil
+ const ioGH * const iogh = (const ioGH *)CCTK_GHExtension (cgh, "IO");
+ assert (iogh);
+
+ // Create the output directory
+ const char* myoutdir = GetStringParameter("out3D_dir", out_dir);
+ if (CCTK_MyProc(cgh)==0) {
+ CCTK_CreateDirectory (0755, myoutdir);
+ }
+
+ // Invent a file name
+ const char* extension = 0;
+ if (CCTK_Equals(out3D_format, "IEEE")) {
+ extension = ".raw";
+#ifdef HDF4
+ } else if (CCTK_Equals(out3D_format, "HDF4")) {
+ extension = ".hdf";
+#endif
+#ifdef HDF5
+ } else if (CCTK_Equals(out3D_format, "HDF5")) {
+ extension = ".h5";
+#endif
+ } else {
+ assert (0);
+ }
+ extension = GetStringParameter ("out3D_extension", extension);
+
+ ostringstream filenamebuf;
+ filenamebuf << myoutdir << "/" << alias << extension;
+ string filenamestr = filenamebuf.str();
+ const char * const filename = filenamestr.c_str();
+
+ IObase* writer = 0;
+ AMRwriter* amrwriter = 0;
+
+ // Write the file only on the root processor
+ if (CCTK_MyProc(cgh)==0) {
+
+ // If this is the first time, then create and truncate the file
+ if (do_truncate[n]) {
+ struct stat fileinfo;
+ if (! iogh->recovered
+ || stat(filename, &fileinfo)!=0) {
+ writer = 0;
+ if (CCTK_Equals(out3D_format, "IEEE")) {
+ writer = new IEEEIO(filename, IObase::Create);
+#ifdef HDF4
+ } else if (CCTK_Equals(out3D_format, "HDF4")) {
+ writer = new HDFIO(filename, IObase::Create);
+#endif
+#ifdef HDF5
+ } else if (CCTK_Equals(out3D_format, "HDF5")) {
+ writer = new H5IO(filename, IObase::Create);
+#endif
+ } else {
+ assert (0);
+ }
+ delete writer;
+ writer = 0;
+ }
+ }
+
+ // Open the file
+ if (CCTK_Equals(out3D_format, "IEEE")) {
+ writer = new IEEEIO(filename, IObase::Append);
+#ifdef HDF4
+ } else if (CCTK_Equals(out3D_format, "HDF4")) {
+ writer = new HDFIO(filename, IObase::Append);
+#endif
+#ifdef HDF5
+ } else if (CCTK_Equals(out3D_format, "HDF5")) {
+ writer = new H5IO(filename, IObase::Append);
+#endif
+ } else {
+ assert (0);
+ }
+ assert (writer->isValid());
+
+
+ amrwriter = new AMRwriter(*writer);
+ }
+
+
+ WriteGF(cgh,writer,amrwriter,request,0);
+
+ // Close the file
+ if (CCTK_MyProc(cgh)==0) {
+ delete amrwriter;
+ amrwriter = 0;
+ delete writer;
+ writer = 0;
+ }
+
+ // Don't truncate again
+ do_truncate[n] = false;
+
+ return 0;
+ }
+
+
+
+ int TimeToOutput (const cGH* const cgh, const int vindex) {
+ DECLARE_CCTK_PARAMETERS;
+
+ assert (vindex>=0 && vindex<(int)last_output[reflevel].size());
+
+ const int myoutevery = GetIntParameter("out3D_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("out3D_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;
+ }
+
+
+
+ int 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;
+ }
+
+
+
+ int ReadGF (const cGH* const cgh, IObase* reader, AmrGridReader* amrreader,int currdataset) {
+
+ /* this functions reads in a variable on the current reflevel from an already open file. At
+ some point it should be called from InputVarAs */
+
+
+ DECLARE_CCTK_PARAMETERS;
+
+ int tl = -1;
+ int mglevel = -1;
+ int rl = -1;
+ int comp = -1;
+ int myproc = CCTK_MyProc (cgh);
+ int rank;
+ int dims[dim];
+ int nbytes;
+
+ char* varname;
+ char warnstring[256];
+ int asize,i;
+ IObase::DataType datatype;
+ int group,varindex;
+ CCTK_REAL cctk_time;
+
+ if(myproc==0) {
+ // read the name of the variable
+ i = reader->readAttributeInfo ("name", datatype, asize);
+ if (i >= 0 && datatype == FLEXIO_CHAR && asize > 0)
+ {
+ varname = (char*) malloc(sizeof(char)*asize+1);
+ reader->readAttribute (i, varname);
+ }
+ else
+ {
+ CCTK_WARN (0, "Something is wrong! Can't read dataset names!!!");
+ }
+
+ varindex = CCTK_VarIndex(varname);
+
+ assert(varindex > -1);
+
+ group = CCTK_GroupIndexFromVarI(varindex);
+ assert(group > -1);
+
+ // Check for storage
+ if (! CCTK_QueryGroupStorageI(cgh, group)) {
+ CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Cannot recover variable \"%s\" because it has no storage",
+ varname);
+ return 0;
+ }
+
+
+ // read reflevel
+ i = reader->readAttributeInfo ("reflevel", datatype, asize);
+ if (i >= 0 && datatype == FLEXIO_INT && asize > 0)
+ {
+ reader->readAttribute (i, &rl);
+ }
+ else
+ {
+ CCTK_WARN (0, "Something is wrong! Can't read refinement level!!!");
+ }
+
+ i = reader->readAttributeInfo ("component", datatype, asize);
+ if (i >= 0 && datatype == FLEXIO_INT && asize > 0)
+ {
+ reader->readAttribute (i, &comp);
+ }
+ else
+ {
+ CCTK_WARN (0, "Something is wrong! Can't read component!!!");
+ }
+
+ i = reader->readAttributeInfo ("timelevel", datatype, asize);
+ if (i >= 0 && datatype == FLEXIO_INT && asize > 0)
+ {
+ reader->readAttribute (i, &tl);
+ }
+ else
+ {
+ CCTK_WARN (0, "Something is wrong! Can't read timelevel!!!");
+ }
+
+ i = reader->readAttributeInfo ("mglevel", datatype, asize);
+ if (i >= 0 && datatype == FLEXIO_INT && asize > 0)
+ {
+ reader->readAttribute (i, &mglevel);
+ }
+ else
+ {
+ CCTK_WARN (0, "Something is wrong! Can't read multi group level!!!");
+ }
+
+ i = reader->readAttributeInfo ("cctk_time", datatype, asize);
+ if (i >= 0 && datatype == FLEXIO_REAL && asize > 0)
+ {
+ reader->readAttribute (i, &cctk_time);
+ }
+ else
+ {
+ CCTK_WARN (0, "Something is wrong! Can't read coordinate time!!!");
+ }
+
+
+
+ // Read information about dataset
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Reading dataset info");
+ reader->readInfo (datatype, rank, dims);
+ nbytes = IObase::nBytes(datatype,rank,dims);
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "type=%d rank=%d dims=[%d,%d,%d] nbytes=%d", (int)datatype, rank, dims[0], dims[1], dims[2], nbytes);
+
+ int gpdim = CCTK_GroupDimI(group);
+ const int grouptype = CCTK_GroupTypeI(group);
+
+ // need gpdim=1 if scalar (flexio wants this)
+ if(gpdim == 0)
+ if(grouptype == CCTK_SCALAR)
+ gpdim = 1;
+ else
+ CCTK_WARN(0,"Non-scalar variable with dimension 0!!!");
+
+
+ CCTK_VInfo(CCTK_THORNSTRING,"Recovering varindex: %d grouptype: %d varname: %s tl: %d, rl: %d, c: %d",varindex,grouptype,varname,tl,rl,comp);
+
+ free(varname);
+
+ // Check rank
+ assert (rank==gpdim);
+ }
+
+
+ // Broadcast varindex,group,rank, dimensions, and nbytes,rl,tl,mglevel
+ MPI_Bcast (&varindex, 1, MPI_INT, 0, dist::comm);
+ assert (varindex>=0);
+ MPI_Bcast (&group, 1, MPI_INT, 0, dist::comm);
+ assert (group>=0);
+ MPI_Bcast (&rank, 1, MPI_INT, 0, dist::comm);
+ assert (rank>=1);
+ MPI_Bcast (&dims, rank, MPI_INT, 0, dist::comm);
+ for (int d=0; d<rank; ++d) assert (dims[d]>=0);
+ MPI_Bcast (&nbytes, 1, MPI_INT, 0, dist::comm);
+ assert (nbytes>=0);
+
+ MPI_Bcast (&rl, 1, MPI_INT, 0, dist::comm);
+ MPI_Bcast (&tl, 1, MPI_INT, 0, dist::comm);
+ MPI_Bcast (&mglevel, 1, MPI_INT, 0, dist::comm);
+ MPI_Bcast (&comp, 1, MPI_INT, 0, dist::comm);
+
+
+ int gpdim = CCTK_GroupDimI(group);
+ const int grouptype = CCTK_GroupTypeI(group);
+
+ cGroup cgdata;
+ int ierr = CCTK_GroupData(group,&cgdata);
+ assert(ierr==0);
+
+ // Read grid
+ AmrGrid* amrgrid = 0;
+ int amr_origin[dim];
+ int amr_dims[dim];
+ if (myproc==0) {
+
+ // Read data
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Reading AMR data");
+
+ amrgrid = amrreader->getGrid(currdataset);
+ assert (amrgrid!=0);
+ assert (amrgrid->data!=0);
+
+ IObase::DataType atype;
+ int alength;
+ // If iorigin attribute is absent, assume file has unigrid
+ // data. Initialize iorigin to 0.
+ if (reader->readAttributeInfo("iorigin", atype, alength) < 0) {
+ for (int d=0; d<gpdim; ++d) {
+ amrgrid->iorigin[d] = 0;
+ }
+ }
+ for (int d=0; d<gpdim; ++d) {
+ amr_origin[d] = amrgrid->iorigin[d];
+ // fprintf(stderr,"\n amr_origin[%d]=%d",d,amr_origin[d]);
+ amr_dims[d] = amrgrid->dims[d];
+ //fprintf(stderr,"\n amr_dims[%d]=%d",d,amr_dims[d]);
+ }
+ for (int d=gpdim; d<dim; ++d) {
+ amr_origin[d] = 0;
+ amr_dims[d] = 1;
+ }
+
+ } // MyProc == 0
+
+
+ MPI_Bcast (amr_origin, dim, MPI_INT, 0, dist::comm);
+ MPI_Bcast (amr_dims, dim, MPI_INT, 0, dist::comm);
+
+ const int n0 = CCTK_FirstVarIndexI(group);
+ assert (n0>=0 && n0<CCTK_NumVars());
+ const int var = varindex - n0;
+ assert (var>=0 && var<CCTK_NumVars());
+
+
+ // Traverse all components on this refinement and multigrid
+ // level
+
+
+ // fprintf(stderr,"\n bogus! reflevel:%d mglevel:%d\n",reflevel,mglevel);
+ //fprintf(stderr,"blahblah: rank: %d dims[0,1,2]: %d,%d,%d\n",rank,dims[0],dims[1],dims[2]);
+
+ // cout << "var " << varindex << " has " << CCTK_NumTimeLevelsFromVarI(varindex) << " timelevels" << endl;
+
+ // BEGIN_COMPONENT_LOOP(cgh, grouptype) {
+
+ //cout << "compontents " << hh->components(rl) << endl;
+
+ //cout << "myproc: " << CCTK_MyProc(cgh) << endl;
+ // fprintf(stderr,"%d amr_dims: %d,%d,%d\n",CCTK_MyProc(cgh),amr_dims[0],amr_dims[1],amr_dims[2]);
+ //fprintf(stderr,"%d amr_origin: %d,%d,%d\n",CCTK_MyProc(cgh),amr_origin[0],amr_origin[1],amr_origin[2]);
+
+ for(int c=0;c<hh->components(rl);c++) {
+
+ ggf<dim>* ff = 0;
+
+ assert (var < (int)arrdata[group].data.size());
+ ff = (ggf<dim>*)arrdata[group].data[var];
+
+ gdata<dim>* const data = (*ff) (tl, rl, c, mglevel);
+
+
+ // Create temporary data storage on processor 0
+ const int reflevelfact_local=ipow(reffact,rl);
+ vect<int,dim> str = vect<int,dim>(maxreflevelfact/reflevelfact_local);
+
+ if(grouptype == CCTK_SCALAR || grouptype == CCTK_ARRAY)
+ str = vect<int,dim> (1);
+
+
+ vect<int,dim> lb = vect<int,dim>(amr_origin) * str;
+ vect<int,dim> ub = lb + (vect<int,dim>(amr_dims) - 1) * str;
+
+
+
+ gdata<dim>* const tmp = data->make_typed (varindex);
+
+
+ // Copy into grid function
+
+ if (cgdata.disttype == CCTK_DISTRIB_CONSTANT) {
+ assert(grouptype == CCTK_ARRAY || grouptype == CCTK_SCALAR);
+ if (grouptype == CCTK_SCALAR) {
+ lb[0] = hh->processors.at(rl).at(c);
+ ub[0] = hh->processors.at(rl).at(c);
+ } else {
+ lb[dim-1] = lb[dim-1] + (ub[dim-1]-lb[dim-1]+1)*hh->processors.at(rl).at(c);
+ ub[dim-1] = ub[dim-1] + (ub[dim-1]-lb[dim-1]+1)*hh->processors.at(rl).at(c);
+ }
+ }
+
+ const bbox<int,dim> ext(lb,ub,str);
+
+
+ if (myproc==0) {
+ tmp->allocate (ext, 0, amrgrid->data);
+ } else {
+ tmp->allocate (ext, 0, 0);
+ }
+
+ for (comm_state<dim> state; !state.done(); state.step()) {
+ data->copy_from (state, tmp, ext & data->extent() );
+ }
+
+
+ // Delete temporary copy
+ delete tmp;
+
+ } // manual component loop
+
+ if (myproc==0) {
+ free (amrgrid->data);
+ free (amrgrid);
+ amrgrid = 0;
+ }
+
+
+ return 0;
+ }
+
+
+
+ int InputGH (const cGH* const cgh) {
+ int retval = 0;
+ for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) {
+ if (CheckForVariable(cgh, GetStringParameter("in3D_vars",""), vindex)) {
+ char* varname = CCTK_FullName(vindex);
+ retval = InputVarAs (cgh, varname, CCTK_VarName(vindex));
+ free (varname);
+ if (retval != 0) return retval;
+ }
+ }
+ return retval;
+ }
+
+
+
+ int InputVarAs (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; // CCTK_VInfo (CCTK_THORNSTRING, "boguscheck reflevel,component,mglevel %d,%d,%d",reflevel,component,mglevel);
+
+
+ // Check for storage
+ if (! CCTK_QueryGroupStorageI(cgh, group)) {
+ CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Cannot input variable \"%s\" because it has no storage",
+ varname);
+ return 0;
+ }
+
+ const int grouptype = CCTK_GroupTypeI(group);
+ const int rl = grouptype==CCTK_GF ? reflevel : 0;
+
+ // Find the input directory
+ const char* myindir = GetStringParameter("indir3D", "");
+
+ // Invent a file name
+ const char* extension = 0;
+ if (CCTK_Equals(in3D_format, "IEEE")) {
+ extension = ".raw";
+#ifdef HDF4
+ } else if (CCTK_Equals(in3D_format, "HDF4")) {
+ extension = ".hdf";
+#endif
+#ifdef HDF5
+ } else if (CCTK_Equals(in3D_format, "HDF5")) {
+ extension = ".h5";
+#endif
+ } else {
+ assert (0);
+ }
+ extension = GetStringParameter ("in3D_extension", extension);
+
+ ostringstream filenamebuf;
+ filenamebuf << myindir << "/" << alias << extension;
+ string filenamestr = filenamebuf.str();
+ const char * const filename = filenamestr.c_str();
+
+ IObase* reader = 0;
+ AmrGridReader* amrreader = 0;
+ int ndatasets = -1;
+
+ const int gpdim = CCTK_GroupDimI(group);
+
+ int rank;
+ int dims[dim];
+ int nbytes;
+
+ // Read the file only on the root processor
+ if (CCTK_MyProc(cgh)==0) {
+
+ // Open the file
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Opening file \"%s\"", filename);
+ if (CCTK_Equals(in3D_format, "IEEE")) {
+ reader = new IEEEIO(filename, IObase::Read);
+#ifdef HDF4
+ } else if (CCTK_Equals(in3D_format, "HDF4")) {
+ reader = new HDFIO(filename, IObase::Read);
+#endif
+#ifdef HDF5
+ } else if (CCTK_Equals(in3D_format, "HDF5")) {
+ reader = new H5IO(filename, IObase::Read);
+#endif
+ } else {
+ assert (0);
+ }
+ if (!reader->isValid()) {
+ CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+ "Could not open file \"%s\" for reading", filename);
+ }
+ assert (reader->isValid());
+
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Reading AMR info");
+ amrreader = new AmrGridReader(*reader);
+
+ // Read information about dataset
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Reading dataset info");
+ IObase::DataType numbertype;
+ reader->readInfo (numbertype, rank, dims);
+ nbytes = IObase::nBytes(numbertype,rank,dims);
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "type=%d rank=%d dims=[%d,%d,%d] nbytes=%d", (int)numbertype, rank, dims[0], dims[1], dims[2], nbytes);
+
+ // Check rank
+ assert (rank==gpdim);
+
+ // Check datatype
+ // TODO: Check datatype correctly
+ assert (CCTK_VarTypeI(n) == CCTK_VARIABLE_REAL8
+ || (sizeof(CCTK_REAL) == sizeof(CCTK_REAL8)
+ && CCTK_VarTypeI(n) == CCTK_VARIABLE_REAL));
+
+ // TODO: check grid spacing
+
+ // Number of datasets
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Reading number of datasets");
+ ndatasets = reader->nDatasets();
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "ndatasets=%d", ndatasets);
+ assert (ndatasets>=0);
+ }
+
+ // Broadcast rank, dimensions, and nbytes
+ MPI_Bcast (&rank, 1, MPI_INT, 0, dist::comm);
+ assert (rank>=1);
+ MPI_Bcast (&dims, rank, MPI_INT, 0, dist::comm);
+ for (int d=0; d<rank; ++d) assert (dims[d]>=0);
+ MPI_Bcast (&nbytes, 1, MPI_INT, 0, dist::comm);
+ assert (nbytes>=0);
+
+ // Broadcast number of datasets
+ MPI_Bcast (&ndatasets, 1, MPI_INT, 0, dist::comm);
+ assert (ndatasets>=0);
+
+ // Read all datasets
+ // TODO: read only some datasets
+ for (int dataset=0; dataset<ndatasets; ++dataset) {
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Handling dataset #%d", dataset);
+
+ // Read grid
+ AmrGrid* amrgrid = 0;
+ int amr_origin[dim];
+ int amr_dims[dim];
+
+ if (CCTK_MyProc(cgh)==0) {
+
+ // Read data
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Reading AMR data");
+ amrgrid = amrreader->getGrid(dataset);
+ assert (amrgrid!=0);
+ assert (amrgrid->data!=0);
+
+ // If iorigin attribute is absent, assume file has unigrid
+ // data. Initialize iorigin to 0.
+ IObase::DataType atype;
+ int alength;
+ if (reader->readAttributeInfo("iorigin", atype, alength) < 0) {
+ for (int d=0; d<gpdim; ++d) {
+ amrgrid->iorigin[d] = 0;
+ }
+ }
+
+ for (int d=0; d<gpdim; ++d) {
+ amr_origin[d] = amrgrid->iorigin[d];
+ amr_dims[d] = amrgrid->dims[d];
+ }
+ for (int d=gpdim; d<dim; ++d) {
+ amr_origin[d] = 0;
+ amr_dims[d] = 1;
+ }
+
+ } // MyProc == 0
+ MPI_Bcast (amr_origin, dim, MPI_INT, 0, dist::comm);
+ MPI_Bcast (amr_dims, dim, MPI_INT, 0, dist::comm);
+
+ // Traverse all components on this refinement and multigrid
+ // level
+ BEGIN_COMPONENT_LOOP(cgh, grouptype) {
+
+ ggf<dim>* ff = 0;
+
+ assert (var < (int)arrdata[group].data.size());
+ ff = (ggf<dim>*)arrdata[group].data[var];
+
+ gdata<dim>* const data = (*ff) (tl, rl, component, mglevel);
+
+ // Create temporary data storage on processor 0
+ const vect<int,dim> str = vect<int,dim>(reflevelfact);
+ const vect<int,dim> lb = vect<int,dim>(amr_origin) * str;
+ const vect<int,dim> ub
+ = lb + (vect<int,dim>(amr_dims) - 1) * str;
+ const bbox<int,dim> ext(lb,ub,str);
+ gdata<dim>* const tmp = data->make_typed (n);
+
+ if (CCTK_MyProc(cgh)==0) {
+ tmp->allocate (ext, 0, amrgrid->data);
+ } else {
+ tmp->allocate (ext, 0, 0);
+ }
+
+ // Copy into grid function
+ for (comm_state<dim> state; !state.done(); state.step()) {
+ data->copy_from (state, tmp, ext);
+ }
+
+ // Delete temporary copy
+ delete tmp;
+
+ } END_COMPONENT_LOOP;
+
+ if (CCTK_MyProc(cgh)==0) {
+ free (amrgrid->data);
+ free (amrgrid);
+ amrgrid = 0;
+ }
+
+ } // loop over datasets
+
+ // Close the file
+ if (CCTK_MyProc(cgh)==0) {
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Deleting AMR info");
+ delete amrreader;
+ amrreader = 0;
+ if (verbose) CCTK_VInfo (CCTK_THORNSTRING, "Closing file");
+ delete reader;
+ reader = 0;
+ }
+
+ return 0;
+ }
+
+
+
+ int CarpetIOFlexIO_ReadData (CCTK_ARGUMENTS)
+ {
+ DECLARE_CCTK_ARGUMENTS;
+ return InputGH(cctkGH);
+ }
+
+
+
+ const char* GetStringParameter (const char* const parametername,
+ const char* const fallback)
+ {
+ 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;
+ }
+
+
+
+ int GetIntParameter (const char* const parametername, int fallback)
+ {
+ 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);
+
+ vector<bool> flags(numvars);
+
+ CCTK_TraverseString (varlist, SetFlag, &flags, CCTK_GROUP_OR_VAR);
+
+ return flags[vindex];
+ }
+
+ void SetFlag (int index, const char* optstring, void* arg)
+ {
+ vector<bool>& flags = *(vector<bool>*)arg;
+ flags[index] = true;
+ }
+
+
+
+} // namespace CarpetIOFlexIO
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-12 19:53:29 +0000