eschnett/cactusamrex - Change 7SFTOTPQMTRL2WMRCWHJGLWYFM42I2VTMMBOKAYCHIN7NEY2U52AC

Weyl: Implement tetrad, NP spin coefficients, Ricci/Weyl scalars

Created by Erik Schnetter on April 17, 2020

7SFTOTPQMTRL2WMRCWHJGLWYFM42I2VTMMBOKAYCHIN7NEY2U52AC

Dependencies

In channels

main

master

Change contents

Insertion in Weyl/interface.ccl at line 54 [3.308]

[3.1730]


CCTK_REAL tetrad_l TYPE=gf TAGS='index={0 0 0}'
{
  lt lx ly lz
} "Tetrad vector l^a"
CCTK_REAL tetrad_n TYPE=gf TAGS='index={0 0 0}'
{
  nt nx ny nz
} "Tetrad vector n^a"
# TODO: Combine these, represent as complex number
CCTK_REAL tetrad_mre TYPE=gf TAGS='index={0 0 0}'
{
  mret mrex mrey mrez
} "Tetrad vector real(m^a)"
CCTK_REAL tetrad_mim TYPE=gf TAGS='index={0 0 0}'
{
  mimt mimx mimy mimz
} "Tetrad vector imag(m^a)"
CCTK_REAL ricci_scalars TYPE=gf TAGS='index={0 0 0}'
{
  Lambda
  Phi00
  Phi11
  Phi22
  Phi10re Phi10im
  Phi20re Phi20im
  Phi21re Phi21im
} "Ricci scalars"
CCTK_REAL weyl_scalars TYPE=gf TAGS='index={0 0 0}'
{
  Psi0re Psi0im
  Psi1re Psi1im
  Psi2re Psi2im
  Psi3re Psi3im
  Psi4re Psi4im
} "Weyl scalars"
CCTK_REAL spin_coefficients TYPE=gf TAGS='index={0 0 0}'
{
  npkappare   npkappaim
  npsigmare   npsigmaim
  nprhore     nprhoim
  nptaure     nptauim
  npepsilonre npepsilonim
  npbetare    npbetaim
  npalphare   npalphaim
  npgammare   npgammaim
  nppire      nppiim
  npmure      npmuim
  nplambdare  nplambdaim
  npnure      npnuim
} "Newman-Penrose spin coefficients"

File move: brill-lindquist.par → weyl-brill-lindquist.par
BF:BFD[2.118] → [2.119:140]
[2.118]
[2.141]
Replacement in Weyl/par/weyl-brill-lindquist.par at line 58 [2.141]
B:BD[2.1216] → [2.1216:1234]
```
IO::out_every = 0
```
[2.1216]
[2.1234]
```
IO::out_every = 1
```

Replacement in Weyl/par/weyl-brill-lindquist.par at line 63 [2.141]

B:BD[2.1312] → [2.1312:1340]

CarpetX::out_silo_vars = ""

[2.1312]

[2.1340]

CarpetX::out_silo_vars = "
    Weyl::metric4
    Weyl::ricci4
    Weyl::weyl4
    Weyl::tetrad_l
    Weyl::tetrad_n
    Weyl::tetrad_mre
    Weyl::tetrad_mim
    Weyl::ricci_scalars
    Weyl::weyl_scalars
    Weyl::spin_coefficients
"

Insertion in Weyl/schedule.ccl at line 9 [3.1799]
[3.1943]
[3.1943]
```
STORAGE: tetrad_l tetrad_n tetrad_mre tetrad_mim
```

Insertion in Weyl/schedule.ccl at line 39 [3.1799]

[3.2612]

  WRITES: tetrad_l(interior)
  WRITES: tetrad_n(interior)
  WRITES: tetrad_mre(interior)
  WRITES: tetrad_mim(interior)
  WRITES: ricci_scalars(interior)
  WRITES: weyl_scalars(interior)
  WRITES: spin_coefficients(interior)

Insertion in Weyl/src/physics.hxx at line 5 [3.11226]
[3.11290]
[3.11290]
```
#include <dual.hxx>
#include <cmath>
```
Insertion in Weyl/src/physics.hxx at line 11 [3.11226]
[3.11308]
[3.11308]
```
using namespace std;
```

Insertion in Weyl/src/physics.hxx at line 130 [3.11226]

[3.15915]

}
template <typename T>
constexpr vec4<T, UP> raise(const mat4<T, UP, UP> &gu, const vec4<T, DN> &vl) {
  return vec4<T, UP>(
      [&](int a) { return sum41([&](int x) { return gu(a, x) * vl(x); }); });
}
template <typename T>
constexpr vec4<T, DN> lower(const mat4<T, DN, DN> &g, const vec4<T, UP> &v) {
  return vec4<T, DN>(
      [&](int a) { return sum41([&](int x) { return g(a, x) * v(x); }); });
}
template <typename T>
constexpr T dot(const vec4<T, DN> &vl, const vec4<T, UP> &v) {
  return sum41([&](int x) { return vl(x) * v(x); });
}
template <typename T>
constexpr vec4<T, UP> normalized(const mat4<T, DN, DN> &g,
                                 const vec4<T, UP> &v) {
  return v / sqrt(dot(lower(g, v), v));
}
template <typename T>
constexpr vec4<T, UP> projected(const mat4<T, DN, DN> &g, const vec4<T, UP> &v,
                                const vec4<T, UP> &w) {
  const auto z = zero<T>()();
  const auto wl = lower(g, w);
  const auto wlen2 = dot(wl, w);
  if (wlen2 < pow(1.0e-12, 2))
    return vec4<T, UP>([&](int a) { return z; });
  return dot(wl, v) / wlen2 * w;
}
template <typename T>
constexpr vec4<T, UP> projected1(const vec4<T, UP> &v, const vec4<T, DN> &wl,
                                 const vec4<T, UP> &w) {
  // assuming dot(wl, w) == 1
  return dot(wl, v) * w;
}
template <typename T>
constexpr vec4<T, UP> rejected(const mat4<T, DN, DN> &g, const vec4<T, UP> &v,
                               const vec4<T, UP> &w) {
  return v - projected(g, v, w);
}
template <typename T>
constexpr vec4<T, UP> rejected1(const vec4<T, UP> &v, const vec4<T, DN> &wl,
                                const vec4<T, UP> &w) {
  return v - projected1(v, wl, w);
}
template <typename T> constexpr vec4<T, UP> calc_et(const mat4<T, UP, UP> &gu) {
  const auto z = zero<T>()();
  const auto e = one<T>()();
  const vec4<T, DN> etl([&](int a) { return a == 0 ? e : z; });
  const auto et = raise(gu, etl);
  const auto etlen = sqrt(-dot(etl, et));
  assert(!isnan(etlen) && etlen >= 1.0e-12);
  return et / etlen;
}
template <typename T>
constexpr vec4<T, UP> calc_ephi(const vec4<T, UP> &x,
                                const mat4<T, DN, DN> &g) {
  const auto z = zero<T>()();
  vec4<T, UP> ephi(z, -x(2), x(1), z);
  const auto ephil = lower(g, ephi);
  const auto ephi_len2 = dot(ephil, ephi);
  if (ephi_len2 < pow(1.0e-12, 2))
    return vec4<T, UP>([&](int a) { return z; });
  ephi /= sqrt(ephi_len2);
  return ephi;

Insertion in Weyl/src/physics.hxx at line 208 [3.11226]

[3.15918]

template <typename T>
constexpr vec4<T, UP> calc_etheta(const vec4<T, UP> &x,
                                  const mat4<T, DN, DN> &g,
                                  const vec4<T, UP> &ephi) {
  const auto z = zero<T>()();
  const T rho2 = pow(x(1), 2) + pow(x(2), 2);
  vec4<T, UP> etheta(z, x(1) * x(3), x(2) * x(3), -rho2);
  const auto ethetal = lower(g, etheta);
  const auto etheta_len2 = dot(ethetal, etheta);
  if (etheta_len2 < pow(1.0e-12, 2))
    return vec4<T, UP>([&](int a) { return z; });
  etheta /= sqrt(etheta_len2); // to improve accuracy
  etheta = rejected(g, etheta, ephi);
  etheta = normalized(g, etheta);
  return etheta;
}
template <typename T>
constexpr vec4<T, UP> calc_er(const vec4<T, UP> &x, const mat4<T, DN, DN> &g,
                              const vec4<T, UP> &etheta,
                              const vec4<T, UP> &ephi) {
  const auto z = zero<T>()();
  vec4<T, UP> er(z, x(1), x(2), x(3));
  const auto erl = lower(g, er);
  const auto er_len2 = dot(erl, er);
  if (er_len2 < pow(1.0e-12, 2))
    return vec4<T, UP>([&](int a) { return z; });
  er /= sqrt(er_len2); // to improve accuracy
  er = rejected(g, er, etheta);
  er = rejected(g, er, ephi);
  er = normalized(g, er);
  return er;
}
template <typename T>
constexpr vec4<vec4<T, DN>, UP>
calc_det(const mat4<T, UP, UP> &gu, const mat4<vec4<T, DN>, UP, UP> &dgu,
         const vec4<T, UP> &et, const vec4<mat4<T, DN, DN>, UP> &Gamma) {
  typedef vec4<T, DN> DT;
  typedef dual<T, DT> TDT;
  const mat4<TDT, UP, UP> gu1(
      [&](int a, int b) { return TDT(gu(a, b), dgu(a, b)); });
  const auto det1 = calc_et(gu1);
  const vec4<vec4<T, DN>, UP> det([&](int a) {
    return vec4<T, DN>([&](int b) {
      return det1(a).eps(b) +
             sum41([&](int x) { return Gamma(a)(b, x) * et(x); });
    });
  });
  for (int a = 0; a < 4; ++a)
    for (int b = 0; b < 4; ++b)
      assert(!isnan(det(a)(b)));
  return det;
}
template <typename T>
constexpr vec4<vec4<T, DN>, UP>
calc_dephi(const vec4<T, UP> &x, const mat4<T, DN, DN> &g,
           const mat4<vec4<T, DN>, DN, DN> &dg, const vec4<T, UP> &ephi,
           const vec4<mat4<T, DN, DN>, UP> &Gamma) {
  typedef vec4<T, DN> DT;
  typedef dual<T, DT> TDT;
  const vec4<TDT, UP> x1(
      [&](int a) { return TDT(x(a), DT([&](int b) { return T(a == b); })); });
  const mat4<TDT, DN, DN> g1(
      [&](int a, int b) { return TDT(g(a, b), dg(a, b)); });
  const auto dephi1 = calc_ephi(x1, g1);
  const vec4<vec4<T, DN>, UP> dephi([&](int a) {
    return vec4<T, DN>([&](int b) {
      return dephi1(a).eps(b) +
             sum41([&](int x) { return Gamma(a)(b, x) * ephi(x); });
    });
  });
  for (int a = 0; a < 4; ++a)
    for (int b = 0; b < 4; ++b)
      assert(!isnan(dephi(a)(b)));
  return dephi;
}
template <typename T>
constexpr vec4<vec4<T, DN>, UP>
calc_detheta(const vec4<T, UP> &x, const mat4<T, DN, DN> &g,
             const mat4<vec4<T, DN>, DN, DN> &dg, const vec4<T, UP> &ephi,
             const vec4<vec4<T, DN>, UP> &dephi, const vec4<T, UP> &etheta,
             const vec4<mat4<T, DN, DN>, UP> &Gamma) {
  typedef vec4<T, DN> DT;
  typedef dual<T, DT> TDT;
  const vec4<TDT, UP> x1(
      [&](int a) { return TDT(x(a), DT([&](int b) { return T(a == b); })); });
  const mat4<TDT, DN, DN> g1(
      [&](int a, int b) { return TDT(g(a, b), dg(a, b)); });
  const vec4<TDT, UP> ephi1([&](int a) { return TDT(ephi(a), dephi(a)); });
  const auto detheta1 = calc_etheta(x1, g1, ephi1);
  const vec4<vec4<T, DN>, UP> detheta([&](int a) {
    return vec4<T, DN>([&](int b) {
      return detheta1(a).eps(b) +
             sum41([&](int x) { return Gamma(a)(b, x) * etheta(x); });
    });
  });
  for (int a = 0; a < 4; ++a)
    for (int b = 0; b < 4; ++b)
      assert(!isnan(detheta(a)(b)));
  return detheta;
}
template <typename T>
constexpr vec4<vec4<T, DN>, UP>
calc_der(const vec4<T, UP> &x, const mat4<T, DN, DN> &g,
         const mat4<vec4<T, DN>, DN, DN> &dg, const vec4<T, UP> &ephi,
         const vec4<vec4<T, DN>, UP> &dephi, const vec4<T, UP> &etheta,
         const vec4<vec4<T, DN>, UP> &detheta, const vec4<T, UP> &er,
         const vec4<mat4<T, DN, DN>, UP> &Gamma) {
  typedef vec4<T, DN> DT;
  typedef dual<T, DT> TDT;
  const vec4<TDT, UP> x1(
      [&](int a) { return TDT(x(a), DT([&](int b) { return T(a == b); })); });
  const mat4<TDT, DN, DN> g1(
      [&](int a, int b) { return TDT(g(a, b), dg(a, b)); });
  const vec4<TDT, UP> ephi1([&](int a) { return TDT(ephi(a), dephi(a)); });
  const vec4<TDT, UP> etheta1(
      [&](int a) { return TDT(etheta(a), detheta(a)); });
  const auto der1 = calc_er(x1, g1, ephi1, etheta1);
  const vec4<vec4<T, DN>, UP> der([&](int a) {
    return vec4<T, DN>([&](int b) {
      return der1(a).eps(b) +
             sum41([&](int x) { return Gamma(a)(b, x) * er(x); });
    });
  });
  for (int a = 0; a < 4; ++a)
    for (int b = 0; b < 4; ++b)
      assert(!isnan(der(a)(b)));
  return der;
}

Insertion in Weyl/src/tensor.hxx at line 4 [3.15983]
[3.16023]
[3.16023]
```
#include <dual.hxx>
```
Insertion in Weyl/src/tensor.hxx at line 6 [3.15983]
[3.16043]
[3.16043]
```
#include <vect.hxx>
```

Insertion in Weyl/src/tensor.hxx at line 640 [3.15983]

[3.35838]

  }
  friend constexpr vec4<T, dnup> operator*(const vec4<T, dnup> &x, const T &a) {
    return {x.elts * a};
  }
  friend constexpr vec4<T, dnup> operator/(const vec4<T, dnup> &x, const T &a) {
    return {x.elts / a};

Insertion in Weyl/src/tensor.hxx at line 647 [3.15983]

[3.35842]


  vec4<T, dnup> &operator*=(const T &a) { return *this = *this * a; }
  vec4<T, dnup> &operator/=(const T &a) { return *this = *this / a; }

Insertion in Weyl/src/tensor.hxx at line 660 [3.15983]

[3.36166]

  constexpr vec4<bool, dnup> isnan() const {
    return vec4<bool, dnup>(std::isnan(elts));
  }
  constexpr bool any() const { return elts.any(); }

Replacement in Weyl/src/tensor.hxx at line 885 [3.15983]

B:BD[3.43720] → [3.43720:43876]

    gf_Axx_(I) = A(1, 0);
    gf_Axy_(I) = A(1, 1);
    gf_Axz_(I) = A(1, 2);
    gf_Ayy_(I) = A(2, 1);
    gf_Ayz_(I) = A(2, 2);
    gf_Azz_(I) = A(3, 2);

[3.43720]

[3.43876]

    gf_Axx_(I) = A(1, 1);
    gf_Axy_(I) = A(1, 2);
    gf_Axz_(I) = A(1, 3);
    gf_Ayy_(I) = A(2, 2);
    gf_Ayz_(I) = A(2, 3);
    gf_Azz_(I) = A(3, 3);

Replacement in Weyl/src/tensor.hxx at line 1275 [3.15983]

B:BD[3.56393] → [3.56393:56471]

    gf_Axy_(I) = A(1, 1);
    gf_Axz_(I) = A(1, 2);
    gf_Ayz_(I) = A(2, 2);

[3.56393]

[3.56471]

    gf_Axy_(I) = A(1, 2);
    gf_Axz_(I) = A(1, 3);
    gf_Ayz_(I) = A(2, 3);

Insertion in Weyl/src/tensor.hxx at line 1457 [3.15983]

[3.62589]

[2.10081]

template <typename F, typename R = remove_cv_t<remove_reference_t<
                          result_of_t<F(int, int, int, int)> > > >
constexpr R sum44(const F &f) {
  R s = zero<R>()();
  for (int x = 0; x < 4; ++x)
    for (int y = 0; y < 4; ++y)
      for (int z = 0; z < 4; ++z)
        for (int w = 0; w < 4; ++w)
          s += f(x, y, z, w);
  return s;
}

Insertion in Weyl/src/weyl.cxx at line 231 [3.62651]

[3.73244]


  //
  const GF3D<CCTK_REAL, 0, 0, 0> gf_lt_(cctkGH, lt);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_lx_(cctkGH, lx);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_ly_(cctkGH, ly);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_lz_(cctkGH, lz);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nt_(cctkGH, nt);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nx_(cctkGH, nx);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_ny_(cctkGH, ny);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nz_(cctkGH, nz);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_mret_(cctkGH, mret);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_mrex_(cctkGH, mrex);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_mrey_(cctkGH, mrey);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_mrez_(cctkGH, mrez);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_mimt_(cctkGH, mimt);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_mimx_(cctkGH, mimx);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_mimy_(cctkGH, mimy);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_mimz_(cctkGH, mimz);
  //
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Lambda_(cctkGH, Lambda);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Phi00_(cctkGH, Phi00);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Phi11_(cctkGH, Phi11);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Phi22_(cctkGH, Phi22);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Phi10re_(cctkGH, Phi10re);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Phi10im_(cctkGH, Phi10im);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Phi20re_(cctkGH, Phi20re);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Phi20im_(cctkGH, Phi20im);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Phi21re_(cctkGH, Phi21re);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Phi21im_(cctkGH, Phi21im);
  //
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi0re_(cctkGH, Psi0re);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi0im_(cctkGH, Psi0im);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi1re_(cctkGH, Psi1re);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi1im_(cctkGH, Psi1im);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi2re_(cctkGH, Psi2re);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi2im_(cctkGH, Psi2im);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi3re_(cctkGH, Psi3re);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi3im_(cctkGH, Psi3im);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi4re_(cctkGH, Psi4re);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_Psi4im_(cctkGH, Psi4im);
  //
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npkappare_(cctkGH, npkappare);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npkappaim_(cctkGH, npkappaim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npsigmare_(cctkGH, npsigmare);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npsigmaim_(cctkGH, npsigmaim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nprhore_(cctkGH, nprhore);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nprhoim_(cctkGH, nprhoim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nptaure_(cctkGH, nptaure);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nptauim_(cctkGH, nptauim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npepsilonre_(cctkGH, npepsilonre);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npepsilonim_(cctkGH, npepsilonim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npbetare_(cctkGH, npbetare);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npbetaim_(cctkGH, npbetaim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npalphare_(cctkGH, npalphare);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npalphaim_(cctkGH, npalphaim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npgammare_(cctkGH, npgammare);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npgammaim_(cctkGH, npgammaim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nppire_(cctkGH, nppire);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nppiim_(cctkGH, nppiim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npmure_(cctkGH, npmure);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npmuim_(cctkGH, npmuim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nplambdare_(cctkGH, nplambdare);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_nplambdaim_(cctkGH, nplambdaim);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npnure_(cctkGH, npnure);
  const GF3D<CCTK_REAL, 0, 0, 0> gf_npnuim_(cctkGH, npnuim);

Insertion in Weyl/src/weyl.cxx at line 312 [3.62651]
[3.73371]
[3.73371]
```
        const vec3<CCTK_REAL, UP> coord3{p.x, p.y, p.z};
```

Insertion in Weyl/src/weyl.cxx at line 316 [3.62651]

[3.73413]

            cctk_time, coord3,                                 //

Insertion in Weyl/src/weyl.cxx at line 430 [3.62651]

[3.78149]


        vars.l.store(gf_lt_, gf_lx_, gf_ly_, gf_lz_, p.I);
        vars.n.store(gf_nt_, gf_nx_, gf_ny_, gf_nz_, p.I);
        gf_mret_(p.I) = real(vars.m(0));
        gf_mrex_(p.I) = real(vars.m(1));
        gf_mrey_(p.I) = real(vars.m(2));
        gf_mrez_(p.I) = real(vars.m(3));
        gf_mimt_(p.I) = imag(vars.m(0));
        gf_mimx_(p.I) = imag(vars.m(1));
        gf_mimy_(p.I) = imag(vars.m(2));
        gf_mimz_(p.I) = imag(vars.m(3));
        gf_Lambda_(p.I) = vars.Lambda;
        gf_Phi00_(p.I) = vars.Phi00;
        gf_Phi11_(p.I) = vars.Phi11;
        gf_Phi22_(p.I) = vars.Phi22;
        gf_Phi10re_(p.I) = real(vars.Phi10);
        gf_Phi10im_(p.I) = imag(vars.Phi10);
        gf_Phi20re_(p.I) = real(vars.Phi20);
        gf_Phi20im_(p.I) = imag(vars.Phi20);
        gf_Phi21re_(p.I) = real(vars.Phi21);
        gf_Phi21im_(p.I) = imag(vars.Phi21);
        gf_Psi0re_(p.I) = real(vars.Psi0);
        gf_Psi0im_(p.I) = imag(vars.Psi0);
        gf_Psi1re_(p.I) = real(vars.Psi1);
        gf_Psi1im_(p.I) = imag(vars.Psi1);
        gf_Psi2re_(p.I) = real(vars.Psi2);
        gf_Psi2im_(p.I) = imag(vars.Psi2);
        gf_Psi3re_(p.I) = real(vars.Psi3);
        gf_Psi3im_(p.I) = imag(vars.Psi3);
        gf_Psi4re_(p.I) = real(vars.Psi4);
        gf_Psi4im_(p.I) = imag(vars.Psi4);
        gf_npkappare_(p.I) = real(vars.npkappa);
        gf_npkappaim_(p.I) = imag(vars.npkappa);
        gf_npsigmare_(p.I) = real(vars.npsigma);
        gf_npsigmaim_(p.I) = imag(vars.npsigma);
        gf_nprhore_(p.I) = real(vars.nprho);
        gf_nprhoim_(p.I) = imag(vars.nprho);
        gf_nptaure_(p.I) = real(vars.nptau);
        gf_nptauim_(p.I) = imag(vars.nptau);
        gf_npepsilonre_(p.I) = real(vars.npepsilon);
        gf_npepsilonim_(p.I) = imag(vars.npepsilon);
        gf_npbetare_(p.I) = real(vars.npbeta);
        gf_npbetaim_(p.I) = imag(vars.npbeta);
        gf_npalphare_(p.I) = real(vars.npalpha);
        gf_npalphaim_(p.I) = imag(vars.npalpha);
        gf_npgammare_(p.I) = real(vars.npgamma);
        gf_npgammaim_(p.I) = imag(vars.npgamma);
        gf_nppire_(p.I) = real(vars.nppi);
        gf_nppiim_(p.I) = imag(vars.nppi);
        gf_npmure_(p.I) = real(vars.npmu);
        gf_npmuim_(p.I) = imag(vars.npmu);
        gf_nplambdare_(p.I) = real(vars.nplambda);
        gf_nplambdaim_(p.I) = imag(vars.nplambda);
        gf_npnure_(p.I) = real(vars.npnu);
        gf_npnuim_(p.I) = imag(vars.npnu);

Insertion in Weyl/src/weyl_vars.hxx at line 9 [3.78199]
[3.78320]
[3.78320]
```
#include <complex>
```
Insertion in Weyl/src/weyl_vars.hxx at line 14 [3.78199]
[3.78389]
[3.78389]
```
  // Position
  const vec4<T, UP> coord;
```

Replacement in Weyl/src/weyl_vars.hxx at line 33 [3.78199]

B:BD[3.78650] → [3.78650:78717]

  weyl_vars_noderivs(const mat3<T, DN, DN> &gamma, const T &alpha,

[3.78650]

[3.78717]

  // Tetrad
  const vec4<T, UP> et, ephi, etheta, er;
  const vec4<T, UP> l, n;
  const vec4<complex<T>, UP> m;
  weyl_vars_noderivs(const T time, const vec3<T, UP> &coord3,
                     const mat3<T, DN, DN> &gamma, const T &alpha,

Replacement in Weyl/src/weyl_vars.hxx at line 41 [3.78199]

B:BD[3.78763] → [3.78763:78811]

      : gamma(gamma), alpha(alpha), beta(beta),

[3.78763]

[3.78811]

      : coord(time, coord3(0), coord3(1), coord3(2)),
        //
        gamma(gamma), alpha(alpha), beta(beta),

Replacement in Weyl/src/weyl_vars.hxx at line 51 [3.78199]
B:BD[3.79003] → [3.79003:79029]
```
            return alpha;
```
[3.79003]
[3.79029]
```
            return -pow(alpha, 2);
```

Replacement in Weyl/src/weyl_vars.hxx at line 58 [3.78199]

B:BD[3.79173] → [3.79173:79223]

        detg(g.det()), //
        gu(g.inv(detg))

[3.79173]

[3.79223]

        detg(g.det()),                       //
        gu(g.inv(detg)),                     //
        et(calc_et(gu)),                     //
        ephi(calc_ephi(coord, g)),           //
        etheta(calc_etheta(coord, g, ephi)), //
        er(calc_er(coord, g, etheta, ephi)), //
        l([&](int a) { return (et(a) + er(a)) / sqrt(T(2)); }),
        n([&](int a) { return (et(a) - er(a)) / sqrt(T(2)); }),
        m([&](int a) { return complex<T>(etheta(a), ephi(a)) / sqrt(T(2)); })

Insertion in Weyl/src/weyl_vars.hxx at line 75 [3.78199]
[3.80507]
[3.80507]
```
  // Position
  using weyl_vars_noderivs<T>::coord;
```

Insertion in Weyl/src/weyl_vars.hxx at line 92 [3.78199]

[3.80851]


  // Tetrad
  using weyl_vars_noderivs<T>::et;
  using weyl_vars_noderivs<T>::ephi;
  using weyl_vars_noderivs<T>::etheta;
  using weyl_vars_noderivs<T>::er;
  using weyl_vars_noderivs<T>::l;
  using weyl_vars_noderivs<T>::n;
  using weyl_vars_noderivs<T>::m;

Insertion in Weyl/src/weyl_vars.hxx at line 156 [3.78199]

[3.82452]


  // Ricci and Weyl scalars
  const T Lambda;
  const T Phi00, Phi11, Phi22;
  const complex<T> Phi10, Phi20, Phi21;
  const complex<T> Psi0, Psi1, Psi2, Psi3, Psi4;
  // Gradient of tetrad
  const vec4<vec4<T, DN>, UP> det, dephi, detheta, der;
  const vec4<vec4<T, DN>, UP> dl, dn;
  const vec4<vec4<complex<T>, DN>, UP> dm;
  // Newman-Penrose spin coefficients
  const complex<T> npkappa, npsigma, nprho, nptau, npepsilon, npbeta, npalpha,
      npgamma, nppi, npmu, nplambda, npnu;

Insertion in Weyl/src/weyl_vars.hxx at line 173 [3.78199]
[3.82466]
[3.82466]
```
      const T time, const vec3<T, UP> &coord3,
      //
```

Replacement in Weyl/src/weyl_vars.hxx at line 189 [3.78199]

B:BD[3.83109] → [3.83109:83160]

      : weyl_vars_noderivs<T>(gamma, alpha, beta),

[3.83109]

[3.83160]

      : weyl_vars_noderivs<T>(time, coord3, gamma, alpha, beta),

Replacement in Weyl/src/weyl_vars.hxx at line 356 [3.78199]

∅:D[2.17784] → [3.89412:89448]

B:BD[3.89412] → [3.89412:89448]

        C(calc_weyl(g, Rm, R, Rsc))

[2.17784]

[3.94176]

        C(calc_weyl(g, Rm, R, Rsc)),
        //
        // Badri Krishnan's PhD thesis, appendix A
        Lambda(Rsc / 24), //
        Phi00([&] {
          return sum42([&](int a, int b) { return R(a, b) * l(a) * l(b) / 2; });
        }()),
        Phi11([&] {
          return sum42([&](int a, int b) {
            return R(a, b) * (l(a) * n(b) + real(m(a) * conj(m(b)))) / 4;
          });
        }()),
        Phi22([&] {
          return sum42([&](int a, int b) { return R(a, b) * n(a) * n(b) / 2; });
        }()),
        Phi10([&] {
          return sum42(
              [&](int a, int b) { return R(a, b) * l(a) * conj(m(b)) / T(2); });
        }()),
        Phi20([&] {
          return sum42([&](int a, int b) {
            return R(a, b) * conj(m(a)) * conj(m(b)) / T(2);
          });
        }()),
        Phi21([&] {
          return sum42(
              [&](int a, int b) { return R(a, b) * conj(m(a)) * n(b) / T(2); });
        }()),
        Psi0([&] {
          return sum44([&](int a, int b, int c, int d) {
            return C(a, b)(c, d) * l(a) * m(b) * l(c) * m(d);
          });
        }()),
        Psi1([&] {
          return sum44([&](int a, int b, int c, int d) {
            return C(a, b)(c, d) * l(a) * m(b) * l(c) * n(d);
          });
        }()),
        Psi2([&] {
          return sum44([&](int a, int b, int c, int d) {
            return C(a, b)(c, d) * l(a) * m(b) * conj(m(c)) * n(d);
          });
        }()),
        Psi3([&] {
          return sum44([&](int a, int b, int c, int d) {
            return C(a, b)(c, d) * l(a) * n(b) * conj(m(c)) * n(d);
          });
        }()),
        Psi4([&] {
          return sum44([&](int a, int b, int c, int d) {
            return C(a, b)(c, d) * conj(m(a)) * n(b) * conj(m(c)) * n(d);
          });
        }()),
        det(calc_det(gu, dgu, et, Gamma)),                                    //
        dephi(calc_dephi(coord, g, dg, ephi, Gamma)),                         //
        detheta(calc_detheta(coord, g, dg, ephi, dephi, etheta, Gamma)),      //
        der(calc_der(coord, g, dg, ephi, dephi, etheta, detheta, er, Gamma)), //
        dl([&](int a) { return (det(a) + der(a)) / sqrt(T(2)); }),
        dn([&](int a) { return (det(a) - der(a)) / sqrt(T(2)); }),
        dm([&](int a) {
          return vec4<complex<T>, DN>([&](int b) {
            return complex<T>(detheta(a)(b), dephi(a)(b)) / sqrt(T(2));
          });
        }),
        npkappa([&] {
          return sum42([&](int a, int b) { return -m(a) * l(b) * dl(a)(b); });
        }()),
        npsigma([&] {
          return sum42([&](int a, int b) { return -m(a) * m(b) * dl(a)(b); });
        }()),
        nprho([&] {
          return sum42(
              [&](int a, int b) { return -m(a) * conj(m(b)) * dl(a)(b); });
        }()),
        nptau([&] {
          return sum42([&](int a, int b) { return -m(a) * n(b) * dl(a)(b); });
        }()),
        npepsilon([&] {
          return sum42([&](int a, int b) {
            return (conj(m(a)) * l(b) * dm(a)(b) - n(a) * l(b) * dl(a)(b)) /
                   T(2);
          });
        }()),
        npbeta([&] {
          return sum42([&](int a, int b) {
            return (conj(m(a)) * m(b) * dm(a)(b) - n(a) * m(b) * dl(a)(b)) /
                   T(2);
          });
        }()),
        npalpha([&] {
          return sum42([&](int a, int b) {
            return (conj(m(a)) * conj(m(b)) * dm(a)(b) -
                    n(a) * conj(m(b)) * dl(a)(b)) /
                   T(2);
          });
        }()),
        npgamma([&] {
          return sum42([&](int a, int b) {
            return (conj(m(a)) * n(b) * dm(a)(b) - n(a) * n(b) * dl(a)(b)) /
                   T(2);
          });
        }()),
        nppi([&] {
          return sum42(
              [&](int a, int b) { return conj(m(a)) * l(b) * dn(a)(b); });
        }()),
        npmu([&] {
          return sum42(
              [&](int a, int b) { return conj(m(a)) * m(b) * dn(a)(b); });
        }()),
        nplambda([&] {
          return sum42(
              [&](int a, int b) { return conj(m(a)) * conj(m(b)) * dn(a)(b); });
        }()),
        npnu([&] {
          return sum42(
              [&](int a, int b) { return conj(m(a)) * n(b) * dn(a)(b); });
        }())

Weyl: Implement tetrad, NP spin coefficients, Ricci/Weyl scalars

Dependencies

In channels

Change contents

Insertion in Weyl/interface.ccl at line 54 [3.308]

File move: brill-lindquist.par → weyl-brill-lindquist.par

Replacement in Weyl/par/weyl-brill-lindquist.par at line 58 [2.141]

Replacement in Weyl/par/weyl-brill-lindquist.par at line 63 [2.141]

Insertion in Weyl/schedule.ccl at line 9 [3.1799]

Insertion in Weyl/schedule.ccl at line 39 [3.1799]

Insertion in Weyl/src/physics.hxx at line 5 [3.11226]

Insertion in Weyl/src/physics.hxx at line 11 [3.11226]

Insertion in Weyl/src/physics.hxx at line 130 [3.11226]

Insertion in Weyl/src/physics.hxx at line 208 [3.11226]

Insertion in Weyl/src/tensor.hxx at line 4 [3.15983]

Insertion in Weyl/src/tensor.hxx at line 6 [3.15983]