EJWFGXTWIBKJ62LO6QM62LM7MGTYEHOCZDKA3UVR5UYWQZ2O54DQC
6D6MZBKTLXUK4ZCIFUPNV7XT64XSFFG4CL6HN6ONK5ONM6MR64AAC
GKKJ75HX2ERLVBZVE2CUB6T3J2SUT7R3UKEKTEYNOG43ZKX6X5MQC
BZ7T742UYYUIK2X3Y6UWJ2G33ER4PJJBA7YLCVAQKAJQUUECHVRQC
RTNZAS3UPI6GG3KY4Z5WVXJ4R2YF5427BB6WAV3GHRS5W7XPOSUQC
722HZ7UFINNE3YKSYKP2NHZ5XEG5QQLQHSKC7PREJZR3EX6RDYUAC
RZPXXQUUAG7JIADVCLYQOY3I4CQDR23PF57L6XDGHNCT4EBPSO6QC
NUOLOGCKMF5UOBGBYEOX4O7NQ5AEVVLCH6KRBQRJQXIRDNJ2C2ZQC
5XGIB7XMEZNBLA5ZLQTXRTC3AZ5CTRGMXWBPVMWXG4DPHKWDF4ZAC
RCLGQ2LZMFVPBPTU2G55DJ6HZPOGGTPZRZCY54VGP6YLHANJ2LAQC
2MTD37PDBF7KZKC6LXG2IVDYMMUMNNS557WMFVSCZM3SAKOGDC5QC
FEMASUBNU32NSG4DNXZX54CGCA57PVRGYO46L3A6F2EJ4BCSJ3SAC
BVR7DVINVPQG7PA6Z7QYVYNQ43YZL7XCC6AOMSMWMGAAB2Q43STAC
MSBBCXVGD3GRLE5KAI6BKAFRV7SQUWI2SNN43AJAUD3ISRCEXY6QC
ActiveThorns = "
ADMBase
BrillLindquist
CarpetX
Formaline
IOUtil
ODESolvers
Z4c
"
CarpetX::verbose = yes
Cactus::cctk_show_schedule = yes
Cactus::terminate = "time"
Cactus::cctk_final_time = 0.0
CarpetX::xmin = -1.0
CarpetX::ymin = -1.0
CarpetX::zmin = -1.0
CarpetX::xmax = +1.0
CarpetX::ymax = +1.0
CarpetX::zmax = +1.0
CarpetX::ncells_x = 64
CarpetX::ncells_y = 64
CarpetX::ncells_z = 64
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::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 = yes
CarpetX::ghost_size = 3
CarpetX::max_num_levels = 1
CarpetX::regrid_every = 0
CarpetX::prolongation_type = "ddf"
CarpetX::prolongation_order = 5
CarpetX::interpolation_order = 3
# Initial data are set up by the file reader
ADMBase::initial_data = "none"
ADMBase::initial_lapse = "none"
ADMBase::initial_shift = "none"
BrillLindquist::x0 = 0.0
BrillLindquist::mass = 1.0
IO::out_mode = "np"
IO::out_proc_every = 1
CarpetX::in_dir = "/Users/eschnett/simulations/brill-lindquist-write/output-0000/brill-lindquist-write"
CarpetX::in_file = "brill-lindquist-write"
CarpetX::in_silo_vars = "
ADMBase::metric
ADMBase::curv
ADMBase::lapse
ADMBase::dtlapse
ADMBase::shift
ADMBase::dtshift
"
IO::out_dir = $parfile
IO::out_every = 1
CarpetX::out_plotfile_groups = ""
CarpetX::out_silo_vars = "all"
CarpetX::out_tsv = no
ActiveThorns = "
ADMBase
BrillLindquist
CarpetX
Formaline
IOUtil
ODESolvers
Z4c
"
CarpetX::verbose = yes
Cactus::cctk_show_schedule = yes
Cactus::terminate = "time"
Cactus::cctk_final_time = 0.0
CarpetX::xmin = -1.0
CarpetX::ymin = -1.0
CarpetX::zmin = -1.0
CarpetX::xmax = +1.0
CarpetX::ymax = +1.0
CarpetX::zmax = +1.0
CarpetX::ncells_x = 64
CarpetX::ncells_y = 64
CarpetX::ncells_z = 64
CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::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 = yes
CarpetX::ghost_size = 3
CarpetX::max_num_levels = 1
CarpetX::regrid_every = 0
CarpetX::prolongation_type = "ddf"
CarpetX::prolongation_order = 5
CarpetX::interpolation_order = 3
ADMBase::initial_data = "Brill-Lindquist"
ADMBase::initial_lapse = "Brill-Lindquist"
BrillLindquist::x0 = 0.0
BrillLindquist::mass = 1.0
IO::out_dir = $parfile
IO::out_every = 1
IO::out_mode = "np"
IO::out_proc_every = 1
CarpetX::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 available
InputSilo(cctkGH);
#endif
////////////////////////////////////////////////////////////////////////////////
#if 0
// Input
int 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 remembered
static 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());
else
CCTK_VINFO("Recovering variables from file \"%s\"", projectname.c_str());
else
CCTK_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 read
const 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" tag
const 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 0
io_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);
#endif
return 0; // no error
}
////////////////////////////////////////////////////////////////////////////////
// Recovering
extern "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 timers
static 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 library
DBShowErrors(DB_ALL_AND_DRVR, nullptr);
// DBSetAllowEmptyObjects(1);
DBSetCompression("METHOD=GZIP");
DBSetEnableChecksums(1);
// Determine input file name
const 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 groups
const 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_HDF5
file = DB::make(DBOpen(filename.c_str(), DB_HDF5, DB_READ));
assert(file);
}
// Loop over levels
for (const auto &leveldata : ghext->leveldata) {
if (io_verbose)
CCTK_VINFO("Reading level %d", leveldata.level);
// Loop over groups
set<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); // exterior
array<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 1
static Timer timer_mpi("InputSilo.mpi");
auto interval_mpi = make_unique<Interval>(timer_mpi);
const int mpi_tag = 22901; // randomly chosen
vector<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 variable
if (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 2
static 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 data
interval_data = nullptr;
}
////////////////////////////////////////////////////////////////////////////////