SpECTRE Documentation Coverage Report
Current view: top level - PointwiseFunctions/AnalyticSolutions/GeneralRelativity - Tov.cpp Hit Total Coverage
Commit: 2db722c93a8e9b106e406b439b79c8e05c2057fb Lines: 0 4 0.0 %
Date: 2021-03-03 22:01:00
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #include "PointwiseFunctions/AnalyticSolutions/GeneralRelativity/Tov.hpp"
       5             : 
       6             : // Need Boost MultiArray because it is used internally by ODEINT
       7             : #include "DataStructures/BoostMultiArray.hpp"  // IWYU pragma: keep
       8             : 
       9             : #include <algorithm>
      10             : #include <array>
      11             : #include <boost/numeric/odeint.hpp>  // IWYU pragma: keep
      12             : #include <cmath>
      13             : #include <cstddef>
      14             : #include <functional>
      15             : #include <ostream>
      16             : #include <pup.h>
      17             : #include <vector>
      18             : 
      19             : #include "DataStructures/DataVector.hpp"
      20             : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
      21             : #include "Utilities/ConstantExpressions.hpp"
      22             : #include "Utilities/ContainerHelpers.hpp"
      23             : #include "Utilities/ErrorHandling/Assert.hpp"
      24             : #include "Utilities/Gsl.hpp"
      25             : #include "Utilities/MakeWithValue.hpp"
      26             : 
      27             : // IWYU pragma: no_include <boost/numeric/odeint/integrate/integrate_adaptive.hpp>
      28             : // IWYU pragma: no_include <boost/numeric/odeint/stepper/controlled_runge_kutta.hpp>
      29             : // IWYU pragma: no_include <boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp>
      30             : // IWYU pragma: no_include <boost/numeric/odeint/stepper/generation/make_dense_output.hpp>
      31             : // IWYU pragma: no_include <boost/numeric/odeint/stepper/runge_kutta_dopri5.hpp>
      32             : // IWYU pragma: no_include <complex>
      33             : 
      34             : // IWYU pragma: no_forward_declare boost::numeric::odeint::controlled_runge_kutta
      35             : // IWYU pragma: no_forward_declare EquationsOfState::EquationOfState
      36             : // IWYU pragma: no_forward_declare Tensor
      37             : 
      38           0 : namespace {
      39             : 
      40             : void lindblom_rhs(const gsl::not_null<std::array<double, 2>*> dvars,
      41             :                   const std::array<double, 2>& vars, const double log_enthalpy,
      42             :                   const EquationsOfState::EquationOfState<true, 1>&
      43             :                       equation_of_state) noexcept {
      44             :   const double& radius_squared = vars[0];
      45             :   const double& mass_over_radius = vars[1];
      46             :   double& d_radius_squared = (*dvars)[0];
      47             :   double& d_mass_over_radius = (*dvars)[1];
      48             :   const double specific_enthalpy = std::exp(log_enthalpy);
      49             :   const double rest_mass_density =
      50             :       get(equation_of_state.rest_mass_density_from_enthalpy(
      51             :           Scalar<double>{specific_enthalpy}));
      52             :   const double pressure = get(equation_of_state.pressure_from_density(
      53             :       Scalar<double>{rest_mass_density}));
      54             :   const double energy_density =
      55             :       specific_enthalpy * rest_mass_density - pressure;
      56             : 
      57             :   // At the center of the star: (u,v) = (0,0)
      58             :   if (UNLIKELY((radius_squared == 0.0) and (mass_over_radius == 0.0))) {
      59             :     d_radius_squared = -3.0 / (2.0 * M_PI * (energy_density + 3.0 * pressure));
      60             :     d_mass_over_radius =
      61             :         -2.0 * energy_density / (energy_density + 3.0 * pressure);
      62             :   } else {
      63             :     const double common_factor =
      64             :         (1.0 - 2.0 * mass_over_radius) /
      65             :         (4.0 * M_PI * radius_squared * pressure + mass_over_radius);
      66             :     d_radius_squared = -2.0 * radius_squared * common_factor;
      67             :     d_mass_over_radius =
      68             :         -(4.0 * M_PI * radius_squared * energy_density - mass_over_radius) *
      69             :         common_factor;
      70             :   }
      71             : }
      72             : 
      73             : class Observer {
      74             :  public:
      75             :   void operator()(const std::array<double, 2>& vars,
      76             :                   const double current_log_enthalpy) noexcept {
      77             :     radius.push_back(std::sqrt(vars[0]));
      78             :     mass_over_radius.push_back(vars[1]);
      79             :     log_enthalpy.push_back(current_log_enthalpy);
      80             :   }
      81             :   std::vector<double> radius;
      82             :   std::vector<double> mass_over_radius;
      83             :   std::vector<double> log_enthalpy;
      84             : };
      85             : 
      86             : template <typename DataType>
      87             : RelativisticEuler::Solutions::TovStar<
      88             :     gr::Solutions::TovSolution>::RadialVariables<DataType>
      89             : interior_solution(
      90             :     const EquationsOfState::EquationOfState<true, 1>& equation_of_state,
      91             :     const DataType& radius, const DataType& mass_over_radius,
      92             :     const DataType& log_specific_enthalpy,
      93             :     const double log_lapse_at_outer_radius) noexcept {
      94             :   RelativisticEuler::Solutions::TovStar<
      95             :       gr::Solutions::TovSolution>::RadialVariables<DataType>
      96             :       result(radius);
      97             :   result.specific_enthalpy = Scalar<DataType>{exp(log_specific_enthalpy)};
      98             :   result.rest_mass_density = equation_of_state.rest_mass_density_from_enthalpy(
      99             :       result.specific_enthalpy);
     100             :   result.pressure =
     101             :       equation_of_state.pressure_from_density(result.rest_mass_density);
     102             :   result.specific_internal_energy =
     103             :       equation_of_state.specific_internal_energy_from_density(
     104             :           result.rest_mass_density);
     105             :   result.metric_time_potential =
     106             :       log_lapse_at_outer_radius - log_specific_enthalpy;
     107             :   result.dr_metric_time_potential =
     108             :       (mass_over_radius / radius + 4.0 * M_PI * get(result.pressure) * radius) /
     109             :       (1.0 - 2.0 * mass_over_radius);
     110             :   result.metric_radial_potential = -0.5 * log(1.0 - 2.0 * mass_over_radius);
     111             :   result.dr_metric_radial_potential =
     112             :       (4.0 * M_PI * radius *
     113             :            (get(result.specific_enthalpy) * get(result.rest_mass_density) -
     114             :             get(result.pressure)) -
     115             :        mass_over_radius / radius) /
     116             :       (1.0 - 2.0 * mass_over_radius);
     117             :   result.metric_angular_potential = make_with_value<DataType>(radius, 0.0);
     118             :   result.dr_metric_angular_potential = make_with_value<DataType>(radius, 0.0);
     119             :   return result;
     120             : }
     121             : 
     122             : template <typename DataType>
     123             : RelativisticEuler::Solutions::TovStar<
     124             :     gr::Solutions::TovSolution>::RadialVariables<DataType>
     125             : vacuum_solution(const DataType& radius, const double total_mass) noexcept {
     126             :   RelativisticEuler::Solutions::TovStar<
     127             :       gr::Solutions::TovSolution>::RadialVariables<DataType>
     128             :       result(radius);
     129             :   result.specific_enthalpy = make_with_value<Scalar<DataType>>(radius, 1.0);
     130             :   result.rest_mass_density = make_with_value<Scalar<DataType>>(radius, 0.0);
     131             :   result.pressure = make_with_value<Scalar<DataType>>(radius, 0.0);
     132             :   result.specific_internal_energy =
     133             :       make_with_value<Scalar<DataType>>(radius, 0.0);
     134             :   const DataType one_minus_two_m_over_r = 1.0 - 2.0 * total_mass / radius;
     135             :   result.metric_time_potential = 0.5 * log(one_minus_two_m_over_r);
     136             :   result.dr_metric_time_potential =
     137             :       total_mass / square(radius) / one_minus_two_m_over_r;
     138             :   result.metric_radial_potential = -result.metric_time_potential;
     139             :   result.dr_metric_radial_potential = -result.dr_metric_time_potential;
     140             :   result.metric_angular_potential = make_with_value<DataType>(radius, 0.0);
     141             :   result.dr_metric_angular_potential = make_with_value<DataType>(radius, 0.0);
     142             :   return result;
     143             : }
     144             : 
     145             : }  // namespace
     146             : 
     147             : namespace gr::Solutions {
     148             : 
     149             : TovSolution::TovSolution(
     150             :     const EquationsOfState::EquationOfState<true, 1>& equation_of_state,
     151             :     const double central_mass_density,
     152             :     const double log_enthalpy_at_outer_radius, const double absolute_tolerance,
     153             :     const double relative_tolerance) {
     154             :   std::array<double, 2> u_and_v = {{0.0, 0.0}};
     155             :   std::array<double, 2> dudh_and_dvdh{};
     156             :   const double central_log_enthalpy =
     157             :       std::log(get(equation_of_state.specific_enthalpy_from_density(
     158             :           Scalar<double>{central_mass_density})));
     159             :   lindblom_rhs(&dudh_and_dvdh, u_and_v, central_log_enthalpy,
     160             :                equation_of_state);
     161             :   const double initial_step = -std::min(std::abs(1.0 / dudh_and_dvdh[0]),
     162             :                                         std::abs(1.0 / dudh_and_dvdh[1]));
     163             :   using StateDopri5 =
     164             :       boost::numeric::odeint::runge_kutta_dopri5<std::array<double, 2>>;
     165             :   boost::numeric::odeint::dense_output_runge_kutta<
     166             :       boost::numeric::odeint::controlled_runge_kutta<StateDopri5>>
     167             :       dopri5 = make_dense_output(absolute_tolerance, relative_tolerance,
     168             :                                  StateDopri5{});
     169             :   Observer observer{};
     170             :   boost::numeric::odeint::integrate_adaptive(
     171             :       dopri5,
     172             :       [&equation_of_state](const std::array<double, 2>& lindblom_u_and_v,
     173             :                            std::array<double, 2>& lindblom_dudh_and_dvdh,
     174             :                            const double lindblom_enthalpy) noexcept {
     175             :         return lindblom_rhs(&lindblom_dudh_and_dvdh, lindblom_u_and_v,
     176             :                             lindblom_enthalpy, equation_of_state);
     177             :       },
     178             :       u_and_v, central_log_enthalpy, log_enthalpy_at_outer_radius, initial_step,
     179             :       std::ref(observer));
     180             :   outer_radius_ = observer.radius.back();
     181             :   const double total_mass_over_radius = observer.mass_over_radius.back();
     182             :   total_mass_ = total_mass_over_radius * outer_radius_;
     183             :   log_lapse_at_outer_radius_ = 0.5 * log(1.0 - 2.0 * total_mass_over_radius);
     184             :   mass_over_radius_interpolant_ =
     185             :       intrp::BarycentricRational(observer.radius, observer.mass_over_radius, 5);
     186             :   // log_enthalpy(radius) is almost linear so an interpolant of order 3
     187             :   // maximizes precision
     188             :   log_enthalpy_interpolant_ =
     189             :       intrp::BarycentricRational(observer.radius, observer.log_enthalpy, 3);
     190             : }
     191             : 
     192             : double TovSolution::outer_radius() const noexcept { return outer_radius_; }
     193             : 
     194             : double TovSolution::mass_over_radius(const double r) const noexcept {
     195             :   ASSERT(r >= 0.0 and r <= outer_radius_,
     196             :          "Invalid radius: " << r << " not in [0.0, " << outer_radius_ << "]\n");
     197             :   return mass_over_radius_interpolant_(r);
     198             : }
     199             : 
     200             : double TovSolution::mass(const double r) const noexcept {
     201             :   return mass_over_radius(r) * r;
     202             : }
     203             : 
     204             : double TovSolution::log_specific_enthalpy(const double r) const noexcept {
     205             :   ASSERT(r >= 0.0 and r <= outer_radius_,
     206             :          "Invalid radius: " << r << " not in [0.0, " << outer_radius_ << "]\n");
     207             :   return log_enthalpy_interpolant_(r);
     208             : }
     209             : 
     210             : template <>
     211             : RelativisticEuler::Solutions::TovStar<TovSolution>::RadialVariables<double>
     212           0 : TovSolution::radial_variables(
     213             :     const EquationsOfState::EquationOfState<true, 1>& equation_of_state,
     214             :     const tnsr::I<double, 3>& x) const noexcept {
     215             :   // add small number to avoid FPEs at origin
     216             :   const double radius = get(magnitude(x)) + 1.e-30 * outer_radius_;
     217             :   if (radius >= outer_radius_) {
     218             :     return vacuum_solution(radius, total_mass_);
     219             :   }
     220             :   return interior_solution(equation_of_state, radius, mass_over_radius(radius),
     221             :                            log_specific_enthalpy(radius),
     222             :                            log_lapse_at_outer_radius_);
     223             : }
     224             : 
     225             : template <>
     226             : RelativisticEuler::Solutions::TovStar<TovSolution>::RadialVariables<DataVector>
     227           0 : TovSolution::radial_variables(
     228             :     const EquationsOfState::EquationOfState<true, 1>& equation_of_state,
     229             :     const tnsr::I<DataVector, 3>& x) const noexcept {
     230             :   // add small number to avoid FPEs at origin
     231             :   const DataVector radius = get(magnitude(x)) + 1.e-30 * outer_radius_;
     232             :   if (min(radius) >= outer_radius_) {
     233             :     return vacuum_solution(radius, total_mass_);
     234             :   }
     235             :   if (max(radius) <= outer_radius_) {
     236             :     DataVector mass_over_radius_data(radius.size());
     237             :     DataVector log_of_specific_enthalpy(radius.size());
     238             :     for (size_t i = 0; i < get_size(radius); i++) {
     239             :       const double r = get_element(radius, i);
     240             :       get_element(mass_over_radius_data, i) = mass_over_radius(r);
     241             :       get_element(log_of_specific_enthalpy, i) = log_specific_enthalpy(r);
     242             :     }
     243             :     return interior_solution(equation_of_state, radius, mass_over_radius_data,
     244             :                              log_of_specific_enthalpy,
     245             :                              log_lapse_at_outer_radius_);
     246             :   }
     247             :   RelativisticEuler::Solutions::TovStar<TovSolution>::RadialVariables<
     248             :       DataVector>
     249             :       result(radius);
     250             :   for (size_t i = 0; i < radius.size(); i++) {
     251             :     const double r = radius[i];
     252             :     auto radial_vars_at_r =
     253             :         (r <= outer_radius_
     254             :              ? interior_solution(equation_of_state, r, mass_over_radius(r),
     255             :                                  log_specific_enthalpy(r),
     256             :                                  log_lapse_at_outer_radius_)
     257             :              : vacuum_solution(r, total_mass_));
     258             :     get(result.rest_mass_density)[i] = get(radial_vars_at_r.rest_mass_density);
     259             :     get(result.pressure)[i] = get(radial_vars_at_r.pressure);
     260             :     get(result.specific_internal_energy)[i] =
     261             :         get(radial_vars_at_r.specific_internal_energy);
     262             :     get(result.specific_enthalpy)[i] = get(radial_vars_at_r.specific_enthalpy);
     263             :     result.metric_time_potential[i] = radial_vars_at_r.metric_time_potential;
     264             :     result.dr_metric_time_potential[i] =
     265             :         radial_vars_at_r.dr_metric_time_potential;
     266             :     result.metric_radial_potential[i] =
     267             :         radial_vars_at_r.metric_radial_potential;
     268             :     result.dr_metric_radial_potential[i] =
     269             :         radial_vars_at_r.dr_metric_radial_potential;
     270             :   }
     271             :   result.metric_angular_potential = make_with_value<DataVector>(radius, 0.0);
     272             :   result.dr_metric_angular_potential = make_with_value<DataVector>(radius, 0.0);
     273             :   return result;
     274             : }
     275             : 
     276             : void TovSolution::pup(PUP::er& p) noexcept {  // NOLINT
     277             :   p | outer_radius_;
     278             :   p | total_mass_;
     279             :   p | log_lapse_at_outer_radius_;
     280             :   p | mass_over_radius_interpolant_;
     281             :   p | log_enthalpy_interpolant_;
     282             : }
     283             : 
     284             : }  // namespace gr::Solutions

Generated by: LCOV version 1.14