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 "Utilities/Gsl.hpp"
      15             : #include "Utilities/Serialization/CharmPupable.hpp"
      16             : 
      17             : /// \cond
      18             : class DataVector;
      19             : /// \endcond
      20             : 
      21             : namespace Particles::MonteCarlo {
      22             : 
      23             : /// Class responsible for reading neutrino-matter interaction
      24             : /// tables.
      25             : template <size_t EnergyBins, size_t NeutrinoSpecies>
      26           1 : class NeutrinoInteractionTable : public PUP::able {
      27             :  public:
      28             :   /// Read table from disk and stores interaction rates.
      29           1 :   explicit NeutrinoInteractionTable(const std::string& filename);
      30             : 
      31             :   /// Explicit instantiation from table values, for tests
      32           1 :   NeutrinoInteractionTable(
      33             :       std::vector<double> table_data_,
      34             :       const std::array<double, EnergyBins>& table_neutrino_energies_,
      35             :       std::vector<double> table_log_density_,
      36             :       std::vector<double> table_log_temperature_,
      37             :       std::vector<double> table_electron_fraction_);
      38             : 
      39           0 :   explicit NeutrinoInteractionTable(CkMigrateMessage* msg) : PUP::able(msg) {}
      40             : 
      41             :   using PUP::able::register_constructor;
      42           0 :   void pup(PUP::er& p) override;
      43           0 :   WRAPPED_PUPable_decl_template(NeutrinoInteractionTable);
      44             : 
      45             :   /// Interpolate interaction rates to given values of density,
      46             :   /// temperature and electron fraction.
      47           1 :   void get_neutrino_matter_interactions(
      48             :       gsl::not_null<
      49             :           std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
      50             :           emissivity_in_cell,
      51             :       gsl::not_null<
      52             :           std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
      53             :           absorption_opacity,
      54             :       gsl::not_null<
      55             :           std::array<std::array<DataVector, EnergyBins>, NeutrinoSpecies>*>
      56             :           scattering_opacity,
      57             :       const Scalar<DataVector>& electron_fraction,
      58             :       const Scalar<DataVector>& rest_mass_density,
      59             :       const Scalar<DataVector>& temperature,
      60             :       const double& minimum_temperature) const;
      61             : 
      62           0 :   const std::array<double, EnergyBins>& get_neutrino_energies() const {
      63             :     return table_neutrino_energies;
      64             :   }
      65             : 
      66             :  private:
      67           0 :   void initialize_interpolator();
      68             : 
      69             :   // Stores emissivity, absorption_opacity, scattering_opacity
      70             :   // For each quantities, there are NeutrinoSpecies * EnergyBins
      71             :   // variables store as a function of log(density), log(temperature)
      72             :   // and electron fraction.
      73             :   // The indexing varies fastest in EnergyBins, then Species, then
      74             :   // log(density), then log(temperature), and finally Ye.
      75           0 :   std::vector<double> table_data{};
      76             :   // Central energy of each bin
      77           0 :   std::array<double, EnergyBins> table_neutrino_energies;
      78             :   // Table discretization
      79           0 :   std::vector<double> table_log_density{};
      80           0 :   std::vector<double> table_log_temperature{};
      81           0 :   std::vector<double> table_electron_fraction{};
      82             : 
      83             :   intrp::UniformMultiLinearSpanInterpolation<3,
      84             :                                              3 * EnergyBins * NeutrinoSpecies>
      85           0 :       interpolator_{};
      86             : 
      87           0 :   const double min_kappa = 1.e-70;
      88           0 :   const double max_kappa = 1.e70;
      89             : };
      90             : 
      91             : }  // namespace Particles::MonteCarlo