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 : #include <optional>
9 :
10 : #include "DataStructures/CachedTempBuffer.hpp"
11 : #include "DataStructures/TaggedTuple.hpp"
12 : #include "Evolution/NumericInitialData.hpp"
13 : #include "Options/String.hpp"
14 : #include "PointwiseFunctions/AnalyticData/GrMhd/AnalyticData.hpp"
15 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
16 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
17 : #include "PointwiseFunctions/Hydro/Tags.hpp"
18 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
19 : #include "Utilities/Serialization/CharmPupable.hpp"
20 : #include "Utilities/TMPL.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
|
