Cactus Code Thorn Maxwell
Author(s)    : Erik Schnetter <schnetter@gmail.com>
Maintainer(s): Erik Schnetter <schnetter@gmail.com>
Licence      : LGPL
--------------------------------------------------------------------------
1. Purpose
Solve the Maxwell equations
2. Formulation (see <https://en.wikipedia.org/wiki/Maxwell%27s_equations>)
Electric potential: phi
Magnetic potential: A
Lorenz gauge:       div A = - dt phi
Electric field: E = - grad phi - dt A
Magnetic field: B = curl A
Homogeneous electric equation: curl E + dt B = 0
Homogeneous magnetic equation: div B = 0
Inhomogeneous electric equation: div E = rho
Inhomogeneous magnetic equation: curl B - dt E = J
Charge density:  rho
Current density: J
Charge conservation: dt rho + div J = 0
State vector:        phi
                     A
                     E
                     B
Equations of motion: dt phi = - div A
                     dt A   = - grad phi - E
                     dt E   = curl B - J
                     dt B   = - curl E
Constraints:         div E = rho
                     div B = 0
                     B     = curl A
Ansatz for initial conditions:
phi = aphi cos (omega t - k x)
A = aA cos (omega t - k x)
Lorenz gauge: omega = aA k / aphi
zero charge: aphi = aA k omega / k^2

# Configuration definitions for thorn Maxwell
REQUIRES Vectors

# Interface definition for thorn Maxwell
IMPLEMENTS: Maxwell
USES INCLUDE HEADER: loop.hxx
USES INCLUDE HEADER: vectors.h
void FUNCTION GetTileExtent(
  CCTK_POINTER_TO_CONST IN cctkGH,
  CCTK_INT ARRAY OUT tile_min,
  CCTK_INT ARRAY OUT tile_max)
REQUIRES FUNCTION GetTileExtent
# State vector
# TODO: phi needs no ghost zones
CCTK_REAL phi TYPE=gf TAGS='index={0 0 0} rhs="dtphi"' "Electric potential"
CCTK_REAL ax TYPE=gf TAGS='index={1 0 0} rhs="dtax"' "Magnetic potential"
CCTK_REAL ay TYPE=gf TAGS='index={0 1 0} rhs="dtay"' "Magnetic potential"
CCTK_REAL az TYPE=gf TAGS='index={0 0 1} rhs="dtaz"' "Magnetic potential"
CCTK_REAL ex TYPE=gf TAGS='index={1 0 0} rhs="dtex"' "Electric field"
CCTK_REAL ey TYPE=gf TAGS='index={0 1 0} rhs="dtey"' "Electric field"
CCTK_REAL ez TYPE=gf TAGS='index={0 0 1} rhs="dtez"' "Electric field"
CCTK_REAL byz TYPE=gf TAGS='index={0 1 1} rhs="dtbyz"' "Magnetic field"
CCTK_REAL bzx TYPE=gf TAGS='index={1 0 1} rhs="dtbzx"' "Magnetic field"
CCTK_REAL bxy TYPE=gf TAGS='index={1 1 0} rhs="dtbxy"' "Magnetic field"
# RHS
CCTK_REAL dtphi TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "Electric potential"
CCTK_REAL dtax TYPE=gf TAGS='index={1 0 0} checkpoint="no"' "Magnetic potential"
CCTK_REAL dtay TYPE=gf TAGS='index={0 1 0} checkpoint="no"' "Magnetic potential"
CCTK_REAL dtaz TYPE=gf TAGS='index={0 0 1} checkpoint="no"' "Magnetic potential"
CCTK_REAL dtex TYPE=gf TAGS='index={1 0 0} checkpoint="no"' "Electric field"
CCTK_REAL dtey TYPE=gf TAGS='index={0 1 0} checkpoint="no"' "Electric field"
CCTK_REAL dtez TYPE=gf TAGS='index={0 0 1} checkpoint="no"' "Electric field"
CCTK_REAL dtbyz TYPE=gf TAGS='index={0 1 1} checkpoint="no"' "Magnetic field"
CCTK_REAL dtbzx TYPE=gf TAGS='index={1 0 1} checkpoint="no"' "Magnetic field"
CCTK_REAL dtbxy TYPE=gf TAGS='index={1 1 0} checkpoint="no"' "Magnetic field"
# Constraints
CCTK_REAL curlax TYPE=gf TAGS='index={0 1 1} checkpoint="no"' "Magnetic potential constraint"
CCTK_REAL curlay TYPE=gf TAGS='index={1 0 1} checkpoint="no"' "Magnetic potential constraint"
CCTK_REAL curlaz TYPE=gf TAGS='index={1 1 0} checkpoint="no"' "Magnetic potential constraint"
CCTK_REAL dive TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "Electric constraint"
CCTK_REAL divb TYPE=gf TAGS='index={1 1 1} checkpoint="no"' "Magnetic constraint"

ActiveThorns = "
    CarpetX
    Formaline
    IOUtil
    Maxwell
    ODESolvers
    SystemTopology
    TimerReport
"
 
$nlevels = 1   #TODO 3
$ncells = 32   #TODO 32
Cactus::cctk_show_schedule = no
# Cactus::terminate = "iteration"
# Cactus::cctk_itlast = 0
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = yes   #TODO no
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 = $ncells
CarpetX::ncells_y = $ncells
CarpetX::ncells_z = $ncells
CarpetX::max_num_levels = $nlevels
CarpetX::regrid_every = 16
CarpetX::regrid_error_threshold = 0.01
CarpetX::prolongation_type = "conservative"
CarpetX::do_reflux = no
CarpetX::do_restrict = no
CarpetX::restrict_during_sync = yes
CarpetX::dtfac = 0.25
Maxwell::setup = "plane wave"   #TODO "Gaussian wave"
Maxwell::spatial_frequency_x = 1.0
Maxwell::spatial_frequency_y = 0.0
Maxwell::spatial_frequency_z = 0.0
Maxwell::amplitude_x = 1.0
Maxwell::amplitude_y = 0.0
Maxwell::amplitude_z = 0.0
Maxwell::width = 0.1
Maxwell::output_every = 1
IO::out_dir = $parfile
IO::out_every = $ncells * 2 ** ($nlevels - 1) / 32
CarpetX::out_tsv = no
TimerReport::out_every = $ncells * 2 ** ($nlevels - 1) / 32
TimerReport::out_filename = "TimerReport"
TimerReport::output_all_timers_together = yes
TimerReport::output_all_timers_readable = yes
TimerReport::n_top_timers = 50

# Parameter definitions for thorn Maxwell
KEYWORD setup "Initial condition"
{
  "plane wave" :: ""
  "Gaussian wave" :: ""
} "Gaussian wave"
CCTK_REAL spatial_frequency_x "spatial frequency"
{
  *:* :: ""
} 0.5
CCTK_REAL spatial_frequency_y "spatial frequency"
{
  *:* :: ""
} 0.5
CCTK_REAL spatial_frequency_z "spatial frequency"
{
  *:* :: ""
} 0.5
CCTK_REAL amplitude_x "amplitude"
{
  *:* :: ""
} 1.0
CCTK_REAL amplitude_y "amplitude"
{
  *:* :: ""
} 1.0
CCTK_REAL amplitude_z "amplitude"
{
  *:* :: ""
} 1.0
CCTK_REAL width "width"
{
  (0.0:* :: ""
} 0.1
CCTK_INT output_every "Calculate output quantities every that many iterations" STEERABLE=recover
{
  0:* :: ""
} 1

# Schedule definitions for thorn Maxwell
STORAGE: phi
STORAGE: ax ay az
STORAGE: ex ey ez
STORAGE: byz bzx bxy
STORAGE: dtphi
STORAGE: dtax dtay dtaz
STORAGE: dtex dtey dtez
STORAGE: dtbyz dtbzx dtbxy
STORAGE: dive
STORAGE: divb
# Initial conditions
SCHEDULE Maxwell_Initial AT initial
{
  LANG: C
  WRITES: phi(interior)
  WRITES: ax(interior) ay(interior) az(interior)
  WRITES: ex(interior) ey(interior) ez(interior)
  WRITES: byz(interior) bzx(interior) bxy(interior)
  SYNC: phi
  SYNC: ax ay az
  SYNC: ex ey ez
  SYNC: byz bzx bxy
} "Set up hydro initial conditions"
# Regridding
SCHEDULE Maxwell_Boundaries AT postregrid
{
  LANG: C
  SYNC: phi
  SYNC: ax ay az
  SYNC: ex ey ez
  SYNC: byz bzx bxy
} "Apply boundary conditions"
# Time evolution
SCHEDULE Maxwell_RHS IN ODESolvers_RHS
{
  LANG: C
  READS: phi(interior)
  READS: ax(everywhere) ay(everywhere) az(everywhere)
  READS: ex(interior) ey(interior) ez(interior)
  READS: byz(everywhere) bzx(everywhere) bxy(everywhere)
  WRITES: dtphi(interior)
  WRITES: dtax(interior) dtay(interior) dtaz(interior)
  WRITES: dtex(interior) dtey(interior) dtez(interior)
  WRITES: dtbyz(interior) dtbzx(interior) dtbxy(interior)
  SYNC: dtphi
  SYNC: dtax dtay dtaz
  SYNC: dtex dtey dtez
  SYNC: dtbyz dtbzx dtbxy       # we could calculate B everywhere
} "Calculate RHS"
# Analysis
SCHEDULE Maxwell_Constraints AT poststep
{
  LANG: C
  READS: ax(interior) ay(interior) az(interior)
  READS: ex(everywhere) ey(everywhere) ez(everywhere)
  READS: byz(interior) bzx(interior) bxy(interior)
  WRITES: curlax(interior) curlay(interior) curlaz(interior)
  WRITES: dive(interior)
  WRITES: divb(interior)
} "Calculate constraints"

#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
namespace Maxwell {
extern "C" void Maxwell_Boundaries(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Maxwell_Boundaries;
  DECLARE_CCTK_PARAMETERS;
  // do nothing
}
} // namespace Maxwell

#include <loop.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
namespace Maxwell {
extern "C" void Maxwell_Constraints(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Maxwell_Constraints;
  DECLARE_CCTK_PARAMETERS;
  const auto DI = Loop::vect<int, Loop::dim>::unit(0);
  const auto DJ = Loop::vect<int, Loop::dim>::unit(1);
  const auto DK = Loop::vect<int, Loop::dim>::unit(2);
  const CCTK_REAL dx = CCTK_DELTA_SPACE(0);
  const CCTK_REAL dy = CCTK_DELTA_SPACE(1);
  const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
  const auto dxm{[&](const auto &u, const auto &p) {
    return (u(p.I) - u(p.I - DI)) / dx;
  }};
  const auto dym{[&](const auto &u, const auto &p) {
    return (u(p.I) - u(p.I - DJ)) / dy;
  }};
  const auto dzm{[&](const auto &u, const auto &p) {
    return (u(p.I) - u(p.I - DK)) / dz;
  }};
  const auto dxp{[&](const auto &u, const auto &p) {
    return (u(p.I + DI) - u(p.I)) / dx;
  }};
  const auto dyp{[&](const auto &u, const auto &p) {
    return (u(p.I + DJ) - u(p.I)) / dy;
  }};
  const auto dzp{[&](const auto &u, const auto &p) {
    return (u(p.I + DK) - u(p.I)) / dz;
  }};
  const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ax_(cctkGH, ax);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ay_(cctkGH, ay);
  const Loop::GF3D<const CCTK_REAL, 0, 0, 1> az_(cctkGH, az);
  const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ex_(cctkGH, ex);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ey_(cctkGH, ey);
  const Loop::GF3D<const CCTK_REAL, 0, 0, 1> ez_(cctkGH, ez);
  const Loop::GF3D<CCTK_REAL, 0, 1, 1> byz_(cctkGH, byz);
  const Loop::GF3D<CCTK_REAL, 1, 0, 1> bzx_(cctkGH, bzx);
  const Loop::GF3D<CCTK_REAL, 1, 1, 0> bxy_(cctkGH, bxy);
  const Loop::GF3D<CCTK_REAL, 0, 1, 1> curlax_(cctkGH, curlax);
  const Loop::GF3D<CCTK_REAL, 1, 0, 1> curlay_(cctkGH, curlay);
  const Loop::GF3D<CCTK_REAL, 1, 1, 0> curlaz_(cctkGH, curlaz);
  const Loop::GF3D<CCTK_REAL, 0, 0, 0> dive_(cctkGH, dive);
  const Loop::GF3D<CCTK_REAL, 1, 1, 1> divb_(cctkGH, divb);
  Loop::loop_int<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    curlax_(p.I) = byz_(p.I) - (dyp(az_, p) - dzp(ay_, p));
  });
  Loop::loop_int<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    curlay_(p.I) = bzx_(p.I) - (dzp(ax_, p) - dxp(az_, p));
  });
  Loop::loop_int<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
    curlaz_(p.I) = bxy_(p.I) - (dxp(ay_, p) - dyp(ax_, p));
  });
  Loop::loop_int<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
    dive_(p.I) = dxm(ex_, p) + dym(ey_, p) + dzm(ez_, p);
  });
  Loop::loop_int<1, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    divb_(p.I) = dxp(byz_, p) + dyp(bzx_, p) + dzp(bxy_, p);
  });
}
} // namespace Maxwell

#include "dual.hxx"
#include <functional>
namespace Maxwell {
using namespace std;
void TestDual() {
  typedef dual<CCTK_REAL> DREAL;
  constexpr equal_to<CCTK_REAL> eq;
  constexpr equal_to<DREAL> eqd;
  static_assert(eq(DREAL().val, 0));
  static_assert(eq(DREAL().eps, 0));
  static_assert(eq(DREAL(1).val, 1));
  static_assert(eq(DREAL(1).eps, 0));
  static_assert(eq(DREAL(1, 2).val, 1));
  static_assert(eq(DREAL(1, 2).eps, 2));
  static_assert(eqd(DREAL(1, 2), DREAL(1, 2)));
  static_assert(!eqd(DREAL(1, 2), DREAL(2, 3)));
  static_assert(eqd(+DREAL(1, 2), DREAL(1, 2)));
  static_assert(eqd(-DREAL(1, 2), DREAL(-1, -2)));
  static_assert(eqd(DREAL(1, 2) + DREAL(3, 4), DREAL(4, 6)));
  static_assert(eqd(DREAL(1, 2) - DREAL(3, 4), DREAL(-2, -2)));
  static_assert(eqd(2 * DREAL(3, 4), DREAL(6, 8)));
  static_assert(eqd(DREAL(3, 4) * 2, DREAL(6, 8)));
  static_assert(eqd(DREAL(3, 4) / 2, DREAL(1.5, 2)));
  static_assert(eqd(DREAL(2, 3) * DREAL(4, 5), DREAL(8, 22)));
  static_assert(eqd(sqrt(DREAL(4, 3)), DREAL(2, 0.75)));
}
} // namespace Maxwell

#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
#include <cassert>
#include <cmath>
#include <functional>
#include <ostream>
namespace Maxwell {
using namespace std;
template <typename T> struct dual {
  T val, eps;
  dual(const dual &) = default;
  dual(dual &&) = default;
  dual &operator=(const dual &) = default;
  dual &operator=(dual &&) = default;
  constexpr dual() : val(), eps() {}
  constexpr dual(const T &x) : val(x), eps() {}
  constexpr dual(const T &x, const T &y) : val(x), eps(y) {}
  friend constexpr dual operator+(const dual &x) { return {+x.val, +x.eps}; }
  friend constexpr dual operator-(const dual &x) { return {-x.val, -x.eps}; }
  friend constexpr dual operator+(const dual &x, const dual &y) {
    return {x.val + y.val, x.eps + y.eps};
  }
  friend constexpr dual operator-(const dual &x, const dual &y) {
    return {x.val - y.val, x.eps - y.eps};
  }
  friend constexpr dual operator+(const dual &x, const T &y) {
    return {x.val + y, x.eps};
  }
  friend constexpr dual operator-(const dual &x, const T &y) {
    return {x.val - y, x.eps};
  }
  friend constexpr dual operator*(const dual &x, const T &y) {
    return {x.val * y, x.eps * y};
  }
  friend constexpr dual operator*(const T &x, const dual &y) {
    return {x * y.val, x * y.eps};
  }
  friend constexpr dual operator/(const dual &x, const T &y) {
    return {x.val / y, x.eps / y};
  }
  friend constexpr dual operator/(const dual &x, const dual &y) {
    assert(y.eps == 0);
    return x / y.val;
  }
  friend constexpr dual operator*(const dual &x, const dual &y) {
    return {x.val * y.val, x.val * y.eps + x.eps * y.val};
  }
  dual &operator+=(const dual &x) { return *this = *this + x; }
  dual &operator-=(const dual &x) { return *this = *this - x; }
  dual &operator*=(const T &x) { return *this = *this * x; }
  dual &operator/=(const T &x) { return *this = *this / x; }
  dual &operator*=(const dual &x) { return *this = *this * x; }
  friend constexpr bool operator==(const dual &x, const dual &y) {
    return x.val == y.val;
  };
  friend constexpr bool operator<(const dual &x, const dual &y) {
    return x.val < y.val;
  };
  friend constexpr bool operator!=(const dual &x, const dual &y) {
    return !(x == y);
  }
  friend constexpr bool operator>(const dual &x, const dual &y) {
    return y < x;
  }
  friend constexpr bool operator<=(const dual &x, const dual &y) {
    return !(x > y);
  }
  friend constexpr bool operator>=(const dual &x, const dual &y) {
    return !(x < y);
  }
  friend ostream &operator<<(ostream &os, const dual &x) {
    return os << x.val << "+eps*" << x.val;
  }
};
} // namespace Maxwell
namespace std {
template <typename T> using dual = Maxwell::dual<T>;
template <typename T> struct equal_to<dual<T> > {
  constexpr bool operator()(const dual<T> &x, const dual<T> &y) const {
    return equal_to<T>()(x.val, y.val) && equal_to<T>()(x.eps, y.eps);
  }
};
template <typename T> struct less<dual<T> > {
  constexpr bool operator()(const dual<T> &x, const dual<T> &y) const {
    return less<T>(x.val, y.val) ||
           (equal_to<T>(x.val, y.val) && less<T>(x.eps, y.eps));
  }
};
template <typename T> constexpr dual<T> cos(const dual<T> &x);
template <typename T> constexpr dual<T> exp(const dual<T> &x);
template <typename T> constexpr dual<T> sin(const dual<T> &x);
template <typename T> constexpr dual<T> sqrt(const dual<T> &x);
template <typename T> constexpr dual<T> cos(const dual<T> &x) {
  return {cos(x.val), -sin(x.val) * x.eps};
}
template <typename T> constexpr dual<T> exp(const dual<T> &x) {
  auto r = exp(x.val);
  return {r, r * x.eps};
}
template <typename T> constexpr dual<T> sin(const dual<T> &x) {
  return {sin(x.val), cos(x.val) * x.eps};
}
template <typename T> constexpr dual<T> sqrt(const dual<T> &x) {
  auto r = sqrt(x.val);
  return {r, x.eps / (2 * r)};
}
} // namespace std

#include "dual.hxx"
#include <loop.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
#include <cassert>
#include <cmath>
#include <complex>
#include <ostream>
namespace Maxwell {
using namespace std;
template <typename T> T pow2(T x) { return x * x; }
////////////////////////////////////////////////////////////////////////////////
template <typename T> struct potential { T phi, ax, ay, az; };
// template <typename F, typename T>
// potential<T> calc_dtp(const F &f, T t, T x, T y, T z, T dt) {
//   auto fm = f(t, x, y, z);
//   auto fp = f(t + dt, x, y, z);
//   return {
//       .phi = (fp.phi - fm.phi) / dt,
//       .ax = (fp.ax - fm.ax) / dt,
//       .ay = (fp.ay - fm.ay) / dt,
//       .az = (fp.az - fm.az) / dt,
//   };
// }
template <typename F, typename T>
potential<T> calc_dxp(const F &f, T t, T x, T y, T z, T dx) {
  auto fm = f(t, x, y, z);
  auto fp = f(t, x + dx, y, z);
  return {
      .phi = (fp.phi - fm.phi) / dx,
      .ax = (fp.ax - fm.ax) / dx,
      .ay = (fp.ay - fm.ay) / dx,
      .az = (fp.az - fm.az) / dx,
  };
}
template <typename F, typename T>
potential<T> calc_dyp(const F &f, T t, T x, T y, T z, T dy) {
  auto fm = f(t, x, y, z);
  auto fp = f(t, x, y + dy, z);
  return {
      .phi = (fp.phi - fm.phi) / dy,
      .ax = (fp.ax - fm.ax) / dy,
      .ay = (fp.ay - fm.ay) / dy,
      .az = (fp.az - fm.az) / dy,
  };
}
template <typename F, typename T>
potential<T> calc_dzp(const F &f, T t, T x, T y, T z, T dz) {
  auto fm = f(t, x, y, z);
  auto fp = f(t, x, y, z + dz);
  return {
      .phi = (fp.phi - fm.phi) / dz,
      .ax = (fp.ax - fm.ax) / dz,
      .ay = (fp.ay - fm.ay) / dz,
      .az = (fp.az - fm.az) / dz,
  };
}
// template <typename F, typename T>
// potential<T> calc_dtm(const F &f, T t, T x, T y, T z, T dt) {
//   auto fm = f(t - dt, x, y, z);
//   auto fp = f(t, x, y, z);
//   return {
//       .phi = (fp.phi - fm.phi) / dt,
//       .ax = (fp.ax - fm.ax) / dt,
//       .ay = (fp.ay - fm.ay) / dt,
//       .az = (fp.az - fm.az) / dt,
//   };
// }
//
// template <typename F, typename T>
// potential<T> calc_dxm(const F &f, T t, T x, T y, T z, T dx) {
//   auto fm = f(t, x - dx, y, z);
//   auto fp = f(t, x, y, z);
//   return {
//       .phi = (fp.phi - fm.phi) / dx,
//       .ax = (fp.ax - fm.ax) / dx,
//       .ay = (fp.ay - fm.ay) / dx,
//       .az = (fp.az - fm.az) / dx,
//   };
// }
//
// template <typename F, typename T>
// potential<T> calc_dym(const F &f, T t, T x, T y, T z, T dy) {
//   auto fm = f(t, x, y - dy, z);
//   auto fp = f(t, x, y, z);
//   return {
//       .phi = (fp.phi - fm.phi) / dy,
//       .ax = (fp.ax - fm.ax) / dy,
//       .ay = (fp.ay - fm.ay) / dy,
//       .az = (fp.az - fm.az) / dy,
//   };
// }
//
// template <typename F, typename T>
// potential<T> calc_dzm(const F &f, T t, T x, T y, T z, T dz) {
//   auto fm = f(t, x, y, z - dz);
//   auto fp = f(t, x, y, z);
//   return {
//       .phi = (fp.phi - fm.phi) / dz,
//       .ax = (fp.ax - fm.ax) / dz,
//       .ay = (fp.ay - fm.ay) / dz,
//       .az = (fp.az - fm.az) / dz,
//   };
// }
template <typename F, typename T>
potential<T> calc_dt(const F &f, T t, T x, T y, T z) {
  auto ff = f(dual<T>(t, 1), dual<T>(x), dual<T>(y), dual<T>(z));
  return {
      ff.phi.eps,
      ff.ax.eps,
      ff.ay.eps,
      ff.az.eps,
  };
}
// template <typename F, typename T>
// potential<T> calc_dx(const F &f, T t, T x, T y, T z) {
//   auto ff = f(dual<T>(t), dual<T>(x, 1), dual<T>(y), dual<T>(z));
//   return {
//       ff.phi.eps,
//       ff.ax.eps,
//       ff.ay.eps,
//       ff.az.eps,
//   };
// }
//
// template <typename F, typename T>
// potential<T> calc_dy(const F &f, T t, T x, T y, T z) {
//   auto ff = f(dual<T>(t), dual<T>(x), dual<T>(y, 1), dual<T>(z));
//   return {
//       ff.phi.eps,
//       ff.ax.eps,
//       ff.ay.eps,
//       ff.az.eps,
//   };
// }
//
// template <typename F, typename T>
// potential<T> calc_dz(const F &f, T t, T x, T y, T z) {
//   auto ff = f(dual<T>(t), dual<T>(x), dual<T>(y), dual<T>(z, 1));
//   return {
//       ff.phi.eps,
//       ff.ax.eps,
//       ff.ay.eps,
//       ff.az.eps,
//   };
// }
////////////////////////////////////////////////////////////////////////////////
// Plane wave
template <typename T>
potential<T> plane_wave(T t, T x, T y, T z, T dx, T dy, T dz) {
  DECLARE_CCTK_PARAMETERS;
  // wave number
  T kx = M_PI * spatial_frequency_x;
  T ky = M_PI * spatial_frequency_y;
  T kz = M_PI * spatial_frequency_z;
  // choose omega to ensure Lorenz gauge
  // continuum
  T omega = sqrt(pow2(kx) + pow2(ky) + pow2(kz));
  // note: we don't have a discrete choice for omega
  T hx = amplitude_x;
  T hy = amplitude_y;
  T hz = amplitude_z;
  // choose ht to ensure div E = 0
  // continuum
  // T ht =
  //     omega * (hx * kx + hy * ky + hz * kz) / (pow2(kx) + pow2(ky) +
  //     pow2(kz));
  // discrete
  typedef complex<T> CT;
  CT ht =
      omega *
      (hx / dx * CT(sin(dx * kx), 2 * pow2(sin(dx * kx / 2))) +
       hy / dy * CT(sin(dy * ky), 2 * pow2(sin(dy * ky / 2))) +
       hz / dz * CT(sin(dz * kz), 2 * pow2(sin(dz * kz / 2)))) /
      CT(4 * ((pow2(sin(dx * kx / 2) / dx)) + (pow2(sin(dy * ky / 2) / dy)) +
              (pow2(sin(dz * kz / 2) / dz))));
  CT u = CT(cos(omega * t - kx * x - ky * y - kz * z),
            sin(omega * t - kx * x - ky * y - kz * z));
  return {
      .phi = real(ht * u),
      .ax = real(hx * u),
      .ay = real(hy * u),
      .az = real(hz * u),
  };
}
// Plane wave with Gaussian profile
template <typename T> potential<T> gaussian_wave(T t, T x, T y, T z) {
  DECLARE_CCTK_PARAMETERS;
  T kx = M_PI * spatial_frequency_x;
  T ky = M_PI * spatial_frequency_y;
  T kz = M_PI * spatial_frequency_z;
  T omega = sqrt(pow2(kx) + pow2(ky) + pow2(kz));
  T hx = amplitude_x;
  T hy = amplitude_y;
  T hz = amplitude_z;
  T ht =
      omega * (hx * kx + hy * ky + hz * kz) / (pow2(kx) + pow2(ky) + pow2(kz));
  T u = exp(-pow2(sin(omega * t - kx * x - ky * y - kz * z) / width) / 2);
  return {
      .phi = ht * u,
      .ax = hx * u,
      .ay = hy * u,
      .az = hz * u,
  };
}
////////////////////////////////////////////////////////////////////////////////
extern "C" void Maxwell_Initial(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Maxwell_Initial;
  DECLARE_CCTK_PARAMETERS;
  const CCTK_REAL t = cctk_time;
  // const CCTK_REAL dt = CCTK_DELTA_TIME;
  const CCTK_REAL dx = CCTK_DELTA_SPACE(0);
  const CCTK_REAL dy = CCTK_DELTA_SPACE(1);
  const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
  const Loop::GF3D<CCTK_REAL, 0, 0, 0> phi_(cctkGH, phi);
  const Loop::GF3D<CCTK_REAL, 1, 0, 0> ax_(cctkGH, ax);
  const Loop::GF3D<CCTK_REAL, 0, 1, 0> ay_(cctkGH, ay);
  const Loop::GF3D<CCTK_REAL, 0, 0, 1> az_(cctkGH, az);
  const Loop::GF3D<CCTK_REAL, 1, 0, 0> ex_(cctkGH, ex);
  const Loop::GF3D<CCTK_REAL, 0, 1, 0> ey_(cctkGH, ey);
  const Loop::GF3D<CCTK_REAL, 0, 0, 1> ez_(cctkGH, ez);
  const Loop::GF3D<CCTK_REAL, 0, 1, 1> byz_(cctkGH, byz);
  const Loop::GF3D<CCTK_REAL, 1, 0, 1> bzx_(cctkGH, bzx);
  const Loop::GF3D<CCTK_REAL, 1, 1, 0> bxy_(cctkGH, bxy);
  const auto loop_setup{[&](const auto &f4) {
    const auto f{[&](const auto &p) { return f4(t, p.x, p.y, p.z); }};
    const auto dtf{[&](const auto &p) {
      return calc_dt(
          [&](auto t, auto x, auto y, auto z) { return f4(t, x, y, z); }, t,
          p.x, p.y, p.z);
    }};
    const auto dxpf{[&](const auto &p) {
      return calc_dxp(
          [&](auto t, auto x, auto y, auto z) { return f4(t, x, y, z); }, t,
          p.x, p.y, p.z, dx);
    }};
    const auto dypf{[&](const auto &p) {
      return calc_dyp(
          [&](auto t, auto x, auto y, auto z) { return f4(t, x, y, z); }, t,
          p.x, p.y, p.z, dy);
    }};
    const auto dzpf{[&](const auto &p) {
      return calc_dzp(
          [&](auto t, auto x, auto y, auto z) { return f4(t, x, y, z); }, t,
          p.x, p.y, p.z, dz);
    }};
    Loop::loop_int<0, 0, 0>(
        cctkGH, [&](const Loop::PointDesc &p) { phi_(p.I) = f(p).phi; });
    Loop::loop_int<1, 0, 0>(
        cctkGH, [&](const Loop::PointDesc &p) { ax_(p.I) = f(p).ax; });
    Loop::loop_int<0, 1, 0>(
        cctkGH, [&](const Loop::PointDesc &p) { ay_(p.I) = f(p).ay; });
    Loop::loop_int<0, 0, 1>(
        cctkGH, [&](const Loop::PointDesc &p) { az_(p.I) = f(p).az; });
    Loop::loop_int<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
      ex_(p.I) = -dxpf(p).phi - dtf(p).ax;
    });
    Loop::loop_int<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
      ey_(p.I) = -dypf(p).phi - dtf(p).ay;
    });
    Loop::loop_int<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
      ez_(p.I) = -dzpf(p).phi - dtf(p).az;
    });
    Loop::loop_int<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
      byz_(p.I) = dypf(p).az - dzpf(p).ay;
    });
    Loop::loop_int<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
      bzx_(p.I) = dzpf(p).az - dxpf(p).az;
    });
    Loop::loop_int<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
      bxy_(p.I) = dxpf(p).ay - dypf(p).ax;
    });
  }};
  if (CCTK_EQUALS(setup, "plane wave")) {
    loop_setup([&](auto t, auto x, auto y, auto z) {
      typedef decltype(t) T;
      return plane_wave(t, x, y, z, T(dx), T(dy), T(dz));
    });
    // } else if (CCTK_EQUALS(setup, "Gaussian wave")) {
    //   loop_setup([&](auto t, auto x, auto y, auto z) {
    //     return gaussian_wave(t, x, y, z);
    //   });
  } else {
    assert(0);
  }
}
} // namespace Maxwell

# Main make.code.defn file for thorn Maxwell
# Source files in this directory
SRCS = boundaries.cxx constraints.cxx dual.cxx initial.cxx rhs.cxx
# Subdirectories containing source files
SUBDIRS =

#include <loop.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
#include <cctk_Parameters.h>
namespace Maxwell {
extern "C" void Maxwell_RHS(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Maxwell_RHS;
  DECLARE_CCTK_PARAMETERS;
  const auto DI = Loop::vect<int, Loop::dim>::unit(0);
  const auto DJ = Loop::vect<int, Loop::dim>::unit(1);
  const auto DK = Loop::vect<int, Loop::dim>::unit(2);
  const CCTK_REAL dx = CCTK_DELTA_SPACE(0);
  const CCTK_REAL dy = CCTK_DELTA_SPACE(1);
  const CCTK_REAL dz = CCTK_DELTA_SPACE(2);
  const auto dxm{[&](const auto &u, const auto &p) {
    return (u(p.I) - u(p.I - DI)) / dx;
  }};
  const auto dym{[&](const auto &u, const auto &p) {
    return (u(p.I) - u(p.I - DJ)) / dy;
  }};
  const auto dzm{[&](const auto &u, const auto &p) {
    return (u(p.I) - u(p.I - DK)) / dz;
  }};
  const auto dxp{[&](const auto &u, const auto &p) {
    return (u(p.I + DI) - u(p.I)) / dx;
  }};
  const auto dyp{[&](const auto &u, const auto &p) {
    return (u(p.I + DJ) - u(p.I)) / dy;
  }};
  const auto dzp{[&](const auto &u, const auto &p) {
    return (u(p.I + DK) - u(p.I)) / dz;
  }};
  const Loop::GF3D<const CCTK_REAL, 0, 0, 0> phi_(cctkGH, phi);
  const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ax_(cctkGH, ax);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ay_(cctkGH, ay);
  const Loop::GF3D<const CCTK_REAL, 0, 0, 1> az_(cctkGH, az);
  const Loop::GF3D<const CCTK_REAL, 1, 0, 0> ex_(cctkGH, ex);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 0> ey_(cctkGH, ey);
  const Loop::GF3D<const CCTK_REAL, 0, 0, 1> ez_(cctkGH, ez);
  const Loop::GF3D<const CCTK_REAL, 0, 1, 1> byz_(cctkGH, byz);
  const Loop::GF3D<const CCTK_REAL, 1, 0, 1> bzx_(cctkGH, bzx);
  const Loop::GF3D<const CCTK_REAL, 1, 1, 0> bxy_(cctkGH, bxy);
  const Loop::GF3D<CCTK_REAL, 0, 0, 0> dtphi_(cctkGH, dtphi);
  const Loop::GF3D<CCTK_REAL, 1, 0, 0> dtax_(cctkGH, dtax);
  const Loop::GF3D<CCTK_REAL, 0, 1, 0> dtay_(cctkGH, dtay);
  const Loop::GF3D<CCTK_REAL, 0, 0, 1> dtaz_(cctkGH, dtaz);
  const Loop::GF3D<CCTK_REAL, 1, 0, 0> dtex_(cctkGH, dtex);
  const Loop::GF3D<CCTK_REAL, 0, 1, 0> dtey_(cctkGH, dtey);
  const Loop::GF3D<CCTK_REAL, 0, 0, 1> dtez_(cctkGH, dtez);
  const Loop::GF3D<CCTK_REAL, 0, 1, 1> dtbyz_(cctkGH, dtbyz);
  const Loop::GF3D<CCTK_REAL, 1, 0, 1> dtbzx_(cctkGH, dtbzx);
  const Loop::GF3D<CCTK_REAL, 1, 1, 0> dtbxy_(cctkGH, dtbxy);
  Loop::loop_int<0, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
    dtphi_(p.I) = -(dxm(ax_, p) + dym(ay_, p) + dzm(az_, p));
  });
  Loop::loop_int<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
    dtax_(p.I) = -dxp(phi_, p) - ex_(p.I);
  });
  Loop::loop_int<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
    dtay_(p.I) = -dyp(phi_, p) - ey_(p.I);
  });
  Loop::loop_int<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    dtaz_(p.I) = -dzp(phi_, p) - ez_(p.I);
  });
  Loop::loop_int<1, 0, 0>(cctkGH, [&](const Loop::PointDesc &p) {
    dtex_(p.I) = -(dzm(bzx_, p) - dym(bxy_, p));
  });
  Loop::loop_int<0, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
    dtey_(p.I) = -(dxm(bxy_, p) - dzm(byz_, p));
  });
  Loop::loop_int<0, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    dtez_(p.I) = -(dym(byz_, p) - dxm(bzx_, p));
  });
  Loop::loop_int<0, 1, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    dtbyz_(p.i, p.j, p.k) = dzp(ey_, p) - dyp(ez_, p);
  });
  Loop::loop_int<1, 0, 1>(cctkGH, [&](const Loop::PointDesc &p) {
    dtbzx_(p.i, p.j, p.k) = dxp(ez_, p) - dzp(ex_, p);
  });
  Loop::loop_int<1, 1, 0>(cctkGH, [&](const Loop::PointDesc &p) {
    dtbxy_(p.i, p.j, p.k) = dyp(ex_, p) - dxp(ey_, p);
  });
}
} // namespace Maxwell