67ABKXMOSVKKVRSSD6QVF4BLBQLZCEAST2UMLLKOBV7WSVTQRSEQC
const mat3<CCTK_REAL, DN, DN> gammat_old(
gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,
gf_gammatyz_, gf_gammatzz_, p.I);
const mat3<CCTK_REAL, DN, DN> At_old(gf_Atxx_, gf_Atxy_, gf_Atxz_,
gf_Atyy_, gf_Atyz_, gf_Atzz_, p.I);
const mat3<CCTK_REAL, DN, DN> gammat_old(gf_gammatxx_, gf_gammatxy_,
gf_gammatxz_, gf_gammatyy_,
gf_gammatyz_, gf_gammatzz_, p.I);
const mat3<CCTK_REAL, DN, DN> At_old(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_,
gf_Atyz_, gf_Atzz_, p.I);
const CCTK_REAL detgammat_old = gammat_old.det();
const CCTK_REAL chi1_old = 1 / cbrt(detgammat_old);
const mat3<CCTK_REAL, DN, DN> gammat(
[&](int a, int b) { return chi1_old * gammat_old(a, b); });
const CCTK_REAL detgammat_old = gammat_old.det();
const CCTK_REAL chi1_old = 1 / cbrt(detgammat_old);
const mat3<CCTK_REAL, DN, DN> gammat(
[&](int a, int b) { return chi1_old * gammat_old(a, b); });
const CCTK_REAL detgammat = gammat.det();
const CCTK_REAL gammat_norm = gammat.maxabs();
const CCTK_REAL gammat_scale = gammat_norm;
if (!(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale)) {
ostringstream buf;
buf << "det gammat is not one: gammat=" << gammat
<< " det(gammat)=" << detgammat;
CCTK_VERROR("%s", buf.str().c_str());
}
assert(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale);
const CCTK_REAL detgammat = gammat.det();
const CCTK_REAL gammat_norm = gammat.maxabs();
const CCTK_REAL gammat_scale = gammat_norm;
if (!(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale)) {
ostringstream buf;
buf << "det gammat is not one: gammat=" << gammat
<< " det(gammat)=" << detgammat;
CCTK_VERROR("%s", buf.str().c_str());
}
assert(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale);
const CCTK_REAL traceAt_old =
sum2([&](int x, int y) { return gammatu(x, y) * At_old(x, y); });
const mat3<CCTK_REAL, DN, DN> At([&](int a, int b) {
return At_old(a, b) - traceAt_old / 3 * gammat(a, b);
});
const CCTK_REAL traceAt_old =
sum2([&](int x, int y) { return gammatu(x, y) * At_old(x, y); });
const mat3<CCTK_REAL, DN, DN> At([&](int a, int b) {
return At_old(a, b) - traceAt_old / 3 * gammat(a, b);
});
const CCTK_REAL traceAt =
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) <= 1.0e-12 * At_scale)) {
ostringstream buf;
buf << "tr At: At=" << At << " tr(At)=" << traceAt;
CCTK_VERROR("%s", buf.str().c_str());
}
assert(fabs(traceAt) <= 1.0e-12 * At_scale);
const CCTK_REAL traceAt =
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) <= 1.0e-12 * At_scale)) {
ostringstream buf;
buf << "tr At: At=" << At << " tr(At)=" << traceAt;
CCTK_VERROR("%s", buf.str().c_str());
}
assert(fabs(traceAt) <= 1.0e-12 * At_scale);
// Store
gf_chi_(p.I) = chi;
gammat.store(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,
gf_gammatyz_, gf_gammatzz_, p.I);
At.store(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_, gf_Atzz_,
p.I);
gf_alphaG_(p.I) = alphaG;
});
// Store
gf_chi_(p.I) = chi;
gammat.store(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,
gf_gammatyz_, gf_gammatzz_, p.I);
At.store(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_, gf_Atzz_, p.I);
gf_alphaG_(p.I) = alphaG;
});
vars.gammat_rhs.store(gf_gammatxx_rhs_, gf_gammatxy_rhs_,
gf_gammatxz_rhs_, gf_gammatyy_rhs_,
gf_gammatyz_rhs_, gf_gammatzz_rhs_, p.I);
vars.gammat_rhs.store(gf_gammatxx_rhs_, gf_gammatxy_rhs_, gf_gammatxz_rhs_,
gf_gammatyy_rhs_, gf_gammatyz_rhs_, gf_gammatzz_rhs_,
p.I);