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

STRING in_dir "Input directory" STEERABLE=always
{
  ".*" :: ""
} ""

STRING in_file "Input file name" STEERABLE=always
{
  ".*" :: ""
} ""

STRING in_silo_vars "Variables to input in Silo format" STEERABLE=always
{
  ".*" :: ""
} ""

# USES STRING filereader_ID_dir
# USES STRING filereader_ID_files
# 
# USES STRING checkpoint_dir
# USES STRING recover_dir

#include <CactusBase/IOUtil/src/ioGH.h>
#include <CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h>

#include <util_Table.h>

#include <cassert>

////////////////////////////////////////////////////////////////////////////////
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

  // TODO: This should be the number of variables input
  return 0;
}

////////////////////////////////////////////////////////////////////////////////
#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

int InputGH(const cGH *cctkGH);

#include <algorithm>

////////////////////////////////////////////////////////////////////////////////
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;
}
////////////////////////////////////////////////////////////////////////////////

  // Create output file

  // Determine output file name

  const vector<bool> group_enabled = [&]() {

  const vector<bool> group_enabled = [&] {

    // int dummy = 0;
    // MPI_Bcast(&dummy, 1, MPI_INT, 0, mpi_comm);

          // TODO: Check whether data are valid

void OutputSilo(const cGH *restrict const cctkGH);

void InputSilo(const cGH *cctkGH);
void OutputSilo(const cGH *cctkGH);

}

} // namespace CarpetX

#include "io.hxx"

      InputGH(cctkGH);