REQUIRES Arith

USES INCLUDE HEADER: vect.hxx

CCTK_REAL rhs TYPE=gf TAGS='index={0 0 0}' "rhs"

CCTK_REAL rhs TYPE=gf TAGS='index={0 0 0}' "Right hand side"

CCTK_REAL phi TYPE=gf TAGS='index={0 0 0}' "phi"
CCTK_REAL res TYPE=gf TAGS='index={0 0 0}' "res"

CCTK_REAL phi TYPE=gf TAGS='index={0 0 0}' "Potential"
CCTK_REAL res TYPE=gf TAGS='index={0 0 0}' "Residual (as reported by AMReX)"

CCTK_REAL ires TYPE=gf TAGS='index={0 0 0}' "ires"

CCTK_REAL ires TYPE=gf TAGS='index={0 0 0}' "Independently evaluated residual"

CCTK_REAL asol TYPE=gf TAGS='index={0 0 0}' "asol"
CCTK_REAL aerr TYPE=gf TAGS='index={0 0 0}' "aerr"

CCTK_REAL asol TYPE=gf TAGS='index={0 0 0}' "Analytic solution"
CCTK_REAL aerr TYPE=gf TAGS='index={0 0 0}' "Error wrt. analytic solution"

$nlevels = 1
$ncells = 8

$nlevels = 8
$ncells = 128

CarpetX::out_tsv = yes

CarpetX::out_tsv = no

# SCHEDULE Poisson_fixup AT initial AFTER  Poisson_solve
# {
#   LANG: C
#   WRITES: phi(boundary)
# } "Apply boundary conditions again"

SCHEDULE Poisson_residual AT initial AFTER (Poisson_solve, Poisson_fixup)

SCHEDULE Poisson_residual AT initial AFTER Poisson_solve

SCHEDULE Poisson_error AT initial AFTER (Poisson_solve, Poisson_fixup)

SCHEDULE Poisson_error AT initial AFTER Poisson_solve

#include <vect.hxx>

// AMReX uses a weird Laplace stencil. It takes a second derivative in
// one direction, and smoothes with weights [1, 4, 1] / 6 in the other
// directions.
template <typename T>
T laplace(const GF3D<const T, 0, 0, 0> &u, const PointDesc &p) {
  T r = 0;
  for (int dir = 0; dir < dim; ++dir) {
    vect<T, dim> lap{1, -2, 1};
    lap /= pow(p.DX[dir], 2);
    vect<T, dim> smo{1, 4, 1};
    smo /= 6;
    for (int k = -1; k <= 1; ++k) {
      for (int j = -1; j <= 1; ++j) {
        for (int i = -1; i <= 1; ++i) {
          const vect<int, dim> di{i, j, k};
          T w = 1;
          for (int d = 0; d < dim; ++d)
            if (d == dir)
              w *= lap[1 + di[d]];
            else
              w *= smo[1 + di[d]];
          auto I = p.I;
          for (int d = 0; d < dim; ++d)
            I += di[d] * p.DI(d);
          r += w * u(I);
        }
      }
    }
  }
  return r;
}

// extern "C" void Poisson_fixup(CCTK_ARGUMENTS) {
//   DECLARE_CCTK_ARGUMENTS_Poisson_fixup;
//   DECLARE_CCTK_PARAMETERS;
//
//   const GF3D<CCTK_REAL, 0, 0, 0> phi_(cctkGH, phi);
//
//   // Set boundary conditions
//   loop_bnd<0, 0, 0>(
//       cctkGH, [&](const PointDesc &p) { phi_(p.I) = fbnd(p.x, p.y, p.z); });
// }

      ires_(p.I) = (phi_(p.I - p.DI(0)) - 2 * phi_(p.I) + phi_(p.I + p.DI(0))) /
                       pow(p.DX[0], 2) +
                   (phi_(p.I - p.DI(1)) - 2 * phi_(p.I) + phi_(p.I + p.DI(1))) /
                       pow(p.DX[1], 2) +
                   (phi_(p.I - p.DI(2)) - 2 * phi_(p.I) + phi_(p.I + p.DI(2))) /
                       pow(p.DX[2], 2) -
                   frhs(p.x, p.y, p.z);

      ires_(p.I) = laplace(phi_, p) - frhs(p.x, p.y, p.z);

  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) {
    // ires_(p.I) = phi_(p.I) - fbnd(p.x, p.y, p.z);
    ires_(p.I) = 0.0;
  });

  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { ires_(p.I) = 0.0; });