Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <limits>
       7             : #include <memory>
       8             : #include <pup.h>
       9             : 
      10             : #include "DataStructures/DataBox/Prefixes.hpp"  // IWYU pragma: keep
      11             : #include "DataStructures/Tensor/Tensor.hpp"     // IWYU pragma: keep
      12             : #include "Elliptic/Systems/Xcts/Tags.hpp"       // IWYU pragma: keep
      13             : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
      14             : #include "Options/Context.hpp"
      15             : #include "Options/String.hpp"
      16             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      17             : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
      18             : #include "Utilities/Serialization/CharmPupable.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : #include "Utilities/TaggedTuple.hpp"
      21             : 
      22             : // IWYU pragma: no_forward_declare Tensor
      23             : 
      24             : /// \cond
      25             : namespace PUP {
      26             : class er;
      27             : }  // namespace PUP
      28             : /// \endcond
      29             : 
      30             : namespace Xcts::Solutions {
      31             : 
      32             : /*!
      33             :  * \brief A constant density star in general relativity
      34             :  *
      35             :  * \details This solution describes a star with constant density \f$\rho_0\f$
      36             :  * that extends to a (conformal) radius \f$R\f$. It solves the XCTS Hamiltonian
      37             :  * constraint that reduces to the non-linear elliptic equation
      38             :  * \f[
      39             :  * \Delta^2\psi+2\pi\rho\psi^5=0
      40             :  * \f]
      41             :  * for the conformal factor \f$\psi\f$ (see `Xcts`) under the following
      42             :  * assumptions \cite Baumgarte2006ug :
      43             :  *
      44             :  * - Time-symmetry \f$K_{ij}=0\f$
      45             :  * - Conformal flatness \f$\overline{\gamma}=\delta\f$, so \f$\Delta\f$ is the
      46             :  * flat-space Laplacian
      47             :  * - Spherical symmetry
      48             :  *
      49             :  * Imposing boundary conditions
      50             :  * \f[
      51             :  * \frac{\partial\psi}{\partial r}=0 \quad \text{for} \quad r=0\\
      52             :  * \psi\rightarrow 1 \quad \text{for} \quad r\rightarrow\infty
      53             :  * \f]
      54             :  * and considering the energy density
      55             :  * \f[
      56             :  * \rho(r\leq R)=\rho_0 \quad \text{and} \quad \rho(r>R)=0
      57             :  * \f]
      58             :  * of the star the authors of \cite Baumgarte2006ug find the solution
      59             :  * \f[
      60             :  * \psi(r\leq R)=C u_\alpha(r) \quad \text{and}
      61             :  * \quad \psi(r>R)=\frac{\beta}{r} + 1
      62             :  * \f]
      63             :  * with \f$C=(2\pi\rho_0/3)^{-1/4}\f$, the Sobolev functions
      64             :  * \f[
      65             :  * u_\alpha(r)=\sqrt{\frac{\alpha R}{r^2+(\alpha R)^2}}
      66             :  * \f]
      67             :  * and real parameters \f$\alpha\f$ and \f$\beta\f$ that are determined by
      68             :  * the following relations:
      69             :  * \f[
      70             :  * \rho_0 R^2=\frac{3}{2\pi}f^2(\alpha) \quad \text{with}
      71             :  * \quad f(\alpha)=\frac{\alpha^5}{(1+\alpha^2)^3} \\
      72             :  * \frac{\beta}{R} + 1 = C u_\alpha(R)
      73             :  * \f]
      74             :  *
      75             :  * This solution is described in detail in \cite Baumgarte2006ug , and also in
      76             :  * Exercise 3.8 in \cite BaumgarteShapiro , since it
      77             :  * exhibits the non-uniqueness properties that are typical for the XCTS system.
      78             :  * In the simple case of the constant-density star the non-uniqueness is
      79             :  * apparent from the function \f$f(\alpha)\f$, which has two solutions for any
      80             :  * \f$\rho_0\f$ smaller than a critical density
      81             :  * \f[
      82             :  * \rho_\mathrm{crit}=\frac{3}{2\pi R^2}\frac{5^2}{6^6}
      83             :  * \approx\frac{0.0320}{R^2} \text{,}
      84             :  * \f]
      85             :  * a unique solution for \f$\rho_0=\rho_\mathrm{crit}\f$ and no solutions
      86             :  * above the critical density \cite Baumgarte2006ug . The authors identify the
      87             :  * \f$\alpha < \alpha_\mathrm{crit}=\sqrt{5}\f$ and \f$\alpha >
      88             :  * \alpha_\mathrm{crit}\f$ branches of solutions with the strong-field and
      89             :  * weak-field regimes, respectively (see \cite Baumgarte2006ug for details).
      90             :  * In this implementation we compute the weak-field solution by choosing the
      91             :  * \f$\alpha > \alpha_\mathrm{crit}\f$ that corresponds to the density
      92             :  * \f$\rho_0\f$ of the star. Therefore, we supply initial data
      93             :  * \f$\psi_\mathrm{init}=1\f$ so that a nonlinear iterative numerical solver
      94             :  * will converge to the same weak-field solution.
      95             :  */
      96           1 : class ConstantDensityStar : public elliptic::analytic_data::AnalyticSolution {
      97             :  public:
      98           0 :   struct Density {
      99           0 :     using type = double;
     100           0 :     static constexpr Options::String help{
     101             :         "The constant density within the star"};
     102           0 :     static double lower_bound() { return 0.; }
     103             :   };
     104           0 :   struct Radius {
     105           0 :     using type = double;
     106           0 :     static constexpr Options::String help{"The conformal radius of the star"};
     107           0 :     static double lower_bound() { return 0.; }
     108             :   };
     109             : 
     110           0 :   using options = tmpl::list<Density, Radius>;
     111           0 :   static constexpr Options::String help{
     112             :       "A constant density star in general relativity"};
     113             : 
     114           0 :   ConstantDensityStar() = default;
     115           0 :   ConstantDensityStar(const ConstantDensityStar&) = default;
     116           0 :   ConstantDensityStar& operator=(const ConstantDensityStar&) = default;
     117           0 :   ConstantDensityStar(ConstantDensityStar&&) = default;
     118           0 :   ConstantDensityStar& operator=(ConstantDensityStar&&) = default;
     119           0 :   ~ConstantDensityStar() = default;
     120           0 :   std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
     121             :       const override {
     122             :     return std::make_unique<ConstantDensityStar>(*this);
     123             :   }
     124             : 
     125             :   /// \cond
     126             :   explicit ConstantDensityStar(CkMigrateMessage* m)
     127             :       : elliptic::analytic_data::AnalyticSolution(m) {}
     128             :   using PUP::able::register_constructor;
     129             :   WRAPPED_PUPable_decl_template(ConstantDensityStar);
     130             :   /// \endcond
     131             : 
     132           0 :   ConstantDensityStar(double density, double radius,
     133             :                       const Options::Context& context = {});
     134             : 
     135           0 :   double density() const { return density_; }
     136           0 :   double radius() const { return radius_; }
     137             : 
     138             :   /// @{
     139             :   /// Retrieve variable at coordinates `x`
     140             :   template <typename DataType>
     141           1 :   auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
     142             :                  tmpl::list<Xcts::Tags::ConformalFactor<DataType>> /*meta*/)
     143             :       const -> tuples::TaggedTuple<Xcts::Tags::ConformalFactor<DataType>>;
     144             : 
     145             :   template <typename DataType>
     146           1 :   auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
     147             :                  tmpl::list<::Tags::Initial<
     148             :                      Xcts::Tags::ConformalFactor<DataType>>> /*meta*/) const
     149             :       -> tuples::TaggedTuple<
     150             :           ::Tags::Initial<Xcts::Tags::ConformalFactor<DataType>>>;
     151             : 
     152             :   template <typename DataType>
     153           1 :   auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
     154             :                  tmpl::list<::Tags::Initial<
     155             :                      ::Tags::deriv<Xcts::Tags::ConformalFactor<DataType>,
     156             :                                    tmpl::size_t<3>, Frame::Inertial>>> /*meta*/)
     157             :       const -> tuples::TaggedTuple<
     158             :           ::Tags::Initial<::Tags::deriv<Xcts::Tags::ConformalFactor<DataType>,
     159             :                                         tmpl::size_t<3>, Frame::Inertial>>>;
     160             : 
     161             :   template <typename DataType>
     162           1 :   auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
     163             :                  tmpl::list<::Tags::FixedSource<
     164             :                      Xcts::Tags::ConformalFactor<DataType>>> /*meta*/) const
     165             :       -> tuples::TaggedTuple<
     166             :           ::Tags::FixedSource<Xcts::Tags::ConformalFactor<DataType>>>;
     167             : 
     168             :   template <typename DataType>
     169           1 :   auto variables(const tnsr::I<DataType, 3, Frame::Inertial>& x,
     170             :                  tmpl::list<gr::Tags::EnergyDensity<DataType>> /*meta*/) const
     171             :       -> tuples::TaggedTuple<gr::Tags::EnergyDensity<DataType>>;
     172             :   /// @}
     173             : 
     174             :   /// Retrieve a collection of variables at coordinates `x`
     175             :   template <typename DataType, typename... Tags>
     176           1 :   tuples::TaggedTuple<Tags...> variables(
     177             :       const tnsr::I<DataType, 3, Frame::Inertial>& x,
     178             :       tmpl::list<Tags...> /*meta*/) const {
     179             :     static_assert(sizeof...(Tags) > 1,
     180             :                   "The generic template will recurse infinitely if only one "
     181             :                   "tag is being retrieved.");
     182             :     return {tuples::get<Tags>(variables(x, tmpl::list<Tags>{}))...};
     183             :   }
     184             : 
     185             :   // NOLINTNEXTLINE(google-runtime-references)
     186           0 :   void pup(PUP::er& p) override;
     187             : 
     188             :  private:
     189           0 :   double density_ = std::numeric_limits<double>::signaling_NaN();
     190           0 :   double radius_ = std::numeric_limits<double>::signaling_NaN();
     191           0 :   double alpha_ = std::numeric_limits<double>::signaling_NaN();
     192             : };
     193             : 
     194           0 : bool operator==(const ConstantDensityStar& /*lhs*/,
     195             :                 const ConstantDensityStar& /*rhs*/);
     196             : 
     197           0 : bool operator!=(const ConstantDensityStar& lhs, const ConstantDensityStar& rhs);
     198             : 
     199             : }  // namespace Xcts::Solutions