EJWFGXTWIBKJ62LO6QM62LM7MGTYEHOCZDKA3UVR5UYWQZ2O54DQC 6D6MZBKTLXUK4ZCIFUPNV7XT64XSFFG4CL6HN6ONK5ONM6MR64AAC GKKJ75HX2ERLVBZVE2CUB6T3J2SUT7R3UKEKTEYNOG43ZKX6X5MQC BZ7T742UYYUIK2X3Y6UWJ2G33ER4PJJBA7YLCVAQKAJQUUECHVRQC RTNZAS3UPI6GG3KY4Z5WVXJ4R2YF5427BB6WAV3GHRS5W7XPOSUQC 722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC RZPXXQUUAG7JIADVCLYQOY3I4CQDR23PF57L6XDGHNCT4EBPSO6QC NUOLOGCKMF5UOBGBYEOX4O7NQ5AEVVLCH6KRBQRJQXIRDNJ2C2ZQC 5XGIB7XMEZNBLA5ZLQTXRTC3AZ5CTRGMXWBPVMWXG4DPHKWDF4ZAC RCLGQ2LZMFVPBPTU2G55DJ6HZPOGGTPZRZCY54VGP6YLHANJ2LAQC 2MTD37PDBF7KZKC6LXG2IVDYMMUMNNS557WMFVSCZM3SAKOGDC5QC FEMASUBNU32NSG4DNXZX54CGCA57PVRGYO46L3A6F2EJ4BCSJ3SAC BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC MSBBCXVGD3GRLE5KAI6BKAFRV7SQUWI2SNN43AJAUD3ISRCEXY6QC ActiveThorns = "ADMBaseBrillLindquistCarpetXFormalineIOUtilODESolversZ4c"CarpetX::verbose = yesCactus::cctk_show_schedule = yesCactus::terminate = "time"Cactus::cctk_final_time = 0.0CarpetX::xmin = -1.0CarpetX::ymin = -1.0CarpetX::zmin = -1.0CarpetX::xmax = +1.0CarpetX::ymax = +1.0CarpetX::zmax = +1.0CarpetX::ncells_x = 64CarpetX::ncells_y = 64CarpetX::ncells_z = 64CarpetX::periodic_x = noCarpetX::periodic_y = noCarpetX::periodic_z = noCarpetX::periodic = no# CarpetX::reflection_x = yes# CarpetX::reflection_y = yes# CarpetX::reflection_z = yes# CarpetX::reflection_upper_x = yes# CarpetX::reflection_upper_y = yes# CarpetX::reflection_upper_z = yesCarpetX::ghost_size = 3CarpetX::max_num_levels = 1CarpetX::regrid_every = 0CarpetX::prolongation_type = "ddf"CarpetX::prolongation_order = 5CarpetX::interpolation_order = 3# Initial data are set up by the file readerADMBase::initial_data = "none"ADMBase::initial_lapse = "none"ADMBase::initial_shift = "none"BrillLindquist::x0 = 0.0BrillLindquist::mass = 1.0IO::out_mode = "np"IO::out_proc_every = 1CarpetX::in_dir = "/Users/eschnett/simulations/brill-lindquist-write/output-0000/brill-lindquist-write"CarpetX::in_file = "brill-lindquist-write"CarpetX::in_silo_vars = "ADMBase::metricADMBase::curvADMBase::lapseADMBase::dtlapseADMBase::shiftADMBase::dtshift"IO::out_dir = $parfileIO::out_every = 1CarpetX::out_plotfile_groups = ""CarpetX::out_silo_vars = "all"CarpetX::out_tsv = no
ActiveThorns = "ADMBaseBrillLindquistCarpetXFormalineIOUtilODESolversZ4c"CarpetX::verbose = yesCactus::cctk_show_schedule = yesCactus::terminate = "time"Cactus::cctk_final_time = 0.0CarpetX::xmin = -1.0CarpetX::ymin = -1.0CarpetX::zmin = -1.0CarpetX::xmax = +1.0CarpetX::ymax = +1.0CarpetX::zmax = +1.0CarpetX::ncells_x = 64CarpetX::ncells_y = 64CarpetX::ncells_z = 64CarpetX::periodic_x = noCarpetX::periodic_y = noCarpetX::periodic_z = noCarpetX::periodic = no# CarpetX::reflection_x = yes# CarpetX::reflection_y = yes# CarpetX::reflection_z = yes# CarpetX::reflection_upper_x = yes# CarpetX::reflection_upper_y = yes# CarpetX::reflection_upper_z = yesCarpetX::ghost_size = 3CarpetX::max_num_levels = 1CarpetX::regrid_every = 0CarpetX::prolongation_type = "ddf"CarpetX::prolongation_order = 5CarpetX::interpolation_order = 3ADMBase::initial_data = "Brill-Lindquist"ADMBase::initial_lapse = "Brill-Lindquist"BrillLindquist::x0 = 0.0BrillLindquist::mass = 1.0IO::out_dir = $parfileIO::out_every = 1IO::out_mode = "np"IO::out_proc_every = 1CarpetX::out_plotfile_groups = ""CarpetX::out_silo_vars = "all"CarpetX::out_tsv = no
# USES STRING filereader_ID_dir# USES STRING filereader_ID_files## USES STRING checkpoint_dir# USES STRING recover_dir
////////////////////////////////////////////////////////////////////////////////int InputGH(const cGH *restrict cctkGH) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;static Timer timer("InputGH");Interval interval(timer);const bool is_root = CCTK_MyProc(nullptr) == 0;if (is_root)cout << "InputGH: iteration " << cctk_iteration << ", time " << cctk_time<< "\n";#ifdef HAVE_CAPABILITY_Silo// TODO: Stop at paramcheck time when Silo input parameters are// set, but Silo is not availableInputSilo(cctkGH);#endif
////////////////////////////////////////////////////////////////////////////////#if 0// Inputint Input(cGH *const cctkGH, const char *const basefilename,const int called_from) {DECLARE_CCTK_PARAMETERS;assert(called_from == CP_RECOVER_PARAMETERS ||called_from == CP_RECOVER_DATA || called_from == FILEREADER_DATA);const bool do_recover =called_from == CP_RECOVER_PARAMETERS || called_from == CP_RECOVER_DATA;const bool read_parameters = called_from == CP_RECOVER_PARAMETERS;static Timer timer("InputGH");Interval interval(timer);const string projectname = [&]() -> string {if (do_recover) {// basename is only passed for CP_RECOVER_PARAMETERS, and needs to// be rememberedstatic string saved_basefilename;if (read_parameters)saved_basefilename = basefilename;assert(!saved_basefilename.empty());return saved_basefilename;} else {return basefilename;}}();if (do_recover)if (read_parameters)CCTK_VINFO("Recovering parameters from file \"%s\"",projectname.c_str());elseCCTK_VINFO("Recovering variables from file \"%s\"", projectname.c_str());elseCCTK_VINFO("Reading variables from file \"%s\"", projectname.c_str());if (do_recover && !read_parameters) {// Set global Cactus variables// CCTK_SetMainLoopIndex(main_loop_index);#warning "TODO"// cctkGH->cctk_iteration = iteration;cctkGH->cctk_iteration = 0;}
// Determine which variables to readconst auto ioUtilGH =static_cast<const ioGH *>(CCTK_GHExtension(cctkGH, "IO"));const vector<bool> groups_enabled = [&] {vector<bool> enabled(CCTK_NumGroups(), false);if (do_recover) {for (int gi = 0; gi < CCTK_NumGroups(); ++gi) {if (!CCTK_QueryGroupStorageI(cctkGH, gi))continue;cGroup gdata;int ierr = CCTK_GroupData(gi, &gdata);assert(!ierr);// Do not recover groups with a "checkpoint=no" tagconst int len =Util_TableGetString(gdata.tagstable, 0, nullptr, "checkpoint");if (len > 0) {array<char, 10> buf;Util_TableGetString(gdata.tagstable, buf.size(), buf.data(),"checkpoint");if (CCTK_EQUALS(buf.data(), "no"))continue;assert(CCTK_EQUALS(buf.data(), "yes"));}enabled.at(gi) = true;}} else {for (int gi = 0; gi < CCTK_NumGroups(); ++gi) {const int v0 = CCTK_FirstVarIndexI(gi);assert(v0 >= 0);const int nv = CCTK_NumVarsInGroupI(gi);assert(nv >= 0);for (int vi = 0; vi < nv; ++vi)if (!ioUtilGH->do_inVars || ioUtilGH->do_inVars[v0 + vi])enabled.at(gi) = true;}}return enabled;}();#warning "CONTINUE HERE"#if 0io_dir_t io_dir = do_recover ? io_dir_t::recover : io_dir_t::input;bool did_read_parameters = false;bool did_read_grid_structure = false;if (not input_file_hdf5_ptrs.count(projectname)) {static HighResTimer::HighResTimer timer1("SimulationIO::read_file_hdf5");auto timer1_clock = timer1.start();input_file_hdf5_ptrs[projectname] = make_unique<input_file_t>(io_dir, projectname, file_format::hdf5, iteration, -1, -1);timer1_clock.stop(0);}const auto &input_file_ptr = input_file_hdf5_ptrs.at(projectname);if (read_parameters) {static HighResTimer::HighResTimer timer1("SimulationIO::read_parameters_hdf5");auto timer1_clock = timer1.start();input_file_ptr->read_params();did_read_parameters = true;timer1_clock.stop(0);} else {if (not did_read_grid_structure) {static HighResTimer::HighResTimer timer1("SimulationIO::read_grid_structure_hdf5");auto timer1_clock = timer1.start();input_file_ptr->read_grid_structure(cctkGH);did_read_grid_structure = true;timer1_clock.stop(0);}static HighResTimer::HighResTimer timer1("SimulationIO::read_variables_hdf5");auto timer1_clock = timer1.start();input_file_ptr->read_vars(input_vars, -1, -1);timer1_clock.stop(0);}if (read_parameters)return did_read_parameters ? 1 : 0;assert(did_read_grid_structure);#endifreturn 0; // no error}////////////////////////////////////////////////////////////////////////////////// Recoveringextern "C" int CarpetX_RecoverParameters() {DECLARE_CCTK_PARAMETERS;const char *const out_extension = ".silo";const int iret = IOUtil_RecoverParameters(Input, out_extension, "CarpetX");return iret;}#endif
////////////////////////////////////////////////////////////////////////////////void InputSilo(const cGH *restrict const cctkGH) {DECLARE_CCTK_ARGUMENTS;DECLARE_CCTK_PARAMETERS;int ierr;// Set up timersstatic Timer timer("InputSilo");Interval interval(timer);static Timer timer_setup("InputSilo.setup");auto interval_setup = make_unique<Interval>(timer_setup);const MPI_Comm mpi_comm = MPI_COMM_WORLD;const int myproc = CCTK_MyProc(cctkGH);const int nprocs = CCTK_nProcs(cctkGH);const int ioproc_every = [&]() {if (CCTK_EQUALS(out_mode, "proc"))return 1;if (CCTK_EQUALS(out_mode, "np"))return out_proc_every;if (CCTK_EQUALS(out_mode, "onefile"))return nprocs;assert(0);}();assert(ioproc_every > 0);const int myioproc = myproc / ioproc_every * ioproc_every;assert(myioproc <= myproc);assert(myproc < myioproc + ioproc_every);const bool read_file = myproc % ioproc_every == 0;assert((myioproc == myproc) == read_file);const int metafile_ioproc = 0;const bool read_metafile = myproc == metafile_ioproc;// Configure Silo libraryDBShowErrors(DB_ALL_AND_DRVR, nullptr);// DBSetAllowEmptyObjects(1);DBSetCompression("METHOD=GZIP");DBSetEnableChecksums(1);// Determine input file nameconst string parfilename = [&]() {string buf(1024, '\0');const int len =CCTK_ParameterFilename(buf.length(), const_cast<char *>(buf.data()));buf.resize(len);return buf;}();// TODO: Check wether this parameter contains multiple file names// const string simulation_name = filereader_ID_files;const string simulation_name = in_file;// TODO: directories instead of carefully chosen names// Find input groupsconst vector<bool> group_enabled = [&]() {vector<bool> enabled(CCTK_NumGroups(), false);const auto callback{[](const int index, const char *const optstring, void *const arg) {vector<bool> &enabled = *static_cast<vector<bool> *>(arg);enabled.at(CCTK_GroupIndexFromVarI(index)) = true;}};CCTK_TraverseString(in_silo_vars, callback, &enabled, CCTK_GROUP_OR_VAR);if (io_verbose) {CCTK_VINFO("Silo input for groups:");for (int gi = 0; gi < CCTK_NumGroups(); ++gi) {if (group_enabled.at(gi)) {char *const groupname = CCTK_GroupName(gi);CCTK_VINFO(" %s", groupname);free(groupname);}}}return enabled;}();const auto num_in_vars =count(group_enabled.begin(), group_enabled.end(), true);if (num_in_vars == 0)return;
constexpr int ndims = dim;interval_setup = nullptr;static Timer timer_data("InputSilo.data");auto interval_data = make_unique<Interval>(timer_data);// Read data{const string subdirname = make_subdirname(simulation_name, cctk_iteration);// TODO// const string pathname = string(filereader_ID_dir) + "/" + subdirname;const string pathname = string(in_dir) + "/" + subdirname;DB::ptr<DBfile> file;if (read_file) {const string filename =pathname + "/" +make_filename(simulation_name, cctk_iteration, myproc / ioproc_every);// We could use DB_UNKNOWN instead of DB_HDF5file = DB::make(DBOpen(filename.c_str(), DB_HDF5, DB_READ));assert(file);}// Loop over levelsfor (const auto &leveldata : ghext->leveldata) {if (io_verbose)CCTK_VINFO("Reading level %d", leveldata.level);// Loop over groupsset<mesh_props_t> have_meshes;for (int gi = 0; gi < CCTK_NumGroups(); ++gi) {if (!group_enabled.at(gi))continue;if (CCTK_GroupTypeI(gi) != CCTK_GF)continue;if (io_verbose) {char *const groupname = CCTK_GroupName(gi);CCTK_VINFO(" Reading group %s", groupname);free(groupname);}auto &groupdata = *leveldata.groupdata.at(gi);const int numvars = groupdata.numvars;const int tl = 0;MultiFab &mfab = *groupdata.mfab[tl];const IndexType &indextype = mfab.ixType();const IntVect &ngrow = mfab.nGrowVect();const DistributionMapping &dm = mfab.DistributionMap();const mesh_props_t mesh_props{ngrow};const bool have_mesh = have_meshes.count(mesh_props);// Loop over components (AMReX boxes)const int nfabs = dm.size();for (int component = 0; component < nfabs; ++component) {if (io_verbose)CCTK_VINFO(" Reading component %d", component);const int proc = dm[component];const int ioproc = proc / ioproc_every * ioproc_every;const bool recv_this_fab = proc == myproc;const bool read_this_fab = ioproc == myproc;if (!(recv_this_fab || read_this_fab))continue;const Box &fabbox = mfab.fabbox(component); // exteriorarray<int, ndims> dims;for (int d = 0; d < ndims; ++d)dims[d] = fabbox.length(d);ptrdiff_t zonecount = 1;for (int d = 0; d < ndims; ++d)zonecount *= dims[d];assert(zonecount >= 0 && zonecount <= INT_MAX);// Communicate variable, part 1static Timer timer_mpi("InputSilo.mpi");auto interval_mpi = make_unique<Interval>(timer_mpi);const int mpi_tag = 22901; // randomly chosenvector<double> buffer;MPI_Request mpi_req;double *data = nullptr;if (recv_this_fab && read_this_fab) {FArrayBox &fab = mfab[component];data = fab.dataPtr();} else if (recv_this_fab) {FArrayBox &fab = mfab[component];assert(numvars * zonecount <= INT_MAX);MPI_Irecv(fab.dataPtr(), numvars * zonecount, MPI_DOUBLE, ioproc,mpi_tag, mpi_comm, &mpi_req);} else {buffer.resize(numvars * zonecount);assert(numvars * zonecount <= INT_MAX);data = buffer.data();}interval_mpi = nullptr;// Read variableif (read_file) {static Timer timer_var("InputSilo.var");Interval interval_var(timer_var);const string meshname = make_meshname(leveldata.level, component);const int centering = [&]() {if (indextype.nodeCentered())return DB_NODECENT;if (indextype.cellCentered())return DB_ZONECENT;assert(0);}();for (int vi = 0; vi < numvars; ++vi) {const string varname =make_varname(gi, vi, leveldata.level, component);if (io_verbose)CCTK_VINFO(" Reading variable %s", varname.c_str());const DB::ptr<DBquadvar> quadvar =DB::make(DBGetQuadvar(file.get(), varname.c_str()));assert(quadvar);assert(quadvar->ndims == ndims);assert(ndims <= 3);for (int d = 0; d < ndims; ++d)assert(quadvar->dims[d] == dims[d]);assert(quadvar->datatype == DB_DOUBLE);assert(quadvar->centering == centering);assert(quadvar->nvals == 1);// TODO: check DBOPT_COORDSYS: int cartesian = DB_CARTESIAN;const int column_major = 0;assert(quadvar->major_order == column_major);const void *const read_ptr = quadvar->vals[0];void *const data_ptr = data + vi * zonecount;memcpy(data_ptr, read_ptr, zonecount * sizeof(double));} // for vi} // if read_file// Communicate variable, part 2static Timer timer_wait("InputSilo.wait");auto interval_wait = make_unique<Interval>(timer_wait);if (recv_this_fab && read_this_fab) {// do nothing} else if (recv_this_fab) {MPI_Wait(&mpi_req, MPI_STATUS_IGNORE);} else {buffer.resize(numvars * zonecount);assert(numvars * zonecount <= INT_MAX);MPI_Send(buffer.data(), numvars * zonecount, MPI_DOUBLE, proc,mpi_tag, mpi_comm);}if (recv_this_fab)for (int vi = 0; vi < numvars; ++vi)groupdata.valid.at(tl).at(vi).set(valid_t(true), [] { return "read from file"; });interval_wait = nullptr;} // for component} // for gi} // for leveldata} // write datainterval_data = nullptr;}////////////////////////////////////////////////////////////////////////////////