Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : #include "DataStructures/DataBox/Prefixes.hpp"
       9             : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
      10             : #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
      11             : #include "Elliptic/Protocols/FirstOrderSystem.hpp"
      12             : #include "Elliptic/Systems/Xcts/FluxesAndSources.hpp"
      13             : #include "Elliptic/Systems/Xcts/Geometry.hpp"
      14             : #include "Elliptic/Systems/Xcts/Tags.hpp"
      15             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      16             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      17             : #include "PointwiseFunctions/GeneralRelativity/Tags/Conformal.hpp"
      18             : #include "Utilities/ProtocolHelpers.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : 
      21             : namespace Xcts {
      22             : 
      23             : /*!
      24             :  * \brief The Extended Conformal Thin Sandwich (XCTS) decomposition of the
      25             :  * Einstein constraint equations, formulated as a set of coupled first-order
      26             :  * partial differential equations
      27             :  *
      28             :  * See \ref Xcts for details on the XCTS equations.
      29             :  * The system can be formulated in terms of these fluxes and sources (see
      30             :  * `elliptic::protocols::FirstOrderSystem`):
      31             :  *
      32             :  * \f{align}
      33             :  * F^i_\psi ={} &\bar{\gamma}^{ij} \partial_j \psi \\
      34             :  * S_\psi ={} &-\bar{\Gamma}^i_{ij} F^j_\psi
      35             :  * + \frac{1}{8}\psi\bar{R} + \frac{1}{12}\psi^5 K^2
      36             :  * - \frac{1}{8}\psi^{-7}\bar{A}^2 - 2\pi\psi^5\rho
      37             :  * \f}
      38             :  *
      39             :  * for the Hamiltonian constraint,
      40             :  *
      41             :  * \f{align}
      42             :  * F^i_{\alpha\psi} ={} &\bar{\gamma}^{ij} \partial_j \alpha\psi \\
      43             :  * S_{\alpha\psi} ={} &-\bar{\Gamma}^i_{ij} F^j_{\alpha\psi}
      44             :  * + \alpha\psi \left(\frac{7}{8}\psi^{-8} \bar{A}^2
      45             :  * + \frac{5}{12} \psi^4 K^2 + \frac{1}{8}\bar{R}
      46             :  * + 2\pi\psi^4\left(\rho + 2S\right) \right) \\
      47             :  * &- \psi^5\partial_t K + \psi^5\left(\beta^i\bar{D}_i K
      48             :  * + \beta_\mathrm{background}^i\bar{D}_i K\right)
      49             :  * \f}
      50             :  *
      51             :  * for the lapse equation, and
      52             :  *
      53             :  * \f{align}
      54             :  * F^{ij}_\beta ={} &\left(\bar{L}\beta\right)^{ij} = \bar{\nabla}^i \beta^j
      55             :  * + \bar{\nabla}^j \beta^i
      56             :  * - \frac{2}{3} \bar{\gamma}^{ij} \bar{\nabla}_k \beta^k \\
      57             :  * S^i_\beta ={} &-\bar{\Gamma}^j_{jk} F^{ik}_\beta
      58             :  * - \bar{\Gamma}^i_{jk} F^{jk}_\beta
      59             :  * + \left(F^{ij}_\beta
      60             :  * + \left(\bar{L}\beta_\mathrm{background}\right)^{ij} - \bar{u}^{ij}\right)
      61             :  * \bar{\gamma}_{jk} \left(\frac{F^k_{\alpha\psi}}{\alpha\psi}
      62             :  * - 7 \frac{F^k_\psi}{\psi}\right) \\
      63             :  * &- \bar{D}_j\left(\left(\bar{L}\beta_\mathrm{background}\right)^{ij}
      64             :  * - \bar{u}^{ij}\right)
      65             :  * + \frac{4}{3}\frac{\alpha\psi}{\psi}\bar{D}^i K
      66             :  * + 16\pi\left(\alpha\psi\right)\psi^3 S^i
      67             :  * \f}
      68             :  *
      69             :  * for the momentum constraint, with
      70             :  *
      71             :  * \f{align}
      72             :  * \bar{A}^{ij} ={} &\frac{\psi^7}{2\alpha\psi}\left(
      73             :  * \left(\bar{L}\beta\right)^{ij} +
      74             :  * \left(\bar{L}\beta_\mathrm{background}\right)^{ij} - \bar{u}^{ij} \right)
      75             :  * \f}
      76             :  *
      77             :  * and all \f$f_\alpha=0\f$.
      78             :  *
      79             :  * Note that the symbol \f$\beta\f$ in the equations above means
      80             :  * \f$\beta_\mathrm{excess}\f$. The full shift is \f$\beta_\mathrm{excess} +
      81             :  * \beta_\mathrm{background}\f$. See `Xcts::Tags::ShiftBackground` and
      82             :  * `Xcts::Tags::ShiftExcess` for details on this split. Also note that the
      83             :  * background shift is degenerate with \f$\bar{u}\f$ so we treat the quantity
      84             :  * \f$\left(\bar{L}\beta_\mathrm{background}\right)^{ij} - \bar{u}^{ij}\f$ as a
      85             :  * single background field (see
      86             :  * `Xcts::Tags::LongitudinalShiftBackgroundMinusDtConformalMetric`). The
      87             :  * covariant divergence of this quantity w.r.t. the conformal metric is also a
      88             :  * background field.
      89             :  *
      90             :  * \par Solving a subset of equations:
      91             :  * This system allows you to select a subset of `Xcts::Equations` so you don't
      92             :  * have to solve for all variables if some are analytically known. Specify the
      93             :  * set of enabled equations as the first template parameter. The set of required
      94             :  * background fields depends on your choice of equations.
      95             :  *
      96             :  * \par Conformal background geometry:
      97             :  * The equations simplify significantly if the conformal metric is flat
      98             :  * ("conformal flatness") and in Cartesian coordinates. In this case you can
      99             :  * specify `Xcts::Geometry::FlatCartesian` as the second template parameter so
     100             :  * computations are optimized for a flat background geometry and you don't have
     101             :  * to supply geometric background fields. Else, specify
     102             :  * `Xcts::Geometry::Curved`.
     103             :  *
     104             :  * \par Conformal matter scale:
     105             :  * The matter source terms in the XCTS equations have the known defect that they
     106             :  * can spoil uniqueness of the solutions. See e.g. \cite Baumgarte2006ug for a
     107             :  * detailed study. To cure this defect one can conformally re-scale the matter
     108             :  * source terms as \f$\bar{\rho}=\psi^n\rho\f$, \f$\bar{S}=\psi^n S\f$ and
     109             :  * \f$\bar{S^i}=\psi^n S^i\f$ and treat the re-scaled fields as
     110             :  * freely-specifyable background data for the XCTS equations. You can select the
     111             :  * `ConformalMatterScale` \f$n\f$ as the third template parameter. Common
     112             :  * choices are \f$n=0\f$ for vacuum systems where the matter sources are
     113             :  * irrelevant, \f$n=6\f$ as suggested in \cite Foucart2008qt or \f$n=8\f$ as
     114             :  * suggested in \cite Baumgarte2006ug.
     115             :  */
     116             : template <Equations EnabledEquations, Geometry ConformalGeometry,
     117             :           int ConformalMatterScale>
     118           1 : struct FirstOrderSystem
     119             :     : tt::ConformsTo<elliptic::protocols::FirstOrderSystem> {
     120           0 :   static constexpr Equations enabled_equations = EnabledEquations;
     121           0 :   static constexpr Geometry conformal_geometry = ConformalGeometry;
     122           0 :   static constexpr int conformal_matter_scale = ConformalMatterScale;
     123           0 :   static constexpr size_t volume_dim = 3;
     124             : 
     125           0 :   using primal_fields = tmpl::flatten<tmpl::list<
     126             :       Tags::ConformalFactorMinusOne<DataVector>,
     127             :       tmpl::conditional_t<
     128             :           EnabledEquations == Equations::HamiltonianAndLapse or
     129             :               EnabledEquations == Equations::HamiltonianLapseAndShift,
     130             :           Tags::LapseTimesConformalFactorMinusOne<DataVector>, tmpl::list<>>,
     131             :       tmpl::conditional_t<
     132             :           EnabledEquations == Equations::HamiltonianLapseAndShift,
     133             :           Tags::ShiftExcess<DataVector, 3, Frame::Inertial>, tmpl::list<>>>>;
     134             : 
     135             :   // As fluxes we use the gradients with raised indices for the Hamiltonian and
     136             :   // lapse equation, and the longitudinal shift excess for the momentum
     137             :   // constraint. The gradient fluxes don't have symmetries and no particular
     138             :   // meaning so we use the standard `Flux` tags, but for the symmetric
     139             :   // longitudinal shift we use the corresponding symmetric tag.
     140           0 :   using primal_fluxes = tmpl::flatten<tmpl::list<
     141             :       ::Tags::Flux<Tags::ConformalFactorMinusOne<DataVector>, tmpl::size_t<3>,
     142             :                    Frame::Inertial>,
     143             :       tmpl::conditional_t<
     144             :           EnabledEquations == Equations::HamiltonianAndLapse or
     145             :               EnabledEquations == Equations::HamiltonianLapseAndShift,
     146             :           ::Tags::Flux<Tags::LapseTimesConformalFactorMinusOne<DataVector>,
     147             :                        tmpl::size_t<3>, Frame::Inertial>,
     148             :           tmpl::list<>>,
     149             :       tmpl::conditional_t<
     150             :           EnabledEquations == Equations::HamiltonianLapseAndShift,
     151             :           Tags::LongitudinalShiftExcess<DataVector, 3, Frame::Inertial>,
     152             :           tmpl::list<>>>>;
     153             : 
     154           0 :   using background_fields = tmpl::flatten<tmpl::list<
     155             :       // Quantities for Hamiltonian constraint
     156             :       gr::Tags::Conformal<gr::Tags::EnergyDensity<DataVector>,
     157             :                           ConformalMatterScale>,
     158             :       gr::Tags::TraceExtrinsicCurvature<DataVector>,
     159             :       tmpl::conditional_t<
     160             :           ConformalGeometry == Geometry::Curved,
     161             :           tmpl::list<
     162             :               Tags::InverseConformalMetric<DataVector, 3, Frame::Inertial>,
     163             :               Tags::ConformalRicciScalar<DataVector>,
     164             :               Tags::ConformalChristoffelContracted<DataVector, 3,
     165             :                                                    Frame::Inertial>,
     166             :               ::Tags::deriv<
     167             :                   Tags::ConformalMetric<DataVector, 3, Frame::Inertial>,
     168             :                   tmpl::size_t<3>, Frame::Inertial>>,
     169             :           tmpl::list<>>,
     170             :       tmpl::conditional_t<
     171             :           EnabledEquations ==
     172             :               Equations::Hamiltonian,
     173             :           Tags::LongitudinalShiftMinusDtConformalMetricOverLapseSquare<
     174             :               DataVector>,
     175             :           tmpl::list<>>,
     176             :       // Additional quantities for lapse equation
     177             :       tmpl::conditional_t<
     178             :           EnabledEquations == Equations::HamiltonianAndLapse or
     179             :               EnabledEquations ==
     180             :                   Equations::HamiltonianLapseAndShift,
     181             :           tmpl::list<gr::Tags::Conformal<gr::Tags::StressTrace<DataVector>,
     182             :                                          ConformalMatterScale>,
     183             :                      ::Tags::dt<gr::Tags::TraceExtrinsicCurvature<DataVector>>>,
     184             :           tmpl::list<>>,
     185             :       tmpl::conditional_t<
     186             :           EnabledEquations ==
     187             :               Equations::HamiltonianAndLapse,
     188             :           tmpl::list<
     189             :               Tags::LongitudinalShiftMinusDtConformalMetricSquare<DataVector>,
     190             :               Tags::ShiftDotDerivExtrinsicCurvatureTrace<DataVector>>,
     191             :           tmpl::list<>>,
     192             :       // Additional quantities for momentum constraint
     193             :       tmpl::conditional_t<
     194             :           EnabledEquations ==
     195             :               Equations::HamiltonianLapseAndShift,
     196             :           tmpl::list<
     197             :               gr::Tags::Conformal<gr::Tags::MomentumDensity<DataVector, 3>,
     198             :                                   ConformalMatterScale>,
     199             :               ::Tags::deriv<gr::Tags::TraceExtrinsicCurvature<DataVector>,
     200             :                             tmpl::size_t<3>, Frame::Inertial>,
     201             :               Tags::ShiftBackground<DataVector, 3, Frame::Inertial>,
     202             :               Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<
     203             :                   DataVector, 3, Frame::Inertial>,
     204             :               // Note that this is the plain divergence, i.e. with no
     205             :               // Christoffel symbol terms added
     206             :               ::Tags::div<
     207             :                   Tags::LongitudinalShiftBackgroundMinusDtConformalMetric<
     208             :                       DataVector, 3, Frame::Inertial>>,
     209             :               tmpl::conditional_t<
     210             :                   ConformalGeometry == Geometry::Curved,
     211             :                   tmpl::list<
     212             :                       Tags::ConformalMetric<DataVector, 3, Frame::Inertial>,
     213             :                       Tags::ConformalChristoffelFirstKind<DataVector, 3,
     214             :                                                           Frame::Inertial>,
     215             :                       Tags::ConformalChristoffelSecondKind<DataVector, 3,
     216             :                                                            Frame::Inertial>>,
     217             :                   tmpl::list<>>>,
     218             :           tmpl::list<>>>>;
     219           0 :   using inv_metric_tag = tmpl::conditional_t<
     220             :       ConformalGeometry == Geometry::FlatCartesian, void,
     221             :       Tags::InverseConformalMetric<DataVector, 3, Frame::Inertial>>;
     222             : 
     223           0 :   using fluxes_computer = Fluxes<EnabledEquations, ConformalGeometry>;
     224           0 :   using sources_computer =
     225             :       Sources<EnabledEquations, ConformalGeometry, ConformalMatterScale>;
     226           0 :   using sources_computer_linearized =
     227             :       LinearizedSources<EnabledEquations, ConformalGeometry,
     228             :                         ConformalMatterScale>;
     229             : 
     230           0 :   using boundary_conditions_base =
     231             :       elliptic::BoundaryConditions::BoundaryCondition<3>;
     232             : };
     233             : 
     234             : }  // namespace Xcts