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             : #include <limits>
       8             : 
       9             : #include "DataStructures/Tensor/TypeAliases.hpp"
      10             : #include "Options/Options.hpp"
      11             : #include "Options/String.hpp"
      12             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      13             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      14             : #include "PointwiseFunctions/Hydro/TagsDeclarations.hpp"
      15             : #include "Utilities/Gsl.hpp"
      16             : #include "Utilities/TMPL.hpp"
      17             : 
      18             : /// \cond
      19             : class DataVector;
      20             : 
      21             : namespace PUP {
      22             : class er;
      23             : }  // namespace PUP
      24             : /// \endcond
      25             : 
      26             : namespace VariableFixing {
      27             : 
      28             : /*!
      29             :  * \ingroup VariableFixingGroup
      30             :  * \brief Adjust the electron fraction (Ye) based on rest mass density (rho).
      31             :  *
      32             :  * Based on \cite Liebendorfer2005 spherically symmetric Boltzmann calculations,
      33             :  * during the collapse phase just before a core-collapse supernova, the electron
      34             :  * fraction naturally follows the rest mass density.  Intuitively, the higher
      35             :  * the density, the more electrons are captured onto protons, leading to
      36             :  * neutrino production.  As these neutrinos diffuse out of the center of the
      37             :  * collapsing star, the total number of leptons near center of the CCSN
      38             :  * drops---the process of deleptonization.
      39             :  *
      40             :  * https://iopscience.iop.org/article/10.1086/466517
      41             :  */
      42           1 : class ParameterizedDeleptonization {
      43             :  public:
      44             :   /// \brief Use an analytic expression vs rest mass density profile
      45             :   ///
      46             :   /// The option to choose between analytic and tabulated will be added in the
      47             :   /// future.  If analytic is chosen (currently the only option), use the below
      48             :   /// parameters.  If tabulated is chosen, a Ye (rho) profile will need to be
      49             :   /// provided.
      50             :   //   struct DeleptonizationFormat {
      51             :   //     using type = double;
      52             :   //     static type lower_bound() { return 0.0; }
      53             :   //     static constexpr Options::String help = {"Choose an 'analytic' or
      54             :   //     'tabulated' expression to express electron fraction as a function of
      55             :   //     rest mass density."};
      56             :   //   };
      57             : 
      58             :   /// \brief Enable parameterized deleptonizations
      59           1 :   struct Enable {
      60           0 :     using type = bool;
      61           0 :     static constexpr Options::String help = {
      62             :         "Whether or not to activate parameterized deleptonization for "
      63             :         "supernovae."};
      64             :   };
      65             : 
      66             :   /// \brief Density near the center of the supernova at bounce, above which
      67             :   /// the central Ye is assumed
      68             :   ///
      69             :   /// In practice: rho(high) ~ 2x10^13 g / cm^3
      70           1 :   struct HighDensityScale {
      71           0 :     using type = double;
      72           0 :     static type lower_bound() { return 0.0; }
      73           0 :     static constexpr Options::String help = {
      74             :         "High end of density scale for parameterized deleptonization."};
      75             :   };
      76             :   /// \brief Density near the Silicon-Oxygen interface, below which the lower
      77             :   /// Ye is assumed
      78             :   ///
      79             :   /// In practice: rho(low) ~ 2x10^7 g / cm^3
      80           1 :   struct LowDensityScale {
      81           0 :     using type = double;
      82           0 :     static type lower_bound() { return 0.0; }
      83           0 :     static constexpr Options::String help = {
      84             :         "Low end of density scale for parameterized deleptonization."};
      85             :   };
      86             :   /// \brief Electron fraction of material when the rest mass density is above
      87             :   /// HighDensityScale
      88             :   ///
      89             :   /// In practice: Y_e (high density) ~ 0.28.
      90           1 :   struct ElectronFractionAtHighDensity {
      91           0 :     using type = double;
      92           0 :     static type lower_bound() { return 0.0; }
      93           0 :     static type upper_bound() { return 0.5; }
      94           0 :     static constexpr Options::String help = {
      95             :         "For densities above HighDensityScale, the electron fraction will "
      96             :         "take this value."};
      97             :   };
      98             :   /// \brief Electron fraction of material when the rest mass density is below
      99             :   /// LowDensityScale
     100             :   ///
     101             :   /// In practice: Y_e (low density) ~ 0.5
     102           1 :   struct ElectronFractionAtLowDensity {
     103           0 :     using type = double;
     104           0 :     static type lower_bound() { return 0.0; }
     105           0 :     static type upper_bound() { return 0.5; }
     106           0 :     static constexpr Options::String help = {
     107             :         "For densities below LowDensityScale, the electron fraction will "
     108             :         "take this value."};
     109             :   };
     110             : 
     111             :   /// \brief Electron fraction correction term.  The larger this value, the
     112             :   /// higher the Ye of matter at densities between LowDensityScale and
     113             :   /// HighDensityScale.
     114             :   ///
     115             :   /// In practice: Y_e (correction) ~ 0.035
     116           1 :   struct ElectronFractionCorrectionScale {
     117           0 :     using type = double;
     118           0 :     static type lower_bound() { return 0.0; }
     119           0 :     static constexpr Options::String help = {
     120             :         "For densities between low and high limits, a higher value of "
     121             :         "ElectronFractionCorrectionScale will increase the value of Ye."};
     122             :   };
     123             : 
     124           0 :   using options =
     125             :       tmpl::list<Enable, HighDensityScale, LowDensityScale,
     126             :                  ElectronFractionAtHighDensity, ElectronFractionAtLowDensity,
     127             :                  ElectronFractionCorrectionScale>;
     128           0 :   static constexpr Options::String help = {
     129             :       "Set electron fraction based on rest mass density.  "
     130             :       "(Low/High)DensityScale sets the limits of the density, beyond which "
     131             :       "the ElectronFractionAt(Low/High)Density is assumed.  At intermediate "
     132             :       "densities, the higher the ElectronFractionCorrectionScale, the higher "
     133             :       "the Ye."};
     134             : 
     135           0 :   ParameterizedDeleptonization(bool enable, double high_density_scale,
     136             :                                double low_density_scale,
     137             :                                double electron_fraction_at_high_density,
     138             :                                double electron_fraction_at_low_density,
     139             :                                double electron_fraction_correction_scale,
     140             :                                const Options::Context& context = {});
     141             : 
     142           0 :   ParameterizedDeleptonization() = default;
     143           0 :   ParameterizedDeleptonization(const ParameterizedDeleptonization& /*rhs*/) =
     144             :       default;
     145           0 :   ParameterizedDeleptonization& operator=(
     146             :       const ParameterizedDeleptonization& /*rhs*/) = default;
     147           0 :   ParameterizedDeleptonization(ParameterizedDeleptonization&& /*rhs*/) =
     148             :       default;
     149           0 :   ParameterizedDeleptonization& operator=(
     150             :       ParameterizedDeleptonization&& /*rhs*/) = default;
     151           0 :   ~ParameterizedDeleptonization() = default;
     152             : 
     153             :   // NOLINTNEXTLINE(google-runtime-references)
     154           0 :   void pup(PUP::er& p);
     155             : 
     156             :   // Quantities being mutated
     157           0 :   using return_tags =
     158             :       tmpl::list<hydro::Tags::SpecificInternalEnergy<DataVector>,
     159             :                  hydro::Tags::ElectronFraction<DataVector>,
     160             :                  hydro::Tags::Pressure<DataVector>,
     161             :                  hydro::Tags::Temperature<DataVector>>;
     162             : 
     163             :   // Things you want from DataBox that won't be change and are passed in as
     164             :   // const-refs
     165           0 :   using argument_tags = tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
     166             :                                    hydro::Tags::EquationOfStateBase>;
     167             : 
     168             :   // for use in `db::mutate_apply`
     169             :   template <size_t ThermodynamicDim>
     170           0 :   void operator()(
     171             :       gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
     172             :       gsl::not_null<Scalar<DataVector>*> electron_fraction,
     173             :       gsl::not_null<Scalar<DataVector>*> pressure,
     174             :       gsl::not_null<Scalar<DataVector>*> temperature,
     175             :       const Scalar<DataVector>& rest_mass_density,
     176             :       const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
     177             :           equation_of_state) const;
     178             : 
     179             :  private:
     180             :   template <size_t ThermodynamicDim>
     181           0 :   void correct_electron_fraction(
     182             :       gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
     183             :       gsl::not_null<Scalar<DataVector>*> electron_fraction,
     184             :       gsl::not_null<Scalar<DataVector>*> pressure,
     185             :       gsl::not_null<Scalar<DataVector>*> temperature,
     186             :       const Scalar<DataVector>& rest_mass_density,
     187             :       const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
     188             :           equation_of_state,
     189             :       size_t grid_index) const;
     190             : 
     191             :   // NOLINTNEXTLINE(readability-redundant-declaration)
     192           0 :   friend bool operator==(const ParameterizedDeleptonization& lhs,
     193             :                          const ParameterizedDeleptonization& rhs);
     194             : 
     195           0 :   bool enable_{false};
     196           0 :   double high_density_scale_{std::numeric_limits<double>::signaling_NaN()};
     197           0 :   double low_density_scale_{std::numeric_limits<double>::signaling_NaN()};
     198           0 :   double electron_fraction_at_high_density_{
     199             :       std::numeric_limits<double>::signaling_NaN()};
     200           0 :   double electron_fraction_at_low_density_{
     201             :       std::numeric_limits<double>::signaling_NaN()};
     202           0 :   double electron_fraction_half_sum_{
     203             :       std::numeric_limits<double>::signaling_NaN()};
     204           0 :   double electron_fraction_half_difference_{
     205             :       std::numeric_limits<double>::signaling_NaN()};
     206           0 :   double electron_fraction_correction_scale_{
     207             :       std::numeric_limits<double>::signaling_NaN()};
     208             : };
     209           0 : bool operator!=(const ParameterizedDeleptonization& lhs,
     210             :                 const ParameterizedDeleptonization& rhs);
     211             : 
     212             : }  // namespace VariableFixing