Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <Exporter.hpp>  // The SpEC Exporter
       7             : #include <optional>
       8             : #include <memory>
       9             : 
      10             : #include "DataStructures/CachedTempBuffer.hpp"
      11             : #include "Evolution/NumericInitialData.hpp"
      12             : #include "Options/String.hpp"
      13             : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
      14             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      15             : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
      16             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      17             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      18             : #include "Utilities/Serialization/CharmPupable.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : #include "Utilities/TaggedTuple.hpp"
      21             : 
      22             : /// \cond
      23             : namespace PUP {
      24             : class er;
      25             : }  // namespace PUP
      26             : /// \endcond
      27             : 
      28             : namespace grmhd::AnalyticData {
      29             : 
      30             : /*!
      31             :  * \brief Hydro initial data generated by SpEC.
      32             :  *
      33             :  * This class loads numerical data written out by the SpEC initial data solver.
      34             :  * It uses the `spec::Exporter` linked in from SpEC to interpolate to arbitrary
      35             :  * grid points. The coordinates are assumed to be in SpEC's "grid" frame.
      36             :  * We interpolate the following quantities:
      37             :  *
      38             :  * - "g": spatial metric
      39             :  * - "K": (lower) extrinsic curvature
      40             :  * - "Lapse": lapse
      41             :  * - "Shift": (upper) shift
      42             :  * - "BaryonDensity": rest mass density
      43             :  * - "u_i": lower spatial four-velocity
      44             :  *
      45             :  * The remaining hydro quantities are computed from the interpolated data and
      46             :  * the equation of state.
      47             :  * The magnetic field is set to zero and the electron fraction is set to a
      48             :  * constant read from the input file.
      49             :  */
      50             : template <size_t ThermodynamicDim>
      51           1 : class SpecInitialData : public evolution::initial_data::InitialData,
      52             :                         public evolution::NumericInitialData,
      53             :                         public AnalyticDataBase {
      54             :  private:
      55           0 :   struct FromEos {};
      56             : 
      57             :  public:
      58           0 :   using equation_of_state_type =
      59             :       EquationsOfState::EquationOfState<true, ThermodynamicDim>;
      60             : 
      61             :   template <typename DataType>
      62           0 :   using tags = tmpl::append<
      63             :       tmpl::list<gr::Tags::SpatialMetric<DataType, 3>,
      64             :                  gr::Tags::ExtrinsicCurvature<DataType, 3>,
      65             :                  gr::Tags::Lapse<DataType>, gr::Tags::Shift<DataType, 3>>,
      66             :       hydro::grmhd_tags<DataType>>;
      67             : 
      68           0 :   static std::string name() {
      69             :     return "SpecInitialData" + std::to_string(ThermodynamicDim) + "dEos";
      70             :   }
      71             : 
      72           0 :   struct DataDirectory {
      73           0 :     using type = std::string;
      74           0 :     static constexpr Options::String help = {
      75             :         "Path to a directory of data produced by SpEC. The directory is "
      76             :         "expected to contain 'GrDomain.input' and 'Vars*.h5' files for all the "
      77             :         "subdomains in GrDomain.input."};
      78             :   };
      79             : 
      80           0 :   struct DensityCutoff {
      81           0 :     using type = double;
      82           0 :     static constexpr Options::String help =
      83             :         "Where the density is below this cutoff the fluid variables are set to "
      84             :         "vacuum (atmosphere density, atmosphere pressure, atmosphere energy "
      85             :         "density/temperature and velocity, unit Lorentz factor and enthalpy). "
      86             :         "During the evolution, atmosphere treatment will typically kick in and "
      87             :         "fix the value of the fluid variables in these regions. Therefore, "
      88             :         "it makes sense to set this density cutoff to the same value as the "
      89             :         "atmosphere density cutoff.";
      90             :   };
      91             : 
      92           0 :   struct AtmosphereDensity {
      93           0 :     using type = double;
      94           0 :     static constexpr Options::String help =
      95             :         "The value to set the density to in atmosphere. This must match what "
      96             :         "is done during the evolution.";
      97             :   };
      98             : 
      99           0 :   struct ElectronFraction {
     100           0 :     using type = Options::Auto<double, FromEos>;
     101           0 :     static constexpr Options::String help = {
     102             :         "Constant electron fraction.\n"
     103             :         "To calculate electron fraction from the Eos, set to 'FromEos'."};
     104             :   };
     105             : 
     106           0 :   using options = tmpl::list<
     107             :       DataDirectory,
     108             :       hydro::OptionTags::InitialDataEquationOfState<true, ThermodynamicDim>,
     109             :       DensityCutoff, AtmosphereDensity, ElectronFraction>;
     110             : 
     111           0 :   static constexpr Options::String help = {"Initial data generated by SpEC"};
     112             : 
     113           0 :   SpecInitialData() = default;
     114           0 :   SpecInitialData(const SpecInitialData& rhs);
     115           0 :   SpecInitialData& operator=(const SpecInitialData& rhs);
     116           0 :   SpecInitialData(SpecInitialData&& /*rhs*/) = default;
     117           0 :   SpecInitialData& operator=(SpecInitialData&& /*rhs*/) = default;
     118           0 :   ~SpecInitialData() override = default;
     119             : 
     120           0 :   SpecInitialData(std::string data_directory,
     121             :                   std::unique_ptr<equation_of_state_type> equation_of_state,
     122             :                   double density_cutoff, double atmosphere_density,
     123             :                   std::optional<double> electron_fraction);
     124             : 
     125           0 :   auto get_clone() const
     126             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
     127             : 
     128             :   /// \cond
     129             :   explicit SpecInitialData(CkMigrateMessage* msg);
     130             :   using PUP::able::register_constructor;
     131             :   WRAPPED_PUPable_decl_template(SpecInitialData);
     132             :   /// \endcond
     133             : 
     134           0 :   const equation_of_state_type& equation_of_state() const {
     135             :     return *equation_of_state_;
     136             :   }
     137             : 
     138             :   template <typename DataType, typename... Tags>
     139           0 :   tuples::TaggedTuple<Tags...> variables(const tnsr::I<DataType, 3>& x,
     140             :                                          tmpl::list<Tags...> /*meta*/) const {
     141             :     auto interpolated_vars = interpolate_from_spec(x);
     142             :     using requested_tags = tmpl::list<Tags...>;
     143             :     using requested_derived_tags =
     144             :         tmpl::list_difference<requested_tags, interpolated_tags<DataType>>;
     145             :     using requested_interpolated_tags =
     146             :         tmpl::list_difference<requested_tags, requested_derived_tags>;
     147             :     tuples::TaggedTuple<Tags...> result{};
     148             :     // First, compute derived quantities from interpolated data
     149             :     if constexpr (tmpl::size<requested_derived_tags>::value > 0) {
     150             :       VariablesCache<DataType> cache{get_size(get<0>(x))};
     151             :       const VariablesComputer<DataType> computer{
     152             :           interpolated_vars, *equation_of_state_, density_cutoff_,
     153             :           atmosphere_density_, electron_fraction_};
     154             :       tmpl::for_each<requested_derived_tags>(
     155             :           [&result, &cache, &computer](const auto tag_v) {
     156             :             using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
     157             :             get<tag>(result) = cache.get_var(computer, tag{});
     158             :           });
     159             :     }
     160             :     // Then, move interpolated data into result buffer
     161             :     tmpl::for_each<requested_interpolated_tags>(
     162             :         [&result, &interpolated_vars](const auto tag_v) {
     163             :           using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
     164             :           get<tag>(result) = std::move(get<tag>(interpolated_vars));
     165             :         });
     166             :     return result;
     167             :   }
     168             : 
     169             :   // NOLINTNEXTLINE(google-runtime-references)
     170           0 :   void pup(PUP::er& /*p*/) override;
     171             : 
     172             :  private:
     173             :   /// These quantities are supported for interpolation from SpEC
     174             :   template <typename DataType>
     175           1 :   using interpolated_tags = tmpl::list<
     176             :       // GR quantities
     177             :       gr::Tags::SpatialMetric<DataType, 3>,
     178             :       gr::Tags::ExtrinsicCurvature<DataType, 3>, gr::Tags::Lapse<DataType>,
     179             :       gr::Tags::Shift<DataType, 3>,
     180             :       // Hydro quantities
     181             :       hydro::Tags::RestMassDensity<DataType>,
     182             :       hydro::Tags::LowerSpatialFourVelocity<DataType, 3>>;
     183             : 
     184             :   /// These are the names in SpEC datasets corresponding to the quantities above
     185           1 :   static const inline std::vector<std::string> vars_to_interpolate_{
     186             :       // GR quantities
     187             :       "g", "K", "Lapse", "Shift",
     188             :       // Hydro quantities
     189             :       "BaryonDensity", "u_i"};
     190             : 
     191             :   template <typename DataType>
     192             :   tuples::tagged_tuple_from_typelist<interpolated_tags<DataType>>
     193           0 :   interpolate_from_spec(const tnsr::I<DataType, 3>& x) const;
     194             : 
     195             :   /// This cache computes all derived quantities from the interpolated
     196             :   /// quantities on demand
     197             :   template <typename DataType>
     198           1 :   using VariablesCache = cached_temp_buffer_from_typelist<tmpl::push_back<
     199             :       tmpl::list_difference<tags<DataType>, interpolated_tags<DataType>>,
     200             :       hydro::Tags::LorentzFactorTimesSpatialVelocity<DataType, 3>,
     201             :       gr::Tags::InverseSpatialMetric<DataType, 3>>>;
     202             : 
     203             :   /// This is the computer for the `VariablesCache`
     204             :   template <typename DataType>
     205           1 :   struct VariablesComputer {
     206           0 :     using Cache = VariablesCache<DataType>;
     207             : 
     208             :     const tuples::tagged_tuple_from_typelist<interpolated_tags<DataType>>
     209           0 :         interpolated_data;
     210           0 :     const equation_of_state_type& eos;
     211           0 :     const double density_cutoff;
     212           0 :     const double atmosphere_density;
     213           0 :     const std::optional<double>& electron_fraction_value;
     214             : 
     215           0 :     void operator()(
     216             :         gsl::not_null<Scalar<DataType>*> specific_internal_energy,
     217             :         gsl::not_null<Cache*> cache,
     218             :         hydro::Tags::SpecificInternalEnergy<DataType> /*meta*/) const;
     219           0 :     void operator()(gsl::not_null<Scalar<DataType>*> pressure,
     220             :                     gsl::not_null<Cache*> cache,
     221             :                     hydro::Tags::Pressure<DataType> /*meta*/) const;
     222           0 :     void operator()(gsl::not_null<Scalar<DataType>*> specific_enthalpy,
     223             :                     gsl::not_null<Cache*> cache,
     224             :                     hydro::Tags::SpecificEnthalpy<DataType> /*meta*/) const;
     225           0 :     void operator()(gsl::not_null<Scalar<DataType>*> temperature,
     226             :                     gsl::not_null<Cache*> cache,
     227             :                     hydro::Tags::Temperature<DataType> /*meta*/) const;
     228           0 :     void operator()(gsl::not_null<tnsr::II<DataType, 3>*> inv_spatial_metric,
     229             :                     gsl::not_null<Cache*> cache,
     230             :                     gr::Tags::InverseSpatialMetric<DataType, 3> /*meta*/) const;
     231           0 :     void operator()(
     232             :         gsl::not_null<tnsr::I<DataType, 3>*>
     233             :             lorentz_factor_times_spatial_velocity,
     234             :         gsl::not_null<Cache*> cache,
     235             :         hydro::Tags::LorentzFactorTimesSpatialVelocity<DataType, 3> /*meta*/)
     236             :         const;
     237           0 :     void operator()(gsl::not_null<tnsr::I<DataType, 3>*> spatial_velocity,
     238             :                     gsl::not_null<Cache*> cache,
     239             :                     hydro::Tags::SpatialVelocity<DataType, 3> /*meta*/) const;
     240           0 :     void operator()(gsl::not_null<Scalar<DataType>*> lorentz_factor,
     241             :                     gsl::not_null<Cache*> cache,
     242             :                     hydro::Tags::LorentzFactor<DataType> /*meta*/) const;
     243           0 :     void operator()(gsl::not_null<Scalar<DataType>*> electron_fraction,
     244             :                     gsl::not_null<Cache*> cache,
     245             :                     hydro::Tags::ElectronFraction<DataType> /*meta*/) const;
     246           0 :     void operator()(gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
     247             :                     gsl::not_null<Cache*> cache,
     248             :                     hydro::Tags::MagneticField<DataType, 3> /*meta*/) const;
     249           0 :     void operator()(
     250             :         gsl::not_null<Scalar<DataType>*> div_cleaning_field,
     251             :         gsl::not_null<Cache*> cache,
     252             :         hydro::Tags::DivergenceCleaningField<DataType> /*meta*/) const;
     253             :   };
     254             : 
     255           0 :   std::string data_directory_{};
     256           0 :   std::unique_ptr<equation_of_state_type> equation_of_state_{nullptr};
     257           0 :   double density_cutoff_ = std::numeric_limits<double>::signaling_NaN();
     258           0 :   double atmosphere_density_ = std::numeric_limits<double>::signaling_NaN();
     259           0 :   std::optional<double> electron_fraction_ = std::nullopt;
     260             : 
     261           0 :   std::unique_ptr<spec::Exporter> spec_exporter_{nullptr};
     262             : };
     263             : 
     264             : }  // namespace grmhd::AnalyticData