GP6UYE5IVVG25EHEMT4EPKFUUWI2QUZSJCRSHOIVTM4AX4IVXBLAC Cactus Code Thorn Z4cAuthor(s) : Erik Schnetter <schnetter@gmail.com>Maintainer(s): Erik Schnetter <schnetter@gmail.com>Licence : LGPL--------------------------------------------------------------------------1. PurposeSolve the Einstein equations in the Z4c formulationDavid Hilditch, Sebastiano Bernuzzi, Marcus Thierfelder, Zhoujian Cao,Wolfgang Tichy, Bernd Brügmann: Compact binary evolutions with the Z4cformulation, Phys. Rev. D 88, 084057 (2013), DOI:10.1103/PhysRevD.88.084057, arXiv:1212.2901 [gr-qc].
# Interface definition for thorn Z4cIMPLEMENTS: Z4cINHERITS: ADMBaseUSES INCLUDE HEADER: loop.hxxCCTK_REAL chi TYPE=gf TAGS='index={0 0 0} rhs="chi_rhs"' "chi"CCTK_REAL gamma_tilde TYPE=gf TAGS='index={0 0 0} rhs="gamma_tilde_rhs"' { gammatxx gammatxy gammatxz gammatyy gammatyz gammatzz } "gamma-tilde"CCTK_REAL K_hat TYPE=gf TAGS='index={0 0 0} rhs="K_hat_rhs"' { Kh } "K-hat"CCTK_REAL A_tilde TYPE=gf TAGS='index={0 0 0} rhs="A_tilde_rhs"' { Atxx Atxy Atxz Atyy Atyz Atzz } "A-tilde"CCTK_REAL Gam_tilde TYPE=gf TAGS='index={0 0 0} rhs="Gam_tilde_rhs"' { Gamtx Gamty Gamtz } "Gamma-tilde"CCTK_REAL Theta TYPE=gf TAGS='index={0 0 0} rhs="Theta_rhs"' "Theta"CCTK_REAL alphaG TYPE=gf TAGS='index={0 0 0} rhs="alphaG_rhs"' "alpha"CCTK_REAL betaG TYPE=gf TAGS='index={0 0 0} rhs="betaG_rhs"' { betaGx betaGy betaGz } "beta"CCTK_REAL ZtC TYPE=gf TAGS='index={0 0 0} checkpoint="no"' { ZtCx ZtCy ZtCz } "Z-tilde"CCTK_REAL HC TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "H"CCTK_REAL MtC TYPE=gf TAGS='index={0 0 0} checkpoint="no"' { MtCx MtCy MtCz } "M-tilde"CCTK_REAL chi_rhs TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "chi"CCTK_REAL gamma_tilde_rhs TYPE=gf TAGS='index={0 0 0} checkpoint="no"' { gammatxx_rhs gammatxy_rhs gammatxz_rhs gammatyy_rhs gammatyz_rhs gammatzz_rhs } "gamma-tilde"CCTK_REAL K_hat_rhs TYPE=gf TAGS='index={0 0 0} checkpoint="no"' { Kh_rhs } "K-hat"CCTK_REAL A_tilde_rhs TYPE=gf TAGS='index={0 0 0} checkpoint="no"' { Atxx_rhs Atxy_rhs Atxz_rhs Atyy_rhs Atyz_rhs Atzz_rhs } "A-tilde"CCTK_REAL Gam_tilde_rhs TYPE=gf TAGS='index={0 0 0} checkpoint="no"' { Gamtx_rhs Gamty_rhs Gamtz_rhs } "Gamma-tilde"CCTK_REAL Theta_rhs TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "Theta"CCTK_REAL alphaG_rhs TYPE=gf TAGS='index={0 0 0} checkpoint="no"' "alpha"CCTK_REAL betaG_rhs TYPE=gf TAGS='index={0 0 0} checkpoint="no"' { betaGx_rhs betaGy_rhs betaGz_rhs } "beta"
ActiveThorns = "ADMBaseBaikalXBrillLindquistCarpetXErrorEstimatorFormalineIOUtilNewRadODESolvers"$nlevels = 1 # 2 #TODO 8$ncells = 128 #TODO 256 #TODO 32Cactus::cctk_show_schedule = yesCactus::terminate = "time"Cactus::cctk_final_time = 1.0CarpetX::verbose = yesCarpetX::xmin = -16.0CarpetX::ymin = -16.0CarpetX::zmin = -16.0CarpetX::xmax = +16.0CarpetX::ymax = +16.0CarpetX::zmax = +16.0CarpetX::ncells_x = $ncellsCarpetX::ncells_y = $ncellsCarpetX::ncells_z = $ncellsCarpetX::periodic_x = noCarpetX::periodic_y = noCarpetX::periodic_z = noCarpetX::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 = yesCarpetX::ghost_size = 4CarpetX::max_num_levels = $nlevelsCarpetX::regrid_every = 16CarpetX::regrid_error_threshold = 1.0 #TODO 1.0 / 16.0ErrorEstimator::region_shape = "cube"ErrorEstimator::scale_by_resolution = no #TODO yesCarpetX::prolongation_type = "ddf"CarpetX::prolongation_order = 5ODESolvers::method = "RK4"CarpetX::dtfac = 0.25ADMBase::initial_data = "Brill-Lindquist"ADMBase::initial_lapse = "Brill-Lindquist"IO::out_dir = $parfileIO::out_every = 1 #TODO $ncells * 2 ** ($nlevels - 1) / 32CarpetX::out_tsv = no
ActiveThorns = "ADMBaseCarpetXErrorEstimatorFormalineIOUtilODESolversZ4c"$nlevels = 1$ncells = 32Cactus::cctk_show_schedule = yesCactus::terminate = "time"Cactus::cctk_final_time = 0.0 #TODO 1.0CarpetX::verbose = yesCarpetX::xmin = 0.0CarpetX::ymin = 0.0CarpetX::zmin = 0.0CarpetX::xmax = 1.0CarpetX::ymax = 1.0CarpetX::zmax = 1.0CarpetX::ncells_x = $ncellsCarpetX::ncells_y = $ncellsCarpetX::ncells_z = $ncellsCarpetX::periodic_x = noCarpetX::periodic_y = noCarpetX::periodic_z = noCarpetX::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 = yesCarpetX::ghost_size = 1CarpetX::max_num_levels = $nlevelsCarpetX::regrid_every = 16CarpetX::regrid_error_threshold = 1.0 #TODO 1.0 / 16.0ErrorEstimator::region_shape = "cube"ErrorEstimator::scale_by_resolution = no #TODO yesCarpetX::prolongation_type = "ddf"CarpetX::prolongation_order = 1ODESolvers::method = "RK2"CarpetX::dtfac = 0.25ADMBase::initial_data = "Cartesian Minkowski"ADMBase::initial_lapse = "one"ADMBase::initial_shift = "zero"IO::out_dir = $parfileIO::out_every = 1 #TODO $ncells * 2 ** ($nlevels - 1) / 32CarpetX::out_tsv = no
# Parameter definitions for thorn Z4cCCTK_REAL kappa1 "kappa1" STEERABLE=always{*:* :: ""} 0.02CCTK_REAL kappa2 "kappa2" STEERABLE=always{*:* :: ""} 0.0CCTK_REAL f_mu_L "mu_L = f_mu_L / alpha" STEERABLE=always{*:* :: ""} 2.0CCTK_REAL f_mu_S "mu_S = f_mu_S / alpha^2" STEERABLE=always{*:* :: ""} 1.0CCTK_REAL eta "eta" STEERABLE=always{*:* :: ""} 2.0
# Schedule definitions for thorn Z4cSTORAGE: chiSTORAGE: gamma_tildeSTORAGE: K_hatSTORAGE: A_tildeSTORAGE: Gam_tildeSTORAGE: ThetaSTORAGE: alphaGSTORAGE: betaGSTORAGE: ZtCSTORAGE: HCSTORAGE: MtCSTORAGE: chi_rhsSTORAGE: gamma_tilde_rhsSTORAGE: K_hat_rhsSTORAGE: A_tilde_rhsSTORAGE: Gam_tilde_rhsSTORAGE: Theta_rhsSTORAGE: alphaG_rhsSTORAGE: betaG_rhsSCHEDULE Z4c_Initial1 AT initial AFTER ADMBase_PostInitial{LANG: CREADS: ADMBase::metric(everywhere)READS: ADMBase::curv(everywhere)READS: ADMBase::lapse(everywhere)READS: ADMBase::shift(everywhere)WRITES: chi(everywhere)WRITES: gamma_tilde(everywhere)WRITES: K_hat(everywhere)WRITES: A_tilde(everywhere)WRITES: Theta(everywhere)WRITES: alphaG(everywhere)WRITES: betaG(everywhere)SYNC: chiSYNC: gamma_tildeSYNC: K_hatSYNC: A_tildeSYNC: ThetaSYNC: alphaGSYNC: betaG} "Convert ADM to Z4c variables, part 1"SCHEDULE Z4c_Initial2 AT initial AFTER Z4c_Initial1{LANG: CREADS: gamma_tilde(everywhere)WRITES: Gam_tilde(interior)SYNC: Gam_tilde} "Convert ADM to Z4c variables, part 2"SCHEDULE Z4c_Boundaries AT initial AFTER Z4c_Initial2{LANG: CWRITES: chi(boundary)WRITES: gamma_tilde(boundary)WRITES: K_hat(boundary)WRITES: A_tilde(boundary)WRITES: Gam_tilde(boundary)WRITES: Theta(boundary)WRITES: alphaG(boundary)WRITES: betaG(boundary)} "Apply boundary conditions to Z4c variables"SCHEDULE Z4c_Enforce AT initial AFTER Z4c_Boundaries{LANG: CREADS: gamma_tilde(everywhere)READS: A_tilde(everywhere)WRITES: gamma_tilde(everywhere)WRITES: A_tilde(everywhere)} "Enforce algebraic Z4c constraints"SCHEDULE Z4c_Constraints AT poststep{LANG: CREADS: chi(everywhere)READS: gamma_tilde(everywhere)READS: K_hat(everywhere)READS: A_tilde(everywhere)READS: Gam_tilde(everywhere)READS: Theta(everywhere)WRITES: ZtC(interior)WRITES: HC(interior)WRITES: MtC(interior)} "Calculate Z4c constraints"SCHEDULE Z4c_ADM AT poststep{LANG: CREADS: 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{LANG: CREADS: chi(everywhere)READS: gamma_tilde(everywhere)READS: K_hat(everywhere)READS: A_tilde(everywhere)READS: Gam_tilde(everywhere)READS: Theta(everywhere)READS: alphaG(everywhere)READS: betaG(everywhere)WRITES: chi_rhs(interior)WRITES: gamma_tilde_rhs(interior)WRITES: K_hat_rhs(interior)WRITES: A_tilde_rhs(interior)WRITES: Gam_tilde_rhs(interior)WRITES: Theta_rhs(interior)WRITES: alphaG_rhs(interior)WRITES: betaG_rhs(interior)SYNC: chi_rhsSYNC: gamma_tilde_rhsSYNC: K_hat_rhsSYNC: A_tilde_rhsSYNC: Gam_tilde_rhsSYNC: Theta_rhsSYNC: alphaG_rhsSYNC: betaG_rhs} "Calculate Z4c RHS"
#include "tensor.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cmath>namespace Z4c {using namespace Loop;using namespace std;extern "C" void Z4c_ADM(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Z4c_ADM;const GF3D<const CCTK_REAL, 0, 0, 0> gf_chi_(cctkGH, chi);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxx_(cctkGH, gammatxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxy_(cctkGH, gammatxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxz_(cctkGH, gammatxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyy_(cctkGH, gammatyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyz_(cctkGH, gammatyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatzz_(cctkGH, gammatzz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Kh_(cctkGH, Kh);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxx_(cctkGH, Atxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxy_(cctkGH, Atxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxz_(cctkGH, Atxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyy_(cctkGH, Atyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyz_(cctkGH, Atyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atzz_(cctkGH, Atzz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Theta_(cctkGH, Theta);const GF3D<const CCTK_REAL, 0, 0, 0> gf_alphaG_(cctkGH, alphaG);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaGx_(cctkGH, betaGx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaGy_(cctkGH, betaGy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaGz_(cctkGH, betaGz);const GF3D<CCTK_REAL, 0, 0, 0> gf_gxx_(cctkGH, gxx);const GF3D<CCTK_REAL, 0, 0, 0> gf_gxy_(cctkGH, gxy);const GF3D<CCTK_REAL, 0, 0, 0> gf_gxz_(cctkGH, gxz);const GF3D<CCTK_REAL, 0, 0, 0> gf_gyy_(cctkGH, gyy);const GF3D<CCTK_REAL, 0, 0, 0> gf_gyz_(cctkGH, gyz);const GF3D<CCTK_REAL, 0, 0, 0> gf_gzz_(cctkGH, gzz);const GF3D<CCTK_REAL, 0, 0, 0> gf_kxx_(cctkGH, kxx);const GF3D<CCTK_REAL, 0, 0, 0> gf_kxy_(cctkGH, kxy);const GF3D<CCTK_REAL, 0, 0, 0> gf_kxz_(cctkGH, kxz);const GF3D<CCTK_REAL, 0, 0, 0> gf_kyy_(cctkGH, kyy);const GF3D<CCTK_REAL, 0, 0, 0> gf_kyz_(cctkGH, kyz);const GF3D<CCTK_REAL, 0, 0, 0> gf_kzz_(cctkGH, kzz);const GF3D<CCTK_REAL, 0, 0, 0> gf_alp_(cctkGH, alp);const GF3D<CCTK_REAL, 0, 0, 0> gf_betax_(cctkGH, betax);const GF3D<CCTK_REAL, 0, 0, 0> gf_betay_(cctkGH, betay);const GF3D<CCTK_REAL, 0, 0, 0> gf_betaz_(cctkGH, betaz);const GF3D<CCTK_REAL, 0, 0, 0> gf_dtbetax_(cctkGH, dtbetax);const GF3D<CCTK_REAL, 0, 0, 0> gf_dtbetay_(cctkGH, dtbetay);const GF3D<CCTK_REAL, 0, 0, 0> gf_dtbetaz_(cctkGH, dtbetaz);loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) {// Loadconst 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 variablesmat3<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);// Storeg.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);});loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_dtbetax_(p.I) = 0; });loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_dtbetay_(p.I) = 0; });loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_dtbetaz_(p.I) = 0; });}} // namespace Z4c
#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>namespace Z4c {using namespace Loop;extern "C" void Z4c_Boundaries(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Z4c_Boundaries;const GF3D<CCTK_REAL, 0, 0, 0> gf_chi_(cctkGH, chi);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxx_(cctkGH, gammatxx);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxy_(cctkGH, gammatxy);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxz_(cctkGH, gammatxz);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatyy_(cctkGH, gammatyy);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatyz_(cctkGH, gammatyz);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatzz_(cctkGH, gammatzz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Kh_(cctkGH, Kh);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxx_(cctkGH, Atxx);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxy_(cctkGH, Atxy);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxz_(cctkGH, Atxz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atyy_(cctkGH, Atyy);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atyz_(cctkGH, Atyz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atzz_(cctkGH, Atzz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Gamtx_(cctkGH, Gamtx);const GF3D<CCTK_REAL, 0, 0, 0> gf_Gamty_(cctkGH, Gamty);const GF3D<CCTK_REAL, 0, 0, 0> gf_Gamtz_(cctkGH, Gamtz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Theta_(cctkGH, Theta);const GF3D<CCTK_REAL, 0, 0, 0> gf_alphaG_(cctkGH, alphaG);const GF3D<CCTK_REAL, 0, 0, 0> gf_betaGx_(cctkGH, betaGx);const GF3D<CCTK_REAL, 0, 0, 0> gf_betaGy_(cctkGH, betaGy);const GF3D<CCTK_REAL, 0, 0, 0> gf_betaGz_(cctkGH, betaGz);loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_chi_(p.I) = 1; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_gammatxx_(p.I) = 1; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_gammatxy_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_gammatxz_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_gammatyy_(p.I) = 1; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_gammatyz_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_gammatzz_(p.I) = 1; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Kh_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atxx_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atxy_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atxz_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atyy_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atyz_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Atzz_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Gamtx_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Gamty_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Gamtz_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_Theta_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_alphaG_(p.I) = 1; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_betaGx_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_betaGy_(p.I) = 0; });loop_bnd<0, 0, 0>(cctkGH, [&](const PointDesc &p) { gf_betaGz_(p.I) = 0; });}} // namespace Z4c
#include "physics.hxx"#include "tensor.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cmath>namespace Z4c {using namespace Loop;using namespace std;extern "C" void Z4c_Constraints(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Z4c_Constraints;const vec3<CCTK_REAL> dx{CCTK_DELTA_SPACE(0),CCTK_DELTA_SPACE(1),CCTK_DELTA_SPACE(2),};const GF3D<const CCTK_REAL, 0, 0, 0> gf_chi_(cctkGH, chi);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxx_(cctkGH, gammatxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxy_(cctkGH, gammatxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxz_(cctkGH, gammatxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyy_(cctkGH, gammatyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyz_(cctkGH, gammatyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatzz_(cctkGH, gammatzz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Kh_(cctkGH, Kh);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxx_(cctkGH, Atxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxy_(cctkGH, Atxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxz_(cctkGH, Atxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyy_(cctkGH, Atyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyz_(cctkGH, Atyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atzz_(cctkGH, Atzz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamtx_(cctkGH, Gamtx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamty_(cctkGH, Gamty);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamtz_(cctkGH, Gamtz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Theta_(cctkGH, Theta);const GF3D<CCTK_REAL, 0, 0, 0> gf_ZtCx_(cctkGH, ZtCx);const GF3D<CCTK_REAL, 0, 0, 0> gf_ZtCy_(cctkGH, ZtCy);const GF3D<CCTK_REAL, 0, 0, 0> gf_ZtCz_(cctkGH, ZtCz);const GF3D<CCTK_REAL, 0, 0, 0> gf_HC_(cctkGH, HC);const GF3D<CCTK_REAL, 0, 0, 0> gf_MtCx_(cctkGH, MtCx);const GF3D<CCTK_REAL, 0, 0, 0> gf_MtCy_(cctkGH, MtCy);const GF3D<CCTK_REAL, 0, 0, 0> gf_MtCz_(cctkGH, MtCz);loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {// Loadconst 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;// (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;}// StoreZtC.store(gf_ZtCx_, gf_ZtCy_, gf_ZtCz_, p);gf_HC_(p.I) = HC;MtC.store(gf_MtCx_, gf_MtCy_, gf_MtCz_, p);});}} // namespace Z4c
#include "physics.hxx"#include "tensor.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cmath>namespace Z4c {using namespace Loop;using namespace std;extern "C" void Z4c_Enforce(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Z4c_Enforce;const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxx_(cctkGH, gammatxx);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxy_(cctkGH, gammatxy);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxz_(cctkGH, gammatxz);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatyy_(cctkGH, gammatyy);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatyz_(cctkGH, gammatyz);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatzz_(cctkGH, gammatzz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxx_(cctkGH, Atxx);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxy_(cctkGH, Atxy);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxz_(cctkGH, Atxz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atyy_(cctkGH, Atyy);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atyz_(cctkGH, Atyz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atzz_(cctkGH, Atzz);loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {// Loadmat3<CCTK_REAL> gammat(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_,gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p);mat3<CCTK_REAL> At(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_,gf_Atzz_, p);// Enforce algebraic constraints// See arXiv:1212.2901 [gr-qc].const CCTK_REAL det_gammat = gammat.det();const CCTK_REAL psi1 = 1 / cbrt(det_gammat);for (int a = 0; a < 3; ++a)for (int b = a; b < 3; ++b)gammat(a, b) *= psi1;const mat3<CCTK_REAL> gammatu = gammat.inv(1);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);for (int a = 0; a < 3; ++a)for (int b = a; b < 3; ++b)At(a, b) -= 1 / 3 * traceAt * gammat(a, b);// Storegammat.store(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, p);At.store(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_, gf_Atzz_, p);});}} // namespace Z4c
#include "tensor.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cmath>namespace Z4c {using namespace Loop;using namespace std;extern "C" void Z4c_Initial1(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Z4c_Initial1;const GF3D<const CCTK_REAL, 0, 0, 0> gf_gxx_(cctkGH, gxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gxy_(cctkGH, gxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gxz_(cctkGH, gxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gyy_(cctkGH, gyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gyz_(cctkGH, gyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gzz_(cctkGH, gzz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_kxx_(cctkGH, kxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_kxy_(cctkGH, kxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_kxz_(cctkGH, kxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_kyy_(cctkGH, kyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_kyz_(cctkGH, kyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_kzz_(cctkGH, kzz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_alp_(cctkGH, alp);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betax_(cctkGH, betax);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betay_(cctkGH, betay);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaz_(cctkGH, betaz);const GF3D<CCTK_REAL, 0, 0, 0> gf_chi_(cctkGH, chi);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxx_(cctkGH, gammatxx);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxy_(cctkGH, gammatxy);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatxz_(cctkGH, gammatxz);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatyy_(cctkGH, gammatyy);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatyz_(cctkGH, gammatyz);const GF3D<CCTK_REAL, 0, 0, 0> gf_gammatzz_(cctkGH, gammatzz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Kh_(cctkGH, Kh);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxx_(cctkGH, Atxx);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxy_(cctkGH, Atxy);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atxz_(cctkGH, Atxz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atyy_(cctkGH, Atyy);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atyz_(cctkGH, Atyz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Atzz_(cctkGH, Atzz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Theta_(cctkGH, Theta);const GF3D<CCTK_REAL, 0, 0, 0> gf_alphaG_(cctkGH, alphaG);const GF3D<CCTK_REAL, 0, 0, 0> gf_betaGx_(cctkGH, betaGx);const GF3D<CCTK_REAL, 0, 0, 0> gf_betaGy_(cctkGH, betaGy);const GF3D<CCTK_REAL, 0, 0, 0> gf_betaGz_(cctkGH, betaGz);loop_all<0, 0, 0>(cctkGH, [&](const PointDesc &p) {// Loadconst mat3<CCTK_REAL> g(gf_gxx_, gf_gxy_, gf_gxz_, gf_gyy_, gf_gyz_,gf_gzz_, p);const mat3<CCTK_REAL> k(gf_kxx_, gf_kxy_, gf_kxz_, gf_kyy_, gf_kyz_,gf_kzz_, p);const CCTK_REAL alp = gf_alp_(p.I);const vec3<CCTK_REAL> beta(gf_betax_, gf_betay_, gf_betaz_, p);// Calculate Z4c variables (all except Gammat)const CCTK_REAL detg = g.det();const mat3<CCTK_REAL> gu = g.inv(detg);const CCTK_REAL chi = cbrt(detg);const CCTK_REAL chi1 = 1 / chi;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);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 Theta = 0;const CCTK_REAL Kh = K - 2 * Theta;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 CCTK_REAL alphaG = alp;vec3<CCTK_REAL> betaG;for (int a = 0; a < 3; ++a)betaG(a) = beta(a);// Storegf_chi_(p.I) = chi;gammat.store(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, gf_gammatyy_,gf_gammatyz_, gf_gammatzz_, p);gf_Kh_(p.I) = Kh;At.store(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_, gf_Atzz_, p);gf_Theta_(p.I) = Theta;gf_alphaG_(p.I) = alphaG;betaG.store(gf_betaGx_, gf_betaGy_, gf_betaGz_, p);});}} // namespace Z4c
#include "physics.hxx"#include "tensor.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cmath>namespace Z4c {using namespace Loop;using namespace std;extern "C" void Z4c_Initial2(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Z4c_Initial2;const vec3<CCTK_REAL> dx{CCTK_DELTA_SPACE(0),CCTK_DELTA_SPACE(1),CCTK_DELTA_SPACE(2),};const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxx_(cctkGH, gammatxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxy_(cctkGH, gammatxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxz_(cctkGH, gammatxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyy_(cctkGH, gammatyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyz_(cctkGH, gammatyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatzz_(cctkGH, gammatzz);const GF3D<CCTK_REAL, 0, 0, 0> gf_Gamtx_(cctkGH, Gamtx);const GF3D<CCTK_REAL, 0, 0, 0> gf_Gamty_(cctkGH, Gamty);const GF3D<CCTK_REAL, 0, 0, 0> gf_Gamtz_(cctkGH, Gamtz);loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {// Loadconst mat3<CCTK_REAL> gammat(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_,gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, p);// Calculate Z4c variables (only Gamt)const mat3<CCTK_REAL> gammatu = gammat.inv(1);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),};vec3<mat3<CCTK_REAL> > Gammatl;vec3<mat3<CCTK_REAL> > Gammat;vec3<CCTK_REAL> Gamt;calc_gamma(gammat, gammatu, dgammat, Gammatl, Gammat, Gamt);// StoreGamt.store(gf_Gamtx_, gf_Gamty_, gf_Gamtz_, p);});}} // namespace Z4c
# Main make.code.defn file for thorn Maxwell# Source files in this directorySRCS = \adm.cxx \boundaries.cxx \constraints.cxx \enforce.cxx \initial1.cxx \initial2.cxx \rhs.cxx# Subdirectories containing source filesSUBDIRS =
#ifndef PHYSICS_HXX#define PHYSICS_HXX#include "tensor.hxx"namespace Z4c {template <typename T>mat3<vec3<T> > calc_dgu(const mat3<T> &gu, const mat3<vec3<T> > &dg) {// g_xy = g_xa g_yb g^ab// g_xy,c = (g_xa g_yb g^ab),c// = g_xa,c g_yb g^ab + g_xa g_yb,c g^ab + g_xa g_yb g^ab,c// g_xa g_yb g^ab,c = g_xy,c - g_xa,c g_yb g^ab - g_xa g_yb,c g^ab// = g_xy,c - g_xy,c - g_xy,c// = - g_xy,c// g^ab,c = - g^ax g^by g_xy,cmat3<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;}template <typename T>mat3<vec3<T> > calc_dAu(const mat3<T> &gu, const mat3<vec3<T> > &dgu,const mat3<T> &A, const mat3<vec3<T> > &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,cmat3<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;}template <typename T>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) {// Gammal_abcfor (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_bcfor (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^afor (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;}}template <typename T>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;}} // namespace Z4c#endif // #ifndef PHYSICS_HXX
#warning "TODO"#include <iostream>#include "physics.hxx"#include "tensor.hxx"#include <loop.hxx>#include <cctk.h>#include <cctk_Arguments_Checked.h>#include <cctk_Parameters.h>#include <cmath>namespace Z4c {using namespace Loop;using namespace std;extern "C" void Z4c_RHS(CCTK_ARGUMENTS) {DECLARE_CCTK_ARGUMENTS_Z4c_RHS;DECLARE_CCTK_PARAMETERS;const vec3<CCTK_REAL> dx{CCTK_DELTA_SPACE(0),CCTK_DELTA_SPACE(1),CCTK_DELTA_SPACE(2),};const GF3D<const CCTK_REAL, 0, 0, 0> gf_chi_(cctkGH, chi);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxx_(cctkGH, gammatxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxy_(cctkGH, gammatxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatxz_(cctkGH, gammatxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyy_(cctkGH, gammatyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatyz_(cctkGH, gammatyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_gammatzz_(cctkGH, gammatzz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Kh_(cctkGH, Kh);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxx_(cctkGH, Atxx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxy_(cctkGH, Atxy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atxz_(cctkGH, Atxz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyy_(cctkGH, Atyy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atyz_(cctkGH, Atyz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Atzz_(cctkGH, Atzz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamtx_(cctkGH, Gamtx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamty_(cctkGH, Gamty);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Gamtz_(cctkGH, Gamtz);const GF3D<const CCTK_REAL, 0, 0, 0> gf_Theta_(cctkGH, Theta);const GF3D<const CCTK_REAL, 0, 0, 0> gf_alphaG_(cctkGH, alphaG);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaGx_(cctkGH, betaGx);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaGy_(cctkGH, betaGy);const GF3D<const CCTK_REAL, 0, 0, 0> gf_betaGz_(cctkGH, betaGz);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_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) {// Loadconst 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);const CCTK_REAL alphaG = gf_alphaG_(p.I);const vec3<CCTK_REAL> betaG(gf_betaGx_, gf_betaGy_, gf_betaGz_, p);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};// Calculate RHS// 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),};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);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 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),};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);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),};// (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;// (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;}// (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);// (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;}// Storegf_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);});}} // namespace Z4c
#ifndef TENSOR_HXX#define TENSOR_HXX#include <loop.hxx>#include <algorithm>#include <cmath>#include <initializer_list>#include <ostream>namespace Z4c {using namespace Loop;using namespace std;template <typename T> constexpr T pow2(const T x) { return x * x; }////////////////////////////////////////////////////////////////////////////////template <typename T> struct nan {constexpr T operator()() const { return 0.0 / 0.0; }};template <typename T, int D> struct nan<vect<T, D> > {constexpr vect<T, D> operator()() const {return vect<T, D>::pure(nan<T>()());}};////////////////////////////////////////////////////////////////////////////////// 3-vectortemplate <typename T> class vec3 {vect<T, 3> elts;static constexpr int ind(const int n) {assert(n >= 0 && n < 3);return n;}public:explicit constexpr vec3() : elts{nan<vect<T, 3> >()()} {}constexpr vec3(initializer_list<T> v) : elts(v) {}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)} {}vec3(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): vec3{gf_vx_(p.I), gf_vy_(p.I), gf_vz_(p.I)} {}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 {const auto &v = *this;#ifdef CCTK_DEBUGassert(!CCTK_isnan(v(0)));assert(!CCTK_isnan(v(1)));assert(!CCTK_isnan(v(2)));#endifgf_vx_(p.I) = v(0);gf_vy_(p.I) = v(1);gf_vz_(p.I) = v(2);}const T &operator()(int i) const { return elts[ind(i)]; }T &operator()(int i) { return elts[ind(i)]; }friend ostream &operator<<(ostream &os, const vec3<T> &v) {return os << "[" << v(0) << "," << v(1) << "," << v(2) << "]";}};template <typename T> struct nan<vec3<T> > {constexpr vec3<T> operator()() const { return vec3<T>(); }};////////////////////////////////////////////////////////////////////////////////// Symmetric 3-matrixtemplate <typename T> class mat3 {vect<T, 6> elts;static constexpr int symind(const int i, const int j) {assert(i >= 0 && i <= j && j < 3);const int n = 3 * i - i * (i + 1) / 2 + j;assert(n >= 0 && n < 6);return n;}static constexpr int ind(const int i, const int j) {return symind(min(i, j), max(i, j));}static_assert(symind(0, 0) == 0, "");static_assert(symind(0, 1) == 1, "");static_assert(symind(0, 2) == 2, "");static_assert(symind(1, 1) == 3, "");static_assert(symind(1, 2) == 4, "");static_assert(symind(2, 2) == 5, "");static_assert(ind(1, 0) == ind(0, 1), "");static_assert(ind(2, 0) == ind(0, 2), "");static_assert(ind(2, 1) == ind(1, 2), "");public:explicit constexpr mat3() : elts{nan<vect<T, 6> >()()} {}constexpr mat3(initializer_list<T> A) : elts(A) {}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)} {}mat3(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): 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)} {}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 auto &A = *this;#ifdef CCTK_DEBUGassert(!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)));#endifgf_Axx_(p.I) = A(0, 0);gf_Axy_(p.I) = A(0, 1);gf_Axz_(p.I) = A(0, 2);gf_Ayy_(p.I) = A(1, 1);gf_Ayz_(p.I) = A(1, 2);gf_Azz_(p.I) = A(2, 2);}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)]; }constexpr T det() const {const auto &A = *this;return A(0, 0) * (A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1)) -A(1, 0) * (A(0, 1) * A(2, 2) - A(0, 2) * A(2, 1)) +A(2, 0) * (A(0, 1) * A(1, 2) - A(0, 2) * A(1, 1));}constexpr mat3 inv(const T detA) const {const auto &A = *this;const T detA1 = 1 / detA;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)),detA1 * (A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0))};}constexpr T trace() const {const auto &A = *this;return A(0, 0) + A(1, 1) + A(2, 2);}friend ostream &operator<<(ostream &os, const mat3<T> &A) {return os << "[[" << A(0, 0) << "," << A(0, 1) << "," << A(0, 2) << "],["<< A(1, 0) << "," << A(1, 1) << "," << A(1, 2) << "],[" << A(2, 0)<< "," << A(2, 1) << "," << A(2, 2) << "]]";}};template <typename T> struct nan<mat3<T> > {constexpr mat3<T> operator()() const { return mat3<T>(); }};template <typename T>constexpr mat3<T> mul(const mat3<T> &A, const mat3<T> &B) {mat3<T> C;// C[a,b] = A[a,c] B[c,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;}////////////////////////////////////////////////////////////////////////////////template <typename T>T deriv(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 + DI[dir]) - gf_(I - DI[dir])) / (2 * dx(dir));}template <typename T>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>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) {constexpr vect<vect<int, dim>, dim> DI{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};if (dir1 == dir2)return (gf_(I + DI[dir1]) - 2 * gf_(I) + gf_(I - DI[dir1])) /(dx(dir1) * dx(dir2));elsereturn (deriv(gf_, I + DI[dir2], dx, dir1) -deriv(gf_, I - DI[dir2], dx, dir1)) /(2 * dx(dir2));}template <typename T>mat3<T> deriv2(const GF3D<const T, 0, 0, 0> &gf_, const vect<int, dim> &I,const vec3<T> &dx) {return {deriv2(gf_, I, dx, 0, 0), deriv2(gf_, I, dx, 0, 1),deriv2(gf_, I, dx, 0, 2), deriv2(gf_, I, dx, 1, 1),deriv2(gf_, I, dx, 1, 2), deriv2(gf_, I, dx, 2, 2),};}} // namespace Z4c#endif // #ifndef TENSOR_HXX