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/TaggedTuple.hpp"
11 : #include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
12 : #include "DataStructures/Tensor/Tensor.hpp"
13 : #include "Evolution/NumericInitialData.hpp"
14 : #include "IO/External/InterpolateFromFuka.hpp"
15 : #include "Options/String.hpp"
16 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
17 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
18 : #include "PointwiseFunctions/Hydro/SpecificEnthalpy.hpp"
19 : #include "PointwiseFunctions/Hydro/Tags.hpp"
20 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
21 : #include "Utilities/Serialization/CharmPupable.hpp"
22 : #include "Utilities/TMPL.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 : if constexpr (tmpl::list_contains_v<
117 : tmpl::list<RequestedTags...>,
118 : hydro::Tags::SpecificEnthalpy<DataVector>>) {
119 : // Compute enthalpy from internal energy and pressure if requested
120 : auto& specific_enthalpy =
121 : get<hydro::Tags::SpecificEnthalpy<DataVector>>(result);
122 : get(specific_enthalpy) = DataVector(num_points);
123 : for (size_t i = 0; i < num_points; ++i) {
124 : const double local_rest_mass_density = get(rest_mass_density)[i];
125 : if (equal_within_roundoff(local_rest_mass_density, 0.)) {
126 : get(specific_enthalpy)[i] = 1.;
127 : } else {
128 : get(specific_enthalpy)[i] = get(hydro::relativistic_specific_enthalpy(
129 : Scalar<double>(local_rest_mass_density),
130 : Scalar<double>(get(specific_internal_energy)[i]),
131 : Scalar<double>(get(pressure)[i])));
132 : }
133 : }
134 : }
135 : // Constant electron fraction specified by input file
136 : auto& electron_fraction =
137 : get<hydro::Tags::ElectronFraction<DataVector>>(result);
138 : get(electron_fraction) = DataVector(num_points, electron_fraction_);
139 : // Zero magnetic field and divergence cleaning field
140 : auto& magnetic_field =
141 : get<hydro::Tags::MagneticField<DataVector, 3>>(result);
142 : get<0>(magnetic_field) = DataVector(num_points, 0.);
143 : get<1>(magnetic_field) = DataVector(num_points, 0.);
144 : get<2>(magnetic_field) = DataVector(num_points, 0.);
145 : auto& div_cleaning_field =
146 : get<hydro::Tags::DivergenceCleaningField<DataVector>>(result);
147 : get(div_cleaning_field) = DataVector(num_points, 0.);
148 : // Compute Lorentz factor from spatial velocity
149 : auto& lorentz_factor = get<hydro::Tags::LorentzFactor<DataVector>>(result);
150 : get(lorentz_factor) =
151 : 1. / sqrt(1. - get(dot_product(spatial_velocity, spatial_velocity,
152 : spatial_metric)));
153 : // Set temperature to zero for now
154 : auto& temperature = get<hydro::Tags::Temperature<DataVector>>(result);
155 : get(temperature) = DataVector(num_points, 0.);
156 : return result;
157 : }
158 :
159 : // NOLINTNEXTLINE(google-runtime-references)
160 0 : void pup(PUP::er& /*p*/) override;
161 :
162 : private:
163 0 : std::string info_filename_{};
164 0 : double electron_fraction_ = std::numeric_limits<double>::signaling_NaN();
165 :
166 : // This lock is used to ensure that only one thread at a time is calling the
167 : // FUKA interpolation routines. We make some assumptions here to guarantee
168 : // thread-safety:
169 : // - This analytic data class exists only once per node (in the global cache).
170 : // This means we don't have to copy or PUP the lock or pass it around
171 : // instances.
172 : // - This also allows the lock to be mutable, which is necessary for the
173 : // const-ness of the `variables` function.
174 0 : mutable std::mutex fuka_lock_{}; // NOLINT(spectre-mutable)
175 : };
176 :
177 : } // namespace grmhd::AnalyticData
|
