INHERITS: IO

BOOLEAN verbose "Verbose output"
{
} "no"

  1:* :: ""

  1:* :: "must be a multiple of the blocking factor"

  1:* :: ""

  1:* :: "must be a multiple of the blocking factor"

{
  1:* :: "must be a multiple of the blocking factor"
} 128
CCTK_INT blocking_factor_x "Blocking factor"
{
  1:* :: ""
} 8
CCTK_INT blocking_factor_y "Blocking factor"

} 128

} 8
CCTK_INT blocking_factor_z "Blocking factor"
{
  1:* :: "must be a power of 2"
} 8
CCTK_INT max_grid_size_x "Maximum grid size"
{
  1:* :: "must be a multiple of the blocking factor"
} 32
CCTK_INT max_grid_size_y "Maximum grid size"
{
  1:* :: "must be a multiple of the blocking factor"
} 32
CCTK_INT max_grid_size_z "Maximum grid size"
{
  1:* :: "must be a multiple of the blocking factor"
} 32
CCTK_INT max_tile_size_x "Maximum tile size"
{
  1:* :: ""
} 1024000
CCTK_INT max_tile_size_y "Maximum tile size"
{
  1:* :: ""
} 16
CCTK_INT max_tile_size_z "Maximum tile size"
{
  1:* :: ""
} 32

CCTK_INT out_every "How often to output" STEERABLE=always

CCTK_REAL dtfac "The standard timestep condition dt = dtfac*min(delta_space)"

  0 :: "no output"
  1:* :: "every that many iterations"
} 1

  *:* :: ""
} 0.5 

CCTK_REAL dtfac "The standard timestep condition dt = dtfac*min(delta_space)"

BOOLEAN out_tsv "Output in TSV (tab-separated values) format?" STEERABLE=always

  *:* :: ""
} 0.5 

} no

SHARES: IO
USES STRING out_dir
USES CCTK_INT out_every

#include <AMReX_ParmParse.H>

  DECLARE_CCTK_PARAMETERS;

  CCTK_VINFO("ErrorEst level %d", level);

  if (verbose)
    CCTK_VINFO("ErrorEst level %d", level);

  auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling({1024000, 16, 32});

  auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
      {max_tile_size_x, max_tile_size_y, max_tile_size_z});

    // Note: This excludes ghosts, it's not the proper Cactus lsh

    const Array4<CCTK_REAL> &ctagarr = groupdata.mfab.at(tl)->array(mfi);
    const CCTK_REAL *restrict ctags = ctagarr.ptr(imin.x, imin.y, imin.z, vi);
    const Array4<char> &atagarr = tags.array(mfi);
    char *restrict atags = atagarr.ptr(imin.x, imin.y, imin.z, vi);

    const Array4<CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
    const CCTK_REAL *restrict ptr = vars.ptr(imin.x, imin.y, imin.z, vi);
    const Array4<char> &tagarr = tags.array(mfi);

          atags[idx] = ctags[idx] == 0.0 ? TagBox::CLEAR : TagBox::SET;

          tagarr(imin.x + i, imin.y + j, imin.z + k) =
              ptr[idx] == 0.0 ? TagBox::CLEAR : TagBox::SET;

  CCTK_VINFO("MakeNewLevelFromScratch level %d", lev);

  DECLARE_CCTK_PARAMETERS;
  if (verbose)
    CCTK_VINFO("MakeNewLevelFromScratch level %d", lev);

  CCTK_VINFO("MakeNewLevelFromCoarse level %d", lev);

  DECLARE_CCTK_PARAMETERS;
  if (verbose)
    CCTK_VINFO("MakeNewLevelFromCoarse level %d", lev);

  CCTK_VINFO("RemakeLevel level %d", lev);

  DECLARE_CCTK_PARAMETERS;
  if (verbose)
    CCTK_VINFO("RemakeLevel level %d", lev);

  CCTK_VINFO("ClearLevel level %d", lev);

  DECLARE_CCTK_PARAMETERS;
  if (verbose)
    CCTK_VINFO("ClearLevel level %d", lev);

  CCTK_VINFO("Startup");

  DECLARE_CCTK_PARAMETERS;
  if (verbose)
    CCTK_VINFO("Startup");

  CCTK_VINFO("SetupGH");

  DECLARE_CCTK_PARAMETERS;
  if (verbose)
    CCTK_VINFO("SetupGH");

  CCTK_VINFO("InitGH");

  if (verbose)
    CCTK_VINFO("InitGH");

  // Set blocking factors via parameter table since AmrMesh needs to
  // know them when its constructor is running, but there are no
  // constructor arguments for them
  ParmParse pp;
  pp.add("amr.blocking_factor_x", blocking_factor_x);
  pp.add("amr.blocking_factor_y", blocking_factor_y);
  pp.add("amr.blocking_factor_z", blocking_factor_z);
  pp.add("amr.max_grid_size_x", max_grid_size_x);
  pp.add("amr.max_grid_size_y", max_grid_size_y);
  pp.add("amr.max_grid_size_z", max_grid_size_z);

#warning "TODO: increase blocking factor"
  const int blocking_factor = 8;
  // const int blocking_factor = 1;
  ghext->amrcore->SetBlockingFactor(blocking_factor);
  const int max_grid_size = 32;
  ghext->amrcore->SetMaxGridSize(max_grid_size);

  int maxnumlevels = ghext->amrcore->maxLevel() + 1;
  for (int level = 0; level < maxnumlevels; ++level) {
    CCTK_VINFO("Geometry level %d:", level);
    cout << ghext->amrcore->Geom(level) << "\n";

  if (verbose) {
    int maxnumlevels = ghext->amrcore->maxLevel() + 1;
    for (int level = 0; level < maxnumlevels; ++level) {
      CCTK_VINFO("Geometry level %d:", level);
      cout << ghext->amrcore->Geom(level) << "\n";
    }

  CCTK_VINFO("CreateRefinedGrid level %d", level);

  DECLARE_CCTK_PARAMETERS;
  if (verbose)
    CCTK_VINFO("CreateRefinedGrid level %d", level);

  cout << "CRG numlevels=" << numlevels << "\n";
  cout << "CRG maxnumlevels=" << maxnumlevels << "\n";

  CCTK_VINFO("SetupLevel level %d", level);

  if (verbose)
    CCTK_VINFO("SetupLevel level %d", level);

  CCTK_VINFO("ScheduleTraverseGH [%d] %s", cctkGH->cctk_iteration, where);

  DECLARE_CCTK_PARAMETERS;
  if (verbose)
    CCTK_VINFO("ScheduleTraverseGH [%d] %s", cctkGH->cctk_iteration, where);

  CCTK_VINFO("Shutdown");

  DECLARE_CCTK_PARAMETERS;
  if (verbose)
    CCTK_VINFO("Shutdown");

#include <AMReX_Orientation.H>

#include <fstream>

namespace {
// Convert a (direction, face) pair to an AMReX Orientation
Orientation orient(int d, int f) {
  return Orientation(d, Orientation::Side(f));
}
} // namespace
template <typename T>
void write_arrays(ostream &os, const cGH *restrict cctkGH, int level,
                  int component, const vector<const T *> ptrs,
                  const int *restrict ash, const int *restrict lbnd,
                  const int *restrict lsh, const int *restrict bbox,
                  const int *restrict nghostzones) {
  DECLARE_CCTK_ARGUMENTS;
  const int levfac = 1 << level;
  CCTK_REAL x0[dim], dx[dim];
  for (int d = 0; d < dim; ++d) {
    dx[d] = cctk_delta_space[d] / levfac;
    x0[d] = cctk_origin_space[d] + 0.5 * dx[d];
  }
  int imin[dim], imax[dim];
  for (int d = 0; d < dim; ++d) {
    imin[d] = cctk_bbox[2 * d] ? 0 : cctk_nghostzones[d];
    imax[d] = cctk_lsh[d] - (cctk_bbox[2 * d + 1] ? 0 : cctk_nghostzones[d]);
  }
  constexpr int di = 1;
  const int dj = di * ash[0];
  const int dk = dj * ash[1];
  for (int k = 0; k < lsh[2]; ++k) {
    for (int j = 0; j < lsh[1]; ++j) {
      for (int i = 0; i < lsh[0]; ++i) {
        const int idx = di * i + dj * j + dk * k;
        const int isghost = !(i >= imin[0] && i < imax[0] && j >= imin[1] &&
                              j < imax[1] && k >= imin[2] && k < imax[2]);
        T x = x0[0] + (lbnd[0] + i) * dx[0];
        T y = x0[1] + (lbnd[1] + j) * dx[1];
        T z = x0[2] + (lbnd[2] + k) * dx[2];
        os << cctk_iteration << "\t" << cctk_time << "\t" << level << "\t"
           << component << "\t" << isghost << "\t" << (lbnd[0] + i) << "\t"
           << (lbnd[1] + j) << "\t" << (lbnd[2] + k) << "\t" << x << "\t" << y
           << "\t" << z;
        for (const auto &ptr : ptrs)
          os << "\t" << ptr[idx];
        os << "\n";
      }
    }
  }
}
void WriteASCII(const cGH *restrict cctkGH, const string &filename, int gi,
                const vector<string> &varnames) {
  ostringstream buf;
  buf << filename << ".tsv";
  ofstream file(buf.str());

  // Output header
  file << "#"
       << "\t"
       << "time"
       << "\t"
       << "level"
       << "\t"
       << "component"
       << "\t"
       << "isghost"
       << "\t"
       << "i"
       << "\t"
       << "j"
       << "\t"
       << "k"
       << "\t"
       << "x"
       << "\t"
       << "y"
       << "\t"
       << "z";
  for (const auto &varname : varnames)
    file << "\t" << varname;
  file << "\n";
  for (const auto &leveldata : ghext->leveldata) {
    const auto &groupdata = leveldata.groupdata.at(gi);
    const int tl = 0;
    const MultiFab &mfab = *groupdata.mfab.at(tl);
    for (MFIter mfi(mfab); mfi.isValid(); ++mfi) {
      const Box &fbx = mfi.fabbox();       // allocated array
      const Box &bx = mfi.tilebox();       // current region (without ghosts)
      const Box &gbx = mfi.growntilebox(); // current region (with ghosts)
      const Dim3 imin = lbound(gbx);
      // Local shape
      int lsh[dim];
      for (int d = 0; d < dim; ++d)
        lsh[d] = gbx[orient(d, 1)] - gbx[orient(d, 0)] + 1;
      // Allocated shape
      int ash[dim];
      for (int d = 0; d < dim; ++d)
        ash[d] = fbx[orient(d, 1)] - fbx[orient(d, 0)] + 1;
      // Local extent
      int lbnd[dim];
      for (int d = 0; d < dim; ++d)
        lbnd[d] = gbx[orient(d, 0)];
      // Boundaries
      int bbox[2 * dim];
      for (int d = 0; d < dim; ++d)
        for (int f = 0; f < 2; ++f)
          bbox[2 * d + f] = bx[orient(d, f)] == gbx[orient(d, f)];
      // Number of ghost zones
      int nghostzones[dim];
      for (int d = 0; d < dim; ++d)
        nghostzones[d] = mfab.fb_nghost[d];
      const Array4<CCTK_REAL> &vars = groupdata.mfab.at(tl)->array(mfi);
      vector<const CCTK_REAL *> ptrs(groupdata.numvars);
      for (int vi = 0; vi < groupdata.numvars; ++vi)
        ptrs.at(vi) = vars.ptr(imin.x, imin.y, imin.z, vi);
      write_arrays(file, cctkGH, leveldata.level, mfi.index(), ptrs, ash, lbnd,
                   lsh, bbox, nghostzones);
    }
  }
  file.close();
}

  const bool is_root = CCTK_MyProc(nullptr) == 0;
  if (is_root)
    cout << "OutputGH: iteration " << cctk_iteration << ", time " << cctk_time
         << ", run time " << CCTK_RunTime() << " s\n";

      buf << "wavetoy/" << groupname;

      buf << out_dir << "/" << groupname;

      if (out_tsv) {
        ostringstream procbuf;
        procbuf << out_dir << "/" << groupname;
        procbuf << ".it" << setw(6) << setfill('0') << cctk_iteration;
        procbuf << ".p" << setw(4) << setfill('0') << CCTK_MyProc(nullptr);
        string procfilename = procbuf.str();

        WriteASCII(cctkGH, procfilename, gi, varnames);
      }

      CCTK_VINFO("maxabs(%s)=%g", CCTK_VarName(groupdata.firstvarindex + vi),
                 red.maxabs);

      if (is_root)
        cout << "  "
             << unique_ptr<char>(CCTK_FullName(groupdata.firstvarindex + vi))
                    .get()
             << ": maxabs=" << red.maxabs << "\n";

#include <cctk_Parameters.h>

reduction<T> reduce_array(const Geometry &geom, const T *restrict var,

reduction<T> reduce_array(const Geometry &geom, const T *restrict ptr,

        red = reduction<T>(red, reduction<T>(dV, var[idx]));

        red = reduction<T>(red, reduction<T>(dV, ptr[idx]));

  DECLARE_CCTK_PARAMETERS;

      auto mfitinfo =
          MFItInfo().SetDynamic(true).EnableTiling({1024000, 16, 32});

      auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
          {max_tile_size_x, max_tile_size_y, max_tile_size_z});

        // Note: This excludes ghosts, it's not the proper Cactus lsh

        const CCTK_REAL *restrict var = vars.ptr(imin.x, imin.y, imin.z, vi);

        const CCTK_REAL *restrict ptr = vars.ptr(imin.x, imin.y, imin.z, vi);

        red += reduce_array(geom, var, ash, lsh);

#warning "TODO: Skip refined regions"
        red += reduce_array(geom, ptr, ash, lsh);

  for (int d = 0; d < dim; ++d)

  for (int d = 0; d < dim; ++d) {

  for (int d = 0; d < dim; ++d)

  }

  // CCTK_REAL mindx = 1.0 / 0.0;
  // const int numlevels = ghext->amrcore->finestLevel() + 1;
  // for (int level = 0; level < numlevels; ++level) {
  //   const Geometry &geom = ghext->amrcore->Geom(level);
  //   const CCTK_REAL *restrict dx = geom.CellSize();
  //   for (int d = 0; d < dim; ++d)
  //     mindx = fmin(mindx, dx[d]);
  // }

  const int numlevels = ghext->amrcore->finestLevel() + 1;
  for (int level = 0; level < numlevels; ++level) {
    const Geometry &geom = ghext->amrcore->Geom(level);
    const CCTK_REAL *restrict dx = geom.CellSize();
    for (int d = 0; d < dim; ++d)
      mindx = fmin(mindx, dx[d]);
  }

  for (int d = 0; d < dim; ++d)
    mindx = fmin(mindx, dx[d]);
  mindx = mindx / (1 << (max_num_levels - 1));

  // TODO: Get this from mfab

  // TODO: Get this from mfab (mfab.fb_ghosts)

  // Offset between this level's and the coarsest level's origin

  // Offset between this level's and the coarsest level's origin as multiple of
  // the grid spacing

    cctkGH->cctk_lsh[d] = bx[orient(d, 1)] - bx[orient(d, 0)] + 1;

    cctkGH->cctk_lsh[d] = gbx[orient(d, 1)] - gbx[orient(d, 0)] + 1;

    cctkGH->cctk_lbnd[d] = bx[orient(d, 0)];
    cctkGH->cctk_ubnd[d] = bx[orient(d, 1)];

    cctkGH->cctk_lbnd[d] = gbx[orient(d, 0)];
    cctkGH->cctk_ubnd[d] = gbx[orient(d, 1)];

      cctkGH->cctk_bbox[2 * d + f] = bx[orient(d, f)] != gbx[orient(d, f)];

      cctkGH->cctk_bbox[2 * d + f] = bx[orient(d, f)] == gbx[orient(d, f)];

  const Dim3 imin = lbound(bx);

  const Dim3 imin = lbound(gbx);

    for (int d = 0; d < dim; ++d)

    for (int d = 0; d < dim; ++d) {

    for (int d = 0; d < dim; ++d)
      x0[d] = cctkGH->cctk_origin_space[d] - 0.0 * dx[d];
    for (int d = 0; d < dim; ++d)

      x0[d] = cctkGH->cctk_origin_space[d];

    }

    CCTK_VINFO("  dx=[%g,%g,%g]", dx[0], dx[1], dx[2]);

    CCTK_VINFO("  base dx=[%g,%g,%g]", dx[0], dx[1], dx[2]);

    CCTK_VINFO("Recovering from checkpoint...");

    CCTK_VINFO("Setting up initial conditions...");

      CCTK_VINFO("Initializing level %d...", current_level);

      cout << "I old_numlevels=" << old_numlevels << "\n";
      cout << "I new_numlevels=" << new_numlevels << "\n";
      cout << "I max_numlevels=" << int(ghext->amrcore->maxLevel() + 1) << "\n";

      if (new_numlevels <= old_numlevels)

      const bool did_create_new_level = new_numlevels > old_numlevels;
      if (!did_create_new_level)

  CCTK_VINFO("Initialized %d levels", int(ghext->leveldata.size()));

  // TODO: Get ghost_size from mfab

  CCTK_VINFO("Starting evolution...");

  CCTK_VINFO("Shutting down...");

  DECLARE_CCTK_PARAMETERS;

  CCTK_VINFO("CallFunction [%d] %s: %s::%s", cctkGH->cctk_iteration,
             attribute->where, attribute->thorn, attribute->routine);

  if (verbose)
    CCTK_VINFO("CallFunction [%d] %s: %s::%s", cctkGH->cctk_iteration,
               attribute->where, attribute->thorn, attribute->routine);

        auto mfitinfo =
            MFItInfo().SetDynamic(true).EnableTiling({1024000, 16, 32});

        auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
            {max_tile_size_x, max_tile_size_y, max_tile_size_z});

#warning "TODO: during evolution, sync only one time level"

        // Refined level: Prolongate from next coarser level, and then
        // copy from adjacent boxes on same level

        // Refined level: Prolongate boundaries from next coarser
        // level, and copy from adjacent boxes on same level

        CellBilinear interp;

        // CellBilinear interp;

        for (int tl = 0; tl < sync_tl; ++tl)
          InterpFromCoarseLevel(
              *groupdata.mfab.at(tl), 0.0, *coarsegroupdata.mfab.at(tl), 0, 0,
              groupdata.numvars, ghext->amrcore->Geom(level - 1),
              ghext->amrcore->Geom(level), cphysbc, 0, fphysbc, 0, reffact,
              &interp, bcs, 0);

        for (int tl = 0; tl < sync_tl; ++tl) {
          // InterpFromCoarseLevel(
          //     *groupdata.mfab.at(tl), 0.0, *coarsegroupdata.mfab.at(tl), 0,
          //     0, groupdata.numvars, ghext->amrcore->Geom(level - 1),
          //     ghext->amrcore->Geom(level), cphysbc, 0, fphysbc, 0, reffact,
          //     &interp, bcs, 0);
#warning "TODO: make copy of fine level"
          FillPatchTwoLevels(
              *groupdata.mfab.at(tl), 0.0, {&*coarsegroupdata.mfab.at(tl)},
              {0.0}, {&*groupdata.mfab.at(tl)}, {0.0}, 0, 0, groupdata.numvars,
              ghext->amrcore->Geom(level - 1), ghext->amrcore->Geom(level),
              cphysbc, 0, fphysbc, 0, reffact, &cell_bilinear_interp, bcs, 0);
        }

#warning "TODO: during evolution, restrict only one time level"

    for (int tl = 0; tl < restrict_tl; ++tl) {

    for (int tl = 0; tl < restrict_tl; ++tl)

    }

2. First-order equations:
Scalar wave equation:
  phi,tt = phi,xx
Introduce new variable:
  psi = phi,t
New set of equations:
  phi,t = psi
  psi,t = phi,tt = phi,xx
Discretization:
  phi_0 = phi_1 + dt psi_1
  psi_0 = psi_1 + dt phi_1,xx
3. Fancy equations:
Scalar wave equation:
  phi,tt = phi,xx
Direct second order discretization:
  phi_0 = 2 phi_1 - phi_2 + dt^2 phi_1,xx
Introduce new variable:
  psi = phi,t
New set of equations:
  phi,t = psi
  psi,t = phi,tt = phi,xx
New, equivalent discretization:
  phi_0 - phi_1 = dt psi_1
  phi_0 - 2 phi_1 + phi_2 = dt^2 phi_1,xx
  dt psi_1 - dt psi_2 = dt^2 phi_1,xx
  psi_1 = psi_2 + dt^2 phi_1,xx
  phi_0 = phi_1 + dt psi_1
  psi_0 = psi_1 + dt phi_0,xx
        = psi_1 + dt phi_1,xx + dt^2 psi_1,xx

  # psi

CCTK_REAL energy TYPE=gf
{
  eps
} "Energy density for wave equation"

  err

  phierr
  # psierr

SCHEDULE WaveToyAMReX_Energy AT analysis
{
  LANG: C
  SYNC: energy
} "Calculate energy density for the wave equation"

#include <iostream>

// Time derivative
template <typename T> auto timederiv(T f(T t, T x, T y, T z), T dt) {
  return [=](T t, T x, T y, T z) {
    // return (f(t + dt, x, y, z) - f(t, x, y, z)) / dt;
    return (f(t, x, y, z) - f(t - dt, x, y, z)) / dt;
  };
}

  T Lx = 2, Ly = 2, Lz = 2;
  T kx = 2 * M_PI / Lx, ky = 2 * M_PI / Ly, kz = 2 * M_PI / Lz;

  T kx = 2 * M_PI / 2, ky = 2 * M_PI / 2, kz = 2 * M_PI / 2;
  // T kx = 2 * M_PI / 2, ky = 0, kz = 0;

  for (int k = 0; k < cctk_lsh[2]; ++k) {
    for (int j = 0; j < cctk_lsh[1]; ++j) {

  int imin[dim], imax[dim];
  for (int d = 0; d < dim; ++d) {
    imin[d] = cctk_bbox[2 * d] ? 0 : cctk_nghostzones[d];
    imax[d] = cctk_lsh[d] - (cctk_bbox[2 * d + 1] ? 0 : cctk_nghostzones[d]);
  }
  for (int k = imin[2]; k < imax[2]; ++k) {
    for (int j = imin[1]; j < imax[1]; ++j) {

      for (int i = 0; i < cctk_lsh[0]; ++i) {

      for (int i = imin[0]; i < imax[0]; ++i) {

        // psi[idx] = timederiv(standing, dt)(t, x, y, z);

        // CCTK_REAL r = fabs(x);

  }
  int imin[dim], imax[dim];
  for (int d = 0; d < dim; ++d) {
    imin[d] = cctk_bbox[2 * d] ? 0 : cctk_nghostzones[d];
    imax[d] = cctk_lsh[d] - (cctk_bbox[2 * d + 1] ? 0 : cctk_nghostzones[d]);

  for (int k = 0; k < cctk_lsh[2]; ++k) {
    for (int j = 0; j < cctk_lsh[1]; ++j) {

  for (int k = imin[2]; k < imax[2]; ++k) {
    for (int j = imin[1]; j < imax[1]; ++j) {

      for (int i = 0; i < cctk_lsh[0]; ++i) {

      for (int i = imin[0]; i < imax[0]; ++i) {

        phi[idx] = -phi_p_p[idx] + 2 * phi_p[idx] +

        // CCTK_REAL ddx_psi =
        //     (psi_p[idx - di] - 2 * psi_p[idx] + psi_p[idx + di]) /
        //     sqr(dx[0]);
        // CCTK_REAL ddy_psi =
        //     (psi_p[idx - dj] - 2 * psi_p[idx] + psi_p[idx + dj]) /
        //     sqr(dx[1]);
        // CCTK_REAL ddz_psi =
        //     (psi_p[idx - dk] - 2 * psi_p[idx] + psi_p[idx + dk]) /
        //     sqr(dx[2]);
        // phi[idx] = phi_p[idx] + dt * psi_p[idx];
        // psi[idx] = psi_p[idx] + dt * (ddx_phi + ddy_phi + ddz_phi) +
        //            sqr(dt) * (ddx_psi + ddy_psi + ddz_psi);
        // phi[idx] = phi_p[idx] + dt * psi_p[idx];
        // psi[idx] = psi_p[idx] + dt * (ddx_phi + ddy_phi + ddz_phi);
        phi[idx] = 2 * phi_p[idx] - phi_p_p[idx] +

      }
    }
  }
}
extern "C" void WaveToyAMReX_Energy(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  const CCTK_REAL dt = CCTK_DELTA_TIME;
  CCTK_REAL dx[dim];
  for (int d = 0; d < dim; ++d)
    dx[d] = CCTK_DELTA_SPACE(d);
  int imin[dim], imax[dim];
  for (int d = 0; d < dim; ++d) {
    imin[d] = cctk_bbox[2 * d] ? 0 : cctk_nghostzones[d];
    imax[d] = cctk_lsh[d] - (cctk_bbox[2 * d + 1] ? 0 : cctk_nghostzones[d]);
  }
  constexpr int di = 1;
  const int dj = di * cctk_ash[0];
  const int dk = dj * cctk_ash[1];
  for (int k = imin[2]; k < imax[2]; ++k) {
    for (int j = imin[1]; j < imax[1]; ++j) {
#pragma omp simd
      for (int i = imin[0]; i < imax[0]; ++i) {
        const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
        CCTK_REAL dt_phi = (phi[idx] - phi_p[idx]) / dt;
        CCTK_REAL dx_phi = (phi[idx + di] - phi[idx - di]) / (2 * dx[0]);
        CCTK_REAL dy_phi = (phi[idx + dj] - phi[idx - dj]) / (2 * dx[1]);
        CCTK_REAL dz_phi = (phi[idx + dk] - phi[idx - dk]) / (2 * dx[2]);
        eps[idx] = (sqr(dt_phi) + sqr(dx_phi) + sqr(dy_phi) + sqr(dz_phi)) / 2;

  // const CCTK_REAL dt = CCTK_DELTA_TIME;

        err[idx] = phi[idx] - standing(t, x, y, z);

        phierr[idx] = phi[idx] - standing(t, x, y, z);
        // psierr[idx] = psi[idx] - timederiv(standing, dt)(t, x, y, z);

ActiveThorns = "
    AMReX
    IOUtil
    WaveToyAMReX
"
$maxlev = 1
$ncells = 128
Cactus::cctk_show_schedule = no
Cactus::cctk_itlast = $ncells * 2 ** ($maxlev - 1) * 2
AMReX::ncells_x = $ncells
AMReX::ncells_y = 1
AMReX::ncells_z = 1
AMReX::blocking_factor_x = 1
AMReX::blocking_factor_y = 1
AMReX::blocking_factor_z = 1
AMReX::max_grid_size_x = 1048576
AMReX::max_grid_size_y = 1048576
AMReX::max_grid_size_z = 1048576
AMReX::max_num_levels = $maxlev
AMReX::dtfac = 0.5 / 4
IO::out_dir = "wavetoy-1"
IO::out_every = $ncells * 2 ** ($maxlev - 1) / 2

ActiveThorns = "
    AMReX
    IOUtil
    WaveToyAMReX
"
$maxlev = 2
$ncells = 128
Cactus::cctk_show_schedule = no
Cactus::cctk_itlast = $ncells * 2 ** ($maxlev - 1) * 2
AMReX::ncells_x = $ncells
AMReX::ncells_y = 1
AMReX::ncells_z = 1
AMReX::blocking_factor_x = 1
AMReX::blocking_factor_y = 1
AMReX::blocking_factor_z = 1
AMReX::max_grid_size_x = 1048576
AMReX::max_grid_size_y = 1048576
AMReX::max_grid_size_z = 1048576
AMReX::max_num_levels = $maxlev
AMReX::dtfac = 0.5 / 2
IO::out_dir = "wavetoy-2"
IO::out_every = $ncells * 2 ** ($maxlev - 1) / 2

ActiveThorns = "
    AMReX
    IOUtil
    WaveToyAMReX
"
$maxlev = 3
$ncells = 128
Cactus::cctk_show_schedule = no
Cactus::cctk_itlast = $ncells * 2 ** ($maxlev - 1) * 2
AMReX::ncells_x = $ncells
AMReX::ncells_y = 1
AMReX::ncells_z = 1
AMReX::blocking_factor_x = 1
AMReX::blocking_factor_y = 1
AMReX::blocking_factor_z = 1
AMReX::max_grid_size_x = 1048576
AMReX::max_grid_size_y = 1048576
AMReX::max_grid_size_z = 1048576
AMReX::max_num_levels = $maxlev
AMReX::dtfac = 0.5
IO::out_dir = "wavetoy-3"
IO::out_every = $ncells * 2 ** ($maxlev - 1) / 2

    IOUtil

 
$nlevels = 3
$ncells = 128

AMReX::ncells_x = 128
AMReX::ncells_y = 128
AMReX::ncells_z = 128

Cactus::cctk_itlast = $ncells * 2 ** ($nlevels - 1) * 2
AMReX::ncells_x = $ncells
AMReX::ncells_y = $ncells
AMReX::ncells_z = $ncells

AMReX::max_num_levels = 3

AMReX::max_num_levels = $nlevels

Cactus::cctk_itlast = 128

AMReX::dtfac = 0.5

AMReX::out_every = 128

IO::out_dir = "wavetoy"
IO::out_every = $ncells * 2 ** ($nlevels - 1) / 2