CCTK_REAL allC TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "constraint monitor"

ActiveThorns = "
    ADMBase
    CarpetX
    ErrorEstimator
    Formaline
    IOUtil
    ODESolvers
    Z4c
"
 
$nlevels = 1
$rho = 1
$ncells = 64 * $rho
Cactus::cctk_show_schedule = yes
Cactus::terminate = "time"
Cactus::cctk_final_time = 100.0
CarpetX::verbose = no
CarpetX::xmin = 0.0
CarpetX::ymin = 0.0
CarpetX::zmin = 0.0
CarpetX::xmax = 1.0
CarpetX::ymax = 5.0 / $ncells
CarpetX::zmax = 5.0 / $ncells
CarpetX::ncells_x = $ncells
CarpetX::ncells_y = 2
CarpetX::ncells_z = 2
CarpetX::blocking_factor_x = 64
CarpetX::blocking_factor_y = 2
CarpetX::blocking_factor_z = 2
# CarpetX::periodic_x = no
# CarpetX::periodic_y = no
# CarpetX::periodic_z = no
# CarpetX::periodic = no
# CarpetX::reflection_x = yes
# CarpetX::reflection_y = yes
# CarpetX::reflection_z = yes
# CarpetX::reflection_upper_x = yes
# CarpetX::reflection_upper_y = yes
# CarpetX::reflection_upper_z = yes
CarpetX::ghost_size = 2
CarpetX::max_num_levels = $nlevels
CarpetX::regrid_every = 16
CarpetX::regrid_error_threshold = 1.0 / 4.0
ErrorEstimator::region_shape = "cube"
ErrorEstimator::scale_by_resolution = yes
CarpetX::prolongation_type = "ddf"
CarpetX::prolongation_order = 1
ODESolvers::method = "RK2"
CarpetX::dtfac = 1.0 / 2
ADMBase::initial_data = "linear wave"
ADMBase::initial_lapse = "one"
ADMBase::initial_shift = "zero"
ADMBase::linear_wave_amplitude = 1.0e-8
ADMBase::linear_wave_wavelength = 1.0
Z4c::epsdiss = 0.0
IO::out_dir = $parfile
IO::out_every = 1   #TODO 2 * $ncells * 2 ** ($nlevels - 1)
CarpetX::out_plotfile_groups = ""
CarpetX::out_tsv = yes   #TODO no

ActiveThorns = "
    ADMBase
    CarpetX
    ErrorEstimator
    Formaline
    IOUtil
    ODESolvers
    Z4c
"
 
$nlevels = 3
$ncells = 32
Cactus::cctk_show_schedule = yes
Cactus::terminate = "time"
Cactus::cctk_final_time = 1.0
CarpetX::verbose = yes
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::periodic_x = no
# CarpetX::periodic_y = no
# CarpetX::periodic_z = no
# CarpetX::periodic = no
# CarpetX::reflection_x = yes
# CarpetX::reflection_y = yes
# CarpetX::reflection_z = yes
# CarpetX::reflection_upper_x = yes
# CarpetX::reflection_upper_y = yes
# CarpetX::reflection_upper_z = yes
CarpetX::ghost_size = 2
CarpetX::max_num_levels = $nlevels
CarpetX::regrid_every = 16
CarpetX::regrid_error_threshold = 1.0 / 4.0
ErrorEstimator::region_shape = "cube"
ErrorEstimator::scale_by_resolution = yes
CarpetX::prolongation_type = "ddf"
CarpetX::prolongation_order = 1
ODESolvers::method = "RK2"
CarpetX::dtfac = 0.25
ADMBase::initial_data = "Cartesian Minkowski"
ADMBase::initial_lapse = "one"
ADMBase::initial_shift = "zero"
ADMBase::noise_amplitude = 1.0e-10
IO::out_dir = $parfile
IO::out_every = 1   #TODO $ncells * 2 ** ($nlevels - 1) / 32
CarpetX::out_tsv = no

$nlevels = 1

$nlevels = 3

Cactus::cctk_final_time = 0.0   #TODO 1.0

Cactus::cctk_final_time = 1.0

CarpetX::xmin = 0.0
CarpetX::ymin = 0.0
CarpetX::zmin = 0.0

CarpetX::xmin = -1.0
CarpetX::ymin = -1.0
CarpetX::zmin = -1.0

CarpetX::xmax = 1.0
CarpetX::ymax = 1.0
CarpetX::zmax = 1.0

CarpetX::xmax = +1.0
CarpetX::ymax = +1.0
CarpetX::zmax = +1.0

CarpetX::periodic_x = no
CarpetX::periodic_y = no
CarpetX::periodic_z = no
CarpetX::periodic = no

# CarpetX::periodic_x = no
# CarpetX::periodic_y = no
# CarpetX::periodic_z = no
# CarpetX::periodic = no

CarpetX::ghost_size = 1

CarpetX::ghost_size = 2

CarpetX::regrid_error_threshold = 1.0   #TODO 1.0 / 16.0

CarpetX::regrid_error_threshold = 1.0 / 4.0

ErrorEstimator::scale_by_resolution = no   #TODO yes

ErrorEstimator::scale_by_resolution = yes

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

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

ActiveThorns = "
    ADMBase
    CarpetX
    ErrorEstimator
    Formaline
    IOUtil
    ODESolvers
    Z4c
"
 
$nlevels = 1
$rho = 1
$ncells = 64 * $rho
Cactus::cctk_show_schedule = no
Cactus::terminate = "time"
Cactus::cctk_final_time = 100.0
CarpetX::verbose = no
CarpetX::xmin = 0.0
CarpetX::ymin = 0.0
CarpetX::zmin = 0.0
CarpetX::xmax = 1.0
CarpetX::ymax = 5.0 / $ncells
CarpetX::zmax = 5.0 / $ncells
CarpetX::ncells_x = $ncells
CarpetX::ncells_y = 2
CarpetX::ncells_z = 2
CarpetX::blocking_factor_x = 64
CarpetX::blocking_factor_y = 2
CarpetX::blocking_factor_z = 2
# CarpetX::periodic_x = no
# CarpetX::periodic_y = no
# CarpetX::periodic_z = no
# CarpetX::periodic = no
# CarpetX::reflection_x = yes
# CarpetX::reflection_y = yes
# CarpetX::reflection_z = yes
# CarpetX::reflection_upper_x = yes
# CarpetX::reflection_upper_y = yes
# CarpetX::reflection_upper_z = yes
CarpetX::ghost_size = 2
CarpetX::max_num_levels = $nlevels
CarpetX::regrid_every = 16
CarpetX::regrid_error_threshold = 1.0 / 4.0
ErrorEstimator::region_shape = "cube"
ErrorEstimator::scale_by_resolution = yes
CarpetX::prolongation_type = "ddf"
CarpetX::prolongation_order = 1
ODESolvers::method = "RK2"
CarpetX::dtfac = 0.25
ADMBase::initial_data = "Cartesian Minkowski"
ADMBase::initial_lapse = "one"
ADMBase::initial_shift = "zero"
ADMBase::noise_amplitude = 1.0e-10 / ($rho*$rho)
IO::out_dir = $parfile
IO::out_every = $ncells * 2 ** ($nlevels - 1)
CarpetX::out_plotfile_groups = ""
CarpetX::out_tsv = no

CCTK_REAL epsdiss "Dissipation coefficient <arXiv:gr-qc/0610128>" STEERABLE=always
{
  0.0:* :: ""
} 0.32

STORAGE: allC

SCHEDULE Z4c_Test AT wragh
{
  LANG: C
  OPTIONS: meta
} "Self-test"

SCHEDULE Z4c_Initial1 AT initial AFTER ADMBase_PostInitial

# We have 4 schedule groups:
# 1. initial: set up core Z4c variables from ADM variables
# 2. poststep: post-process core Z4c variables and calculate other variables
# 3. analysis: calculate constraints etc.
# 4. rhs: calculate RHS of Z4c variables
SCHEDULE Z4c_Initial1 IN Z4c_Initial

  SYNC: chi
  SYNC: gamma_tilde
  SYNC: K_hat
  SYNC: A_tilde
  SYNC: Theta
  SYNC: alphaG
  SYNC: betaG

SCHEDULE Z4c_Initial2 AT initial AFTER Z4c_Initial1

SCHEDULE Z4c_Initial2 IN Z4c_Initial AFTER Z4c_Initial1

  SYNC: Gam_tilde

SCHEDULE Z4c_Boundaries AT initial AFTER Z4c_Initial2

SCHEDULE Z4c_Boundaries IN Z4c_PostStep

  SYNC: chi
  SYNC: gamma_tilde
  SYNC: K_hat
  SYNC: A_tilde
  SYNC: Gam_tilde
  SYNC: Theta
  SYNC: alphaG
  SYNC: betaG

SCHEDULE Z4c_Enforce AT initial AFTER Z4c_Boundaries

SCHEDULE Z4c_Enforce IN Z4c_PostStep AFTER Z4c_Boundaries

SCHEDULE Z4c_ADM IN Z4c_PostStep AFTER Z4c_Enforce
{
  LANG: C
  READS: chi(everywhere)
  READS: gamma_tilde(everywhere)
  READS: K_hat(everywhere)
  READS: A_tilde(everywhere)
  READS: Theta(everywhere)
  READS: alphaG(everywhere)
  READS: betaG(everywhere)
  WRITES: ADMBase::metric(everywhere)
  WRITES: ADMBase::curv(everywhere)
  WRITES: ADMBase::lapse(everywhere)
  WRITES: ADMBase::shift(everywhere)
  WRITES: ADMBase::dtshift(everywhere)
} "Convert Z4c to ADM variables"

SCHEDULE Z4c_Constraints AT poststep

SCHEDULE Z4c_Constraints IN Z4c_Analysis

  WRITES: allC(interior)

SCHEDULE Z4c_ADM AT poststep
{
  LANG: C
  READS: chi(everywhere)
  READS: gamma_tilde(everywhere)
  READS: K_hat(everywhere)
  READS: A_tilde(everywhere)
  READS: Theta(everywhere)
  READS: alphaG(everywhere)
  READS: betaG(everywhere)
  WRITES: ADMBase::metric(everywhere)
  WRITES: ADMBase::curv(everywhere)
  WRITES: ADMBase::lapse(everywhere)
  WRITES: ADMBase::shift(everywhere)
  WRITES: ADMBase::dtshift(everywhere)
} "Convert Z4c to ADM variables"

SCHEDULE Z4c_RHS AT poststep

SCHEDULE Z4c_RHS IN Z4c_RHS

} "Calculate Z4c RHS"
SCHEDULE Z4c_RHSBoundaries IN Z4c_RHS AFTER Z4c_RHS
{
  LANG: C
  WRITES: chi_rhs(boundary)
  WRITES: gamma_tilde_rhs(boundary)
  WRITES: K_hat_rhs(boundary)
  WRITES: A_tilde_rhs(boundary)
  WRITES: Gam_tilde_rhs(boundary)
  WRITES: Theta_rhs(boundary)
  WRITES: alphaG_rhs(boundary)
  WRITES: betaG_rhs(boundary)
} "Apply boundary conditions to Z4c RHS variables"
################################################################################
SCHEDULE GROUP Z4c_Initial AT initial AFTER ADMBase_PostInitial
{
} "Convert ADM to Z4c variables"
SCHEDULE GROUP Z4c_PostStep AT initial AFTER Z4c_Initial
{
} "Post-process Z4c variables"
SCHEDULE GROUP Z4c_PostStep AT postregrid
{
} "Post-process Z4c variables"
SCHEDULE GROUP Z4c_Analysis AT analysis
{
} "Analyse Z4c variables"
SCHEDULE GROUP Z4c_PostStep IN ODESolvers_PostStep
{
} "Post-process Z4c variables"
SCHEDULE GROUP Z4c_RHS IN ODESolvers_RHS
{

SCHEDULE GROUP Z4c_PostSTep AT postrestrict
{
} "Post-process Z4c variables"

#include "z4c_vars.hxx"

#include <cctk_Parameters.h>

  DECLARE_CCTK_PARAMETERS;

    // Load
    const CCTK_REAL chi = gf_chi_(p.I);
    const mat3<CCTK_REAL> gammat(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_,
                                 gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p);
    const CCTK_REAL Kh = gf_Kh_(p.I);
    const mat3<CCTK_REAL> At(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_,
                             gf_Atzz_, p);
    const CCTK_REAL Theta = gf_Theta_(p.I);
    const CCTK_REAL alphaG = gf_alphaG_(p.I);
    const vec3<CCTK_REAL> betaG(gf_betaGx_, gf_betaGy_, gf_betaGz_, p);
    // Calculate ADM variables
    mat3<CCTK_REAL> g;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b)
        g(a, b) = chi * gammat(a, b);
    const CCTK_REAL K = Kh - 2 * Theta;
    mat3<CCTK_REAL> k;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b)
        k(a, b) = chi * At(a, b) + K / 3 * g(a, b);
    const CCTK_REAL alp = alphaG;
    vec3<CCTK_REAL> beta;
    for (int a = 0; a < 3; ++a)
      beta(a) = betaG(a);

    // Load and calculate
    const z4c_vars_noderivs<CCTK_REAL> vars(
        kappa1, kappa2, f_mu_L, f_mu_S, eta, //
        gf_chi_, gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,
        gf_gammatyz_, gf_gammatzz_, gf_Kh_, gf_Atxx_, gf_Atxy_, gf_Atxz_,
        gf_Atyy_, gf_Atyz_, gf_Atzz_,
        //
        // gf_Gamtx_, gf_Gamty_, gf_Gamtz_,
        gf_chi_, gf_chi_, gf_chi_,
        //
        gf_Theta_, gf_alphaG_, gf_betaGx_, gf_betaGy_, gf_betaGz_, //
        p.I);

    g.store(gf_gxx_, gf_gxy_, gf_gxz_, gf_gyy_, gf_gyz_, gf_gzz_, p);
    k.store(gf_kxx_, gf_kxy_, gf_kxz_, gf_kyy_, gf_kyz_, gf_kzz_, p);
    gf_alp_(p.I) = alp;
    beta.store(gf_betax_, gf_betay_, gf_betaz_, p);

    vars.g.store(gf_gxx_, gf_gxy_, gf_gxz_, gf_gyy_, gf_gyz_, gf_gzz_, p);
    vars.k.store(gf_kxx_, gf_kxy_, gf_kxz_, gf_kyy_, gf_kyz_, gf_kzz_, p);
    gf_alp_(p.I) = vars.alp;
    vars.beta.store(gf_betax_, gf_betay_, gf_betaz_, p);

#include "z4c_vars.hxx"

#include <cctk_Parameters.h>

  DECLARE_CCTK_PARAMETERS;
  for (int d = 0; d < 3; ++d)
    if (cctk_nghostzones[d] < 2)
      CCTK_VERROR("Need at least 2 ghost zones");

  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {
    // Load
    const CCTK_REAL chi = gf_chi_(p.I);
    const mat3<CCTK_REAL> gammat(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_,
                                 gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p);
    const CCTK_REAL Kh = gf_Kh_(p.I);
    const mat3<CCTK_REAL> At(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_,
                             gf_Atzz_, p);
    const vec3<CCTK_REAL> Gamt(gf_Gamtx_, gf_Gamty_, gf_Gamtz_, p);
    const CCTK_REAL Theta = gf_Theta_(p.I);
    constexpr CCTK_REAL rho = 0;
    constexpr vec3<CCTK_REAL> S{0, 0, 0};
    // Calculate constraints
    // See arXiv:1212.2901 [gr-qc].
    const CCTK_REAL chi1 = 1 / chi;
    const vec3<CCTK_REAL> dchi = deriv(gf_chi_, p.I, dx);
    const mat3<CCTK_REAL> ddchi = deriv2(gf_chi_, p.I, dx);
    mat3<CCTK_REAL> g;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b)
        g(a, b) = chi * gammat(a, b);
    const mat3<CCTK_REAL> gammatu = gammat.inv(1);
    mat3<CCTK_REAL> gu;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b)
        gu(a, b) = chi1 * gammatu(a, b);
    const mat3<vec3<CCTK_REAL> > dgammat{
        deriv(gf_gammatxx_, p.I, dx), deriv(gf_gammatxy_, p.I, dx),
        deriv(gf_gammatxz_, p.I, dx), deriv(gf_gammatyy_, p.I, dx),
        deriv(gf_gammatyz_, p.I, dx), deriv(gf_gammatzz_, p.I, dx),
    };
    const mat3<mat3<CCTK_REAL> > ddgammat{
        deriv2(gf_gammatxx_, p.I, dx), deriv2(gf_gammatxy_, p.I, dx),
        deriv2(gf_gammatxz_, p.I, dx), deriv2(gf_gammatyy_, p.I, dx),
        deriv2(gf_gammatyz_, p.I, dx), deriv2(gf_gammatzz_, p.I, dx),
    };
    const mat3<vec3<CCTK_REAL> > dgammatu = calc_dgu(gammatu, dgammat);
    const vec3<CCTK_REAL> dKh = deriv(gf_Kh_, p.I, dx);
    mat3<CCTK_REAL> Atu;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b) {
        CCTK_REAL s = 0;
        for (int x = 0; x < 3; ++x)
          for (int y = 0; y < 3; ++y)
            s += gammatu(a, x) * gammatu(b, y) * At(x, y);
        Atu(a, b) = s;
      }
    const mat3<vec3<CCTK_REAL> > dAt{
        deriv(gf_Atxx_, p.I, dx), deriv(gf_Atxy_, p.I, dx),
        deriv(gf_Atxz_, p.I, dx), deriv(gf_Atyy_, p.I, dx),
        deriv(gf_Atyz_, p.I, dx), deriv(gf_Atzz_, p.I, dx),
    };
    const mat3<vec3<CCTK_REAL> > dAtu = calc_dAu(gammatu, dgammatu, At, dAt);
    const vec3<vec3<CCTK_REAL> > dGamt{
        deriv(gf_Gamtx_, p.I, dx),
        deriv(gf_Gamty_, p.I, dx),
        deriv(gf_Gamtz_, p.I, dx),
    };
    vec3<mat3<CCTK_REAL> > Gammatdl;
    vec3<mat3<CCTK_REAL> > Gammatd;
    vec3<CCTK_REAL> Gamtd;
    calc_gamma(gammat, gammatu, dgammat, Gammatdl, Gammatd, Gamtd);
    mat3<CCTK_REAL> DDchi;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b) {
        CCTK_REAL s = ddchi(a, b);
        for (int x = 0; x < 3; ++x)
          s -= Gammatd(x)(a, b) * dchi(x);
        DDchi(a, b) = s;
      }

    const vec3<CCTK_REAL> dTheta = deriv(gf_Theta_, p.I, dx);
    // (13)
    vec3<CCTK_REAL> ZtC;
    for (int a = 0; a < 3; ++a)
      ZtC(a) = (Gamt(a) - Gamtd(a)) / 2;
    const mat3<CCTK_REAL> R =
        calc_ricci(chi, dchi, DDchi, gammat, gammatu, ddgammat, Gammatdl,
                   Gammatd, Gamtd, dGamt);
    CCTK_REAL Rsc = 0;
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        Rsc += gu(x, y) * R(x, y);
    // (14)
    CCTK_REAL HC = 0;
    HC += Rsc;
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        HC += At(x, y) * Atu(x, y);
    HC -= pow2(Kh + 2 * Theta) * 2 / 3;
    HC -= 16 * M_PI * rho;

  const GF3D<CCTK_REAL, 0, 0, 0> gf_allC_(cctkGH, allC);

    // (15)
    vec3<CCTK_REAL> MtC;
    for (int a = 0; a < 3; ++a) {
      CCTK_REAL s = 0;
      for (int x = 0; x < 3; ++x)
        s += dAtu(a, x)(x);
      for (int x = 0; x < 3; ++x)
        for (int y = 0; y < 3; ++y)
          s += Gammatd(a)(x, y) * Atu(x, y);
      for (int x = 0; x < 3; ++x)
        s -= gammatu(a, x) * (dKh(x) + 2 * dTheta(x)) * 2 / 3;
      for (int x = 0; x < 3; ++x)
        s -= Atu(a, x) * dchi(x) / chi * 2 / 3;
      for (int x = 0; x < 3; ++x)
        s -= 8 * M_PI * gammatu(a, x) * S(x);
      MtC(a) = s;
    }

  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {
    // Load and calculate
    const z4c_vars<CCTK_REAL> vars(
        kappa1, kappa2, f_mu_L, f_mu_S, eta, //
        gf_chi_, gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,
        gf_gammatyz_, gf_gammatzz_, gf_Kh_, gf_Atxx_, gf_Atxy_, gf_Atxz_,
        gf_Atyy_, gf_Atyz_, gf_Atzz_, gf_Gamtx_, gf_Gamty_, gf_Gamtz_,
        gf_Theta_,
        //
        // gf_alphaG_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
        gf_chi_, gf_chi_, gf_chi_, gf_chi_,
        //
        p.I, dx);

    ZtC.store(gf_ZtCx_, gf_ZtCy_, gf_ZtCz_, p);
    gf_HC_(p.I) = HC;
    MtC.store(gf_MtCx_, gf_MtCy_, gf_MtCz_, p);

    vars.ZtC.store(gf_ZtCx_, gf_ZtCy_, gf_ZtCz_, p);
    gf_HC_(p.I) = vars.HC;
    vars.MtC.store(gf_MtCx_, gf_MtCy_, gf_MtCz_, p);
    gf_allC_(p.I) = vars.allC;

#include <sstream>

                           gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p);

                           gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p.I);

                       gf_Atzz_, p);

                       gf_Atzz_, p.I);

    const CCTK_REAL det_gammat = gammat.det();
    const CCTK_REAL psi1 = 1 / cbrt(det_gammat);

    const CCTK_REAL detgammat = gammat.det();
    const CCTK_REAL chi_new = 1 / cbrt(detgammat);

        gammat(a, b) *= psi1;

        gammat(a, b) *= chi_new;
#ifdef CCTK_DEBUG
    const CCTK_REAL detgammat_new = gammat.det();
    const CCTK_REAL gammat_norm = gammat.maxabs();
    const CCTK_REAL gammat_scale = gammat_norm;
    if (!(fabs(detgammat_new - 1) <= 1.0e-12 * gammat_scale)) {
      ostringstream buf;
      buf << "det gammat is not one: gammat=" << gammat
          << " det(gammat)=" << detgammat_new;
      CCTK_VERROR("%s", buf.str().c_str());
    }
    assert(fabs(detgammat_new - 1) <= 1.0e-12 * gammat_scale);
#endif

    CCTK_REAL traceAt = 0;
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        traceAt += gammatu(x, y) * At(x, y);

    const CCTK_REAL traceAt =
        sum2([&](int x, int y) { return gammatu(x, y) * At(x, y); });

        At(a, b) -= 1 / 3 * traceAt * gammat(a, b);

        At(a, b) -= traceAt / 3 * gammat(a, b);
#ifdef CCTK_DEBUG
    const CCTK_REAL traceAt_new =
        sum2([&](int x, int y) { return gammatu(x, y) * At(x, y); });
    const CCTK_REAL gammatu_norm = gammatu.maxabs();
    const CCTK_REAL At_norm = At.maxabs();
    const CCTK_REAL At_scale = fmax(fmax(gammat_norm, gammatu_norm), At_norm);
    if (!(fabs(traceAt_new) <= 1.0e-12 * At_scale)) {
      ostringstream buf;
      buf << "tr At: At=" << At << " tr(At)=" << traceAt_new;
      CCTK_VERROR("%s", buf.str().c_str());
    }
    assert(fabs(traceAt_new) <= 1.0e-12 * At_scale);
#endif

#ifndef FIELD_HXX
#define FIELD_HXX
#include <loop.hxx>
#include <cassert>
#include <cstddef>
#include <vector>
namespace Z4c {
using namespace Loop;
using namespace std;
template <typename T> class field {
  const GridDescBase &restrict grid;
  vect<int, dim> imin, imax;
  vect<ptrdiff_t, dim> istr;
  vector<T> elts;
  constexpr ptrdiff_t idx(const vect<T, dim> &I) const {
#ifdef CCTK_DEBUG
    for (int d = 0; d < dim; ++d)
      assert(I[d] >= imin[d] && I[d] < imax[d]);
#endif
    ptrdiff_t n = 0;
    if (dim > 0) {
      n += I[0];
      for (int d = 1; d < dim; ++d)
        n += istr[d] * I[d];
    }
#ifdef CCTK_DEBUG
    assert(n >= 0 && n < ptrdiff_t(elts.size()));
#endif
    return n;
  }
public:
  field(const GridDescBase &grid) : grid(grid) {
    grid.box_all<0, 0, 0>(grid.nghostzones, imin, imax);
    ptrdiff_t str = 1;
    for (int d = 0; d < dim; ++d) {
      istr[d] = str;
      str *= imax[d] - imin[d];
    }
    elts.resize(str);
  }
  const T &operator()(const vect<T, dim> &I) const { return elts[idx(I)]; }
  T &operator()(const vect<T, dim> &I) { return elts[idx(I)]; }
  template <int dir>
  field deriv(const int order, const vect<T, dim> &dx) const {
    static_assert(dir >= 0 && dir < dim, "");
    constexpr auto DI = vect<int, dim>::unit(dir);
    const auto &f = *this;
    field r(grid);
    switch (order) {
    case 2:
      assert(grid.nghostzones[dir] >= 1);
      loop_int(grid.nghostzones, [&](const PointDesc &p) {
        r(p.I) = -1 / T(2) * (f(p.I - DI) - f(p.I + DI)) / dx[dir];
      });
      break;
    case 4:
      assert(grid.nghostzones[dir] >= 2);
      loop_int(grid.nghostzones, [&](const PointDesc &p) {
        r(p.I) = (1 / T(12) * (f(p.I - 2 * DI) - f(p.I + 2 * DI)) +
                  2 / T(3) * (f(p.I - DI) - f(p.I + DI))) /
                 dx[dir];
      });
      break;
    default:
      assert(0);
    }
    return r;
  }
};
} // namespace Z4c
#endif // #ifndef FIELD_HXX

                            gf_gzz_, p);

                            gf_gzz_, p.I);

                            gf_kzz_, p);

                            gf_kzz_, p.I);

    const vec3<CCTK_REAL> beta(gf_betax_, gf_betay_, gf_betaz_, p);

    const vec3<CCTK_REAL> beta(gf_betax_, gf_betay_, gf_betaz_, p.I);

    const CCTK_REAL chi = cbrt(detg);
    const CCTK_REAL chi1 = 1 / chi;

    const CCTK_REAL chi = 1 / cbrt(detg);

    mat3<CCTK_REAL> gammat;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b)
        gammat(a, b) = chi1 * g(a, b);

    const mat3<CCTK_REAL> gammat([&](int a, int b) { return chi * g(a, b); });

    CCTK_REAL K = 0;
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        K += gu(x, y) * k(x, y);

    const CCTK_REAL K = sum2([&](int x, int y) { return gu(x, y) * k(x, y); });

    mat3<CCTK_REAL> At;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b)
        At(a, b) = chi1 * (k(a, b) - K / 3 * g(a, b));

    const mat3<CCTK_REAL> At(
        [&](int a, int b) { return chi * (k(a, b) - K / 3 * g(a, b)); });

    vec3<CCTK_REAL> betaG;
    for (int a = 0; a < 3; ++a)
      betaG(a) = beta(a);

    const vec3<CCTK_REAL> betaG([&](int a) { return beta(a); });

#include "field.hxx"

                                 gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p);

                                 gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p.I);

    vec3<mat3<CCTK_REAL> > Gammatl;
    vec3<mat3<CCTK_REAL> > Gammat;
    vec3<CCTK_REAL> Gamt;
    calc_gamma(gammat, gammatu, dgammat, Gammatl, Gammat, Gamt);

    const vec3<mat3<CCTK_REAL> > Gammatl = calc_gammal(dgammat);
    const vec3<mat3<CCTK_REAL> > Gammat = calc_gamma(gammatu, Gammatl);
    const vec3<CCTK_REAL> Gamt([&](int a) {
      return sum2(
          [&](int x, int y) { return gammatu(x, y) * Gammat(a)(x, y); });
    });

	rhs.cxx

	rhs.cxx					\
	rhsboundaries.cxx			\
	test.cxx

mat3<vec3<T> > calc_dgu(const mat3<T> &gu, const mat3<vec3<T> > &dg) {

/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<vec3<T> >
calc_dgu(const mat3<T> &gu, const mat3<vec3<T> > &dg) {

  mat3<vec3<T> > dgu;
  for (int a = 0; a < 3; ++a)
    for (int b = a; b < 3; ++b)
      for (int c = 0; c < 3; ++c) {
        T s = 0;
        for (int x = 0; x < 3; ++x)
          for (int y = 0; y < 3; ++y)
            s -= gu(a, x) * gu(b, y) * dg(x, y)(c);
        dgu(a, b)(c) = s;
      }
  return dgu;

  return mat3<vec3<T> >([&](int a, int b) {
    return vec3<T>([&](int c) {
      return sum2(
          [&](int x, int y) { return -gu(a, x) * gu(b, y) * dg(x, y)(c); });
    });
  });

mat3<vec3<T> > calc_dAu(const mat3<T> &gu, const mat3<vec3<T> > &dgu,
                        const mat3<T> &A, const mat3<vec3<T> > &dA) {

/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<vec3<T> >
calc_dAu(const mat3<T> &gu, const mat3<vec3<T> > &dgu, const mat3<T> &A,
         const mat3<vec3<T> > &dA) {

  mat3<vec3<T> > dAu;
  for (int a = 0; a < 3; ++a)
    for (int b = a; b < 3; ++b)
      for (int c = 0; c < 3; ++c) {
        T s = 0;
        for (int x = 0; x < 3; ++x)
          for (int y = 0; y < 3; ++y)
            s += dgu(a, x)(c) * gu(b, y) * A(x, y) +
                 gu(a, x) * dgu(b, y)(c) * A(x, y) +
                 gu(a, x) * gu(b, y) * dA(x, y)(c);
        dAu(a, b)(c) = s;
      }
  return dAu;

  return mat3<vec3<T> >([&](int a, int b) {
    return vec3<T>([&](int c) {
      return sum2([&](int x, int y) {
        return dgu(a, x)(c) * gu(b, y) * A(x, y)   //
               + gu(a, x) * dgu(b, y)(c) * A(x, y) //
               + gu(a, x) * gu(b, y) * dA(x, y)(c);
      });
    });
  });

void calc_gamma(const mat3<T> &g, const mat3<T> &gu, const mat3<vec3<T> > dg,
                vec3<mat3<T> > &restrict Gammal, vec3<mat3<T> > &restrict Gamma,
                vec3<T> &restrict Gam) {

/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<mat3<T> >
calc_gammal(const mat3<vec3<T> > dg) {

  for (int a = 0; a < 3; ++a)
    for (int b = 0; b < 3; ++b)
      for (int c = b; c < 3; ++c)
        Gammal(a)(b, c) = (dg(a, b)(c) + dg(a, c)(b) - dg(a, b)(c)) / 2;
  // Gamma^a_bc
  for (int a = 0; a < 3; ++a)
    for (int b = 0; b < 3; ++b)
      for (int c = b; c < 3; ++c) {
        T s = 0;
        for (int x = 0; x < 3; ++x)
          s += gu(a, x) * Gammal(x)(b, c);
        Gamma(a)(b, c) = s;
      }
  // Gam^a
  for (int a = 0; a < 3; ++a) {
    T s = 0;
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        s += gu(x, y) * Gamma(a)(x, y);
    Gam(a) = s;
  }

  return vec3<mat3<T> >([&](int a) {
    return mat3<T>([&](int b, int c) {
      return (dg(a, b)(c) + dg(a, c)(b) - dg(b, c)(a)) / 2;
    });
  });

mat3<T> calc_ricci(const T &chi, const vec3<T> &dchi, const mat3<T> &DDchi,
                   const mat3<T> &gammat, const mat3<T> &gammatu,
                   const mat3<mat3<T> > &ddgammat,
                   const vec3<mat3<T> > &Gammatdl,
                   const vec3<mat3<T> > &Gammatd, const vec3<T> &Gamtd,
                   const vec3<vec3<T> > &dGamt) {
  // (8)
  mat3<T> Rchi;
  for (int a = 0; a < 3; ++a)
    for (int b = 0; b < 3; ++b) {
      T s = 0;
      s += DDchi(a, b) / (2 * chi);
      for (int x = 0; x < 3; ++x)
        for (int y = 0; y < 3; ++y)
          s += gammat(a, b) * gammatu(x, y) * DDchi(x, y);
      s -= dchi(a) * dchi(b) / (4 * chi * chi);
      for (int x = 0; x < 3; ++x)
        for (int y = 0; y < 3; ++y)
          s -= gammat(a, b) * gammatu(x, y) * dchi(x) * dchi(y) * 3 /
               (4 * chi * chi);
      Rchi(a, b) = s;
    }
  // (9)
  mat3<T> Rt;
  for (int a = 0; a < 3; ++a)
    for (int b = 0; b < 3; ++b) {
      T s = 0;
      for (int x = 0; x < 3; ++x)
        for (int y = 0; y < 3; ++y)
          s -= gammatu(x, y) * ddgammat(a, b)(x, y) / 2;
      for (int x = 0; x < 3; ++x)
        s += gammat(x, a) * dGamt(x)(b) + gammat(x, b) * dGamt(x)(a);
      for (int x = 0; x < 3; ++x)
        s += Gamtd(x) * Gammatdl(a)(b, x) + Gamtd(x) * Gammatdl(b)(a, x);
      for (int x = 0; x < 3; ++x)
        for (int y = 0; y < 3; ++y)
          for (int z = 0; z < 3; ++z)
            s += gammatu(y, z) * (2 * Gammatd(x)(y, a) * Gammatdl(b)(x, z) +
                                  2 * Gammatd(x)(y, b) * Gammatdl(a)(x, z) +
                                  Gammatd(x)(a, z) * Gammatdl(x)(y, b) +
                                  Gammatd(x)(b, z) * Gammatdl(x)(y, a));
      Rt(a, b) = s;
    }
  // (7)
  mat3<T> R;
  for (int a = 0; a < 3; ++a)
    for (int b = 0; b < 3; ++b)
      R(a, b) = Rchi(a, b) + Rt(a, b);
  return R;

/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<mat3<T> >
calc_gamma(const mat3<T> &gu, const vec3<mat3<T> > &Gammal) {
  // Gamma^a_bc
  return vec3<mat3<T> >([&](int a) {
    return mat3<T>([&](int b, int c) {
      return sum1([&](int x) { return gu(a, x) * Gammal(x)(b, c); });
    });
  });

#warning "TODO"
#include <iostream>

#include "z4c_vars.hxx"

template <typename T>
void apply_upwind(const cGH *restrict const cctkGH,
                  const GF3D<const T, 0, 0, 0> gf_,
                  const GF3D<const T, 0, 0, 0> gf_betaGx_,
                  const GF3D<const T, 0, 0, 0> gf_betaGy_,
                  const GF3D<const T, 0, 0, 0> gf_betaGz_,
                  const GF3D<T, 0, 0, 0> gf_rhs_) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  const vec3<CCTK_REAL> dx{
      CCTK_DELTA_SPACE(0),
      CCTK_DELTA_SPACE(1),
      CCTK_DELTA_SPACE(2),
  };
  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {
    const vec3<CCTK_REAL> betaG(gf_betaGx_, gf_betaGy_, gf_betaGz_, p.I);
    const vec3<CCTK_REAL> dgf_upwind(deriv_upwind(gf_, p.I, betaG, dx));
    gf_rhs_(p.I) += sum1([&](int x) { return betaG(x) * dgf_upwind(x); });
  });
}
template <typename T>
void apply_diss(const cGH *restrict const cctkGH,
                const GF3D<const T, 0, 0, 0> gf_,
                const GF3D<T, 0, 0, 0> gf_rhs_) {
  DECLARE_CCTK_ARGUMENTS;
  DECLARE_CCTK_PARAMETERS;
  const vec3<CCTK_REAL> dx{
      CCTK_DELTA_SPACE(0),
      CCTK_DELTA_SPACE(1),
      CCTK_DELTA_SPACE(2),
  };
  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {
    gf_rhs_(p.I) += epsdiss * diss(gf_, p.I, dx);
  });
}

  for (int d = 0; d < 3; ++d)
    if (cctk_nghostzones[d] < 2)
      CCTK_VERROR("Need at least 2 ghost zones");

    // Load
    const CCTK_REAL chi = gf_chi_(p.I);
    const mat3<CCTK_REAL> gammat(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_,
                                 gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p);
    const CCTK_REAL Kh = gf_Kh_(p.I);
    const mat3<CCTK_REAL> At(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_,
                             gf_Atzz_, p);
    const vec3<CCTK_REAL> Gamt(gf_Gamtx_, gf_Gamty_, gf_Gamtz_, p);
    const CCTK_REAL Theta = gf_Theta_(p.I);

    // Load and calculate
    const z4c_vars<CCTK_REAL> vars(
        kappa1, kappa2, f_mu_L, f_mu_S, eta, //
        gf_chi_, gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,
        gf_gammatyz_, gf_gammatzz_, gf_Kh_, gf_Atxx_, gf_Atxy_, gf_Atxz_,
        gf_Atyy_, gf_Atyz_, gf_Atzz_, gf_Gamtx_, gf_Gamty_, gf_Gamtz_,
        gf_Theta_, gf_alphaG_, gf_betaGx_, gf_betaGy_, gf_betaGz_, //
        p.I, dx);

    const CCTK_REAL alphaG = gf_alphaG_(p.I);

    // Store
    gf_chi_rhs_(p.I) = vars.chi_rhs;
    vars.gammat_rhs.store(gf_gammatxx_rhs_, gf_gammatxy_rhs_, gf_gammatxz_rhs_,
                          gf_gammatyy_rhs_, gf_gammatyz_rhs_, gf_gammatzz_rhs_,
                          p);
    gf_Kh_rhs_(p.I) = vars.Kh_rhs;
    vars.At_rhs.store(gf_Atxx_rhs_, gf_Atxy_rhs_, gf_Atxz_rhs_, gf_Atyy_rhs_,
                      gf_Atyz_rhs_, gf_Atzz_rhs_, p);
    vars.Gamt_rhs.store(gf_Gamtx_rhs_, gf_Gamty_rhs_, gf_Gamtz_rhs_, p);
    gf_Theta_rhs_(p.I) = vars.Theta_rhs;
    gf_alphaG_rhs_(p.I) = vars.alphaG_rhs;
    vars.betaG_rhs.store(gf_betaGx_rhs_, gf_betaGy_rhs_, gf_betaGz_rhs_, p);
  });

    const vec3<CCTK_REAL> betaG(gf_betaGx_, gf_betaGy_, gf_betaGz_, p);

  // Upwind terms

    constexpr CCTK_REAL rho = 0;
    constexpr vec3<CCTK_REAL> S{0, 0, 0};
    constexpr CCTK_REAL traceT = 0;
    constexpr mat3<CCTK_REAL> T{0, 0, 0, 0, 0, 0};

  apply_upwind(cctkGH, gf_chi_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_chi_rhs_);

    // Calculate RHS
    // See arXiv:1212.2901 [gr-qc].

  apply_upwind(cctkGH, gf_gammatxx_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_gammatxx_rhs_);
  apply_upwind(cctkGH, gf_gammatxy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_gammatxy_rhs_);
  apply_upwind(cctkGH, gf_gammatxz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_gammatxz_rhs_);
  apply_upwind(cctkGH, gf_gammatyy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_gammatyy_rhs_);
  apply_upwind(cctkGH, gf_gammatyz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_gammatyz_rhs_);
  apply_upwind(cctkGH, gf_gammatzz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_gammatzz_rhs_);

    const CCTK_REAL chi1 = 1 / chi;

  apply_upwind(cctkGH, gf_Kh_, gf_betaGx_, gf_betaGy_, gf_betaGz_, gf_Kh_rhs_);

    const vec3<CCTK_REAL> dchi = deriv(gf_chi_, p.I, dx);
    const mat3<CCTK_REAL> ddchi = deriv2(gf_chi_, p.I, dx);

  apply_upwind(cctkGH, gf_Atxx_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Atxx_rhs_);
  apply_upwind(cctkGH, gf_Atxy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Atxy_rhs_);
  apply_upwind(cctkGH, gf_Atxz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Atxz_rhs_);
  apply_upwind(cctkGH, gf_Atyy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Atyy_rhs_);
  apply_upwind(cctkGH, gf_Atyz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Atyz_rhs_);
  apply_upwind(cctkGH, gf_Atzz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Atzz_rhs_);

    mat3<CCTK_REAL> g;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b)
        g(a, b) = chi * gammat(a, b);

  apply_upwind(cctkGH, gf_Gamtx_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Gamtx_rhs_);
  apply_upwind(cctkGH, gf_Gamty_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Gamty_rhs_);
  apply_upwind(cctkGH, gf_Gamtz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Gamtz_rhs_);

    const mat3<CCTK_REAL> gammatu = gammat.inv(1);
    mat3<CCTK_REAL> gu;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b)
        gu(a, b) = chi1 * gammatu(a, b);

  apply_upwind(cctkGH, gf_Theta_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_Theta_rhs_);

    const mat3<vec3<CCTK_REAL> > dgammat{
        deriv(gf_gammatxx_, p.I, dx), deriv(gf_gammatxy_, p.I, dx),
        deriv(gf_gammatxz_, p.I, dx), deriv(gf_gammatyy_, p.I, dx),
        deriv(gf_gammatyz_, p.I, dx), deriv(gf_gammatzz_, p.I, dx),
    };
    mat3<vec3<CCTK_REAL> > dg;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b)
        for (int c = 0; c < 3; ++c)
          dg(a, b)(c) = dchi(c) * g(a, b) + chi * dgammat(a, b)(c);

  apply_upwind(cctkGH, gf_alphaG_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_alphaG_rhs_);

    const mat3<mat3<CCTK_REAL> > ddgammat{
        deriv2(gf_gammatxx_, p.I, dx), deriv2(gf_gammatxy_, p.I, dx),
        deriv2(gf_gammatxz_, p.I, dx), deriv2(gf_gammatyy_, p.I, dx),
        deriv2(gf_gammatyz_, p.I, dx), deriv2(gf_gammatzz_, p.I, dx),
    };

  apply_upwind(cctkGH, gf_betaGx_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_betaGx_rhs_);
  apply_upwind(cctkGH, gf_betaGy_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_betaGy_rhs_);
  apply_upwind(cctkGH, gf_betaGz_, gf_betaGx_, gf_betaGy_, gf_betaGz_,
               gf_betaGz_rhs_);

    const vec3<CCTK_REAL> dKh = deriv(gf_Kh_, p.I, dx);

  // Dissipation

    mat3<CCTK_REAL> Atu;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b) {
        CCTK_REAL s = 0;
        for (int x = 0; x < 3; ++x)
          for (int y = 0; y < 3; ++y)
            s += gammatu(a, x) * gammatu(b, y) * At(x, y);
        Atu(a, b) = s;
      }

  apply_diss(cctkGH, gf_chi_, gf_chi_rhs_);

    const mat3<vec3<CCTK_REAL> > dAt{
        deriv(gf_Atxx_, p.I, dx), deriv(gf_Atxy_, p.I, dx),
        deriv(gf_Atxz_, p.I, dx), deriv(gf_Atyy_, p.I, dx),
        deriv(gf_Atyz_, p.I, dx), deriv(gf_Atzz_, p.I, dx),
    };

  apply_diss(cctkGH, gf_gammatxx_, gf_gammatxx_rhs_);
  apply_diss(cctkGH, gf_gammatxy_, gf_gammatxy_rhs_);
  apply_diss(cctkGH, gf_gammatxz_, gf_gammatxz_rhs_);
  apply_diss(cctkGH, gf_gammatyy_, gf_gammatyy_rhs_);
  apply_diss(cctkGH, gf_gammatyz_, gf_gammatyz_rhs_);
  apply_diss(cctkGH, gf_gammatzz_, gf_gammatzz_rhs_);

    vec3<mat3<CCTK_REAL> > Gammatl;
    vec3<mat3<CCTK_REAL> > Gammat;
    vec3<CCTK_REAL> Gamtd;
    calc_gamma(gammat, gammatu, dgammat, Gammatl, Gammat, Gamtd);
    vec3<mat3<CCTK_REAL> > Gammal;
    vec3<mat3<CCTK_REAL> > Gamma;
    vec3<CCTK_REAL> Gam;
    calc_gamma(g, gu, dg, Gammal, Gamma, Gam);

  apply_diss(cctkGH, gf_Kh_, gf_Kh_rhs_);

    mat3<CCTK_REAL> DDchi;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b) {
        CCTK_REAL s = ddchi(a, b);
        for (int x = 0; x < 3; ++x)
          s -= Gammat(x)(a, b) * dchi(x);
        DDchi(a, b) = s;
      }
    const vec3<vec3<CCTK_REAL> > dGamt{
        deriv(gf_Gamtx_, p.I, dx),
        deriv(gf_Gamty_, p.I, dx),
        deriv(gf_Gamtz_, p.I, dx),
    };
    const mat3<CCTK_REAL> R =
        calc_ricci(chi, dchi, DDchi, gammat, gammatu, ddgammat, Gammatl, Gammat,
                   Gamtd, dGamt);
    CCTK_REAL Rsc = 0;
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        Rsc += gu(x, y) * R(x, y);
    const vec3<CCTK_REAL> dalphaG = deriv(gf_alphaG_, p.I, dx);
    const mat3<CCTK_REAL> ddalphaG = deriv2(gf_alphaG_, p.I, dx);
    const vec3<CCTK_REAL> dTheta = deriv(gf_Theta_, p.I, dx);
    mat3<CCTK_REAL> DDalphaG;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b) {
        CCTK_REAL s = ddalphaG(a, b);
        for (int x = 0; x < 3; ++x)
          s -= Gammat(x)(a, b) * dalphaG(x);
        DDalphaG(a, b) = s;
      }
    const vec3<vec3<CCTK_REAL> > dbetaG{
        deriv(gf_betaGx_, p.I, dx),
        deriv(gf_betaGy_, p.I, dx),
        deriv(gf_betaGz_, p.I, dx),
    };
    const vec3<mat3<CCTK_REAL> > ddbetaG{
        deriv2(gf_betaGx_, p.I, dx),
        deriv2(gf_betaGy_, p.I, dx),
        deriv2(gf_betaGz_, p.I, dx),
    };

  apply_diss(cctkGH, gf_Atxx_, gf_Atxx_rhs_);
  apply_diss(cctkGH, gf_Atxy_, gf_Atxy_rhs_);
  apply_diss(cctkGH, gf_Atxz_, gf_Atxz_rhs_);
  apply_diss(cctkGH, gf_Atyy_, gf_Atyy_rhs_);
  apply_diss(cctkGH, gf_Atyz_, gf_Atyz_rhs_);
  apply_diss(cctkGH, gf_Atzz_, gf_Atzz_rhs_);

    // (1)
    CCTK_REAL chi_rhs = 0;
    chi_rhs += alphaG * (Kh + 2 * Theta);
    for (int x = 0; x < 3; ++x)
      chi_rhs -= dbetaG(x)(x);
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        chi_rhs -= Gamma(x)(x, y) * betaG(y);
    chi_rhs *= chi * 2 / 3;

  apply_diss(cctkGH, gf_Gamtx_, gf_Gamtx_rhs_);
  apply_diss(cctkGH, gf_Gamty_, gf_Gamty_rhs_);
  apply_diss(cctkGH, gf_Gamtz_, gf_Gamtz_rhs_);

    // (2)
    mat3<CCTK_REAL> gammat_rhs;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b) {
        CCTK_REAL s = 0;
        s += -2 * alphaG * At(a, b);
        for (int x = 0; x < 3; ++x)
          s += betaG(x) * dgammat(a, b)(x);
        for (int x = 0; x < 3; ++x)
          s += 2 * (gammat(x, a) * dbetaG(x)(b) + gammat(x, b) * dbetaG(x)(a));
        for (int x = 0; x < 3; ++x)
          s -= gammat(a, b) * dbetaG(x)(x) * 2 / 3;
        gammat_rhs(a, b) = s;
      }

  apply_diss(cctkGH, gf_Theta_, gf_Theta_rhs_);

    // (3)
    CCTK_REAL Kh_rhs = 0;
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        Kh_rhs -= gu(x, y) * DDalphaG(x, y);
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        Kh_rhs += alphaG * (At(x, y) * Atu(x, y) + pow2(Kh + 2 * Theta) / 3);
    Kh_rhs += 4 * M_PI * alphaG * (traceT + rho);
    Kh_rhs += alphaG * kappa1 * (1 - kappa2) * Theta;
    for (int x = 0; x < 3; ++x)
      Kh_rhs += betaG(x) * dKh(x);

  apply_diss(cctkGH, gf_alphaG_, gf_alphaG_rhs_);

    // (4)
    mat3<CCTK_REAL> At_rhs1;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b) {
        CCTK_REAL s = 0;
        s -= DDalphaG(a, b) + alphaG * (R(a, b) - 8 * M_PI * T(a, b));
        At_rhs1(a, b) = s;
      }
    CCTK_REAL trace_At_rhs1 = 0;
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        trace_At_rhs1 += gammatu(x, y) * At_rhs1(x, y);
    mat3<CCTK_REAL> At_rhs;
    for (int a = 0; a < 3; ++a)
      for (int b = a; b < 3; ++b) {
        CCTK_REAL s = 0;
        s += chi * (At_rhs1(a, b) - 1 / 3 * trace_At_rhs1 * gammat(a, b));
        s += alphaG * (Kh + 2 * Theta) * At(a, b);
        for (int x = 0; x < 3; ++x)
          for (int y = 0; y < 3; ++y)
            s -= alphaG * 2 * gammatu(x, y) * At(x, a) * At(y, b);
        for (int x = 0; x < 3; ++x)
          s += betaG(x) * dAt(a, b)(x);
        for (int x = 0; x < 3; ++x)
          s += 2 * (At(x, a) * dbetaG(x)(b) + At(x, b) * dbetaG(x)(a));
        for (int x = 0; x < 3; ++x)
          s -= At(a, b) * dbetaG(x)(x) * 2 / 3;
        At_rhs(a, b) = s;
      }
    // (5)
    vec3<CCTK_REAL> Gamt_rhs;
    for (int a = 0; a < 3; ++a) {
      CCTK_REAL s = 0;
      for (int x = 0; x < 3; ++x)
        s -= 2 * Atu(a, x) * dalphaG(x);
      for (int x = 0; x < 3; ++x)
        for (int y = 0; y < 3; ++y)
          s += 2 * alphaG * Gammat(a)(x, y) * Atu(x, y);
      for (int x = 0; x < 3; ++x)
        s -= 2 * alphaG * Atu(a, x) * dchi(x) / chi * 3 / 2;
      for (int x = 0; x < 3; ++x)
        s -= 2 * alphaG * gammatu(a, x) * (2 * dKh(x) + dTheta(x)) / 3;
      for (int x = 0; x < 3; ++x)
        s -= 2 * alphaG * 8 * M_PI * gammatu(a, x) * S(x);
      for (int x = 0; x < 3; ++x)
        for (int y = 0; y < 3; ++y)
          s += gammatu(x, y) * ddbetaG(a)(x, y);
      for (int x = 0; x < 3; ++x)
        for (int y = 0; y < 3; ++y)
          s += gammatu(a, x) * ddbetaG(y)(x, y) / 3;
      for (int x = 0; x < 3; ++x)
        s += betaG(x) * dGamt(a)(x);
      for (int x = 0; x < 3; ++x)
        s -= Gamtd(x) * dbetaG(a)(x);
      for (int x = 0; x < 3; ++x)
        s += Gamtd(a) * dbetaG(x)(x) * 2 / 3;
      s -= 2 * alphaG * kappa1 * (Gamt(a) - Gamtd(a));
      Gamt_rhs(a) = s;
    }
    // (6)
    CCTK_REAL Theta_rhs = 0;
    Theta_rhs += alphaG * Rsc / 2;
    for (int x = 0; x < 3; ++x)
      for (int y = 0; y < 3; ++y)
        Theta_rhs -= alphaG * At(x, y) * Atu(x, y) / 2;
    Theta_rhs += alphaG * pow2(Kh + 2 * Theta) / 3;
    Theta_rhs -= alphaG * 8 * M_PI * rho;
    Theta_rhs -= alphaG * kappa1 * (2 + kappa2) * Theta;
    for (int x = 0; x < 3; ++x)
      Theta_rhs += betaG(x) * dTheta(x);
    const CCTK_REAL mu_L = f_mu_L / alphaG;
    CCTK_REAL alphaG_rhs = 0;
    alphaG_rhs -= pow2(alphaG) * mu_L * Kh;
    for (int x = 0; x < 3; ++x)
      alphaG_rhs += betaG(x) * dalphaG(x);
    const CCTK_REAL mu_S = f_mu_S / pow2(alphaG);
    vec3<CCTK_REAL> betaG_rhs;
    for (int a = 0; a < 3; ++a) {
      CCTK_REAL s = 0;
      s += pow2(alphaG) * mu_S * Gamt(a) - eta * betaG(a);
      for (int x = 0; x < 3; ++x)
        s += betaG(x) * dbetaG(a)(x);
      betaG_rhs(a) = s;
    }
    // Store
    gf_chi_rhs_(p.I) = chi_rhs;
    gammat_rhs.store(gf_gammatxx_rhs_, gf_gammatxy_rhs_, gf_gammatxz_rhs_,
                     gf_gammatyy_rhs_, gf_gammatyz_rhs_, gf_gammatzz_rhs_, p);
    gf_Kh_rhs_(p.I) = Kh_rhs;
    At_rhs.store(gf_Atxx_rhs_, gf_Atxy_rhs_, gf_Atxz_rhs_, gf_Atyy_rhs_,
                 gf_Atyz_rhs_, gf_Atzz_rhs_, p);
    Gamt_rhs.store(gf_Gamtx_rhs_, gf_Gamty_rhs_, gf_Gamtz_rhs_, p);
    gf_Theta_rhs_(p.I) = Theta_rhs;
    gf_alphaG_rhs_(p.I) = alphaG_rhs;
    betaG_rhs.store(gf_betaGx_rhs_, gf_betaGy_rhs_, gf_betaGz_rhs_, p);
  });

  apply_diss(cctkGH, gf_betaGx_, gf_betaGx_rhs_);
  apply_diss(cctkGH, gf_betaGy_, gf_betaGy_rhs_);
  apply_diss(cctkGH, gf_betaGz_, gf_betaGz_rhs_);

#include <loop.hxx>
#include <cctk.h>
#include <cctk_Arguments_Checked.h>
namespace Z4c {
using namespace Loop;
extern "C" void Z4c_RHSBoundaries(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS_Z4c_RHSBoundaries;
  const GF3D<CCTK_REAL, 0, 0, 0> gf_chi_rhs_(cctkGH, chi_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxx_rhs_(cctkGH, gammatxx_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxy_rhs_(cctkGH, gammatxy_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxz_rhs_(cctkGH, gammatxz_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatyy_rhs_(cctkGH, gammatyy_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatyz_rhs_(cctkGH, gammatyz_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatzz_rhs_(cctkGH, gammatzz_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Kh_rhs_(cctkGH, Kh_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxx_rhs_(cctkGH, Atxx_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxy_rhs_(cctkGH, Atxy_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxz_rhs_(cctkGH, Atxz_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Atyy_rhs_(cctkGH, Atyy_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Atyz_rhs_(cctkGH, Atyz_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Atzz_rhs_(cctkGH, Atzz_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Gamtx_rhs_(cctkGH, Gamtx_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Gamty_rhs_(cctkGH, Gamty_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Gamtz_rhs_(cctkGH, Gamtz_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Theta_rhs_(cctkGH, Theta_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_alphaG_rhs_(cctkGH, alphaG_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_betaGx_rhs_(cctkGH, betaGx_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_betaGy_rhs_(cctkGH, betaGy_rhs);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_betaGz_rhs_(cctkGH, betaGz_rhs);
  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_chi_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_gammatxx_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_gammatxy_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_gammatxz_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_gammatyy_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_gammatyz_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_gammatzz_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Kh_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atxx_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atxy_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atxz_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atyy_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atyz_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atzz_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_Gamtx_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_Gamty_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_Gamtz_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_Theta_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_alphaG_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_betaGx_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_betaGy_rhs_(p.I) = 0; });
  loop_bnd<0, 0, 0>(cctkGH,
                    [&](const PointDesc &p) { gf_betaGz_rhs_(p.I) = 0; });
}
} // namespace Z4c

#include <cctk.h>

#include <sstream>
#include <type_traits>

template <typename T> constexpr T pow2(const T x) { return x * x; }

template <typename T>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T pow2(const T x) {
  return x * x;
}
template <typename T>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T pow3(const T x) {
  return pow2(x) * x;
}
template <typename T>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T pow4(const T x) {
  return pow2(pow2(x));
}

  constexpr T operator()() const { return 0.0 / 0.0; }

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T operator()() const {
    return NAN;
  }

  constexpr vect<T, D> operator()() const {

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vect<T, D> operator()() const {

template <typename T> struct norm1 {
  typedef T result_type;
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_type
  operator()(const T &x) const {
    return abs(x);
  }
};
template <typename T, int D> struct norm1<vect<T, D> > {
  typedef typename norm1<T>::result_type result_type;
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_type
  operator()(const vect<T, D> &xs) const {
    typename norm1<T>::result_type r{0};
    for (int d = 0; d < D; ++d)
      r = max(r, norm1<T>()(xs[d]));
    return r;
  }
};

  static constexpr int ind(const int n) {

  static /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr int ind(const int n) {
#ifdef CCTK_DEBUG

#endif

  explicit constexpr vec3() : elts{nan<vect<T, 3> >()()} {}

  explicit /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3()
      : elts{nan<vect<T, 3> >()()} {}

  constexpr vec3(initializer_list<T> v) : elts(v) {}

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3(const vect<T, 3> &elts)
      : elts(elts) {}
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3(vect<T, 3> &&elts)
      : elts(move(elts)) {}

  vec3(const GF3D<const T, 0, 0, 0> &gf_vx_,
       const GF3D<const T, 0, 0, 0> &gf_vy_,
       const GF3D<const T, 0, 0, 0> &gf_vz_, const PointDesc &p)
      : vec3{gf_vx_(p.I), gf_vy_(p.I), gf_vz_(p.I)} {}

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3(initializer_list<T> v)
      : elts(v) {}
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ vec3(const GF3D<const T, 0, 0, 0> &gf_vx_,
                                        const GF3D<const T, 0, 0, 0> &gf_vy_,
                                        const GF3D<const T, 0, 0, 0> &gf_vz_,
                                        const vect<int, 3> &I)
      : vec3{gf_vx_(I), gf_vy_(I), gf_vz_(I)} {}

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/

       const GF3D<T, 0, 0, 0> &gf_vz_, const PointDesc &p)
      : vec3{gf_vx_(p.I), gf_vy_(p.I), gf_vz_(p.I)} {}

       const GF3D<T, 0, 0, 0> &gf_vz_, const vect<int, 3> &I)
      : vec3{gf_vx_(I), gf_vy_(I), gf_vz_(I)} {}
  template <typename F, typename = result_of_t<F(int)> >
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3(const F &f)
      : elts{f(0), f(1), f(2)} {}

  void store(const GF3D<T, 0, 0, 0> &gf_vx_, const GF3D<T, 0, 0, 0> &gf_vy_,
             const GF3D<T, 0, 0, 0> &gf_vz_, const PointDesc &p) const {

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ void store(const GF3D<T, 0, 0, 0> &gf_vx_,
                                              const GF3D<T, 0, 0, 0> &gf_vy_,
                                              const GF3D<T, 0, 0, 0> &gf_vz_,
                                              const PointDesc &p) const {

    if (!((!CCTK_isnan(v(0))) && (!CCTK_isnan(v(1))) && (!CCTK_isnan(v(2))))) {
      ostringstream buf;
      buf << "v=" << v;
      CCTK_VERROR("nan found: %s", buf.str().c_str());
    }

  const T &operator()(int i) const { return elts[ind(i)]; }

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ const T &operator()(int i) const {
    return elts[ind(i)];
  }

  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<T>
  operator+(const vec3<T> &x) {
    return {+x.elts};
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<T>
  operator-(const vec3<T> &x) {
    return {-x.elts};
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<T>
  operator+(const vec3<T> &x, const vec3<T> &y) {
    return {x.elts + y.elts};
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<T>
  operator-(const vec3<T> &x, const vec3<T> &y) {
    return {x.elts - y.elts};
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<T>
  operator*(const T &a, const vec3<T> &x) {
    return {a * x.elts};
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr bool
  operator==(const vec3<T> &x, const vec3<T> &y) {
    return equal_to<vect<T, 3> >()(x.elts, y.elts);
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr bool
  operator!=(const vec3<T> &x, const vec3<T> &y) {
    return !(x == y);
  }

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T maxabs() const {
    return elts.maxabs();
  }
  friend struct norm1<vec3>;

  constexpr vec3<T> operator()() const { return vec3<T>(); }

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<T> operator()() const {
    return vec3<T>();
  }
};
template <typename T> struct norm1<vec3<T> > {
  typedef typename norm1<vect<T, 3> >::result_type result_type;
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_type
  operator()(const vec3<T> &x) const {
    return norm1<vect<T, 3> >()(x.elts);
  }

  static constexpr int symind(const int i, const int j) {

  static /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr int symind(const int i,
                                                               const int j) {
#ifdef CCTK_DEBUG

#endif

    // i j n
    // 0 0 0
    // 0 1 1
    // 0 2 2
    // 1 1 3
    // 1 2 4
    // 2 2 5
    // const int n = 2 * i + j - (unsigned)i / 2;
#ifdef CCTK_DEBUG

#endif

  static constexpr int ind(const int i, const int j) {

  static /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr int ind(const int i,
                                                            const int j) {

  explicit constexpr mat3() : elts{nan<vect<T, 6> >()()} {}

  explicit /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3()
      : elts{nan<vect<T, 6> >()()} {}

  constexpr mat3(initializer_list<T> A) : elts(A) {}

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3(const vect<T, 6> &elts)
      : elts(elts) {}
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3(vect<T, 6> &&elts)
      : elts(move(elts)) {}

  mat3(const GF3D<const T, 0, 0, 0> &gf_Axx_,
       const GF3D<const T, 0, 0, 0> &gf_Axy_,
       const GF3D<const T, 0, 0, 0> &gf_Axz_,
       const GF3D<const T, 0, 0, 0> &gf_Ayy_,
       const GF3D<const T, 0, 0, 0> &gf_Ayz_,
       const GF3D<const T, 0, 0, 0> &gf_Azz_, const PointDesc &p)
      : mat3{gf_Axx_(p.I), gf_Axy_(p.I), gf_Axz_(p.I),
             gf_Ayy_(p.I), gf_Ayz_(p.I), gf_Azz_(p.I)} {}

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3(initializer_list<T> A)
      : elts(A) {}
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ mat3(const GF3D<const T, 0, 0, 0> &gf_Axx_,
                                        const GF3D<const T, 0, 0, 0> &gf_Axy_,
                                        const GF3D<const T, 0, 0, 0> &gf_Axz_,
                                        const GF3D<const T, 0, 0, 0> &gf_Ayy_,
                                        const GF3D<const T, 0, 0, 0> &gf_Ayz_,
                                        const GF3D<const T, 0, 0, 0> &gf_Azz_,
                                        const vect<int, 3> &I)
      : mat3{gf_Axx_(I), gf_Axy_(I), gf_Axz_(I),
             gf_Ayy_(I), gf_Ayz_(I), gf_Azz_(I)} {}

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/

       const PointDesc &p)
      : mat3{gf_Axx_(p.I), gf_Axy_(p.I), gf_Axz_(p.I),
             gf_Ayy_(p.I), gf_Ayz_(p.I), gf_Azz_(p.I)} {}

       const vect<int, 3> &I)
      : mat3{gf_Axx_(I), gf_Axy_(I), gf_Axz_(I),
             gf_Ayy_(I), gf_Ayz_(I), gf_Azz_(I)} {}
  template <typename F, typename = result_of_t<F(int, int)> >
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3(const F &f)
      : elts{f(0, 0), f(0, 1), f(0, 2), f(1, 1), f(1, 2), f(2, 2)} {
#ifdef CCTK_DEBUG
    // Check symmetry
    const T f10 = f(1, 0);
    const T f20 = f(0, 2);
    const T f21 = f(1, 2);
    const auto is_symmetric{[](const T &fgood, const T &fother) {
      return norm1<T>()(fother - fgood) <=
             1.0e-12 * (1 + norm1<T>()(fgood) + norm1<T>()(fother));
    }};
    if (!(is_symmetric(f(0, 1), f10) && is_symmetric(f(0, 2), f20) &&
          is_symmetric(f(1, 2), f21))) {
      ostringstream buf;
      buf << "f(0,1)=" << f(0, 1) << "\n"
          << "f(1,0)=" << f10 << "\n"
          << "f(0,2)=" << f(0, 2) << "\n"
          << "f(2,0)=" << f20 << "\n"
          << "f(1,2)=" << f(1, 2) << "\n"
          << "f(2,1)=" << f21 << "\n";
      CCTK_VERROR("symmetric matrix is not symmetric:\n%s", buf.str().c_str());
    }
    assert(is_symmetric(f(0, 1), f10));
    assert(is_symmetric(f(0, 2), f20));
    assert(is_symmetric(f(1, 2), f21));
#endif
  }

  void store(const GF3D<T, 0, 0, 0> &gf_Axx_, const GF3D<T, 0, 0, 0> &gf_Axy_,
             const GF3D<T, 0, 0, 0> &gf_Axz_, const GF3D<T, 0, 0, 0> &gf_Ayy_,
             const GF3D<T, 0, 0, 0> &gf_Ayz_, const GF3D<T, 0, 0, 0> &gf_Azz_,
             const PointDesc &p) const {

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ void
  store(const GF3D<T, 0, 0, 0> &gf_Axx_, const GF3D<T, 0, 0, 0> &gf_Axy_,
        const GF3D<T, 0, 0, 0> &gf_Axz_, const GF3D<T, 0, 0, 0> &gf_Ayy_,
        const GF3D<T, 0, 0, 0> &gf_Ayz_, const GF3D<T, 0, 0, 0> &gf_Azz_,
        const PointDesc &p) const {

  const T &operator()(int i, int j) const { return elts[ind(i, j)]; }
  // T &operator()(int i, int j) { return elts[symind(i, j)]; }
  T &operator()(int i, int j) { return elts[ind(i, j)]; }

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ const T &operator()(int i, int j) const {
    return elts[ind(i, j)];
  }
  // /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T &operator()(int i, int j) { return
  // elts[symind(i, j)]; }
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T &operator()(int i, int j) {
    return elts[ind(i, j)];
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T>
  operator+(const mat3<T> &x) {
    return {+x.elts};
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T>
  operator-(const mat3<T> &x) {
    return {-x.elts};
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T>
  operator+(const mat3<T> &x, const mat3<T> &y) {
    return {x.elts + y.elts};
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T>
  operator-(const mat3<T> &x, const mat3<T> &y) {
    return {x.elts - y.elts};
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T>
  operator*(const T &a, const mat3<T> &x) {
    return {a * x.elts};
  }

  constexpr T det() const {

  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr bool
  operator==(const mat3<T> &x, const mat3<T> &y) {
    return equal_to<vect<T, 6> >()(x.elts, y.elts);
  }
  friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr bool
  operator!=(const mat3<T> &x, const mat3<T> &y) {
    return !(x == y);
  }
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T maxabs() const {
    return elts.maxabs();
  }
  friend struct norm1<mat3>;
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T det() const {

  constexpr mat3 inv(const T detA) const {

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3 inv(const T detA) const {

    return mat3{detA1 * (A(1, 1) * A(2, 2) - A(2, 1) * A(2, 1)),
                detA1 * (A(1, 2) * A(2, 0) - A(2, 2) * A(0, 1)),
                detA1 * (A(1, 0) * A(2, 1) - A(2, 0) * A(1, 1)),
                detA1 * (A(2, 2) * A(0, 0) - A(0, 2) * A(2, 0)),
                detA1 * (A(2, 0) * A(0, 1) - A(0, 0) * A(2, 1)),

    return mat3{detA1 * (A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1)),
                detA1 * (A(1, 2) * A(2, 0) - A(1, 0) * A(2, 2)),
                detA1 * (A(1, 0) * A(2, 1) - A(1, 1) * A(2, 0)),
                detA1 * (A(2, 2) * A(0, 0) - A(2, 0) * A(0, 2)),
                detA1 * (A(2, 0) * A(0, 1) - A(2, 1) * A(0, 0)),

  constexpr T trace() const {

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T trace(const mat3 &gu) const {
    const auto &A = *this;
    return sum2([&](int x, int y) { return gu(x, y) * A(x, y); });
  }
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3
  trace_free(const mat3 &g, const mat3 &gu) const {

    return A(0, 0) + A(1, 1) + A(2, 2);

    const T trA = A.trace(gu);
    return mat3([&](int a, int b) { return A(a, b) - trA / 3 * g(a, b); });

  constexpr mat3<T> operator()() const { return mat3<T>(); }

  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T> operator()() const {
    return mat3<T>();
  }
};
template <typename T> struct norm1<mat3<T> > {
  typedef typename norm1<vect<T, 6> >::result_type result_type;
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_type
  operator()(const mat3<T> &x) const {
    return norm1<vect<T, 6> >()(x.elts);
  }

constexpr mat3<T> mul(const mat3<T> &A, const mat3<T> &B) {
  mat3<T> C;

/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T> mul(const mat3<T> &A,
                                                       const mat3<T> &B) {

  for (int a = 0; a < 3; ++a)
    for (int b = a; b < 3; ++b) {
      T s = 0;
      for (int c = 0; c < 3; ++c)
        s += A(a, c) * B(c, b);
      C(a, b) = s;
    }
  return C;

  return mat3<T>([&](int a, int b) {
    return sum1([&](int x) { return A(a, x) * B(x, b); });
  });
}
template <typename T>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T add(const T &x) {
  return x;
}
template <typename T, typename... Ts>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr T add(const T &x, const Ts &... xs) {
  return x + add(xs...);
}
template <typename F>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_of_t<F(int)>
sum1(const F &f) {
  result_of_t<F(int)> s{0};
  for (int x = 0; x < 3; ++x)
    s += f(x);
  return s;
}
template <typename F>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_of_t<F(int, int)>
sum2(const F &f) {
  result_of_t<F(int, int)> s{0};
  for (int x = 0; x < 3; ++x)
    for (int y = 0; y < 3; ++y)
      s += f(x, y);
  return s;
}
template <typename F>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr result_of_t<F(int, int, int)>
sum3(const F &f) {
  result_of_t<F(int, int, int)> s{0};
  for (int x = 0; x < 3; ++x)
    for (int y = 0; y < 3; ++y)
      for (int z = 0; z < 3; ++z)
        s += f(x, y, z);
  return s;

T deriv(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
        const vec3<T> &dx, const int dir) {

/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T deriv(const GF3D<const T, 0, 0, 0> &gf_,
                                         const vect<int, dim> &I,
                                         const vec3<T> &dx, const int dir) {

vec3<T> deriv(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
              const vec3<T> &dx) {
  return {
      deriv(gf_, I, dx, 0),
      deriv(gf_, I, dx, 1),
      deriv(gf_, I, dx, 2),
  };

/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T
deriv_upwind(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
             const bool sign, const vec3<T> &dx, const int dir) {
  constexpr vect<vect<int, dim>, dim> DI{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
  // arXiv:1111.2177 [gr-qc], (71)
  if (sign)
    // +     [ 0   -1   +1    0    0]
    // + 1/2 [+1   -2   +1    0    0]
    //       [+1/2 -2   +3/2  0    0]
    return (+1 / T(2) * gf_(I - 2 * DI[dir]) //
            - 2 * gf_(I - DI[dir])           //
            + 3 / T(2) * gf_(I)) /
           dx(dir);
  else
    // +     [ 0    0   -1   +1    0  ]
    // - 1/2 [ 0    0   +1   -2   +1  ]
    //       [ 0    0   -3/2 +2   -1/2]
    return (-3 / T(2) * gf_(I)     //
            + 2 * gf_(I + DI[dir]) //
            - 1 / T(2) * gf_(I + 2 * DI[dir])) /
           dx(dir);

T deriv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
         const vec3<T> &dx, const int dir1, const int dir2) {

/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T
deriv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
       const vec3<T> &dx, const int dir1, const int dir2) {

    return (gf_(I + DI[dir1]) - 2 * gf_(I) + gf_(I - DI[dir1])) /

    return ((gf_(I - DI[dir1]) + gf_(I + DI[dir1])) - 2 * gf_(I)) /

mat3<T> deriv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
               const vec3<T> &dx) {

/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T deriv4(const GF3D<const T, 0, 0, 0> &gf_,
                                          const vect<int, dim> &I,
                                          const vec3<T> &dx, const int dir) {
  constexpr vect<vect<int, dim>, dim> DI{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
  return ((gf_(I - 2 * DI[dir]) + gf_(I + 2 * DI[dir])) //
          - 4 * (gf_(I - DI[dir]) + gf_(I + DI[dir]))   //
          + 6 * gf_(I)) /
         pow4(dx(dir));
}
template <typename T>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ vec3<T>
deriv(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
      const vec3<T> &dx) {
  return {
      deriv(gf_, I, dx, 0),
      deriv(gf_, I, dx, 1),
      deriv(gf_, I, dx, 2),
  };
}
template <typename T>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ vec3<T>
deriv_upwind(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
             const vec3<T> &dir, const vec3<T> &dx) {
  return {
      deriv_upwind(gf_, I, signbit(dir(0)), dx, 0),
      deriv_upwind(gf_, I, signbit(dir(1)), dx, 1),
      deriv_upwind(gf_, I, signbit(dir(2)), dx, 2),
  };
}
template <typename T>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ mat3<T>
deriv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
       const vec3<T> &dx) {

template <typename T>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T diss(const GF3D<const T, 0, 0, 0> &gf_,
                                        const vect<int, dim> &I,
                                        const vec3<T> &dx) {
  // arXiv:gr-qc/0610128, (63), with r=2
  return -1 / T(16) *
         (pow3(dx(0)) * deriv4(gf_, I, dx, 0)   //
          + pow3(dx(1)) * deriv4(gf_, I, dx, 1) //
          + pow3(dx(2)) * deriv4(gf_, I, dx, 2));
}

#warning "TODO"
#include <iostream>
#include "physics.hxx"
#include "tensor.hxx"
#include <cctk.h>
#include <algorithm>
#include <array>
#include <cassert>
#include <random>
namespace Z4c {
using namespace std;
// TODO: Use GoogleTest instead of assert
extern "C" void Z4c_Test(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  mt19937 engine(42);
  uniform_int_distribution<int> dist(-10, 10);
  const auto rand10{[&]() { return double(dist(engine)); }};
  const auto randmat10{[&]() {
    array<array<double, 3>, 3> arr;
    for (int a = 0; a < 3; ++a)
      for (int b = 0; b < 3; ++b)
        arr[a][b] = rand10();
    return mat3<double>(
        [&](int a, int b) { return arr[min(a, b)][max(a, b)]; });
  }};
  const mat3<double> Z([&](int a, int b) { return double(0); });
  const mat3<double> I([&](int a, int b) { return double(a == b); });
  assert(I != Z);
  for (int n = 0; n < 100; ++n) {
    const mat3<double> A = randmat10();
    const mat3<double> B = randmat10();
    const mat3<double> C = randmat10();
    const double a = rand10();
    const double b = rand10();
    assert((A + B) + C == A + (B + C));
    assert(Z + A == A);
    assert(A + Z == A);
    assert(A + (-A) == Z);
    assert((-A) + A == Z);
    assert(A - B == A + (-B));
    assert(A + B == B + A);
    assert(1 * A == A);
    assert(0 * A == Z);
    assert(-1 * A == -A);
    // assert(mul(a * A, B) == a * mul(A, B));
    assert((a * b) * A == a * (b * A));
    assert(a * (A + B) == a * A + a * B);
    assert((a + b) * A == a * A + b * A);
    // assert(mul(mul(A, B), C) == mul(A, mul(B, C)));
    assert(mul(I, A) == A);
    assert(mul(A, I) == A);
    assert(mul(Z, A) == Z);
    assert(mul(A, Z) == Z);
    assert(Z.det() == 0);
    assert(I.det() == 1);
    assert((a * A).det() == pow3(a) * A.det());
    assert(Z.inv(1) == Z);
    assert(I.inv(1) == I);
    assert(mul(A.inv(1), A) == A.det() * I);
    assert(mul(A, A.inv(1)) == A.det() * I);
    assert((a * A).inv(1) == pow2(a) * A.inv(1));
  }
}
} // namespace Z4c

#ifndef Z4C_HXX
#define Z4C_HXX
#include "physics.hxx"
#include "tensor.hxx"
#include <cmath>
namespace Z4c {
// See arXiv:1212.2901 [gr-qc]
// Note: A_(ij) there means 1/2 (A_ij + A_ji)
template <typename T> struct z4c_vars_noderivs {
  // Parameters
  const T kappa1;
  const T kappa2;
  const T f_mu_L;
  const T f_mu_S;
  const T eta;
  // Tensor densities: TD = (sqrt detg)^W T
  //                   W[sqrt detg] = +1
  // Z4c variables
  const T chi;          // W = -2/3
  const mat3<T> gammat; // W = -2/3
  const T Kh;           // W = 0
  const mat3<T> At;     // W = -2/3
  const vec3<T> Gamt;   // W = +2/3
  const T Theta;        // W = 0
  const T alphaG;       // W = 0
  const vec3<T> betaG;  // W = 0
  // Hydro variables
  const T rho;
  const vec3<T> Si;
  const mat3<T> Sij;
  // ADM variables
  const mat3<T> g;    // W = 0
  const mat3<T> k;    // W = 0
  const T alp;        // W = 0
  const vec3<T> beta; // W = 0
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/
  z4c_vars_noderivs(const T &kappa1, const T &kappa2, const T &f_mu_L,
                    const T &f_mu_S, const T &eta,
                    //
                    const T &chi, const mat3<T> &gammat, const T &Kh,
                    const mat3<T> &At, const vec3<T> &Gamt, const T &Theta,
                    const T &alphaG, const vec3<T> &betaG,
                    //
                    const T &rho, const vec3<T> &Si, const mat3<T> &Sij)
      : kappa1(kappa1), kappa2(kappa2), f_mu_L(f_mu_L), f_mu_S(f_mu_S),
        eta(eta),
        //
        chi(chi), gammat(gammat), Kh(Kh), At(At), Gamt(Gamt), Theta(Theta),
        alphaG(alphaG), betaG(betaG),
        //
        rho(rho), Si(Si), Sij(Sij),
        // ADM variables
        g([&](int a, int b) { return 1 / chi * gammat(a, b); }), //
        k([&](int a, int b) {
          return 1 / chi * (At(a, b) + (Kh + 2 * Theta) / 3 * gammat(a, b));
        }),          //
        alp(alphaG), //
        beta([&](int a) { return betaG(a); })
  //
  {}
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/
  z4c_vars_noderivs(const T &kappa1, const T &kappa2, const T &f_mu_L,
                    const T &f_mu_SL, const T &eta,
                    //
                    const GF3D<const T, 0, 0, 0> &gf_chi_,
                    //
                    const GF3D<const T, 0, 0, 0> &gf_gammatxx_,
                    const GF3D<const T, 0, 0, 0> &gf_gammatxy_,
                    const GF3D<const T, 0, 0, 0> &gf_gammatxz_,
                    const GF3D<const T, 0, 0, 0> &gf_gammatyy_,
                    const GF3D<const T, 0, 0, 0> &gf_gammatyz_,
                    const GF3D<const T, 0, 0, 0> &gf_gammatzz_,
                    //
                    const GF3D<const T, 0, 0, 0> &gf_Kh_,
                    //
                    const GF3D<const T, 0, 0, 0> &gf_Atxx_,
                    const GF3D<const T, 0, 0, 0> &gf_Atxy_,
                    const GF3D<const T, 0, 0, 0> &gf_Atxz_,
                    const GF3D<const T, 0, 0, 0> &gf_Atyy_,
                    const GF3D<const T, 0, 0, 0> &gf_Atyz_,
                    const GF3D<const T, 0, 0, 0> &gf_Atzz_,
                    //
                    const GF3D<const T, 0, 0, 0> &gf_Gamtx_,
                    const GF3D<const T, 0, 0, 0> &gf_Gamty_,
                    const GF3D<const T, 0, 0, 0> &gf_Gamtz_,
                    //
                    const GF3D<const T, 0, 0, 0> &gf_Theta_,
                    //
                    const GF3D<const T, 0, 0, 0> &gf_alphaG_,
                    //
                    const GF3D<const T, 0, 0, 0> &gf_betaGx_,
                    const GF3D<const T, 0, 0, 0> &gf_betaGy_,
                    const GF3D<const T, 0, 0, 0> &gf_betaGz_,
                    //
                    const vect<int, 3> &I)
      : z4c_vars_noderivs<T>(kappa1, kappa2, f_mu_L, f_mu_S, eta,
                             //
                             gf_chi_(I),
                             //
                             mat3<T>(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_,
                                     gf_gammatyy_, gf_gammatyz_, gf_gammatzz_,
                                     I),
                             //
                             gf_Kh_(I),
                             //
                             mat3<T>(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_,
                                     gf_Atyz_, gf_Atzz_, I),
                             //
                             vec3<T>(gf_Gamtx_, gf_Gamty_, gf_Gamtz_, I),
                             //
                             gf_Theta_(I),
                             //
                             gf_alphaG_(I),
                             //
                             vec3<T>(gf_betaGx_, gf_betaGy_, gf_betaGz_, I),
                             //
                             0,
                             //
                             vec3<T>{0, 0, 0},
                             //
                             mat3<T>{0, 0, 0, 0, 0, 0})
  //
  {}
};
template <typename T> struct z4c_vars : z4c_vars_noderivs<T> {
  // C++ is tedious:
  // Parameters
  using z4c_vars_noderivs<T>::kappa1;
  using z4c_vars_noderivs<T>::kappa2;
  using z4c_vars_noderivs<T>::f_mu_L;
  using z4c_vars_noderivs<T>::f_mu_S;
  using z4c_vars_noderivs<T>::eta;
  // Z4c variables
  using z4c_vars_noderivs<T>::chi;
  using z4c_vars_noderivs<T>::gammat;
  using z4c_vars_noderivs<T>::Kh;
  using z4c_vars_noderivs<T>::At;
  using z4c_vars_noderivs<T>::Gamt;
  using z4c_vars_noderivs<T>::Theta;
  using z4c_vars_noderivs<T>::alphaG;
  using z4c_vars_noderivs<T>::betaG;
  // Hydro variables
  using z4c_vars_noderivs<T>::rho;
  using z4c_vars_noderivs<T>::Si;
  using z4c_vars_noderivs<T>::Sij;
  // ADM variables
  using z4c_vars_noderivs<T>::g;
  using z4c_vars_noderivs<T>::k;
  using z4c_vars_noderivs<T>::alp;
  using z4c_vars_noderivs<T>::beta;
  // Derivatives of Z4c variables
  const vec3<T> dchi;
  const mat3<T> ddchi;
  const mat3<vec3<T> > dgammat;
  const mat3<mat3<T> > ddgammat;
  const vec3<T> dKh;
  const mat3<vec3<T> > dAt;
  const vec3<vec3<T> > dGamt;
  const vec3<T> dTheta;
  const vec3<T> dalphaG;
  const mat3<T> ddalphaG;
  const vec3<vec3<T> > dbetaG;
  const vec3<mat3<T> > ddbetaG;
  // Intermediate variables
  const mat3<T> gammatu;
  const mat3<vec3<T> > dgammatu;
  const vec3<mat3<T> > Gammatl;
  const vec3<mat3<T> > Gammat;
  const vec3<T> Gamtd;
  const mat3<T> DDchi;
  const mat3<T> gu;
  const mat3<vec3<T> > dg;
  const vec3<mat3<T> > Gammal;
  const vec3<mat3<T> > Gamma;
  const mat3<T> DDalphaG;
  const mat3<T> Rchi;
  const mat3<T> Rt;
  const mat3<T> R;
  const T Rsc;
  const mat3<T> Atu;
  const mat3<vec3<T> > dAtu;
  const T traceSij;
  // Constraints
  const vec3<T> ZtC;
  const T HC;
  const vec3<T> MtC;
  const T allC;
  // RHS variables
  const T chi_rhs;
  const mat3<T> gammat_rhs;
  const T Kh_rhs;
  const mat3<T> At_rhs;
  const vec3<T> Gamt_rhs;
  const T Theta_rhs;
  const T alphaG_rhs;
  const vec3<T> betaG_rhs;
  // See arXiv:1212.2901 [gr-qc]
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/
  z4c_vars(const T &kappa1, const T &kappa2, const T &f_mu_L, const T &f_mu_S,
           const T &eta,
           //
           const T &chi, const vec3<T> &dchi, const mat3<T> &ddchi, //
           const mat3<T> &gammat, const mat3<vec3<T> > &dgammat,
           const mat3<mat3<T> > &ddgammat,                                   //
           const T &Kh, const vec3<T> &dKh,                                  //
           const mat3<T> &At, const mat3<vec3<T> > &dAt,                     //
           const vec3<T> &Gamt, const vec3<vec3<T> > &dGamt,                 //
           const T &Theta, const vec3<T> &dTheta,                            //
           const T &alphaG, const vec3<T> &dalphaG, const mat3<T> &ddalphaG, //
           const vec3<T> &betaG, const vec3<vec3<T> > &dbetaG,
           const vec3<mat3<T> > &ddbetaG,
           //
           const T &rho, const vec3<T> &Si, const mat3<T> &Sij)
      : z4c_vars_noderivs<T>(kappa1, kappa2, f_mu_L, f_mu_S, eta, //
                             chi, gammat, Kh, At, Gamt, Theta, alphaG, betaG,
                             rho, Si, Sij),
        // Derivatives of Z4c variables
        dchi(dchi), ddchi(ddchi),             //
        dgammat(dgammat), ddgammat(ddgammat), //
        dKh(dKh),                             //
        dAt(dAt),                             //
        dGamt(dGamt),                         //
        dTheta(dTheta),                       //
        dalphaG(dalphaG), ddalphaG(ddalphaG), //
        dbetaG(dbetaG), ddbetaG(ddbetaG),
        // Intermediate variables
        gammatu(gammat.inv(1)),               //
        dgammatu(calc_dgu(gammatu, dgammat)), //
        Gammatl(calc_gammal(dgammat)),        //
        Gammat(calc_gamma(gammatu, Gammatl)), //
        Gamtd([&](int a) {
          return sum2(
              [&](int x, int y) { return gammatu(x, y) * Gammat(a)(x, y); });
        }), //
        DDchi([&](int a, int b) {
          return ddchi(a, b) //
                 - sum1([&](int x) { return Gammat(x)(a, b) * dchi(x); });
        }),                                                    //
        gu([&](int a, int b) { return chi * gammatu(a, b); }), //
        dg([&](int a, int b) {
          return vec3<T>([&](int c) {
            return -dchi(c) / pow2(chi) * gammat(a, b) //
                   + 1 / chi * dgammat(a, b)(c);
          });
        }),                            //
        Gammal(calc_gammal(dg)),       //
        Gamma(calc_gamma(gu, Gammal)), //
        DDalphaG([&](int a, int b) {
          return ddalphaG(a, b) //
                 - sum1([&](int x) { return Gamma(x)(a, b) * dalphaG(x); });
        }),
        // (8)
        Rchi([&](int a, int b) {
          return 1 / (2 * chi) * DDchi(a, b) //
                 + 1 / (2 * chi) * sum2([&](int x, int y) {
                     return gammat(a, b) * gammatu(x, y) * DDchi(x, y);
                   })                                      //
                 - 1 / (4 * pow2(chi)) * dchi(a) * dchi(b) //
                 - 3 / (4 * pow2(chi)) * sum2([&](int x, int y) {
                     return gammat(a, b) * gammatu(x, y) * dchi(x) * dchi(y);
                   });
        }),
        // (9)
        Rt([&](int a, int b) {
          return sum2([&](int x, int y) {
                   return -1 / T(2) * gammatu(x, y) * ddgammat(a, b)(x, y);
                 }) //
                 + sum1([&](int x) {
                     return 1 / T(2) *
                            (gammat(x, a) * dGamt(x)(b) //
                             + gammat(x, b) * dGamt(x)(a));
                   }) //
                 + sum1([&](int x) {
                     return 1 / T(2) *
                            (Gamtd(x) * Gammatl(a)(b, x) //
                             + Gamtd(x) * Gammatl(b)(a, x));
                   }) //
                 + sum3([&](int x, int y, int z) {
                     return gammatu(y, z) *
                            (Gammat(x)(a, y) * Gammatl(b)(x, z)   //
                             + Gammat(x)(b, y) * Gammatl(a)(x, z) //
                             + Gammat(x)(a, y) * Gammatl(x)(b, z));
                   });
        }),
        // (7)
        R([&](int a, int b) { return Rchi(a, b) + Rt(a, b); }), //
        Rsc(R.trace(gu)),                                       //
        Atu([&](int a, int b) {
          return sum2([&](int x, int y) {
            return gammatu(a, x) * gammatu(b, y) * At(x, y);
          });
        }),                                         //
        dAtu(calc_dAu(gammatu, dgammatu, At, dAt)), //
        traceSij(Sij.trace(gu)),
        // Constraints
        // (13)
        ZtC([&](int a) { return (Gamt(a) - Gamtd(a)) / 2; }), //
        // (14)
        HC(Rsc                                                        //
           + sum2([&](int x, int y) { return At(x, y) * Atu(x, y); }) //
           - 2 / T(3) * pow2(Kh + 2 * Theta)                          //
           - 16 * M_PI * rho),
        // (15)
        MtC([&](int a) {
          return sum1([&](int x) { return dAtu(a, x)(x); }) //
                 + sum2([&](int x, int y) {
                     return Gammat(a)(x, y) * Atu(x, y);
                   }) //
                 - sum1([&](int x) {
                     return 2 / T(3) * gammatu(a, x) * (dKh(x) + 2 * dTheta(x));
                   }) //
                 - sum1([&](int x) {
                     return 2 / T(3) * Atu(a, x) * dchi(x) / chi;
                   }) //
                 -
                 sum1([&](int x) { return 8 * M_PI * gammatu(a, x) * Si(x); });
        }),
        // arXiv:1111.2177, (73)
        allC(sqrt(pow2(HC) //
                  + sum2([&](int x, int y) {
                      return gammat(x, y) * MtC(x) * MtC(y);
                    })          //
                  + pow2(Theta) //
                  + sum2([&](int x, int y) {
                      return 2 * gammat(x, y) * ZtC(x) * ZtC(y);
                    }))),
        // RHS
        // (1)
        chi_rhs(2 / T(3) * chi *
                (alphaG * (Kh + 2 * Theta)
                 // chi = detg^(-1/3)
                 // chi^(-3) = detg
                 // chi^(-3/2) = sqrt detg
                 // D_i beta^i = 1/sqrt detg d_i sqrt detg beta^i
                 //            = chi^(3/2) d_i chi^(-3/2) beta^i
                 //            = chi^(3/2) (-3/2 chi^(-5/2) chi,i beta^i
                 //                         + chi^(-3/2) beta^i,i)
                 //            = beta^i,i - 3/2 1/chi beta^i chi,i
                 - sum1([&](int x) { return dbetaG(x)(x); }))),
        // (2)
        gammat_rhs([&](int a, int b) {
          return -2 * alphaG * At(a, b) //
                 + sum1([&](int x) {
                     return gammat(x, a) * dbetaG(x)(b) //
                            + gammat(x, b) * dbetaG(x)(a);
                   }) //
                 - sum1([&](int x) {
                     return 2 / T(3) * gammat(a, b) * dbetaG(x)(x);
                   });
        }),
        // (3)
        Kh_rhs(-sum2([&](int x, int y) { return gu(x, y) * DDalphaG(x, y); }) //
               + alphaG * (sum2([&](int x, int y) {
                             return At(x, y) * Atu(x, y);
                           })                                 //
                           + 1 / T(3) * pow2(Kh + 2 * Theta)) //
               + 4 * M_PI * alphaG * (traceSij + rho)         //
               + alphaG * kappa1 * (1 - kappa2) * Theta),
        // (4)
        At_rhs([&](int a, int b) {
          return chi * ((-DDalphaG(a, b)                      //
                         + alphaG * (R(a, b)                  //
                                     - 8 * M_PI * Sij(a, b))) //
                        - sum2([&](int x, int y) {
                            return 1 / T(3) * g(a, b) * gu(x, y) *
                                   (-DDalphaG(x, y)     //
                                    + alphaG * (R(x, y) //
                                                - 8 * M_PI * Sij(x, y)));
                          }))                            //
                 + alphaG * ((Kh + 2 * Theta) * At(a, b) //
                             - sum2([&](int x, int y) {
                                 return 2 * gammatu(x, y) * At(x, a) * At(y, b);
                               })) //
                 + sum1([&](int x) {
                     return At(x, a) * dbetaG(x)(b) //
                            + At(x, b) * dbetaG(x)(a);
                   }) //
                 - sum1([&](int x) {
                     return 2 / T(3) * At(a, b) * dbetaG(x)(x);
                   }); //
        }),
        // (5)
        Gamt_rhs([&](int a) {
          return -sum1([&](int x) { return 2 * Atu(a, x) * dalphaG(x); }) //
                 + 2 * alphaG *
                       (sum2([&](int x, int y) {
                          return Gammat(a)(x, y) * Atu(x, y);
                        }) //
                        - sum1([&](int x) {
                            return 3 / T(2) * Atu(a, x) * dchi(x) / chi;
                          }) //
                        - sum1([&](int x) {
                            return 1 / T(3) * gammatu(a, x) *
                                   (2 * dKh(x) + dTheta(x));
                          }) //
                        - sum1([&](int x) {
                            return 8 * M_PI * gammatu(a, x) * Si(x);
                          })) //
                 + sum2([&](int x, int y) {
                     return gammatu(x, y) * ddbetaG(a)(x, y);
                   }) //
                 + sum2([&](int x, int y) {
                     return 1 / T(3) * gammatu(a, x) * ddbetaG(y)(x, y);
                   })                                                   //
                 - sum1([&](int x) { return Gamtd(x) * dbetaG(a)(x); }) //
                 + sum1([&](int x) {
                     return 2 / T(3) * Gamtd(a) * dbetaG(x)(x);
                   }) //
                 - 2 * alphaG * kappa1 * (Gamt(a) - Gamtd(a));
        }),
        // (6)
        Theta_rhs(
            1 / T(2) * alphaG *
                (Rsc                                                        //
                 - sum2([&](int x, int y) { return At(x, y) * Atu(x, y); }) //
                 + 2 / T(3) * pow2(Kh + 2 * Theta))                         //
            - alphaG * (8 * M_PI * rho                                      //
                        + kappa1 * (2 + kappa2) * Theta)),
        // (11)
        alphaG_rhs([&] {
          const T mu_L = f_mu_L / alphaG;
          return -pow2(alphaG) * mu_L * Kh;
        }()),
        // (12)
        betaG_rhs([&](int a) {
          const T mu_S = f_mu_S / pow2(alphaG);
          return pow2(alphaG) * mu_S * Gamt(a) //
                 - eta * betaG(a);
        })
  //
  {}
  /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/
  z4c_vars(const T &kappa1, const T &kappa2, const T &f_mu_L, const T &f_mu_SL,
           const T &eta,
           //
           const GF3D<const T, 0, 0, 0> &gf_chi_,
           //
           const GF3D<const T, 0, 0, 0> &gf_gammatxx_,
           const GF3D<const T, 0, 0, 0> &gf_gammatxy_,
           const GF3D<const T, 0, 0, 0> &gf_gammatxz_,
           const GF3D<const T, 0, 0, 0> &gf_gammatyy_,
           const GF3D<const T, 0, 0, 0> &gf_gammatyz_,
           const GF3D<const T, 0, 0, 0> &gf_gammatzz_,
           //
           const GF3D<const T, 0, 0, 0> &gf_Kh_,
           //
           const GF3D<const T, 0, 0, 0> &gf_Atxx_,
           const GF3D<const T, 0, 0, 0> &gf_Atxy_,
           const GF3D<const T, 0, 0, 0> &gf_Atxz_,
           const GF3D<const T, 0, 0, 0> &gf_Atyy_,
           const GF3D<const T, 0, 0, 0> &gf_Atyz_,
           const GF3D<const T, 0, 0, 0> &gf_Atzz_,
           //
           const GF3D<const T, 0, 0, 0> &gf_Gamtx_,
           const GF3D<const T, 0, 0, 0> &gf_Gamty_,
           const GF3D<const T, 0, 0, 0> &gf_Gamtz_,
           //
           const GF3D<const T, 0, 0, 0> &gf_Theta_,
           //
           const GF3D<const T, 0, 0, 0> &gf_alphaG_,
           //
           const GF3D<const T, 0, 0, 0> &gf_betaGx_,
           const GF3D<const T, 0, 0, 0> &gf_betaGy_,
           const GF3D<const T, 0, 0, 0> &gf_betaGz_,
           //
           const vect<int, 3> &I, const vec3<T> &dx)
      : z4c_vars(kappa1, kappa2, f_mu_L, f_mu_S, eta,
                 //
                 gf_chi_(I), deriv(gf_chi_, I, dx), deriv2(gf_chi_, I, dx),
                 //
                 mat3<T>(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,
                         gf_gammatyz_, gf_gammatzz_, I),
                 mat3<vec3<T> >{
                     deriv(gf_gammatxx_, I, dx),
                     deriv(gf_gammatxy_, I, dx),
                     deriv(gf_gammatxz_, I, dx),
                     deriv(gf_gammatyy_, I, dx),
                     deriv(gf_gammatyz_, I, dx),
                     deriv(gf_gammatzz_, I, dx),
                 },
                 mat3<mat3<T> >{
                     deriv2(gf_gammatxx_, I, dx),
                     deriv2(gf_gammatxy_, I, dx),
                     deriv2(gf_gammatxz_, I, dx),
                     deriv2(gf_gammatyy_, I, dx),
                     deriv2(gf_gammatyz_, I, dx),
                     deriv2(gf_gammatzz_, I, dx),
                 },
                 //
                 gf_Kh_(I), deriv(gf_Kh_, I, dx),
                 //
                 mat3<T>(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_,
                         gf_Atzz_, I),
                 mat3<vec3<T> >{
                     deriv(gf_Atxx_, I, dx),
                     deriv(gf_Atxy_, I, dx),
                     deriv(gf_Atxz_, I, dx),
                     deriv(gf_Atyy_, I, dx),
                     deriv(gf_Atyz_, I, dx),
                     deriv(gf_Atzz_, I, dx),
                 },
                 //
                 vec3<T>(gf_Gamtx_, gf_Gamty_, gf_Gamtz_, I),
                 vec3<vec3<T> >{
                     deriv(gf_Gamtx_, I, dx),
                     deriv(gf_Gamty_, I, dx),
                     deriv(gf_Gamtz_, I, dx),
                 },
                 //
                 gf_Theta_(I), deriv(gf_Theta_, I, dx),
                 //
                 gf_alphaG_(I), deriv(gf_alphaG_, I, dx),
                 deriv2(gf_alphaG_, I, dx),
                 //
                 vec3<T>(gf_betaGx_, gf_betaGy_, gf_betaGz_, I),
                 vec3<vec3<T> >{
                     deriv(gf_betaGx_, I, dx),
                     deriv(gf_betaGy_, I, dx),
                     deriv(gf_betaGz_, I, dx),
                 },
                 vec3<mat3<T> >{
                     deriv2(gf_betaGx_, I, dx),
                     deriv2(gf_betaGy_, I, dx),
                     deriv2(gf_betaGz_, I, dx),
                 },
                 //
                 0,
                 //
                 vec3<T>{0, 0, 0},
                 //
                 mat3<T>{0, 0, 0, 0, 0, 0})
  //
  {}
}; // namespace Z4c
} // namespace Z4c
#endif // #ifndef Z4C_HXX