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