Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <pup.h>
       8             : #include <string>
       9             : #include <vector>
      10             : 
      11             : #include "DataStructures/BoostMultiArray.hpp"
      12             : #include "DataStructures/Tensor/TypeAliases.hpp"
      13             : #include "NumericalAlgorithms/Interpolation/MultiLinearSpanInterpolation.hpp"
      14             : #include "Options/String.hpp"
      15             : #include "Utilities/Gsl.hpp"
      16             : #include "Utilities/Serialization/CharmPupable.hpp"
      17             : 
      18             : /// \cond
      19             : class DataVector;
      20             : /// \endcond
      21             : 
      22             : namespace Particles::MonteCarlo {
      23             : 
      24             : /// Class responsible for reading neutrino-matter interaction
      25             : /// tables.
      26             : template <size_t EnergyBins, size_t NeutrinoSpecies>
      27           1 : class NeutrinoInteractionTable : public PUP::able {
      28             :  public:
      29             :   /// Read table from disk and stores interaction rates.
      30           1 :   explicit NeutrinoInteractionTable(const std::string& filename);
      31             : 
      32           0 :   static constexpr Options::String help = {
      33             :       "Tabulated neutrino-matter interactions in NuLib format.\n"
      34             :       "Emissivity, absorption and scattering opacities are    \n"
      35             :       "tabulated as a function of density, eletron fraction   \n"
      36             :       "and temperature."};
      37             : 
      38           0 :   struct TableFilename {
      39           0 :     using type = std::string;
      40           0 :     static constexpr Options::String help{"File name of the NuLib table"};
      41             :   };
      42             : 
      43           0 :   using options = tmpl::list<TableFilename>;
      44             : 
      45             :   /// Explicit instantiation from table values, for tests
      46           1 :   NeutrinoInteractionTable(
      47             :       std::vector<double> table_data_,
      48             :       const std::array<double, EnergyBins>& table_neutrino_energies_,
      49             :       std::vector<double> table_log_density_,
      50             :       std::vector<double> table_log_temperature_,
      51             :       std::vector<double> table_electron_fraction_);
      52             : 
      53           0 :   explicit NeutrinoInteractionTable(CkMigrateMessage* msg) : PUP::able(msg) {}
      54             : 
      55             :   using PUP::able::register_constructor;
      56           0 :   void pup(PUP::er& p) override;
      57           0 :   WRAPPED_PUPable_decl_template(NeutrinoInteractionTable);
      58             : 
      59             :   /// Interpolate interaction rates to given values of density,
      60             :   /// temperature and electron fraction.
      61           1 :   void get_neutrino_matter_interactions(
      62             :       gsl::not_null<
      63             :           std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
      64             :           emissivity_in_cell,
      65             :       gsl::not_null<
      66             :           std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
      67             :           absorption_opacity,
      68             :       gsl::not_null<
      69             :           std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
      70             :           scattering_opacity,
      71             :       const Scalar<DataVector>& electron_fraction,
      72             :       const Scalar<DataVector>& rest_mass_density,
      73             :       const Scalar<DataVector>& temperature,
      74             :       const double& minimum_temperature) const;
      75             : 
      76           0 :   const std::array<double, EnergyBins>& get_neutrino_energies() const {
      77             :     return table_neutrino_energies;
      78             :   }
      79             : 
      80             :  private:
      81           0 :   void initialize_interpolator();
      82             : 
      83             :   // Stores emissivity, absorption_opacity, scattering_opacity
      84             :   // For each quantities, there are NeutrinoSpecies * EnergyBins
      85             :   // variables store as a function of log(density), log(temperature)
      86             :   // and electron fraction.
      87             :   // The indexing varies fastest in EnergyBins, then Species, then
      88             :   // log(density), then log(temperature), and finally Ye.
      89           0 :   std::vector<double> table_data{};
      90             :   // Central energy of each bin
      91           0 :   std::array<double, EnergyBins> table_neutrino_energies;
      92             :   // Table discretization
      93           0 :   std::vector<double> table_log_density{};
      94           0 :   std::vector<double> table_log_temperature{};
      95           0 :   std::vector<double> table_electron_fraction{};
      96             : 
      97             :   intrp::UniformMultiLinearSpanInterpolation<3,
      98             :                                              3 * EnergyBins * NeutrinoSpecies>
      99           0 :       interpolator_{};
     100             : 
     101           0 :   const double min_kappa = 1.e-70;
     102           0 :   const double max_kappa = 1.e70;
     103             : };
     104             : 
     105             : }  // namespace Particles::MonteCarlo