INCLUDES HEADER: loop.hxx IN loop.hxx
void FUNCTION GetTileExtent(
  CCTK_POINTER_TO_CONST IN cctkGH,
  CCTK_INT ARRAY OUT tile_min,
  CCTK_INT ARRAY OUT tile_max)
PROVIDES FUNCTION GetTileExtent WITH AMReX_GetTileExtent LANGUAGE C

$nlevels = 6
$ncells = 16

$nlevels = 5
$ncells = 32

Cactus::cctk_itlast = $ncells * 2 ** ($nlevels - 1) * 2

Cactus::cctk_itlast = $ncells * 2 ** ($nlevels - 1) / 2

AMReX::regrid_error_threshold = 0.01   #TODO 0.1

AMReX::regrid_error_threshold = 0.01

# WaveToyAMReX::spatial_frequency_x = 0.5
# WaveToyAMReX::spatial_frequency_y = 0.0
# WaveToyAMReX::spatial_frequency_z = 0.0

IO::out_every = 16   #TODO $ncells * 2 ** ($nlevels - 1) / 2

IO::out_every = 2   #TODO $ncells * 2 ** ($nlevels - 1) / 32

AMReX::regrid_error_threshold = 0.01   #TODO 0.1

AMReX::regrid_error_threshold = 0.01

ActiveThorns = "
    AMReX
    IOUtil
    WaveToyAMReX
"
 
$nlevels = 7
$ncells = 32
Cactus::cctk_show_schedule = no
# Cactus::terminate = "iteration"
# Cactus::cctk_itlast = 64   #TODO 0
Cactus::terminate = "time"
Cactus::cctk_final_time = 50.0
# AMReX::verbose = yes
AMReX::xmin = -100.0
AMReX::ymin = -100.0
AMReX::zmin = -100.0
AMReX::xmax = +100.0
AMReX::ymax = +100.0
AMReX::zmax = +100.0
AMReX::ncells_x = $ncells
AMReX::ncells_y = $ncells
AMReX::ncells_z = $ncells
# AMReX::periodic_x = no
# AMReX::periodic_y = no
# AMReX::periodic_z = no
AMReX::max_num_levels = $nlevels
AMReX::regrid_every = 16
AMReX::regrid_error_threshold = 0.1
AMReX::dtfac = 0.5
WaveToyAMReX::initial_condition = "central potential"
WaveToyAMReX::boundary_condition = "initial"
WaveToyAMReX::central_point_charge = 1.0
WaveToyAMReX::central_point_radius = 1.0
IO::out_dir = $parfile
IO::out_every = $ncells * 2 ** ($nlevels - 1) / 32

#include <driver.hxx>
#include <io.hxx>
#include <schedule.hxx>

#include "driver.hxx"
#include "io.hxx"
#include "schedule.hxx"

namespace {
// Convert a (direction, face) pair to an AMReX Orientation
Orientation orient(int d, int f) {
  return Orientation(d, Orientation::Side(f));
}
} // namespace

    const Box &fbx = mfi.fabbox();
    const Box &bx = mfi.tilebox();
    const Dim3 imin = lbound(bx);
    int ash[3], lsh[3];
    for (int d = 0; d < dim; ++d)
      ash[d] = fbx[orient(d, 1)] - fbx[orient(d, 0)] + 1;
    // Note: This excludes ghosts, it's not the proper Cactus lsh
    for (int d = 0; d < dim; ++d)
      lsh[d] = bx[orient(d, 1)] - bx[orient(d, 0)] + 1;

    GridPtrDesc grid(leveldata, mfi);

    const CCTK_REAL *restrict err = vars.ptr(imin.x, imin.y, imin.z, vi);
    const Array4<char> &tagarr = tags.array(mfi);
    constexpr int di = 1;
    const int dj = di * ash[0];
    const int dk = dj * ash[1];

    CCTK_REAL *restrict const err = grid.ptr(vars, vi);
    const Array4<char> &tag = tags.array(mfi);

    for (int k = 0; k < lsh[2]; ++k) {
      for (int j = 0; j < lsh[1]; ++j) {
#pragma omp simd
        for (int i = 0; i < lsh[0]; ++i) {
          int idx = di * i + dj * j + dk * k;
          tagarr(imin.x + i, imin.y + j, imin.z + k) =
              err[idx] >= regrid_error_threshold ? TagBox::SET : TagBox::CLEAR;
        }
      }
    }

    grid.loop_int([&](int i, int j, int k, int idx) {
      tag(grid.cactus_offset.x + i, grid.cactus_offset.y + j,
          grid.cactus_offset.z + k) =
          err[idx] >= regrid_error_threshold ? TagBox::SET : TagBox::CLEAR;
    });

    for (int tl = 0; tl < int(groupdata.mfab.size()); ++tl)

    for (int tl = 0; tl < int(groupdata.mfab.size()); ++tl) {

      // Set grid functions to nan
      MultiFab &mfab = *groupdata.mfab.at(tl);
      auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
          {max_tile_size_x, max_tile_size_y, max_tile_size_z});
#pragma omp parallel
      for (MFIter mfi(mfab, mfitinfo); mfi.isValid(); ++mfi) {
        GridPtrDesc grid(leveldata, mfi);
        const Array4<CCTK_REAL> &vars = mfab.array(mfi);
        for (int vi = 0; vi < groupdata.numvars; ++vi) {
          CCTK_REAL *restrict const ptr = grid.ptr(vars, vi);
          grid.loop_all(
              [&](int i, int j, int k, int idx) { ptr[idx] = 0.0 / 0.0; });
        }
      }
    }

        // Set grid functions to nan
        auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
            {max_tile_size_x, max_tile_size_y, max_tile_size_z});
#pragma omp parallel
        for (MFIter mfi(*mfab, mfitinfo); mfi.isValid(); ++mfi) {
          GridPtrDesc grid(leveldata, mfi);
          const Array4<CCTK_REAL> &vars = mfab->array(mfi);
          for (int vi = 0; vi < groupdata.numvars; ++vi) {
            CCTK_REAL *restrict const ptr = grid.ptr(vars, vi);
            grid.loop_all(
                [&](int i, int j, int k, int idx) { ptr[idx] = 0.0 / 0.0; });
          }
        }

        // Set grid functions to nan
        auto mfitinfo = MFItInfo().SetDynamic(true).EnableTiling(
            {max_tile_size_x, max_tile_size_y, max_tile_size_z});
#pragma omp parallel
        for (MFIter mfi(*mfab, mfitinfo); mfi.isValid(); ++mfi) {
          GridPtrDesc grid(leveldata, mfi);
          const Array4<CCTK_REAL> &vars = mfab->array(mfi);
          for (int vi = 0; vi < groupdata.numvars; ++vi) {
            CCTK_REAL *restrict const ptr = grid.ptr(vars, vi);
            grid.loop_all(
                [&](int i, int j, int k, int idx) { ptr[idx] = 0.0 / 0.0; });
          }
        }

  assert(level == int(ghext->leveldata.size()) - 1);

  // assert(level == int(ghext->leveldata.size()) - 1);

  pp.add("amrex.verbose", true);
  pp.add("amrex.signal_handling", false);

  // Don't catch Unix signals. If signals are caught, we don't get
  // core files.
  pp.add("amrex.signal_handling", 0);
  // Throw exceptions for failing AMReX assertions. With exceptions,
  // we get core files.
  pp.add("amrex.throw_exception", 1);

  const Array<int, dim> periodic{periodic_x, periodic_y, periodic_z};

  const Array<int, dim> is_periodic{periodic_x, periodic_y, periodic_z};

      domain, max_num_levels - 1, ncells, coord, reffacts, periodic);

      domain, max_num_levels - 1, ncells, coord, reffacts, is_periodic);

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

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

  // Should we really do this? Cactus's extension handling mechanism
  // becomes inconsistent once extensions have been unregistered.

#include <cctk.h>

#include <cctk.h>

#include <driver.hxx>
#include <io.hxx>
#include <reduction.hxx>

#include "driver.hxx"
#include "io.hxx"
#include "reduction.hxx"
#include "schedule.hxx"

namespace {
// Convert a (direction, face) pair to an AMReX Orientation
Orientation orient(int d, int f) {
  return Orientation(d, Orientation::Side(f));
}
} // namespace

                  const int *restrict ash, const int *restrict lbnd,
                  const int *restrict lsh, const int *restrict bbox,
                  const int *restrict nghostzones) {

                  const GridDesc &grid) {

    imin[d] = cctk_bbox[2 * d] ? 0 : cctk_nghostzones[d];
    imax[d] = cctk_lsh[d] - (cctk_bbox[2 * d + 1] ? 0 : cctk_nghostzones[d]);

    imin[d] = grid.bbox[2 * d] ? 0 : grid.nghostzones[d];
    imax[d] = grid.lsh[d] - (grid.bbox[2 * d + 1] ? 0 : grid.nghostzones[d]);

  const int dj = di * ash[0];
  const int dk = dj * ash[1];

  const int dj = di * grid.ash[0];
  const int dk = dj * grid.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) {

  for (int k = 0; k < grid.lsh[2]; ++k) {
    for (int j = 0; j < grid.lsh[1]; ++j) {
      for (int i = 0; i < grid.lsh[0]; ++i) {

        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];

        T x = x0[0] + (grid.lbnd[0] + i) * dx[0];
        T y = x0[1] + (grid.lbnd[1] + j) * dx[1];
        T z = x0[2] + (grid.lbnd[2] + k) * dx[2];

           << component << "\t" << isghost << "\t" << (lbnd[0] + i) << "\t"
           << (lbnd[1] + j) << "\t" << (lbnd[2] + k) << "\t" << x << "\t" << y
           << "\t" << z;

           << component << "\t" << isghost << "\t" << (grid.lbnd[0] + i) << "\t"
           << (grid.lbnd[1] + j) << "\t" << (grid.lbnd[2] + k) << "\t" << x
           << "\t" << y << "\t" << z;

      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];

      GridPtrDesc grid(leveldata, mfi);

        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);

        ptrs.at(vi) = grid.ptr(vars, vi);
      write_arrays(file, cctkGH, leveldata.level, mfi.index(), ptrs, grid);

#include "schedule.hxx"

#ifndef LOOP_HXX
#define LOOP_HXX
#include <cctk.h>
#include <array>
#include <iostream>
#include <string>
namespace Loop {
using namespace std;
constexpr int dim = 3;
struct GridDescBase {
  array<int, dim> gsh;
  array<int, dim> lbnd, ubnd;
  array<int, dim> lsh;
  array<int, dim> ash;
  array<int, 2 * dim> bbox;
  array<int, dim> nghostzones;
  array<int, dim> tmin, tmax;
  template <typename T, size_t N>
  static void output(ostream &os, const string &str, const array<T, N> &arr) {
    os << str << ":[";
    for (size_t n = 0; n < N; ++n) {
      if (n > 0)
        os << ",";
      os << arr[n];
    }
    os << "]";
  }
  friend ostream &operator<<(ostream &os, const GridDescBase &grid) {
    os << "GridDescBase{";
    output(os, "gsh", grid.gsh);
    output(os, ",lbnd", grid.lbnd);
    output(os, ",ubnd", grid.ubnd);
    output(os, ",lsh", grid.lsh);
    output(os, ",bbox", grid.bbox);
    output(os, ",nghostzones", grid.nghostzones);
    output(os, ",tmin", grid.tmin);
    output(os, ",tmax", grid.tmax);
    os << "}";
    return os;
  }
protected:
  GridDescBase();
public:
  GridDescBase(const cGH *cctkGH);
  // Loop over a given box
  template <typename F>
  void loop_box(const F &f, const array<int, dim> &restrict imin,
                const array<int, dim> &restrict imax) const {
    // cout << *this;
    // output(cout, ",imin", imin);
    // output(cout, ",imax", imax);
    // cout << "\n";
    for (int d = 0; d < dim; ++d)
      if (imin[d] >= imax[d])
        return;
    constexpr int di = 1;
    const int dj = di * ash[0];
    const int dk = dj * 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) {
          int idx = i * di + j * dj + k * dk;
          f(i, j, k, idx);
        }
      }
    }
  }
  // Loop over all points
  template <typename F> void loop_all(const F &f) const {
    array<int, dim> imin, imax;
    for (int d = 0; d < dim; ++d) {
      imin[d] = max(tmin[d], 0);
      imax[d] = min(tmax[d], lsh[d]);
    }
    loop_box(f, imin, imax);
  }
  // Loop over all interior points
  template <typename F> void loop_int(const F &f) const {
    array<int, dim> imin, imax;
    for (int d = 0; d < dim; ++d) {
      imin[d] = max(tmin[d], nghostzones[d]);
      imax[d] = min(tmax[d], lsh[d] - nghostzones[d]);
    }
    loop_box(f, imin, imax);
  }
  // Loop over all outer boundary points. This excludes ghost faces, but
  // includes ghost edges/corners on non-ghost faces.
  template <typename F> void loop_bnd(const F &f) const {
    for (int dir = 0; dir < dim; ++dir) {
      for (int face = 0; face < 2; ++face) {
        if (bbox[2 * dir + face]) {
          array<int, dim> imin, imax;
          for (int d = 0; d < dim; ++d) {
            // by default, include interior and outer boundaries and ghosts
            imin[d] = 0;
            imax[d] = lsh[d];
            // avoid covering edges and corners multiple times
            if (d < dir) {
              if (bbox[2 * d])
                imin[d] = nghostzones[d]; // only interior
              if (bbox[2 * d + 1])
                imax[d] = lsh[d] - nghostzones[d]; // only interior
            }
          }
          // only one face on outer boundary
          if (face == 0)
            imax[dir] = nghostzones[dir];
          else
            imin[dir] = lsh[dir] - nghostzones[dir];
          for (int d = 0; d < dim; ++d) {
            imin[d] = max(tmin[d], imin[d]);
            imax[d] = min(tmax[d], imax[d]);
          }
          loop_box(f, imin, imax);
        }
      }
    }
  }
};
template <typename F> void loop_all(const cGH *cctkGH, const F &f) {
  GridDescBase(cctkGH).loop_all(f);
}
template <typename F> void loop_int(const cGH *cctkGH, const F &f) {
  GridDescBase(cctkGH).loop_int(f);
}
template <typename F> void loop_bnd(const cGH *cctkGH, const F &f) {
  GridDescBase(cctkGH).loop_bnd(f);
}
} // namespace Loop
#endif // #ifndef LOOP_HXX

#include <driver.hxx>
#include <reduction.hxx>

#include "driver.hxx"
#include "reduction.hxx"
#include "schedule.hxx"

namespace {
// Convert a (direction, face) pair to an AMReX Orientation
Orientation orient(int d, int f) {
  return Orientation(d, Orientation::Side(f));
}
} // namespace

                          const int *restrict ash, const int *restrict lsh) {

                          const GridDesc &grid) {

  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) {
#pragma omp simd
      for (int i = 0; i < lsh[0]; ++i) {
        int idx = di * i + dj * j + dk * k;
        red = reduction<T>(red, reduction<T>(dV, ptr[idx]));
      }
    }
  }

  grid.loop_int([&](int i, int j, int k, int idx) {
    red = reduction<T>(red, reduction<T>(dV, ptr[idx]));
  });

        const Box &fbx = mfi.fabbox();
        const Box &bx = mfi.tilebox();

        GridPtrDesc grid(leveldata, mfi);

        const Dim3 imin = lbound(bx);
        int ash[3], lsh[3];
        for (int d = 0; d < dim; ++d)
          ash[d] = fbx[orient(d, 1)] - fbx[orient(d, 0)] + 1;
        // Note: This excludes ghosts, it's not the proper Cactus lsh
        for (int d = 0; d < dim; ++d)
          lsh[d] = bx[orient(d, 1)] - bx[orient(d, 0)] + 1;

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

        const CCTK_REAL *restrict ptr = grid.ptr(vars, vi);

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

        red += reduce_array(geom, ptr, grid);

  T fmax1(T x, T y) {
    if (CCTK_isnan(x))
      return x;
    if (CCTK_isnan(y))
      return y;
    return fmax(x, y);
  }
  T fmin1(T x, T y) {
    if (CCTK_isnan(x))
      return x;
    if (CCTK_isnan(y))
      return y;
    return fmin(x, y);
  }

  T sdv() const { return sqrt(fmax(T(0.0), vol * sum2 - sum * sum)); }

  T sdv() const { return sqrt(fmax1(T(0.0), vol * sum2 - pow(sum, 2))); }

    : min(x), max(x), sum(V * x), sum2((V * x) * (V * x)), vol(V),
      maxabs(fabs(x)), sumabs(fabs(V * x)), sum2abs(fabs(V * x) * fabs(V * x)) {
}

    : min(x), max(x), sum(V * x), sum2(pow(V * x, 2)), vol(V), maxabs(fabs(x)),
      sumabs(fabs(V * x)), sum2abs(pow(fabs(V * x), 2)) {}

    : min(fmin(x.min, y.min)), max(fmax(x.max, y.max)), sum(x.sum + y.sum),

    : min(fmin1(x.min, y.min)), max(fmax1(x.max, y.max)), sum(x.sum + y.sum),

      maxabs(fmax(x.maxabs, y.maxabs)), sumabs(x.sumabs + y.sumabs),

      maxabs(fmax1(x.maxabs, y.maxabs)), sumabs(x.sumabs + y.sumabs),

#include <driver.hxx>
#include <schedule.hxx>

#include "driver.hxx"
#include "schedule.hxx"

#include <memory>
#include <sstream>

// Tile boxes, should be part of cGH
struct TileBox {
  int tile_min[dim];
  int tile_max[dim];
};
vector<TileBox> thread_local_tilebox;

// Set up a GridDesc
} // namespace AMReX
namespace Loop {
GridDescBase::GridDescBase() {}
GridDescBase::GridDescBase(const cGH *restrict cctkGH) {
  for (int d = 0; d < dim; ++d)
    gsh[d] = cctkGH->cctk_gsh[d];
  for (int d = 0; d < dim; ++d) {
    lbnd[d] = cctkGH->cctk_lbnd[d];
    ubnd[d] = cctkGH->cctk_ubnd[d];
  }
  for (int d = 0; d < dim; ++d)
    lsh[d] = cctkGH->cctk_lsh[d];
  for (int d = 0; d < dim; ++d)
    ash[d] = cctkGH->cctk_ash[d];
  for (int d = 0; d < dim; ++d)
    for (int f = 0; f < 2; ++f)
      bbox[2 * d + f] = cctkGH->cctk_bbox[2 * d + f];
  for (int d = 0; d < dim; ++d)
    nghostzones[d] = cctkGH->cctk_nghostzones[d];
  // Check whether we are in local mode
  assert(cctkGH->cctk_bbox[0] != AMReX::undefined);
  int thread_num = omp_get_thread_num();
  const cGH *restrict threadGH = &AMReX::thread_local_cctkGH.at(thread_num);
  // Check whether this is the correct cGH structure
  assert(cctkGH == threadGH);
  const AMReX::TileBox &restrict tilebox =
      AMReX::thread_local_tilebox.at(thread_num);
  for (int d = 0; d < dim; ++d) {
    tmin[d] = tilebox.tile_min[d];
    tmax[d] = tilebox.tile_max[d];
  }
}
} // namespace Loop
namespace AMReX {
GridDesc::GridDesc(const GHExt::LevelData &leveldata, const MFIter &mfi) {
  const Box &fbx = mfi.fabbox();   // allocated array
  const Box &vbx = mfi.validbox(); // interior region (without ghosts)
  // const Box &bx = mfi.tilebox();       // current region (without ghosts)
  const Box &gbx = mfi.growntilebox(); // current region (with ghosts)
  const Box &domain = ghext->amrcore->Geom(leveldata.level).Domain();
  // The number of ghostzones in each direction
  for (int d = 0; d < dim; ++d)
    nghostzones[d] = mfi.theFabArrayBase().nGrowVect()[d];
  // Global shape
  for (int d = 0; d < dim; ++d)
    gsh[d] =
        domain[orient(d, 1)] - domain[orient(d, 0)] + 1 + 2 * nghostzones[d];
  // Local shape
  for (int d = 0; d < dim; ++d)
    lsh[d] = fbx[orient(d, 1)] - fbx[orient(d, 0)] + 1;
  // Allocated shape
  for (int d = 0; d < dim; ++d)
    ash[d] = fbx[orient(d, 1)] - fbx[orient(d, 0)] + 1;
  // Local extent
  for (int d = 0; d < dim; ++d) {
    lbnd[d] = fbx[orient(d, 0)] + nghostzones[d];
    ubnd[d] = fbx[orient(d, 1)] + nghostzones[d];
  }
  // Boundaries
  for (int d = 0; d < dim; ++d)
    for (int f = 0; f < 2; ++f)
      bbox[2 * d + f] = vbx[orient(d, f)] == domain[orient(d, f)];
  // Thread tile box
  for (int d = 0; d < dim; ++d) {
    tmin[d] = gbx[orient(d, 0)] - fbx[orient(d, 0)];
    tmax[d] = gbx[orient(d, 1)] + 1 - fbx[orient(d, 0)];
  }
  // Check constraints
  for (int d = 0; d < dim; ++d) {
    // Domain size
    assert(gsh[d] >= 0);
    // Local size
    assert(lbnd[d] >= 0);
    assert(lsh[d] >= 0);
    assert(lbnd[d] + lsh[d] <= gsh[d]);
    assert(ubnd[d] == lbnd[d] + lsh[d] - 1);

    // Internal representation
    assert(ash[d] >= 0);
    assert(ash[d] >= lsh[d]);
    // Ghost zones
    assert(nghostzones[d] >= 0);
    assert(2 * nghostzones[d] <= lsh[d]);
    // Tiles
    assert(tmin[d] >= 0);
    assert(tmax[d] >= tmin[d]);
    assert(tmax[d] <= lsh[d]);
  }
}
GridPtrDesc::GridPtrDesc(const GHExt::LevelData &leveldata, const MFIter &mfi)
    : GridDesc(leveldata, mfi) {
  const Box &fbx = mfi.fabbox(); // allocated array
  cactus_offset = lbound(fbx);
}

    cctkGH->cctk_origin_space[d] = x0[d] + 0.0 * dx[d];

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

    cctkGH->cctk_gsh[d] = domain[orient(d, 1)] - domain[orient(d, 0)] + 1;

    cctkGH->cctk_gsh[d] = domain[orient(d, 1)] - domain[orient(d, 0)] + 1 +
                          2 * cctkGH->cctk_nghostzones[d];

    cctkGH->cctk_levoff[d] = 1;

    cctkGH->cctk_levoff[d] = 1 - 2 * cctkGH->cctk_nghostzones[d];

void enter_local_mode(cGH *restrict cctkGH,
                      GHExt::LevelData &restrict leveldata, const MFIter &mfi) {
  const Box &fbx = mfi.fabbox(); // allocated array
  // const Box &vbx = mfi.validbox();     // interior region (without ghosts)
  const Box &bx = mfi.tilebox();       // current region (without ghosts)
  const Box &gbx = mfi.growntilebox(); // current region (with ghosts)

void enter_local_mode(cGH *restrict cctkGH, TileBox &restrict tilebox,
                      const GHExt::LevelData &restrict leveldata,
                      const MFIter &mfi) {
  GridPtrDesc grid(leveldata, mfi);

  // Local shape

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

    cctkGH->cctk_lsh[d] = grid.lsh[d];

    cctkGH->cctk_ash[d] = fbx[orient(d, 1)] - fbx[orient(d, 0)] + 1;
  // Local extent

    cctkGH->cctk_ash[d] = grid.ash[d];

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

    cctkGH->cctk_lbnd[d] = grid.lbnd[d];
    cctkGH->cctk_ubnd[d] = grid.ubnd[d];

  // Boundaries

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

      cctkGH->cctk_bbox[2 * d + f] = grid.bbox[2 * d + f];
  for (int d = 0; d < dim; ++d) {
    tilebox.tile_min[d] = grid.tmin[d];
    tilebox.tile_max[d] = grid.tmax[d];
  }

  const Dim3 imin = lbound(gbx);

        cctkGH->data[groupdata.firstvarindex + vi][tl] =
            vars.ptr(imin.x, imin.y, imin.z, vi);

        cctkGH->data[groupdata.firstvarindex + vi][tl] = grid.ptr(vars, vi);

  // Check constraints
  for (int d = 0; d < dim; ++d) {
    // Domain size
    assert(cctkGH->cctk_gsh[d] >= 0);
    // Local size
    assert(cctkGH->cctk_lbnd[d] >= 0);
    assert(cctkGH->cctk_lsh[d] >= 0);
    assert(cctkGH->cctk_lbnd[d] + cctkGH->cctk_lsh[d] <= cctkGH->cctk_gsh[d]);
    assert(cctkGH->cctk_ubnd[d] ==
           cctkGH->cctk_lbnd[d] + cctkGH->cctk_lsh[d] - 1);
    // Internal representation
    assert(cctkGH->cctk_ash[d] >= 0);
    assert(cctkGH->cctk_ash[d] >= cctkGH->cctk_lsh[d]);
    // Ghost zones
    assert(cctkGH->cctk_nghostzones[d] >= 0);
    assert(2 * cctkGH->cctk_nghostzones[d] <= cctkGH->cctk_lsh[d]);
    // Tiles
    assert(tilebox.tile_min[d] >= 0);
    assert(tilebox.tile_max[d] >= tilebox.tile_min[d]);
    assert(tilebox.tile_max[d] <= cctkGH->cctk_lsh[d]);
  }

void leave_local_mode(cGH *restrict cctkGH,

void leave_local_mode(cGH *restrict cctkGH, TileBox &restrict tilebox,

  for (int d = 0; d < dim; ++d) {
    tilebox.tile_min[d] = undefined;
    tilebox.tile_max[d] = undefined;
  }

extern "C" void AMReX_GetTileExtent(const void *restrict cctkGH_,
                                    CCTK_INT *restrict tile_min,
                                    CCTK_INT *restrict tile_max) {
  const cGH *restrict cctkGH = static_cast<const cGH *>(cctkGH_);
  // Check whether we are in local mode
  assert(cctkGH->cctk_bbox[0] != undefined);
  int thread_num = omp_get_thread_num();
  const cGH *restrict threadGH = &thread_local_cctkGH.at(thread_num);
  // Check whether this is the correct cGH structure
  assert(cctkGH == threadGH);

  const TileBox &restrict tilebox = thread_local_tilebox.at(thread_num);
  for (int d = 0; d < dim; ++d) {
    tile_min[d] = tilebox.tile_min[d];
    tile_max[d] = tilebox.tile_max[d];
  }
}

  // enter_level_mode(cctkGH);

  thread_local_tilebox.resize(max_threads);

    // enter_level_mode(threadGH);

      // TODO: Remove "current_level" mechanism
      // current_level = ghext->amrcore->finestLevel();

      CCTK_Traverse(cctkGH, "CCTK_POSTREGRIDINITIAL");

    // current_level = -1;

    if (regrid_every > 0 && cctkGH->cctk_iteration % regrid_every == 0) {

    if (regrid_every > 0 && cctkGH->cctk_iteration % regrid_every == 0 &&
        ghext->amrcore->maxLevel() > 0) {

      CCTK_Traverse(cctkGH, "CCTK_POSTREGRID");

    CCTK_VINFO("CallFunction [%d] %s: %s::%s", cctkGH->cctk_iteration,

    CCTK_VINFO("CallFunction iteration %d %s: %s::%s", cctkGH->cctk_iteration,

      cGH *restrict threadGH = &thread_local_cctkGH.at(omp_get_thread_num());

      int thread_num = omp_get_thread_num();
      cGH *restrict threadGH = &thread_local_cctkGH.at(thread_num);

      TileBox &restrict thread_tilebox = thread_local_tilebox.at(thread_num);

      auto callfunc = [&](auto &restrict leveldata) {

      // Loop over all levels
      // TODO: parallelize this loop
      for (auto &restrict leveldata : ghext->leveldata) {

          enter_local_mode(threadGH, leveldata, mfi);

          enter_local_mode(threadGH, thread_tilebox, leveldata, mfi);

          leave_local_mode(threadGH, leveldata, mfi);

          leave_local_mode(threadGH, thread_tilebox, leveldata, mfi);

      };
      // Loop over all levels
      // TODO: parallelize this loop
      for (auto &restrict leveldata : ghext->leveldata)
        callfunc(leveldata);

      }

  if (verbose) {
    ostringstream buf;
    for (int n = 0; n < numgroups; ++n) {
      if (n != 0)
        buf << ", ";
      buf << unique_ptr<char>(CCTK_GroupName(groups[n])).get();
    }
    CCTK_VINFO("SyncGroups %s", buf.str().c_str());
  }

        for (int tl = 0; tl < sync_tl; ++tl) {
#warning "TODO: make copy of fine level"

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

        }

#include "driver.hxx"
#include "loop.hxx"

#include <algorithm>
#include <array>

using namespace Loop;
using namespace std;
using std::max;
using std::min;

////////////////////////////////////////////////////////////////////////////////
struct GridDesc : GridDescBase {
  GridDesc() = delete;
  GridDesc(const GHExt::LevelData &leveldata, const MFIter &mfi);
  GridDesc(const cGH *cctkGH) : GridDescBase(cctkGH) {}
};
struct GridPtrDesc : GridDesc {
  Dim3 cactus_offset;
  GridPtrDesc() = delete;
  GridPtrDesc(const GHExt::LevelData &leveldata, const MFIter &mfi);
  template <typename T> T *ptr(const Array4<T> &vars, int vi) const {
    return vars.ptr(cactus_offset.x, cactus_offset.y, cactus_offset.z, vi);
  }
  template <typename T>
  T &idx(const Array4<T> &vars, int i, int j, int k, int vi) const {
    return vars(cactus_offset.x + i, cactus_offset.y + i, cactus_offset.z + j,
                vi);
  }
};

3. Scalar particles
Lagrangian:
  S = int 1/2 \eta^ab d_a u d_b u + 1/2 m^2 u^2 + \rho u
equation of motion:
  0 = dS/du = - \eta^ab d_a d_b u + m^2 u + \rho
nonsense: [0 = dS/d\rho = u]
nonsense: [0 = dS/dm = m u^2]
Ansatz:
  \rho = q \delta(x)
  0 = dS/dx = -q (grad u)(x)
  0 = dS/dq = u(x)

USES INCLUDE HEADER: loop.hxx

CCTK_REAL mass "Scalar field mass"
{
  0:* :: ""
} 0.0
CCTK_REAL central_point_charge "Scalar piont charge located at origin"
{
  *:* :: ""
} 0.0
CCTK_REAL central_point_radius "Radius of point charge"
{
  (*:* :: ""
} 1.0

  "Gaussian" :: ""
  "central potential" :: ""

KEYWORD boundary_condition "Boundary condition"
{
  "none" :: ""
  "initial" :: ""
} "none"

CCTK_REAL psi_abs "Typical value of psi to determine absolute error" STEERABLE=always
{
  (0.0:* :: ""
} 1.0

{
  LANG: C
} "Set up initial conditions for the wave equation"
SCHEDULE WaveToyAMReX_Sync AT initial AFTER WaveToyAMReX_Initialize

} "Set up initial conditions for the wave equation"

} "Boundary conditions for the wave equation"
SCHEDULE WaveToyAMReX_Boundaries AT initial AFTER WaveToyAMReX_Sync
{
  LANG: C
} "Boundary conditions for the wave equation"
SCHEDULE WaveToyAMReX_NaNCheck_current AT postinitial
{
  LANG: C
} "Check for nans in the state vector"

SCHEDULE WaveToyAMReX_Sync AT postregridinitial
{
  LANG: C
  SYNC: phi
} "Boundary conditions for the wave equation"
SCHEDULE WaveToyAMReX_Boundaries AT postregridinitial AFTER WaveToyAMReX_Sync
{
  LANG: C
} "Boundary conditions for the wave equation"

SCHEDULE WaveToyAMReX_NaNCheck_current AT postregridinitial AFTER WaveToyAMReX_Boundaries
{
  LANG: C
} "Check for nans in the state vector"
SCHEDULE WaveToyAMReX_Sync AT postregrid
{
  LANG: C
  SYNC: phi
} "Boundary conditions for the wave equation"
SCHEDULE WaveToyAMReX_Boundaries AT postregrid AFTER WaveToyAMReX_Sync
{
  LANG: C
} "Boundary conditions for the wave equation"
SCHEDULE WaveToyAMReX_NaNCheck_current AT postregrid AFTER WaveToyAMReX_Boundaries
{
  LANG: C
} "Check for nans in the state vector"
SCHEDULE WaveToyAMReX_NaNCheck_past AT prestep
{
  LANG: C
} "Check for nans in the state vector"

} "Evolve the wave equation"
SCHEDULE WaveToyAMReX_Sync AT evol AFTER WaveToyAMReX_Evolve
{
  LANG: C

} "Evolve the wave equation"

} "Boundary conditions for the wave equation"
SCHEDULE WaveToyAMReX_Boundaries AT evol AFTER WaveToyAMReX_Sync
{
  LANG: C
} "Boundary conditions for the wave equation"

SCHEDULE WaveToyAMReX_NaNCheck_current AT poststep
{
  LANG: C
} "Check for nans in the state vector"

#include <loop.hxx>

// Square
template <typename T> T sqr(T x) { return x * x; }

// Spline with compact support of radius 1 and volume 1
template <typename T> T spline(T r) {
  if (r >= 1.0)
    return 0.0;
  constexpr CCTK_REAL f =
      dim == 1 ? 1.0
               : dim == 2 ? 24.0 / 7.0 / M_PI : dim == 3 ? 4.0 / M_PI : -1;
  const T r2 = pow(r, 2);
  return f * (r <= 0.5 ? 1 - 2 * r2 : 2 + r * (-4 + 2 * r));
}
// The potential for the spline
template <typename T> T spline_potential(T r) {
  // \Laplace u = 4 \pi \rho
  if (r >= 1.0)
    return -1 / r;
  static_assert(dim == 3, "");
  const T r2 = pow(r, 2);
  return r <= 0.5
             ? -7 / T(3) + r2 * (8 / T(3) - 8 / T(5) * r2)
             : (1 / T(15) +
                r * (-40 / T(15) +
                     r2 * (80 / T(15) + r * (-80 / T(15) + r * 24 / T(15))))) /
                   r;
}

////////////////////////////////////////////////////////////////////////////////

  T omega = sqrt(sqr(kx) + sqr(ky) + sqr(kz));

  T omega = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2));

// Standing wave
template <typename T> T gaussian(T t, T x, T y, T z) {

// Periodic Gaussian
template <typename T> T periodic_gaussian(T t, T x, T y, T z) {

  T omega = sqrt(sqr(kx) + sqr(ky) + sqr(kz));
  return exp(-0.5 * sqr(sin(kx * x + ky * y + kz * z - omega * t) / width));

  T omega = sqrt(pow(kx, 2) + pow(ky, 2) + pow(kz, 2));
  return exp(-0.5 * pow(sin(kx * x + ky * y + kz * z - omega * t) / width, 2));
}
// Gaussian
template <typename T> T gaussian(T t, T x, T y, T z) {
  DECLARE_CCTK_PARAMETERS;
  // u(t,r) = (f(r-t) - f(r+t)) / r
  T r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2));
  auto f = [&](T x) { return exp(-0.5 * pow(x / width, 2)); };
  auto fx = [&](T x) { return -x / pow(width, 2) * f(x); };
  if (r < 1.0e-8)
    // Use L'Hôpital's rule for small r
    return fx(r - t) - fx(r + t);
  else
    return (f(r - t) - f(r + t)) / r;

// Central potential
template <typename T> T central_potential(T t, T x, T y, T z) {
  DECLARE_CCTK_PARAMETERS;
  T r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2));
  return central_point_charge * pow(central_point_radius, -dim) *
         spline_potential(r / central_point_radius);
}

  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) {
#pragma omp simd
        for (int i = imin[0]; i < imax[0]; ++i) {
          const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
          CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
          CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
          CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
          phi[idx] = standing(t, x, y, z);
          // phi_p[idx] = standing(t - dt, x, y, z);
          psi[idx] = timederiv(standing, dt)(t, x, y, z);
        }
      }
    }

    Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {
      CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
      CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
      CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
      phi[idx] = standing(t, x, y, z);
      // phi_p[idx] = standing(t - dt, x, y, z);
      psi[idx] = timederiv(standing, dt)(t, x, y, z);
    });

    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 x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
          CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
          CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
          phi[idx] = gaussian(t, x, y, z);
          // phi_p[idx] = gaussian(t - dt, x, y, z);
          psi[idx] = timederiv(gaussian, dt)(t, x, y, z);
        }
      }
    }

    Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {
      CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
      CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
      CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
      phi[idx] = periodic_gaussian(t, x, y, z);
      // phi_p[idx] = periodic_gaussian(t - dt, x, y, z);
      psi[idx] = timederiv(periodic_gaussian, dt)(t, x, y, z);
    });
  } else if (CCTK_EQUALS(initial_condition, "Gaussian")) {
    Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {
      CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
      CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
      CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
      phi[idx] = gaussian(t, x, y, z);
      // phi_p[idx] = gaussian(t - dt, x, y, z);
      psi[idx] = timederiv(gaussian, dt)(t, x, y, z);
    });
  } else if (CCTK_EQUALS(initial_condition, "central potential")) {
    Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {
      CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
      CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
      CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
      phi[idx] = central_potential(t, x, y, z);
      // phi_p[idx] = central_potential(t - dt, x, y, z);
      psi[idx] = timederiv(central_potential, dt)(t, x, y, z);
    });

  for (int k = 0; k < cctk_lsh[2]; ++k) {
    for (int j = 0; j < cctk_lsh[1]; ++j) {
#pragma omp simd
      for (int i = 0; i < cctk_lsh[0]; ++i) {
        const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
        CCTK_REAL base_phi = fabs(phi[idx]) + fabs(phi_abs);
        CCTK_REAL errx_phi =
            fabs(phi[idx - di] - 2 * phi[idx] + phi[idx + di]) / base_phi;
        CCTK_REAL erry_phi =
            fabs(phi[idx - dj] - 2 * phi[idx] + phi[idx + dj]) / base_phi;
        CCTK_REAL errz_phi =
            fabs(phi[idx - dk] - 2 * phi[idx] + phi[idx + dk]) / base_phi;
        regrid_error[idx] = errx_phi + erry_phi + errz_phi;
      }
    }

  Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {
    CCTK_REAL base_phi = fabs(phi[idx]) + fabs(phi_abs);
    CCTK_REAL errx_phi =
        fabs(phi[idx - di] - 2 * phi[idx] + phi[idx + di]) / base_phi;
    CCTK_REAL erry_phi =
        fabs(phi[idx - dj] - 2 * phi[idx] + phi[idx + dj]) / base_phi;
    CCTK_REAL errz_phi =
        fabs(phi[idx - dk] - 2 * phi[idx] + phi[idx + dk]) / base_phi;
    CCTK_REAL base_psi = fabs(psi[idx]) + fabs(psi_abs);
    CCTK_REAL errx_psi =
        fabs(psi[idx - di] - 2 * psi[idx] + psi[idx + di]) / base_psi;
    CCTK_REAL erry_psi =
        fabs(psi[idx - dj] - 2 * psi[idx] + psi[idx + dj]) / base_psi;
    CCTK_REAL errz_psi =
        fabs(psi[idx - dk] - 2 * psi[idx] + psi[idx + dk]) / base_psi;
    regrid_error[idx] =
        errx_phi + erry_phi + errz_phi + errx_psi + erry_psi + errz_psi;
  });
}
extern "C" void WaveToyAMReX_Evolve(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  const CCTK_REAL dt = CCTK_DELTA_TIME;
  CCTK_REAL x0[dim], dx[dim];
  for (int d = 0; d < dim; ++d) {
    x0[d] = CCTK_ORIGIN_SPACE(d);
    dx[d] = CCTK_DELTA_SPACE(d);

  constexpr int di = 1;
  const int dj = di * cctk_ash[0];
  const int dk = dj * cctk_ash[1];
  Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {
    CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
    CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
    CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
    CCTK_REAL r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2));
    CCTK_REAL ddx_phi =
        (phi_p[idx - di] - 2 * phi_p[idx] + phi_p[idx + di]) / pow(dx[0], 2);
    CCTK_REAL ddy_phi =
        (phi_p[idx - dj] - 2 * phi_p[idx] + phi_p[idx + dj]) / pow(dx[1], 2);
    CCTK_REAL ddz_phi =
        (phi_p[idx - dk] - 2 * phi_p[idx] + phi_p[idx + dk]) / pow(dx[2], 2);
    // phi[idx] = 2 * phi_p[idx] - phi_p_p[idx] +
    //            pow(dt, 2) * (ddx_phi + ddy_phi + ddz_phi);
    psi[idx] = psi_p[idx] +
               dt * (ddx_phi + ddy_phi + ddz_phi - pow(mass, 2) * phi_p[idx] -
                     4 * M_PI * central_point_charge *
                         pow(central_point_radius, -dim) *
                         spline(r / central_point_radius));
    phi[idx] = phi_p[idx] + dt * psi[idx];
  });

extern "C" void WaveToyAMReX_Evolve(CCTK_ARGUMENTS) {

// Miguel Alcubierre, Gabrielle Allen, Gerd Lanfermann
//
// Routines for applying radiation boundary conditions
//
// Taken from
// <Cactus/arrangements/CactusBase/Boundary/src/RadiationBoundary.c>
//
// The radiative boundary condition that is implemented is:
//
//   f  =  f0  +  u(r - v*t) / r  +  h(r + v*t) / r
//
// That is, I assume outgoing radial waves with a 1/r fall off, and
// the correct asymptotic value f0, plus I include the possibility of
// incoming waves as well (these incoming waves should be modeled
// somehow).
//
// The condition above leads to the differential equation:
//
//   (x / r) d f  +  v d f  + v x (f - f0) / r^2  =  v x H / r^2
//     i      t         i        i                      i
//
// where x_i is the normal direction to the given boundaries, and H =
// 2 dh(s)/ds.
//
// So at a given boundary I only worry about derivatives in the normal
// direction. Notice that u(r-v*t) has dissapeared, but we still do
// not know the value of H.
//
// To get H what I do is the following: I evaluate the expression one
// point in from the boundary and solve for H there. We now need a way
// of extrapolation H to the boundary. For this I assume that H falls
// off as a power law:
//
//   H = k/r**n  =>  d H  =  - n H/r
//                    i
//
// The value of n is is defined by the parameter "radpower". If this
// parameter is negative, H is forced to be zero (this corresponds to
// pure outgoing waves and is the default).
//
// The behaviour I have observed is the following: Using H=0 is very
// stable, but has a very bad initial transient. Taking n to be 0 or
// positive improves the initial behaviour considerably, but
// introduces a drift that can kill the evolution at very late times.
// Empirically, the best value I have found is n=2, for which the
// initial behaviour is very nice, and the late time drift is quite
// small.
//
// Another problem with this condition is that it does not use the
// physical characteristic speed, but rather it assumes a wave speed
// of v, so the boundaries should be out in the region where the
// characteristic speed is constant. Notice that this speed does not
// have to be 1. For gauge quantities {alpha, phi, trK} we can have a
// different asymptotic speed, which is why the value of v is passed
// as a parameter.
#if 0
extern "C" void WaveToyAMReX_RadiativeBoundaries(CCTK_ARGUMENTS) {

  constexpr CCTK_REAL phi_inf = 0; // value at infinity
  constexpr CCTK_REAL phi_vel = 1; // propagation speed
  // constexpr CCTK_REAL psi_inf = 0;
  // constexpr CCTK_REAL psi_vel = 1;

  CCTK_REAL dx[dim];

  CCTK_REAL x0[dim], dx[dim];

    x0[d] = CCTK_ORIGIN_SPACE(d);

  // CCTK_VINFO("imin[%d,%d,%d]", imin[0], imin[1], imin[2]);
  // CCTK_VINFO("imax[%d,%d,%d]", imax[0], imax[1], imax[2]);

  int gbox[2 * dim]; // Are the boundaries (if any) ghosts zones?
  for (int d = 0; d < dim; ++d) {
    // #error "NEED lbnd>=0 ALWAYS"
    gbox[2 * d] = cctk_lbnd[d] > 0;
    gbox[2 * d + 1] = cctk_lbnd[d] + cctk_lsh[d] < cctk_gsh[d];
  }
  // CCTK_VINFO("lbnd[%d,%d,%d]", cctk_lbnd[0], cctk_lbnd[1], cctk_lbnd[2]);
  // CCTK_VINFO("lsh[%d,%d,%d]", cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]);
  // CCTK_VINFO("gsh[%d,%d,%d]", cctk_gsh[0], cctk_gsh[1], cctk_gsh[2]);
  // CCTK_VINFO("gbox[%d,%d,%d,%d,%d,%d]", gbox[0], gbox[1], gbox[2], gbox[3],
  //            gbox[4], gbox[5]);

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

  for (int dir = 0; dir < dim; ++dir) {
    for (int face = 0; face < 2; ++face) {
      if (!gbox[2 * dir + face]) {
        int bmin[dim], bmax[dim];
        for (int d = 0; d < dim; ++d) {
          // Skip edges and corners if d < dir
          bmin[d] = d < dir ? imin[d] : 0;
          bmax[d] = d < dir ? imax[d] : cctk_lsh[d];
        }
        if (face == 0)
          bmax[dir] = imin[dir];
        else
          bmin[dir] = imax[dir];
        // CCTK_VINFO("dir=%d face=%d", dir, face);
        // CCTK_VINFO("bmin[%d,%d,%d]", bmin[0], bmin[1], bmin[2]);
        // CCTK_VINFO("bmax[%d,%d,%d]", bmax[0], bmax[1], bmax[2]);
        for (int k = bmin[2]; k < bmax[2]; ++k) {
          for (int j = bmin[1]; j < bmax[1]; ++j) {

      for (int i = imin[0]; i < imax[0]; ++i) {
        const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
        CCTK_REAL ddx_phi =
            (phi_p[idx - di] - 2 * phi_p[idx] + phi_p[idx + di]) / sqr(dx[0]);
        CCTK_REAL ddy_phi =
            (phi_p[idx - dj] - 2 * phi_p[idx] + phi_p[idx + dj]) / sqr(dx[1]);
        CCTK_REAL ddz_phi =
            (phi_p[idx - dk] - 2 * phi_p[idx] + phi_p[idx + dk]) / sqr(dx[2]);
        // phi[idx] = 2 * phi_p[idx] - phi_p_p[idx] +
        //            sqr(dt) * (ddx_phi + ddy_phi + ddz_phi);
        psi[idx] = psi_p[idx] + dt * (ddx_phi + ddy_phi + ddz_phi);
        phi[idx] = phi_p[idx] + dt * psi[idx];

            for (int i = bmin[0]; i < bmax[0]; ++i) {
              const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
              CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
              CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
              CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
              CCTK_REAL r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2));
              int ddir = dir == 0 ? di : dir == 1 ? dj : dk;
              int ndir = 2 * face - 1;
              // CCTK_REAL gradphi =
              //     (phi_p[idx] - phi_p[idx - ndir * ddir]) / dx[dir];
              // phi[idx] =
              //     phi_p[idx] -
              //     dt * phi_vel * (r * gradphi + (phi_p[idx] - phi_inf) / r);
              // CCTK_REAL gradpsi =
              //     (psi_p[idx] - psi_p[idx - ndir * ddir]) / dx[dir];
              // psi[idx] =
              //     psi_p[idx] -
              //     dt * psi_vel * (r * gradpsi + (psi_p[idx] - psi_inf) / r);
              CCTK_REAL gradphi =
                  (phi_p[idx] - phi_p[idx - ndir * ddir]) / dx[dir];
              psi[idx] = -phi_vel * (r * gradphi + (phi_p[idx] - phi_inf) / r);
              phi[idx] = phi_p[idx] + dt * psi[idx];
            }
          }
        }

  }
}
#endif
extern "C" void WaveToyAMReX_Sync(CCTK_ARGUMENTS) {
  // do nothinga
}
extern "C" void WaveToyAMReX_Boundaries(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  const CCTK_REAL t = cctk_time;
  const CCTK_REAL dt = CCTK_DELTA_TIME;
  CCTK_REAL x0[dim], dx[dim];
  for (int d = 0; d < dim; ++d) {
    x0[d] = CCTK_ORIGIN_SPACE(d);
    dx[d] = CCTK_DELTA_SPACE(d);
  }
  if (CCTK_EQUALS(boundary_condition, "none")) {
    // do nothing
  } else if (CCTK_EQUALS(boundary_condition, "initial")) {
    if (CCTK_EQUALS(initial_condition, "central potential")) {
      Loop::loop_bnd(cctkGH, [&](int i, int j, int k, int idx) {
        CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
        CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
        CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
        phi[idx] = central_potential(t, x, y, z);
        psi[idx] = timederiv(central_potential, dt)(t, x, y, z);
      });
    } else {
      assert(0);
    }
  } else {
    assert(0);
  }
}
extern "C" void WaveToyAMReX_NaNCheck_past(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {
    assert(CCTK_isfinite(phi_p[idx]));
    assert(CCTK_isfinite(psi_p[idx]));
  });
}
extern "C" void WaveToyAMReX_NaNCheck_current(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  int levfac = cctk_levfac[0];
  int lev = 0;
  while (levfac > 1) {
    levfac >>= 1;
    lev += 1;

  bool has_error = false;
  Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {
    if (!CCTK_isfinite(phi[idx]) || !CCTK_isfinite(psi[idx])) {
      CCTK_VINFO("level %d idx=[%d,%d,%d] gidx=[%d,%d,%d] (%g,%g)", lev, i, j,
                 k, cctk_lbnd[0] + i, cctk_lbnd[1] + j, cctk_lbnd[2] + k,
                 phi[idx], psi[idx]);
      has_error = true;
    }
  });
  if (has_error)
    CCTK_VERROR("NaNs found");

  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) {
#pragma omp simd
      for (int i = imin[0]; i < imax[0]; ++i) {
        const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
        CCTK_REAL ddx_phi =
            (phi[idx - di] - 2 * phi[idx] + phi[idx + di]) / sqr(dx[0]);
        CCTK_REAL ddy_phi =
            (phi[idx - dj] - 2 * phi[idx] + phi[idx + dj]) / sqr(dx[1]);
        CCTK_REAL ddz_phi =
            (phi[idx - dk] - 2 * phi[idx] + phi[idx + dk]) / sqr(dx[2]);
        CCTK_REAL psi_n = psi[idx] + dt * (ddx_phi + ddy_phi + ddz_phi);
        CCTK_REAL dt_phi_p = psi_n;
        CCTK_REAL dt_phi_m = psi_p[idx];
        CCTK_REAL dx_phi_p = (phi[idx + di] - phi[idx]) / dx[0];
        CCTK_REAL dx_phi_m = (phi[idx] - phi[idx - di]) / dx[0];
        CCTK_REAL dy_phi_p = (phi[idx + dj] - phi[idx]) / dx[1];
        CCTK_REAL dy_phi_m = (phi[idx] - phi[idx - dj]) / dx[1];
        CCTK_REAL dz_phi_p = (phi[idx + dk] - phi[idx]) / dx[2];
        CCTK_REAL dz_phi_m = (phi[idx] - phi[idx - dk]) / dx[2];
        eps[idx] =
            (sqr(dt_phi_p) + sqr(dt_phi_m) + sqr(dx_phi_p) + sqr(dx_phi_m) +
             sqr(dy_phi_p) + sqr(dy_phi_m) + sqr(dz_phi_p) + sqr(dz_phi_m)) /
            4;
      }
    }
  }

  Loop::loop_int(cctkGH, [&](int i, int j, int k, int idx) {
    CCTK_REAL ddx_phi =
        (phi[idx - di] - 2 * phi[idx] + phi[idx + di]) / pow(dx[0], 2);
    CCTK_REAL ddy_phi =
        (phi[idx - dj] - 2 * phi[idx] + phi[idx + dj]) / pow(dx[1], 2);
    CCTK_REAL ddz_phi =
        (phi[idx - dk] - 2 * phi[idx] + phi[idx + dk]) / pow(dx[2], 2);
    CCTK_REAL psi_n = psi[idx] + dt * (ddx_phi + ddy_phi + ddz_phi);
    CCTK_REAL dt_phi_p = psi_n;
    CCTK_REAL dt_phi_m = psi_p[idx];
    CCTK_REAL dx_phi_p = (phi[idx + di] - phi[idx]) / dx[0];
    CCTK_REAL dx_phi_m = (phi[idx] - phi[idx - di]) / dx[0];
    CCTK_REAL dy_phi_p = (phi[idx + dj] - phi[idx]) / dx[1];
    CCTK_REAL dy_phi_m = (phi[idx] - phi[idx - dj]) / dx[1];
    CCTK_REAL dz_phi_p = (phi[idx + dk] - phi[idx]) / dx[2];
    CCTK_REAL dz_phi_m = (phi[idx] - phi[idx - dk]) / dx[2];
    eps[idx] = (pow(dt_phi_p, 2) + pow(dt_phi_m, 2) + pow(dx_phi_p, 2) +
                pow(dx_phi_m, 2) + pow(dy_phi_p, 2) + pow(dy_phi_m, 2) +
                pow(dz_phi_p, 2) + pow(dz_phi_m, 2)) /
               4;
  });

    for (int k = 0; k < cctk_lsh[2]; ++k) {
      for (int j = 0; j < cctk_lsh[1]; ++j) {
#pragma omp simd
        for (int i = 0; i < cctk_lsh[0]; ++i) {
          const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
          CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
          CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
          CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
          phierr[idx] = phi[idx] - standing(t, x, y, z);
          psierr[idx] = psi[idx] - timederiv(standing, dt)(t, x, y, z);
        }
      }
    }

    Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {
      CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
      CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
      CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
      phierr[idx] = phi[idx] - standing(t, x, y, z);
      psierr[idx] = psi[idx] - timederiv(standing, dt)(t, x, y, z);
    });

    for (int k = 0; k < cctk_lsh[2]; ++k) {
      for (int j = 0; j < cctk_lsh[1]; ++j) {
#pragma omp simd
        for (int i = 0; i < cctk_lsh[0]; ++i) {
          const int idx = CCTK_GFINDEX3D(cctkGH, i, j, k);
          CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
          CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
          CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
          phierr[idx] = phi[idx] - gaussian(t, x, y, z);
          psierr[idx] = psi[idx] - timederiv(gaussian, dt)(t, x, y, z);
        }
      }
    }

    Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {
      CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
      CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
      CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
      phierr[idx] = phi[idx] - periodic_gaussian(t, x, y, z);
      psierr[idx] = psi[idx] - timederiv(periodic_gaussian, dt)(t, x, y, z);
    });
  } else if (CCTK_EQUALS(initial_condition, "Gaussian")) {
    Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {
      CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
      CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
      CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
      phierr[idx] = phi[idx] - gaussian(t, x, y, z);
      psierr[idx] = psi[idx] - timederiv(gaussian, dt)(t, x, y, z);
    });
  } else if (CCTK_EQUALS(initial_condition, "central potential")) {
    Loop::loop_all(cctkGH, [&](int i, int j, int k, int idx) {
      CCTK_REAL x = x0[0] + (cctk_lbnd[0] + i) * dx[0];
      CCTK_REAL y = x0[1] + (cctk_lbnd[1] + j) * dx[1];
      CCTK_REAL z = x0[2] + (cctk_lbnd[2] + k) * dx[2];
      phierr[idx] = phi[idx] - central_potential(t, x, y, z);
      psierr[idx] = psi[idx] - timederiv(central_potential, dt)(t, x, y, z);
    });
  } else {
    assert(0);