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