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 <memory>
       8             : #include <mutex>
       9             : 
      10             : #include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
      11             : #include "DataStructures/Tensor/Tensor.hpp"
      12             : #include "Evolution/NumericInitialData.hpp"
      13             : #include "IO/External/InterpolateFromFuka.hpp"
      14             : #include "Options/String.hpp"
      15             : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
      16             : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
      17             : #include "PointwiseFunctions/Hydro/SpecificEnthalpy.hpp"
      18             : #include "PointwiseFunctions/Hydro/Tags.hpp"
      19             : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
      20             : #include "Utilities/Serialization/CharmPupable.hpp"
      21             : #include "Utilities/TMPL.hpp"
      22             : #include "Utilities/TaggedTuple.hpp"
      23             : 
      24             : /// \cond
      25             : namespace PUP {
      26             : class er;
      27             : }  // namespace PUP
      28             : /// \endcond
      29             : 
      30             : namespace grmhd::AnalyticData {
      31             : 
      32             : /*!
      33             :  * \brief Hydro initial data generated by FUKA.
      34             :  *
      35             :  * This class loads numerical data written out by the FUKA initial data solver.
      36             :  *
      37             :  * We choose a constant electron fraction and zero temperature for now because
      38             :  * FUKA doesn't export these quantities. We'll have to improve this later, e.g.
      39             :  * by constructing an EOS consistent with the FUKA data.
      40             :  */
      41           1 : class FukaInitialData : public evolution::initial_data::InitialData,
      42             :                         public evolution::NumericInitialData,
      43             :                         public AnalyticDataBase {
      44             :  public:
      45           0 :   struct InfoFilename {
      46           0 :     using type = std::string;
      47           0 :     static constexpr Options::String help = {
      48             :         "Path to the '.info' file of data produced by FUKA."};
      49             :   };
      50             : 
      51           0 :   struct ElectronFraction {
      52           0 :     using type = double;
      53           0 :     static constexpr Options::String help = {"Constant electron fraction"};
      54             :   };
      55             : 
      56           0 :   using options = tmpl::list<InfoFilename, ElectronFraction>;
      57             : 
      58           0 :   static constexpr Options::String help = {"Initial data generated by FUKA"};
      59             : 
      60           0 :   FukaInitialData() = default;
      61           0 :   FukaInitialData(const FukaInitialData& rhs);
      62           0 :   FukaInitialData& operator=(const FukaInitialData& rhs);
      63           0 :   FukaInitialData(FukaInitialData&& rhs);
      64           0 :   FukaInitialData& operator=(FukaInitialData&& rhs);
      65           0 :   ~FukaInitialData() override = default;
      66             : 
      67           0 :   FukaInitialData(std::string info_filename, double electron_fraction);
      68             : 
      69           0 :   auto get_clone() const
      70             :       -> std::unique_ptr<evolution::initial_data::InitialData> override;
      71             : 
      72             :   /// \cond
      73             :   explicit FukaInitialData(CkMigrateMessage* msg);
      74             :   using PUP::able::register_constructor;
      75             :   WRAPPED_PUPable_decl_template(FukaInitialData);
      76             :   /// \endcond
      77             : 
      78             :   template <typename DataType>
      79           0 :   using tags = tmpl::append<
      80             :       tmpl::list<gr::Tags::SpatialMetric<DataType, 3>,
      81             :                  gr::Tags::ExtrinsicCurvature<DataType, 3>,
      82             :                  gr::Tags::Lapse<DataType>, gr::Tags::Shift<DataType, 3>>,
      83             :       hydro::grmhd_tags<DataType>>;
      84             : 
      85             :   template <typename... RequestedTags>
      86           0 :   tuples::TaggedTuple<RequestedTags...> variables(
      87             :       const tnsr::I<DataVector, 3>& x,
      88             :       tmpl::list<RequestedTags...> /*meta*/) const {
      89             :     auto interpolated_vars = io::interpolate_from_fuka<io::FukaIdType::Bns>(
      90             :         make_not_null(&fuka_lock_), info_filename_, x);
      91             :     tuples::TaggedTuple<RequestedTags...> result{};
      92             :     // Move interpolated data into result buffer
      93             :     tmpl::for_each<
      94             :         tmpl::list<gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>,
      95             :                    gr::Tags::SpatialMetric<DataVector, 3>,
      96             :                    gr::Tags::ExtrinsicCurvature<DataVector, 3>,
      97             :                    hydro::Tags::RestMassDensity<DataVector>,
      98             :                    hydro::Tags::SpecificInternalEnergy<DataVector>,
      99             :                    hydro::Tags::Pressure<DataVector>,
     100             :                    hydro::Tags::SpatialVelocity<DataVector, 3>>>(
     101             :         [&result, &interpolated_vars](const auto tag_v) {
     102             :           using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
     103             :           get<tag>(result) = std::move(get<tag>(interpolated_vars));
     104             :         });
     105             :     // Compute derived quantities from interpolated data
     106             :     const size_t num_points = x.begin()->size();
     107             :     const auto& rest_mass_density =
     108             :         get<hydro::Tags::RestMassDensity<DataVector>>(result);
     109             :     const auto& specific_internal_energy =
     110             :         get<hydro::Tags::SpecificInternalEnergy<DataVector>>(result);
     111             :     const auto& pressure = get<hydro::Tags::Pressure<DataVector>>(result);
     112             :     const auto& spatial_velocity =
     113             :         get<hydro::Tags::SpatialVelocity<DataVector, 3>>(result);
     114             :     const auto& spatial_metric =
     115             :         get<gr::Tags::SpatialMetric<DataVector, 3>>(result);
     116             :     // Compute enthalpy from internal energy and pressure
     117             :     auto& specific_enthalpy =
     118             :         get<hydro::Tags::SpecificEnthalpy<DataVector>>(result);
     119             :     get(specific_enthalpy) = DataVector(num_points);
     120             :     for (size_t i = 0; i < num_points; ++i) {
     121             :       const double local_rest_mass_density = get(rest_mass_density)[i];
     122             :       if (equal_within_roundoff(local_rest_mass_density, 0.)) {
     123             :         get(specific_enthalpy)[i] = 1.;
     124             :       } else {
     125             :         get(specific_enthalpy)[i] = get(hydro::relativistic_specific_enthalpy(
     126             :             Scalar<double>(local_rest_mass_density),
     127             :             Scalar<double>(get(specific_internal_energy)[i]),
     128             :             Scalar<double>(get(pressure)[i])));
     129             :       }
     130             :     }
     131             :     // Constant electron fraction specified by input file
     132             :     auto& electron_fraction =
     133             :         get<hydro::Tags::ElectronFraction<DataVector>>(result);
     134             :     get(electron_fraction) = DataVector(num_points, electron_fraction_);
     135             :     // Zero magnetic field and divergence cleaning field
     136             :     auto& magnetic_field =
     137             :         get<hydro::Tags::MagneticField<DataVector, 3>>(result);
     138             :     get<0>(magnetic_field) = DataVector(num_points, 0.);
     139             :     get<1>(magnetic_field) = DataVector(num_points, 0.);
     140             :     get<2>(magnetic_field) = DataVector(num_points, 0.);
     141             :     auto& div_cleaning_field =
     142             :         get<hydro::Tags::DivergenceCleaningField<DataVector>>(result);
     143             :     get(div_cleaning_field) = DataVector(num_points, 0.);
     144             :     // Compute Lorentz factor from spatial velocity
     145             :     auto& lorentz_factor = get<hydro::Tags::LorentzFactor<DataVector>>(result);
     146             :     get(lorentz_factor) =
     147             :         1. / sqrt(1. - get(dot_product(spatial_velocity, spatial_velocity,
     148             :                                        spatial_metric)));
     149             :     // Set temperature to zero for now
     150             :     auto& temperature = get<hydro::Tags::Temperature<DataVector>>(result);
     151             :     get(temperature) = DataVector(num_points, 0.);
     152             :     return result;
     153             :   }
     154             : 
     155             :   // NOLINTNEXTLINE(google-runtime-references)
     156           0 :   void pup(PUP::er& /*p*/) override;
     157             : 
     158             :  private:
     159           0 :   std::string info_filename_{};
     160           0 :   double electron_fraction_ = std::numeric_limits<double>::signaling_NaN();
     161             : 
     162             :   // This lock is used to ensure that only one thread at a time is calling the
     163             :   // FUKA interpolation routines. We make some assumptions here to guarantee
     164             :   // thread-safety:
     165             :   // - This analytic data class exists only once per node (in the global cache).
     166             :   //   This means we don't have to copy or PUP the lock or pass it around
     167             :   //   instances.
     168             :   // - This also allows the lock to be mutable, which is necessary for the
     169             :   //   const-ness of the `variables` function.
     170           0 :   mutable std::mutex fuka_lock_{};  // NOLINT(spectre-mutable)
     171             : };
     172             : 
     173             : }  // namespace grmhd::AnalyticData