FHE6M7GL6TJSS4VHGV6XY5PQSZBCU6SIXAGZG67Z76NRGKZZAGTQC
ActiveThorns = "
ADMBase
CTGBase
CTGEvolution
CTGGauge
Carpet
CarpetIOASCII
CarpetIOBasic
CarpetIOScalar
CarpetLib
CartGrid3D
CoordBase
CartesianCoordinates
Formaline
GlobalDerivative
IDLinearWaves
IOUtil
MoL
PeriodicCarpet
SphericalSurface
StaticConformal
SummationByParts
SymBase
Time
CarpetReduce
"
Cactus::cctk_show_schedule = yes
# Cactus::terminate = "time"
# Cactus::cctk_final_time = 0.0
Cactus::cctk_itlast = 1
CartGrid3D::type = "coordbase"
CartGrid3D::domain = "full"
CartGrid3D::avoid_origin = no
CoordBase::xmin = 0.0
CoordBase::ymin = 0.0
CoordBase::zmin = 0.0
CoordBase::xmax = 1.0
CoordBase::ymax = 2.0 / 64
CoordBase::zmax = 2.0 / 64
CoordBase::dx = 1.0 / 64
CoordBase::dy = 1.0 / 64
CoordBase::dz = 1.0 / 64
CoordBase::boundary_size_x_lower = 2
CoordBase::boundary_size_y_lower = 2
CoordBase::boundary_size_z_lower = 2
CoordBase::boundary_size_x_upper = 2
CoordBase::boundary_size_y_upper = 2
CoordBase::boundary_size_z_upper = 2
CoordBase::boundary_shiftout_x_lower = 1
CoordBase::boundary_shiftout_y_lower = 1
CoordBase::boundary_shiftout_z_lower = 1
CoordBase::boundary_shiftout_x_upper = 0
CoordBase::boundary_shiftout_y_upper = 0
CoordBase::boundary_shiftout_z_upper = 0
Carpet::domain_from_coordbase = yes
PeriodicCarpet::periodic = yes
driver::ghost_size = 2
ADMBase::initial_data = "sine_planewaves"
ADMBase::initial_lapse = "one"
ADMBase::initial_shift = "zero"
IDLinearWaves::amplitude = 1.0e-8
IDLinearWaves::wavelength = 1.0
IDLinearWaves::wavetheta = 90
ADMBase::evolution_method = "ctgamma"
SummationByParts::order = 2
CTGBase::schedule_rhs_at_analysis = yes
CTGBase::conformal_factor_type = "chi"
CTGBase::evolution_system_type = "Z4c"
CTGEvolution::MaxNumEvolvedVars = 18
CTGEvolution::MaxNumConstrainedVars = 13
CTGEvolution::force_lndetg_zero = yes
CTGEvolution::kappa1 = 0.02
CTGEvolution::kappa2 = 0.0
Time::dtfac = 0.5
MoL::ODE_Method = "Euler"
MoL::MoL_Intermediate_Steps = 1
IO::out_dir = $parfile
IO::out_every = 1
IOASCII::out1d_vars = "
ADMBase::metric
ADMBase::curv
ADMBase::lapse
ADMBase::shift
ADMBase::dtlapse
ADMBase::dtshift
CTGBase::conformal_factor
CTGBase::conformal_metric
CTGBase::curvature_scalar
CTGBase::curvature_scalar_a
CTGBase::curvature_scalar_b
CTGBase::curvature_tensor
CTGBase::Gamma
CTGEvolution::conformal_factor_source
CTGEvolution::conformal_metric_source
CTGEvolution::curvature_scalar_source
CTGEvolution::Theta_source
CTGEvolution::curvature_tensor_source
CTGEvolution::Gamma_source
CTGGauge::shift_aux
CTGGauge::shift_aux_source
"
# 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
mat3<CCTK_REAL> gammat(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_,
gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p.I);
mat3<CCTK_REAL> At(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 = gammat.det();
const CCTK_REAL chi_new = 1 / cbrt(detgammat);
for (int a = 0; a < 3; ++a)
for (int b = a; b < 3; ++b)
gammat(a, b) *= chi_new;
const CCTK_REAL detgammat_old = gammat_old.det();
const CCTK_REAL chi_old = 1 / cbrt(detgammat_old);
const mat3<CCTK_REAL, DN, DN> gammat(
[&](int a, int b) { return chi_old * gammat_old(a, b); });
const CCTK_REAL traceAt =
sum2([&](int x, int y) { return gammatu(x, y) * At(x, y); });
for (int a = 0; a < 3; ++a)
for (int b = a; b < 3; ++b)
At(a, b) -= traceAt / 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 mat3<CCTK_REAL> g(gf_gxx_, gf_gxy_, gf_gxz_, gf_gyy_, gf_gyz_,
gf_gzz_, p.I);
const mat3<CCTK_REAL> k(gf_kxx_, gf_kxy_, gf_kxz_, gf_kyy_, gf_kyz_,
gf_kzz_, p.I);
const mat3<CCTK_REAL, DN, DN> g(gf_gxx_, gf_gxy_, gf_gxz_, gf_gyy_, gf_gyz_,
gf_gzz_, p.I);
const mat3<CCTK_REAL, DN, DN> k(gf_kxx_, gf_kxy_, gf_kxz_, gf_kyy_, gf_kyz_,
gf_kzz_, p.I);
const mat3<CCTK_REAL> gammat(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_,
gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p.I);
const mat3<CCTK_REAL, DN, DN> gammat(gf_gammatxx_, gf_gammatxy_,
gf_gammatxz_, gf_gammatyy_,
gf_gammatyz_, gf_gammatzz_, p.I);
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) {
const vec3<mat3<CCTK_REAL, DN, DN>, DN> Gammatl = calc_gammal(dgammat);
const vec3<mat3<CCTK_REAL, DN, DN>, UP> Gammat =
calc_gamma(gammatu, Gammatl);
const vec3<CCTK_REAL, UP> Gamt([&](int a) {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<vec3<T> >
calc_dgu(const mat3<T> &gu, const mat3<vec3<T> > &dg) {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<vec3<T, DN>, UP, UP>
calc_dgu(const mat3<T, UP, UP> &gu, const mat3<vec3<T, DN>, DN, DN> &dg) {
/*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) {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<vec3<T, DN>, UP, UP>
calc_dAu(const mat3<T, UP, UP> &gu, const mat3<vec3<T, DN>, UP, UP> &dgu,
const mat3<T, DN, DN> &A, const mat3<vec3<T, DN>, DN, DN> &dA) {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<mat3<T> >
calc_gamma(const mat3<T> &gu, const vec3<mat3<T> > &Gammal) {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<mat3<T, DN, DN>, UP>
calc_gamma(const mat3<T, UP, UP> &gu, const vec3<mat3<T, DN, DN>, DN> &Gammal) {
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));
const vec3<CCTK_REAL, UP> betaG(gf_betaGx_, gf_betaGy_, gf_betaGz_, p.I);
const vec3<CCTK_REAL, DN> dgf_upwind(deriv_upwind(gf_, p.I, betaG, dx));
template <typename T> struct nan<vec3<T> > {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<T> operator()() const {
return vec3<T>();
template <typename T, dnup_t dnup> struct nan<vec3<T, dnup> > {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr vec3<T, dnup> operator()() const {
return vec3<T, dnup>();
assert(!CCTK_isnan(A(0, 0)));
assert(!CCTK_isnan(A(0, 1)));
assert(!CCTK_isnan(A(0, 2)));
assert(!CCTK_isnan(A(1, 1)));
assert(!CCTK_isnan(A(1, 2)));
assert(!CCTK_isnan(A(2, 2)));
assert(CCTK_isfinite(A(0, 0)));
assert(CCTK_isfinite(A(0, 1)));
assert(CCTK_isfinite(A(0, 2)));
assert(CCTK_isfinite(A(1, 1)));
assert(CCTK_isfinite(A(1, 2)));
assert(CCTK_isfinite(A(2, 2)));
friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T>
operator+(const mat3<T> &x, const mat3<T> &y) {
friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T, dnup1, dnup2>
operator+(const mat3<T, dnup1, dnup2> &x, const mat3<T, dnup1, dnup2> &y) {
friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T>
operator-(const mat3<T> &x, const mat3<T> &y) {
friend /*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T, dnup1, dnup2>
operator-(const mat3<T, dnup1, dnup2> &x, const mat3<T, dnup1, dnup2> &y) {
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)),
detA1 * (A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0))};
return mat3<T, !dnup1, !dnup2>{
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)),
detA1 * (A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0))};
template <typename T> struct nan<mat3<T> > {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T> operator()() const {
return mat3<T>();
template <typename T, dnup_t dnup1, dnup_t dnup2>
struct nan<mat3<T, dnup1, dnup2> > {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T, dnup1, dnup2>
operator()() const {
return mat3<T, dnup1, dnup2>();
template <typename T>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T> mul(const mat3<T> &A,
const mat3<T> &B) {
template <typename T, dnup_t dnup1, dnup_t dnup2, dnup_t dnup3>
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ constexpr mat3<T, dnup1, dnup2>
mul(const mat3<T, dnup1, dnup3> &A, const mat3<T, !dnup3, dnup2> &B) {
return (deriv(gf_, I + DI[dir2], dx, dir1) -
deriv(gf_, I - DI[dir2], dx, dir1)) /
(2 * dx(dir2));
// return (deriv(gf_, I + DI[dir2], dx, dir1) -
// deriv(gf_, I - DI[dir2], dx, dir1)) /
// (2 * dx(dir2));
return ((gf_(I + DI[dir1] + DI[dir2]) //
- gf_(I - DI[dir1] + DI[dir2])) //
- (gf_(I + DI[dir1] - DI[dir2]) //
- gf_(I - DI[dir1] - DI[dir2]))) //
/ (4 * dx(dir1) * dx(dir2));
/*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) {
/*CCTK_ATTRIBUTE_ALWAYS_INLINE*/ T
deriv4(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,
const vec3<T, UP> &dx, const int dir) {
const mat3<double> Z([&](int a, int b) { return double(0); });
const mat3<double> I([&](int a, int b) { return double(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); });
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
const T chi; // W = -2/3
const mat3<T, DN, DN> gammat; // W = -2/3
const T Kh; // W = 0
const mat3<T, DN, DN> At; // W = -2/3
const vec3<T, UP> Gamt; // W = +2/3
const T Theta; // W = 0
const T alphaG; // W = 0
const vec3<T, UP> betaG; // W = 0
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 &chi, const mat3<T, DN, DN> &gammat, const T &Kh,
const mat3<T, DN, DN> &At, const vec3<T, UP> &Gamt,
const T &Theta, const T &alphaG, const vec3<T, UP> &betaG,
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;
const vec3<T, DN> dchi;
const mat3<T, DN, DN> ddchi;
const mat3<vec3<T, DN>, DN, DN> dgammat;
const mat3<mat3<T, DN, DN>, DN, DN> ddgammat;
const vec3<T, DN> dKh;
const mat3<vec3<T, DN>, DN, DN> dAt;
const vec3<vec3<T, DN>, UP> dGamt;
const vec3<T, DN> dTheta;
const vec3<T, DN> dalphaG;
const mat3<T, DN, DN> ddalphaG;
const vec3<vec3<T, DN>, UP> dbetaG;
const vec3<mat3<T, DN, DN>, UP> ddbetaG;
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 mat3<T, UP, UP> gammatu;
const mat3<vec3<T, DN>, UP, UP> dgammatu;
const vec3<mat3<T, DN, DN>, DN> Gammatl;
const vec3<mat3<T, DN, DN>, UP> Gammat;
const vec3<T, UP> Gamtd;
const mat3<T, DN, DN> DDchi;
const mat3<T, UP, UP> gu;
const mat3<vec3<T, DN>, DN, DN> dg;
const vec3<mat3<T, DN, DN>, DN> Gammal;
const vec3<mat3<T, DN, DN>, UP> Gamma;
const mat3<T, DN, DN> DDalphaG;
const mat3<T, DN, DN> Rchi;
const mat3<T, DN, DN> Rt;
const mat3<T, DN, DN> R;
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 &chi, const vec3<T, DN> &dchi,
const mat3<T, DN, DN> &ddchi, //
const mat3<T, DN, DN> &gammat,
const mat3<vec3<T, DN>, DN, DN> &dgammat,
const mat3<mat3<T, DN, DN>, DN, DN> &ddgammat, //
const T &Kh, const vec3<T, DN> &dKh, //
const mat3<T, DN, DN> &At, const mat3<vec3<T, DN>, DN, DN> &dAt, //
const vec3<T, UP> &Gamt, const vec3<vec3<T, DN>, UP> &dGamt, //
const T &Theta, const vec3<T, DN> &dTheta, //
const T &alphaG, const vec3<T, DN> &dalphaG,
const mat3<T, DN, DN> &ddalphaG, //
const vec3<T, UP> &betaG, const vec3<vec3<T, DN>, UP> &dbetaG,
const vec3<mat3<T, DN, DN>, UP> &ddbetaG,
// return ddalphaG(a, b) //
// - sum1([&](int x) { return Gammat(x)(a, b) * dalphaG(x); })
// //
// // chipsipower = -4
// // oochipsipower = -1
// // df = oochipsipower dchi(a) / chi
// // ddf == oochipsipower ddchi / chi - (-4) df df
// + 1 / (2 * chi) *
// (dchi(a) * dalphaG(b) //
// + dchi(b) * dalphaG(a)) //
// + 1 / (4 * chi) * sum2([&](int x, int y) {
// return gammat(a, b) * gammatu(x, y) * dchi(x) *
// dalphaG(y);
// });