SpECTRE Documentation Coverage Report
Current view: top level - Evolution/Systems/GrMhd/ValenciaDivClean - FixConservatives.hpp Hit Total Coverage
Commit: 664546099c4dbf27a1b708fac45e39c82dd743d2 Lines: 12 65 18.5 %
Date: 2024-04-19 16:28:01
Legend: Lines: hit not hit

          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             : 
       8             : #include "DataStructures/Tensor/TypeAliases.hpp"
       9             : #include "Evolution/Systems/GrMhd/ValenciaDivClean/TagsDeclarations.hpp"  // IWYU pragma: keep
      10             : #include "Options/Context.hpp"
      11             : #include "Options/String.hpp"
      12             : #include "PointwiseFunctions/GeneralRelativity/TagsDeclarations.hpp"  // IWYU pragma: keep
      13             : #include "PointwiseFunctions/Hydro/MagneticFieldTreatment.hpp"
      14             : #include "Utilities/Gsl.hpp"
      15             : #include "Utilities/TMPL.hpp"
      16             : 
      17             : /// \cond
      18             : namespace PUP {
      19             : class er;
      20             : }  // namespace PUP
      21             : class DataVector;
      22             : /// \endcond
      23             : 
      24             : namespace grmhd {
      25             : namespace ValenciaDivClean {
      26             : 
      27             : /*!
      28             :  * \ingroup VariableFixingGroup
      29             :  * \brief Fix conservative variables using method developed by Foucart.
      30             :  *
      31             :  * Adjusts the conservative variables as follows:
      32             :  * - If the electron fraction \f$Y_e\f$ is below the value of the option
      33             :  *   `CutoffYe`, change \f$\tilde{Y}_e\f$ to \f$\tilde{D}\f$ times the value of
      34             :  *   the option `MinimumValueOfYe`.
      35             :  * - Changes \f${\tilde D}\f$, the generalized mass-energy density, such
      36             :  *   that \f$D\f$, the product of the rest mass density \f$\rho\f$ and the
      37             :  *   Lorentz factor \f$W\f$, is set to value of the option `MinimumValueOfD`,
      38             :  *   whenever \f$D\f$ is below the value of the option `CutoffD`.
      39             :  * - Increases \f${\tilde \tau}\f$, the generalized internal energy density,
      40             :  *   such that
      41             :  *   \f${\tilde B}^2 \leq 2 \sqrt{\gamma} (1 - \epsilon_B) {\tilde \tau}\f$,
      42             :  *   where \f${\tilde B}^i\f$ is the generalized magnetic field,
      43             :  *   \f$\gamma\f$ is the determinant of the spatial metric, and
      44             :  *   \f$\epsilon_B\f$ is the option `SafetyFactorForB`.
      45             :  * - Decreases \f${\tilde S}_i\f$, the generalized momentum density, such that
      46             :  *   \f${\tilde S}^2 \leq (1 - \epsilon_S) {\tilde S}^2_{max}\f$, where
      47             :  *   \f$\epsilon_S\f$ is the option `SafetyFactorForS`, and
      48             :  *   \f${\tilde S}^2_{max}\f$ is a complicated function of the conservative
      49             :  *   variables which can only be found through root finding. There are
      50             :  *   sufficient conditions for a set of conservative variables to satisfy the
      51             :  *   inequality, which can be used to avoid root finding at most points.
      52             :  *
      53             :  * \note The routine currently assumes the minimum specific enthalpy is one.
      54             :  *
      55             :  * For more details see Appendix B from the [thesis of Francois
      56             :  * Foucart](https://ecommons.cornell.edu/handle/1813/30652)
      57             :  *
      58             :  * You can plot the function whose root we are finding using:
      59             :  *
      60             :  * \code{.py}
      61             :  *
      62             :  * import numpy as np
      63             :  * import matplotlib.pyplot as plt
      64             :  *
      65             :  * upper_bound = 1.000119047987896748e+00
      66             :  * lower_bound = 1.000000000000000000e+00
      67             :  * s_tilde_squared = 5.513009056734747750e-30
      68             :  * d_tilde = 1.131468709980503465e-12
      69             :  * sqrt_det_g: 1.131468709980503418e+00
      70             :  * tau_tilde = 1.346990732914080573e-16
      71             :  * b_tilde_squared = 3.048155733848927391e-16
      72             :  * b_squared_over_d = 2.380959757934347320e-04
      73             :  * tau_over_d = 1.190479878968363843e-04
      74             :  * normalized_s_dot_b = -9.999999082462245337e-01
      75             :  *
      76             :  *
      77             :  * def function_of_w(lorentz_factor):
      78             :  *     return ((lorentz_factor + b_squared_over_d - tau_over_d - 1.0) *
      79             :  *             (lorentz_factor**2 + b_squared_over_d * normalized_s_dot_b**2 *
      80             :  *              (b_squared_over_d + 2.0 * lorentz_factor)) -
      81             :  *             0.5 * b_squared_over_d -
      82             :  *             0.5 * b_squared_over_d * normalized_s_dot_b**2 *
      83             :  *             (lorentz_factor**2 - 1.0 +
      84             :  *              2.0 * lorentz_factor * b_squared_over_d + b_squared_over_d**2))
      85             :  *
      86             :  *
      87             :  * lorentz_factor = np.linspace(lower_bound, upper_bound, num=10000)
      88             :  *
      89             :  * plt.plot(lorentz_factor, function_of_w(lorentz_factor))
      90             :  * plt.show()
      91             :  *
      92             :  * \endcode
      93             :  */
      94           1 : class FixConservatives {
      95             :  public:
      96             :   /// \brief Minimum value of rest-mass density times lorentz factor
      97           1 :   struct MinimumValueOfD {
      98           0 :     using type = double;
      99           0 :     static type lower_bound() { return 0.0; }
     100           0 :     static constexpr Options::String help = {
     101             :         "Minimum value of rest-mass density times lorentz factor"};
     102             :   };
     103             :   /// \brief Cutoff below which \f$D = \rho W\f$ is set to MinimumValueOfD
     104           1 :   struct CutoffD {
     105           0 :     using type = double;
     106           0 :     static type lower_bound() { return 0.0; }
     107           0 :     static constexpr Options::String help = {
     108             :         "Cutoff below which D is set to MinimumValueOfD"};
     109             :   };
     110             :   /// \brief Minimum value of electron fraction \f$Y_e\f$
     111           1 :   struct MinimumValueOfYe {
     112           0 :     using type = double;
     113           0 :     static type lower_bound() { return 0.0; }
     114           0 :     static constexpr Options::String help = {
     115             :         "Minimum value of electron fraction"};
     116             :   };
     117             :   /// \brief Cutoff below which \f$Y_e\f$ is set to MinimumValueOfYe
     118           1 :   struct CutoffYe {
     119           0 :     using type = double;
     120           0 :     static type lower_bound() { return 0.0; }
     121           0 :     static constexpr Options::String help = {
     122             :         "Cutoff below which Y_e is set to MinimumValueOfYe"};
     123             :   };
     124             :   /// \brief Safety factor \f$\epsilon_B\f$.
     125           1 :   struct SafetyFactorForB {
     126           0 :     using type = double;
     127           0 :     static type lower_bound() { return 0.0; }
     128           0 :     static constexpr Options::String help = {
     129             :         "Safety factor for magnetic field bound."};
     130             :   };
     131             :   /// \brief Safety factor \f$\epsilon_S\f$.
     132           1 :   struct SafetyFactorForS {
     133           0 :     using type = double;
     134           0 :     static type lower_bound() { return 0.0; }
     135           0 :     static constexpr Options::String help = {
     136             :         "Safety factor for momentum density bound above density cutoff."};
     137             :   };
     138             :   /// \brief Cutoff in \f$\rho_0 W\f$ below which we use a stricter safety
     139             :   /// factor for the magnitude of S.
     140           1 :   struct SafetyFactorForSCutoffD {
     141           0 :     using type = double;
     142           0 :     static type lower_bound() { return 0.0; }
     143           0 :     static constexpr Options::String help = {
     144             :         "Below this value of rest mass density time Lorentz factor, limit S "
     145             :         "more agressively."};
     146             :   };
     147             : 
     148             :   /// \brief Below SafetyFactorForSCutoffD, reduce \f$\epsilon_S\f$ by
     149             :   /// SafetyFactorForSSlope times
     150             :   /// \f$\log_{10}(\rho_0 W / SafetyFactorForSCutoffD)\f$
     151           1 :   struct SafetyFactorForSSlope {
     152           0 :     using type = double;
     153           0 :     static type lower_bound() { return 0.0; }
     154           0 :     static constexpr Options::String help = {
     155             :         "Slope of safety factor for momentum density bound below "
     156             :         "SafetyFactorForSCutoffD, express as a function of log10(rho*W)."};
     157             :   };
     158             :   /// Whether or not the limiting is enabled
     159           1 :   struct Enable {
     160           0 :     using type = bool;
     161           0 :     static constexpr Options::String help = {
     162             :         "If true then the limiting is applied."};
     163             :   };
     164             :   /// How to treat the magnetic field
     165           1 :   struct MagneticField {
     166           0 :     using type = hydro::MagneticFieldTreatment;
     167           0 :     static constexpr Options::String help = {
     168             :         "How to treat the magnetic field."};
     169             :   };
     170             : 
     171           0 :   using options =
     172             :       tmpl::list<MinimumValueOfD, CutoffD, MinimumValueOfYe, CutoffYe,
     173             :                  SafetyFactorForB, SafetyFactorForS, SafetyFactorForSCutoffD,
     174             :                  SafetyFactorForSSlope, Enable, MagneticField>;
     175           0 :   static constexpr Options::String help = {
     176             :       "Variable fixing used in Foucart's thesis.\n"};
     177             : 
     178           0 :   FixConservatives(double minimum_rest_mass_density_times_lorentz_factor,
     179             :                    double rest_mass_density_times_lorentz_factor_cutoff,
     180             :                    double minimum_electron_fraction,
     181             :                    double electron_fraction_cutoff,
     182             :                    double safety_factor_for_magnetic_field,
     183             :                    double safety_factor_for_momentum_density,
     184             :                    double safety_factor_for_momentum_density_cutoff_d,
     185             :                    double safety_factor_for_momentum_density_slope, bool enable,
     186             :                    hydro::MagneticFieldTreatment magnetic_field_treatment,
     187             :                    const Options::Context& context = {});
     188             : 
     189           0 :   FixConservatives() = default;
     190           0 :   FixConservatives(const FixConservatives& /*rhs*/) = default;
     191           0 :   FixConservatives& operator=(const FixConservatives& /*rhs*/) = default;
     192           0 :   FixConservatives(FixConservatives&& /*rhs*/) = default;
     193           0 :   FixConservatives& operator=(FixConservatives&& /*rhs*/) = default;
     194           0 :   ~FixConservatives() = default;
     195             : 
     196             :   // NOLINTNEXTLINE(google-runtime-references)
     197           0 :   void pup(PUP::er& p);
     198             : 
     199           0 :   using return_tags = tmpl::list<grmhd::ValenciaDivClean::Tags::TildeD,
     200             :                                  grmhd::ValenciaDivClean::Tags::TildeYe,
     201             :                                  grmhd::ValenciaDivClean::Tags::TildeTau,
     202             :                                  grmhd::ValenciaDivClean::Tags::TildeS<>>;
     203           0 :   using argument_tags =
     204             :       tmpl::list<grmhd::ValenciaDivClean::Tags::TildeB<>,
     205             :                  gr::Tags::SpatialMetric<DataVector, 3>,
     206             :                  gr::Tags::InverseSpatialMetric<DataVector, 3>,
     207             :                  gr::Tags::SqrtDetSpatialMetric<DataVector>>;
     208             : 
     209             :   /// Returns `true` if any variables were fixed.
     210           1 :   bool operator()(
     211             :       gsl::not_null<Scalar<DataVector>*> tilde_d,
     212             :       gsl::not_null<Scalar<DataVector>*> tilde_ye,
     213             :       gsl::not_null<Scalar<DataVector>*> tilde_tau,
     214             :       gsl::not_null<tnsr::i<DataVector, 3, Frame::Inertial>*> tilde_s,
     215             :       const tnsr::I<DataVector, 3, Frame::Inertial>& tilde_b,
     216             :       const tnsr::ii<DataVector, 3, Frame::Inertial>& spatial_metric,
     217             :       const tnsr::II<DataVector, 3, Frame::Inertial>& inv_spatial_metric,
     218             :       const Scalar<DataVector>& sqrt_det_spatial_metric) const;
     219             : 
     220             :  private:
     221           0 :   friend bool operator==(const FixConservatives& lhs,
     222             :                          const FixConservatives& rhs);
     223             : 
     224           0 :   double minimum_rest_mass_density_times_lorentz_factor_{
     225             :       std::numeric_limits<double>::signaling_NaN()};
     226           0 :   double rest_mass_density_times_lorentz_factor_cutoff_{
     227             :       std::numeric_limits<double>::signaling_NaN()};
     228           0 :   double minimum_electron_fraction_{
     229             :       std::numeric_limits<double>::signaling_NaN()};
     230           0 :   double electron_fraction_cutoff_{
     231             :       std::numeric_limits<double>::signaling_NaN()};
     232           0 :   double one_minus_safety_factor_for_magnetic_field_{
     233             :       std::numeric_limits<double>::signaling_NaN()};
     234           0 :   double one_minus_safety_factor_for_momentum_density_{
     235             :       std::numeric_limits<double>::signaling_NaN()};
     236           0 :   double safety_factor_for_momentum_density_cutoff_d_{
     237             :       std::numeric_limits<double>::signaling_NaN()};
     238           0 :   double safety_factor_for_momentum_density_slope_{
     239             :       std::numeric_limits<double>::signaling_NaN()};
     240           0 :   bool enable_{true};
     241           0 :   hydro::MagneticFieldTreatment magnetic_field_treatment_{
     242             :       hydro::MagneticFieldTreatment::AssumeNonZero};
     243             : };
     244             : 
     245           0 : bool operator!=(const FixConservatives& lhs, const FixConservatives& rhs);
     246             : 
     247             : }  // namespace ValenciaDivClean
     248             : }  // namespace grmhd

Generated by: LCOV version 1.14