# Configuration definitions for thorn Weyl
REQUIRES Arith

USES INCLUDE HEADER: vect.hxx

ActiveThorns = "
    ADMBase
    BrillLindquist
    CarpetX
    Formaline
    IOUtil
    ODESolvers
    Weyl
    Z4c
"
CarpetX::verbose = yes
Cactus::cctk_show_schedule = yes
Cactus::terminate = "time"
Cactus::cctk_final_time = 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::ncells_x = 64
CarpetX::ncells_y = 64
CarpetX::ncells_z = 64
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 = 3
CarpetX::max_num_levels = 1
CarpetX::regrid_every = 0
CarpetX::prolongation_type = "ddf"
CarpetX::prolongation_order = 5
CarpetX::interpolation_order = 3
ADMBase::initial_data = "Brill-Lindquist"
ADMBase::initial_lapse = "Brill-Lindquist"
BrillLindquist::x0 = 0.0
BrillLindquist::mass = 1.0
IO::out_dir = $parfile
IO::out_every = 0
IO::out_mode = "np"
IO::out_proc_every = 1
CarpetX::out_plotfile_groups = ""
CarpetX::out_silo_vars = ""
CarpetX::out_tsv = no

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

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

SRCS = weyl.cxx

SRCS = test.cxx weyl.cxx

#if 0

 constexpr mat4<vec4<T, DN>, UP, UP>
calc_dAu(const mat4<T, UP, UP> &gu, const mat4<vec4<T, DN>, UP, UP> &dgu,
         const mat4<T, DN, DN> &A, const mat4<vec4<T, DN>, DN, DN> &dA) {
  // A^ab,c = (g^ax g^by A_xy),c
  //        = g^ax,c g^by A_xy + g^ax g^by,c A_xy + g^ax g^by A_xy,c
  return mat4<vec4<T, DN>, UP, UP>(
      [&](int a, int b)  {
        return vec4<T, DN>([&](int c)  {
          return sum42([&](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);
          });
        });
      });
}
#endif
template <typename T>

calc_riemann(const vec4<mat4<T, DN, DN>, DN> &Gammal,
             const vec4<mat4<T, DN, DN>, UP> &Gamma,
             const vec4<mat4<vec4<T, DN>, DN, DN>, DN> &dGammal) {

calc_riemann(const mat4<T, DN, DN> &g, const vec4<mat4<T, DN, DN>, UP> &Gamma,
             const vec4<mat4<vec4<T, DN>, DN, DN>, UP> &dGamma) {

      return dGammal(a)(d, b)(c)                                              //
             - dGammal(a)(c, b)(d)                                            //
             + sum41([&](int x) { return Gammal(a)(c, x) * Gamma(x)(d, b); }) //
             - sum41([&](int x) { return Gammal(a)(d, x) * Gamma(x)(c, b); });

      return sum41([&](int x) {
        const T rmuxbcd =
            dGamma(x)(d, b)(c)                                              //
            - dGamma(x)(c, b)(d)                                            //
            + sum41([&](int y) { return Gamma(x)(c, y) * Gamma(y)(d, b); }) //
            - sum41([&](int y) { return Gamma(x)(d, y) * Gamma(y)(c, b); });
        return g(a, x) * rmuxbcd;
      });

      return Rm(a, b)(c, d) +
             1 / T{4 - 2} *
                 (R(a, d) * g(b, c) - R(a, c) * g(b, d) + R(b, c) * g(a, d) -
                  R(b, d) * g(a, c)) +
             1 / T{(4 - 1) * (4 - 2)} * Rsc *
                 (g(a, c) * g(b, d) - g(a, d) * g(b, d));

      return Rm(a, b)(c, d) //
             + 1 / T{4 - 2} *
                   (R(a, d) * g(b, c) - R(a, c) * g(b, d)    //
                    + R(b, c) * g(a, d) - R(b, d) * g(a, c)) //
             + 1 / T{(4 - 1) * (4 - 2)} * Rsc *
                   (g(a, c) * g(b, d) - g(a, d) * g(b, c));

using namespace Arith;

    typename norm1<T>::result_type r{0};

    typedef typename norm1<T>::result_type R;
    R r = zero<R>()();

};
} // namespace Weyl
namespace Arith {
template <typename T, Weyl::dnup_t dnup> struct zero<Weyl::vec3<T, dnup> > {
  constexpr Weyl::vec3<T, dnup> operator()() const {
    return Weyl::vec3<T, dnup>(zero<vect<T, 3> >()());
  }

} // namespace Arith
namespace Weyl {

    if (!(is_approx(f(0, 1), f10) && is_approx(f(0, 2), f20) &&
          is_approx(f(1, 2), f21))) {

    if (!(is_approx((*this)(0, 1), f10) && is_approx((*this)(0, 2), f20) &&
          is_approx((*this)(1, 2), f21))) {

      buf << "f(0,1)=" << f(0, 1) << "\n"

      buf << "f(0,1)=" << (*this)(0, 1) << "\n"

          << "f(0,2)=" << f(0, 2) << "\n"

          << "f(0,2)=" << (*this)(0, 2) << "\n"

          << "f(1,2)=" << f(1, 2) << "\n"

          << "f(1,2)=" << (*this)(1, 2) << "\n"

    assert(is_approx(f(0, 1), f10));
    assert(is_approx(f(0, 2), f20));
    assert(is_approx(f(1, 2), f21));

    assert(is_approx(f10, (*this)(0, 1)));
    assert(is_approx(f20, (*this)(0, 2)));
    assert(is_approx(f21, (*this)(1, 2)));

} // namespace Weyl
namespace Arith {
template <typename T, Weyl::dnup_t dnup1, Weyl::dnup_t dnup2>
struct zero<Weyl::mat3<T, dnup1, dnup2> > {
  constexpr Weyl::mat3<T, dnup1, dnup2> operator()() const {
    return Weyl::mat3<T, dnup1, dnup2>(zero<vect<T, 6> >()());
  }
};
} // namespace Arith
namespace Weyl {

}
template <typename F, typename T> constexpr T fold(const F &f, const T &x) {
  return x;
}
template <typename F, typename T, typename... Ts>
constexpr T fold(const F &f, const T &x0, const T &x1, const Ts &... xs) {
  return fold(f, fold(f, x0, x1), xs...);
}
template <typename T> constexpr T add() { return T(0); }
// template <typename T>
//  constexpr  T add(const T &x) {
//   return x;
// }
template <typename T, typename... Ts>
constexpr T add(const T &x, const Ts &... xs) {
  return x + add(xs...);

  R s{0};

  R s = zero<R>()();

  R s{0};

  R s = zero<R>()();

  R s{0};

  R s = zero<R>()();

    assert(n >= 0 && n <= 4);

    assert(n >= 0 && n < 4);

} // namespace Weyl
namespace Arith {
template <typename T, Weyl::dnup_t dnup> struct zero<Weyl::vec4<T, dnup> > {
  constexpr Weyl::vec4<T, dnup> operator()() const {
    return Weyl::vec4<T, dnup>(zero<vect<T, 4> >()());
  }
};
} // namespace Arith
namespace Weyl {

    assert(i >= 0 && i <= j && j <= 4);

    assert(i >= 0 && i <= j && j < 4);

    if (!(is_approx(f(0, 1), f10) && is_approx(f(0, 2), f20) &&
          is_approx(f(0, 3), f30) && is_approx(f(1, 2), f21) &&
          is_approx(f(1, 3), f31) && is_approx(f(2, 3), f32))) {

    if (!(is_approx(f10, (*this)(0, 1)) && is_approx(f20, (*this)(0, 2)) &&
          is_approx(f30, (*this)(0, 3)) && is_approx(f21, (*this)(1, 2)) &&
          is_approx(f31, (*this)(1, 3)) && is_approx(f32, (*this)(2, 3)))) {

      buf << "f(0,1)=" << f(0, 1) << "\n"

      buf << "f(0,1)=" << (*this)(0, 1) << "\n"

          << "f(0,2)=" << f(0, 2) << "\n"

          << "f(0,2)=" << (*this)(0, 2) << "\n"

          << "f(0,3)=" << f(0, 3) << "\n"

          << "f(0,3)=" << (*this)(0, 3) << "\n"

          << "f(1,2)=" << f(1, 2) << "\n"

          << "f(1,2)=" << (*this)(1, 2) << "\n"

          << "f(1,3)=" << f(1, 3) << "\n"

          << "f(1,3)=" << (*this)(1, 3) << "\n"

          << "f(2,3)=" << f(2, 3) << "\n"

          << "f(2,3)=" << (*this)(2, 3) << "\n"

    assert(is_approx(f(0, 1), f10));
    assert(is_approx(f(0, 2), f20));
    assert(is_approx(f(0, 3), f30));
    assert(is_approx(f(1, 2), f21));
    assert(is_approx(f(1, 3), f31));
    assert(is_approx(f(2, 3), f32));

    assert(is_approx((*this)(0, 1), f10));
    assert(is_approx((*this)(0, 2), f20));
    assert(is_approx((*this)(0, 3), f30));
    assert(is_approx((*this)(1, 2), f21));
    assert(is_approx((*this)(1, 3), f31));
    assert(is_approx((*this)(2, 3), f32));

    return pow(A(0, 3), 2) * pow(A(1, 2), 2) -
           2 * A(0, 2) * A(0, 3) * A(1, 2) * A(1, 3) +
           pow(A(0, 2), 2) * pow(A(1, 3), 2) -
           pow(A(0, 3), 2) * A(1, 1) * A(2, 2) +
           2 * A(0, 1) * A(0, 3) * A(1, 3) * A(2, 2) -
           A(0, 0) * pow(A(1, 3), 2) * A(2, 2) +
           2 * A(0, 2) * A(0, 3) * A(1, 1) * A(2, 3) -
           2 * A(0, 1) * A(0, 3) * A(1, 2) * A(2, 3) -
           2 * A(0, 1) * A(0, 2) * A(1, 3) * A(2, 3) +

    return -(A(0, 0) * pow(A(1, 3), 2) * A(2, 2)) +
           pow(A(0, 3), 2) * (pow(A(1, 2), 2) - A(1, 1) * A(2, 2)) +

           pow(A(0, 2), 2) * A(1, 1) * A(3, 3) +
           2 * A(0, 1) * A(0, 2) * A(1, 2) * A(3, 3) -

           2 * A(0, 3) *
               (A(0, 2) * (A(1, 2) * A(1, 3) - A(1, 1) * A(2, 3)) +
                A(0, 1) * (-(A(1, 3) * A(2, 2)) + A(1, 2) * A(2, 3))) -

           A(0, 0) * A(1, 1) * A(2, 2) * A(3, 3);

           A(0, 0) * A(1, 1) * A(2, 2) * A(3, 3) +
           pow(A(0, 2), 2) * (pow(A(1, 3), 2) - A(1, 1) * A(3, 3)) +
           2 * A(0, 1) * A(0, 2) * (-(A(1, 3) * A(2, 3)) + A(1, 2) * A(3, 3));

    return mat4<T, !dnup1, !dnup2>{
        detA1 * (-(pow(A(1, 3), 2) * A(2, 2)) +
                 2 * A(1, 2) * A(1, 3) * A(2, 3) - pow(A(1, 2), 2) * A(3, 3) +
                 A(1, 1) * (-pow(A(2, 3), 2) + A(2, 2) * A(3, 3))),
        detA1 * (A(0, 3) * (A(1, 3) * A(2, 2) - A(1, 2) * A(2, 3)) +
                 A(0, 2) * (-(A(1, 3) * A(2, 3)) + A(1, 2) * A(3, 3)) +
                 A(0, 1) * (pow(A(2, 3), 2) - A(2, 2) * A(3, 3))),
        detA1 * (A(0, 3) * (-(A(1, 2) * A(1, 3)) + A(1, 1) * A(2, 3)) +
                 A(0, 2) * (pow(A(1, 3), 2) - A(1, 1) * A(3, 3)) +
                 A(0, 1) * (-(A(1, 3) * A(2, 3)) + A(1, 2) * A(3, 3))),
        detA1 * (A(0, 3) * (pow(A(1, 2), 2) - A(1, 1) * A(2, 2)) +
                 A(0, 2) * (-(A(1, 2) * A(1, 3)) + A(1, 1) * A(2, 3)) +
                 A(0, 1) * (A(1, 3) * A(2, 2) - A(1, 2) * A(2, 3))),
        detA1 * (-(pow(A(0, 3), 2) * A(2, 2)) +
                 2 * A(0, 2) * A(0, 3) * A(2, 3) - pow(A(0, 2), 2) * A(3, 3) +
                 A(0, 0) * (-pow(A(2, 3), 2) + A(2, 2) * A(3, 3))),
        detA1 * (pow(A(0, 3), 2) * A(1, 2) -
                 A(0, 3) * (A(0, 2) * A(1, 3) + A(0, 1) * A(2, 3)) +
                 A(0, 1) * A(0, 2) * A(3, 3) +
                 A(0, 0) * (A(1, 3) * A(2, 3) - A(1, 2) * A(3, 3))),
        detA1 * (pow(A(0, 2), 2) * A(1, 3) + A(0, 1) * A(0, 3) * A(2, 2) -
                 A(0, 2) * (A(0, 3) * A(1, 2) + A(0, 1) * A(2, 3)) +
                 A(0, 0) * (-(A(1, 3) * A(2, 2)) + A(1, 2) * A(2, 3))),
        detA1 * (-(pow(A(0, 3), 2) * A(1, 1)) +
                 2 * A(0, 1) * A(0, 3) * A(1, 3) - pow(A(0, 1), 2) * A(3, 3) +
                 A(0, 0) * (-pow(A(1, 3), 2) + A(1, 1) * A(3, 3))),
        detA1 * (-(A(0, 1) * A(0, 3) * A(1, 2)) +
                 A(0, 2) * (A(0, 3) * A(1, 1) - A(0, 1) * A(1, 3)) +
                 pow(A(0, 1), 2) * A(2, 3) +
                 A(0, 0) * (A(1, 2) * A(1, 3) - A(1, 1) * A(2, 3))),
        detA1 * (-(pow(A(0, 2), 2) * A(1, 1)) +
                 2 * A(0, 1) * A(0, 2) * A(1, 2) - pow(A(0, 1), 2) * A(2, 2) +
                 A(0, 0) * (-pow(A(1, 2), 2) + A(1, 1) * A(2, 2)))};

  }
};
} // namespace Weyl
namespace Arith {
template <typename T, Weyl::dnup_t dnup1, Weyl::dnup_t dnup2>
struct zero<Weyl::mat4<T, dnup1, dnup2> > {
  constexpr Weyl::mat4<T, dnup1, dnup2> operator()() const {
    return Weyl::mat4<T, dnup1, dnup2>(zero<vect<T, 10> >()());

} // namespace Arith
namespace Weyl {

    assert(i >= 0 && i < j && j <= 4);

    assert(i >= 0 && i < j && j < 4);

  static_assert(sign(1, 0) == -sign(0, 1), "");
  static_assert(sign(2, 0) == -sign(0, 2), "");
  static_assert(sign(3, 0) == -sign(0, 3), "");
  static_assert(sign(2, 1) == -sign(1, 2), "");
  static_assert(sign(3, 1) == -sign(1, 3), "");
  static_assert(sign(3, 2) == -sign(2, 3), "");
  static_assert(sign(0, 0) == 0, "");
  static_assert(sign(1, 1) == 0, "");
  static_assert(sign(2, 2) == 0, "");
  static_assert(sign(3, 3) == 0, "");

    const T f00 = f(0, 0);

    const T f11 = f(1, 1);

    const T f22 = f(2, 2);

    const T f33 = f(3, 3);

    if (!(is_approx(f(0, 0), T{0}) && is_approx(f(0, 1), f10) &&
          is_approx(f(0, 2), f20) && is_approx(f(0, 3), f30) &&
          is_approx(f(1, 1), T{0}) && is_approx(f(1, 2), f21) &&
          is_approx(f(1, 3), f31) && is_approx(f(2, 2), T{0}) &&
          is_approx(f(2, 3), f32) && is_approx(f(3, 3), T{0}))) {

    if (!(is_approx(f00, zero<T>()()) && is_approx(f10, -(*this)(0, 1)) &&
          is_approx(f20, -(*this)(0, 2)) && is_approx(f30, -(*this)(0, 3)) &&
          is_approx(f11, zero<T>()()) && is_approx(f21, -(*this)(1, 2)) &&
          is_approx(f31, -(*this)(1, 3)) && is_approx(f22, zero<T>()()) &&
          is_approx(f32, -(*this)(2, 3)) && is_approx(f33, zero<T>()()))) {

      buf << "f(0,0)=" << f(0, 0) << "\n"
          << "f(0,1)=" << f(0, 1) << "\n"

      buf << "f(0,0)=" << f00 << "\n"
          << "f(0,1)=" << (*this)(0, 1) << "\n"

          << "f(0,2)=" << f(0, 2) << "\n"

          << "f(0,2)=" << (*this)(0, 2) << "\n"

          << "f(0,3)=" << f(0, 3) << "\n"

          << "f(0,3)=" << (*this)(0, 3) << "\n"

          << "f(1,1)=" << f(1, 1) << "\n"
          << "f(1,2)=" << f(1, 2) << "\n"

          << "f(1,1)=" << f11 << "\n"
          << "f(1,2)=" << (*this)(1, 2) << "\n"

          << "f(1,3)=" << f(1, 3) << "\n"

          << "f(1,3)=" << (*this)(1, 3) << "\n"

          << "f(2,2)=" << f(2, 2) << "\n"
          << "f(2,3)=" << f(2, 3) << "\n"

          << "f(2,2)=" << f22 << "\n"
          << "f(2,3)=" << (*this)(2, 3) << "\n"

          << "f(3,3)=" << f(3, 3) << "\n";
      CCTK_VERROR("symmetric matrix is not symmetric:\n%s", buf.str().c_str());

          << "f(3,3)=" << f33 << "\n";
      CCTK_VERROR("antisymmetric matrix is not antisymmetric:\n%s",
                  buf.str().c_str());

    assert(is_approx(f(0, 0), T{0}));
    assert(is_approx(f(0, 1), f10));
    assert(is_approx(f(0, 2), f20));
    assert(is_approx(f(0, 3), f30));
    assert(is_approx(f(1, 1), T{0}));
    assert(is_approx(f(1, 2), f21));
    assert(is_approx(f(1, 3), f31));
    assert(is_approx(f(2, 2), T{0}));
    assert(is_approx(f(2, 3), f32));
    assert(is_approx(f(3, 3), T{0}));

    assert(is_approx(f00, zero<T>()()));
    assert(is_approx(f10, -(*this)(0, 1)));
    assert(is_approx(f20, -(*this)(0, 2)));
    assert(is_approx(f30, -(*this)(0, 3)));
    assert(is_approx(f11, zero<T>()()));
    assert(is_approx(f21, -(*this)(1, 2)));
    assert(is_approx(f31, -(*this)(1, 3)));
    assert(is_approx(f22, zero<T>()()));
    assert(is_approx(f32, -(*this)(2, 3)));
    assert(is_approx(f33, zero<T>()()));

      return T{0};

      return zero<T>()();

            elts[4](I), elts[4](I), elts[5](I)};

            elts[3](I), elts[4](I), elts[5](I)};

    return Power(A(0, 3), 2) * Power(A(1, 2), 2) -

    return pow(A(0, 3), 2) * pow(A(1, 2), 2) -

           Power(A(0, 2), 2) * Power(A(1, 3), 2) +

           pow(A(0, 2), 2) * pow(A(1, 3), 2) +

           Power(A(0, 1), 2) * Power(A(2, 3), 2);

           pow(A(0, 1), 2) * pow(A(2, 3), 2);

  }
};
} // namespace Weyl
namespace Arith {
template <typename T, Weyl::dnup_t dnup1, Weyl::dnup_t dnup2>
struct zero<Weyl::amat4<T, dnup1, dnup2> > {
  constexpr Weyl::amat4<T, dnup1, dnup2> operator()() const {
    return Weyl::amat4<T, dnup1, dnup2>(zero<vect<T, 6> >()());

} // namespace Arith
namespace Weyl {

#if 0
template <typename F, typename T>
constexpr  T fold(const F &f, const T &x) {
  return x;
}
template <typename F, typename T, typename... Ts>
constexpr  T fold(const F &f, const T &x0,
                                              const T &x1, const Ts &... xs) {
  return fold(f, fold(f, x0, x1), xs...);
}

template <typename T> constexpr  T add() {
  return T(0);
}
// template <typename T>
//  constexpr  T add(const T &x) {
//   return x;
// }
template <typename T, typename... Ts>
constexpr  T add(const T &x, const Ts &... xs) {
  return x + add(xs...);
}
#endif

  R s{0};

  R s = zero<R>()();

  R s{0};

  R s = zero<R>()();

  R s{0};

  R s = zero<R>()();

////////////////////////////////////////////////////////////////////////////////
template <typename F, typename T> constexpr T fold(const F &f, const T &x) {
  return x;
}
template <typename F, typename T, typename... Ts>
constexpr T fold(const F &f, const T &x0, const T &x1, const Ts &... xs) {
  return fold(f, fold(f, x0, x1), xs...);
}
template <typename T> constexpr T add() { return T(0); }
// template <typename T>
//  constexpr  T add(const T &x) {
//   return x;
// }
template <typename T, typename... Ts>
constexpr T add(const T &x, const Ts &... xs) {
  return x + add(xs...);
}

#include "derivs.hxx"
#include "physics.hxx"
#include "tensor.hxx"
#include <cctk.h>
#include <algorithm>
#include <array>
#include <cassert>
#include <cmath>
#include <random>
namespace Weyl {
using namespace std;
// TODO: Use GoogleTest instead of assert
extern "C" void Weyl_Test(CCTK_ARGUMENTS) {
  DECLARE_CCTK_ARGUMENTS;
  // Test tensors
  {
    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, DN, DN>(
          [&](int a, int b) { return arr[min(a, b)][max(a, b)]; });
    }};
    const mat3<double, DN, DN> Z([&](int a, int b) { return double(0); });
    const mat3<double, DN, DN> I([&](int a, int b) { return double(a == b); });
    assert(I != Z);
    const mat3<double, UP, UP> Zup([&](int a, int b) { return double(0); });
    const mat3<double, UP, UP> Iup(
        [&](int a, int b) { return double(a == b); });
    assert(Iup != Zup);
    for (int n = 0; n < 100; ++n) {
      const mat3<double, DN, DN> A = randmat10();
      const mat3<double, DN, DN> B = randmat10();
      const mat3<double, DN, DN> 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)));
      // DNUP  assert(mul(I, A) == A);
      // DNUP  assert(mul(A, I) == A);
      // DNUP  assert(mul(Z, A) == Z);
      // DNUP  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) == Zup);
      assert(I.inv(1) == Iup);
      // DNUP assert(mul(A.inv(1), A) == A.det() * Iup);
      // DNUP assert(mul(A, A.inv(1)) == A.det() * Iup);
      assert((a * A).inv(1) == pow2(a) * A.inv(1));
    }
  }
  // Test 4-tensors
  {
    mt19937 engine(42);
    uniform_int_distribution<int> dist(-10, 10);
    const auto rand10{[&]() { return double(dist(engine)); }};
    const auto randmat10{[&]() {
      array<array<double, 4>, 4> arr;
      for (int a = 0; a < 4; ++a)
        for (int b = 0; b < 4; ++b)
          arr[a][b] = rand10();
      return mat4<double, DN, DN>(
          [&](int a, int b) { return arr[min(a, b)][max(a, b)]; });
    }};
    const mat4<double, DN, DN> Z([&](int a, int b) { return double(0); });
    const mat4<double, DN, DN> I([&](int a, int b) { return double(a == b); });
    assert(I != Z);
    const mat4<double, UP, UP> Zup([&](int a, int b) { return double(0); });
    const mat4<double, UP, UP> Iup(
        [&](int a, int b) { return double(a == b); });
    assert(Iup != Zup);
    for (int n = 0; n < 100; ++n) {
      const mat4<double, DN, DN> A = randmat10();
      const mat4<double, DN, DN> B = randmat10();
      const mat4<double, DN, DN> 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)));
      // DNUP  assert(mul(I, A) == A);
      // DNUP  assert(mul(A, I) == A);
      // DNUP  assert(mul(Z, A) == Z);
      // DNUP  assert(mul(A, Z) == Z);
      assert(Z.det() == 0);
      assert(I.det() == 1);
      assert((a * A).det() == pow4(a) * A.det());
      assert(Z.inv(1) == Zup);
      assert(I.inv(1) == Iup);
      // DNUP assert(mul(A.inv(1), A) == A.det() * Iup);
      // DNUP assert(mul(A, A.inv(1)) == A.det() * Iup);
      assert((a * A).inv(1) == pow3(a) * A.inv(1));
    }
  }
  // Test derivatives
  static_assert(deriv_order % 2 == 0, "");
  constexpr int required_ghosts = deriv_order / 2 + 1;
  const double eps = 1.0e-12;
  // deriv
  for (int order = 0; order <= deriv_order; ++order) {
    // CCTK_VINFO("Testing deriv order=%d", order);
    array<double, 2 * required_ghosts + 7> arr;
    for (size_t i = 0; i < arr.size(); ++i)
      arr[i] = NAN;
    double *const var = &arr[arr.size() / 2];
    for (int i = -deriv_order / 2; i <= deriv_order / 2; ++i)
      var[i] = pown(i, order);
    const double expected = order == 1 ? 1 : 0;
    const double found = deriv1d(var, 1, 1.0);
    assert(fabs(found - expected) <= eps);
  }
  // deriv2
  for (int order = 0; order <= deriv_order; ++order) {
    // CCTK_VINFO("Testing deriv2 order=%d", order);
    array<double, 2 * required_ghosts + 7> arr;
    for (size_t i = 0; i < arr.size(); ++i)
      arr[i] = NAN;
    double *const var = &arr[arr.size() / 2];
    for (int i = -deriv_order / 2; i <= deriv_order / 2; ++i)
      var[i] = pown(i, order);
    const double expected = order == 2 ? 2 : 0;
    const double found = deriv2_1d(var, 1, 1.0);
    assert(fabs(found - expected) <= eps);
  }
  // deriv2 (mixed)
  for (int orderj = 0; orderj <= deriv_order; ++orderj) {
    for (int orderi = 0; orderi <= deriv_order; ++orderi) {
      // CCTK_VINFO("Testing deriv2 (mixed) order=%d,%d", orderi, orderj);
      array<array<double, 2 * required_ghosts + 7>, 2 * required_ghosts + 7>
          arr;
      for (size_t j = 0; j < arr.size(); ++j)
        for (size_t i = 0; i < arr[j].size(); ++i)
          arr[j][i] = NAN;
      const int di = 1;
      const int dj = arr[0].size();
      double *const var = &arr[arr.size() / 2][arr[0].size() / 2];
      for (int j = -deriv_order / 2; j <= deriv_order / 2; ++j)
        for (int i = -deriv_order / 2; i <= deriv_order / 2; ++i)
          var[j * dj + i * di] = pown(i, orderi) * pown(j, orderj);
      const double expected = orderi == 1 && orderj == 1 ? 1 : 0;
      const double found = deriv2_2d(var, di, dj, 1.0, 1.0);
      assert(fabs(found - expected) <= eps);
    }
  }
}
} // namespace Weyl

            return betal(b);

            return betal(b - 1);

            return betal(a);
          return gamma(a, b);

            return betal(a - 1);
          return gamma(a - 1, b - 1);

#if 0
  
  weyl_vars_noderivs(const GF3D<const T, 0, 0, 0> &gf_gammaxx_,
                     const GF3D<const T, 0, 0, 0> &gf_gammaxy_,
                     const GF3D<const T, 0, 0, 0> &gf_gammaxz_,
                     const GF3D<const T, 0, 0, 0> &gf_gammayy_,
                     const GF3D<const T, 0, 0, 0> &gf_gammayz_,
                     const GF3D<const T, 0, 0, 0> &gf_gammazz_,
                     //
                     const GF3D<const T, 0, 0, 0> &gf_alpha_,
                     //
                     const GF3D<const T, 0, 0, 0> &gf_betax_,
                     const GF3D<const T, 0, 0, 0> &gf_betay_,
                     const GF3D<const T, 0, 0, 0> &gf_betaz_,
                     //
                     const vect<int, 3> &I)
      : weyl_vars_noderivs<T>(mat3<T, DN, DN>(gf_gammaxx_, gf_gammaxy_,
                                              gf_gammaxz_, gf_gammayy_,
                                              gf_gammayz_, gf_gammazz_, I),
                              //
                              gf_alpha_(I),
                              //
                              vec3<T, UP>(gf_betax_, gf_betay_, gf_betaz_, I))
  // {}
#endif

  // const vec4<mat4<vec4<T, DN>, DN, DN>, UP> dGamma;

  const vec4<mat4<vec4<T, DN>, DN, DN>, UP> dGamma;

        dGammal(calc_dgammal(ddg)), //
        // dGamma(calc_dgamma(gu,dgu,Gammal,dGammal)),   //

        dGammal(calc_dgammal(ddg)),                    //
        dGamma(calc_dgamma(gu, dgu, Gammal, dGammal)), //

        Rm(calc_riemann(Gammal, Gamma, dGammal)), //
        R(calc_ricci(gu, Rm)),                    //
        Rsc(R.trace(gu)),                         //

        Rm(calc_riemann(g, Gamma, dGamma)), //
        R(calc_ricci(gu, Rm)),              //
        Rsc(R.trace(gu)),                   //

  //
  {}
#if 0
  weyl_vars(const GF3D<const T, 0, 0, 0> &gf_gammaxx_,
            const GF3D<const T, 0, 0, 0> &gf_gammaxy_,
            const GF3D<const T, 0, 0, 0> &gf_gammaxz_,
            const GF3D<const T, 0, 0, 0> &gf_gammayy_,
            const GF3D<const T, 0, 0, 0> &gf_gammayz_,
            const GF3D<const T, 0, 0, 0> &gf_gammazz_,
            //
            const GF3D<const T, 0, 0, 0> &gf_alpha_,
            //
            const GF3D<const T, 0, 0, 0> &gf_betax_,
            const GF3D<const T, 0, 0, 0> &gf_betay_,
            const GF3D<const T, 0, 0, 0> &gf_betaz_,
            //
            const GF3D<const T, 0, 0, 0> &gf_kxx_,
            const GF3D<const T, 0, 0, 0> &gf_kxy_,
            const GF3D<const T, 0, 0, 0> &gf_kxz_,
            const GF3D<const T, 0, 0, 0> &gf_kyy_,
            const GF3D<const T, 0, 0, 0> &gf_kyz_,
            const GF3D<const T, 0, 0, 0> &gf_kzz_,
            //
            const GF3D<const T, 0, 0, 0> &gf_dtalpha_,
            //
            const GF3D<const T, 0, 0, 0> &gf_dtbetax_,
            const GF3D<const T, 0, 0, 0> &gf_dtbetay_,
            const GF3D<const T, 0, 0, 0> &gf_dtbetaz_,
            //
            const GF3D<const T, 0, 0, 0> &gf_dtkxx_,
            const GF3D<const T, 0, 0, 0> &gf_dtkxy_,
            const GF3D<const T, 0, 0, 0> &gf_dtkxz_,
            const GF3D<const T, 0, 0, 0> &gf_dtkyy_,
            const GF3D<const T, 0, 0, 0> &gf_dtkyz_,
            const GF3D<const T, 0, 0, 0> &gf_dtkzz_,
            //
            const GF3D<const T, 0, 0, 0> &gf_dt2alpha_,
            //
            const GF3D<const T, 0, 0, 0> &gf_dt2betax_,
            const GF3D<const T, 0, 0, 0> &gf_dt2betay_,
            const GF3D<const T, 0, 0, 0> &gf_dt2betaz_,
            //
            const vect<int, 3> &I, const vec3<T, UP> &dx)
      : weyl_vars(mat3<T, DN, DN>(gf_gammaxx_, gf_gammaxy_, gf_gammaxz_,
                                  gf_gammayy_, gf_gammayz_, gf_gammazz_, I),
                  //
                  gf_alpha_(I),
                  //
                  vec3<T, UP>(gf_betax_, gf_betay_, gf_betaz_, I),
                  //
                  mat3<T, DN, DN>(gf_kxx_, gf_kxy_, gf_kxz_, gf_kyy_, gf_kyz_,
                                  gf_kzz_, I),
                  //
                  gf_dtalpha_(I),
                  //
                  vec3<T, UP>(gf_dtbetax_, gf_dtbetay_, gf_dtbetaz_, I),
                  //
                  mat3<vec3<T, DN>, DN, DN>{
                      deriv(gf_gammaxx_, I, dx),
                      deriv(gf_gammaxy_, I, dx),
                      deriv(gf_gammaxz_, I, dx),
                      deriv(gf_gammayy_, I, dx),
                      deriv(gf_gammayz_, I, dx),
                      deriv(gf_gammazz_, I, dx),
                  },
                  //
                  deriv(gf_alpha_, I, dx),
                  //
                  vec3<vec3<T, DN>, UP>{
                      deriv(gf_betax_, I, dx),
                      deriv(gf_betay_, I, dx),
                      deriv(gf_betaz_, I, dx),
                  },
                  //
                  mat3<T, DN, DN>(gf_dtkxx_, gf_dtkxy_, gf_dtkxz_, gf_dtkyy_,
                                  gf_dtkyz_, gf_dtkzz_, I),
                  //
                  gf_dt2alpha_(I),
                  //
                  vec3<T, UP>(gf_dt2betax_, gf_dt2betay_, gf_dt2betaz_, I),
                  //
                  mat3<vec3<T, DN>, DN, DN>{
                      deriv(gf_kxx_, I, dx),
                      deriv(gf_kxy_, I, dx),
                      deriv(gf_kxz_, I, dx),
                      deriv(gf_kyy_, I, dx),
                      deriv(gf_kyz_, I, dx),
                      deriv(gf_kzz_, I, dx),
                  },
                  //
                  deriv(gf_dtalpha_, I, dx),
                  //
                  vec3<vec3<T, DN>, UP>{
                      deriv(gf_dtbetax_, I, dx),
                      deriv(gf_dtbetay_, I, dx),
                      deriv(gf_dtbetaz_, I, dx),
                  },
                  //
                  mat3<mat3<T, DN, DN>, DN, DN>{
                      deriv2(gf_gammaxx_, I, dx),
                      deriv2(gf_gammaxy_, I, dx),
                      deriv2(gf_gammaxz_, I, dx),
                      deriv2(gf_gammayy_, I, dx),
                      deriv2(gf_gammayz_, I, dx),
                      deriv2(gf_gammazz_, I, dx),
                  },
                  //
                  deriv2(gf_alpha_, I, dx),
                  //
                  vec3<mat3<T, DN, DN>, UP>{
                      deriv2(gf_betax_, I, dx),
                      deriv2(gf_betay_, I, dx),
                      deriv2(gf_betaz_, I, dx),
                  })

#endif