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 <string>
8 : #include <variant>
9 :
10 : #include "DataStructures/DataBox/DataBox.hpp"
11 : #include "DataStructures/TaggedTuple.hpp"
12 : #include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
13 : #include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.hpp"
14 : #include "DataStructures/Tensor/Tensor.hpp"
15 : #include "Domain/Structure/ElementId.hpp"
16 : #include "Domain/Tags.hpp"
17 : #include "IO/Importers/Actions/ReadVolumeData.hpp"
18 : #include "IO/Importers/ElementDataReader.hpp"
19 : #include "IO/Importers/Tags.hpp"
20 : #include "Options/String.hpp"
21 : #include "Parallel/AlgorithmExecution.hpp"
22 : #include "Parallel/GlobalCache.hpp"
23 : #include "Parallel/Invoke.hpp"
24 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
25 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
26 : #include "PointwiseFunctions/Hydro/Tags.hpp"
27 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
28 : #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
29 : #include "Utilities/ErrorHandling/Error.hpp"
30 : #include "Utilities/Gsl.hpp"
31 : #include "Utilities/SetNumberOfGridPoints.hpp"
32 : #include "Utilities/TMPL.hpp"
33 :
34 1 : namespace grmhd::ValenciaDivClean {
35 :
36 : /*!
37 : * \brief Numeric initial data loaded from volume data files
38 : *
39 : * This class can be factory-created in the input file to start an evolution
40 : * from numeric initial data. It selects the hydro variables to load from the
41 : * volume data files and allows to choose constant values for some of them.
42 : *
43 : * Where the density is below the `DensityCutoff` the fluid variables are set to
44 : * vacuum (zero density, pressure, energy and velocity, and unit Lorentz
45 : * factor). To evolve the initial data, an atmosphere treatment is likely
46 : * required to fix the value of the fluid variables in these regions.
47 : */
48 1 : class NumericInitialData : public evolution::initial_data::InitialData {
49 : public:
50 : /// Name of a variable in the volume data file. Can be optional, in which case
51 : /// a constant value can be supplied instead of a dataset name.
52 : template <typename Tag, typename IsRequired>
53 1 : struct VarName {
54 0 : using tag = Tag;
55 0 : static constexpr bool is_required = IsRequired::value;
56 0 : static std::string name() { return db::tag_name<Tag>(); }
57 0 : using type = std::conditional_t<is_required, std::string,
58 : std::variant<double, std::string>>;
59 0 : static constexpr Options::String help =
60 : "Name of the variable in the volume data file. For optional variables "
61 : "you may instead specify a double that is used as a constant value "
62 : "on the entire grid.";
63 : };
64 :
65 : // These are the hydro variables that we support loading from volume
66 : // data files
67 0 : using required_primitive_vars =
68 : tmpl::list<hydro::Tags::RestMassDensity<DataVector>,
69 : hydro::Tags::LowerSpatialFourVelocity<DataVector, 3>>;
70 0 : using optional_primitive_vars =
71 : tmpl::list<hydro::Tags::Temperature<DataVector>,
72 : hydro::Tags::ElectronFraction<DataVector>,
73 : hydro::Tags::MagneticField<DataVector, 3>>;
74 0 : using primitive_vars_option_tags =
75 : tmpl::append<db::wrap_tags_in<VarName, required_primitive_vars,
76 : std::bool_constant<true>>,
77 : db::wrap_tags_in<VarName, optional_primitive_vars,
78 : std::bool_constant<false>>>;
79 0 : struct PrimitiveVars
80 : : tuples::tagged_tuple_from_typelist<primitive_vars_option_tags> {
81 0 : static constexpr Options::String help =
82 : "Primitive hydro variables: 'RestMassDensity' and "
83 : "'LowerSpatialFourVelocity' (which is u_i = W * gamma_ij v^j). ";
84 0 : using options = tags_list;
85 : using TaggedTuple::TaggedTuple;
86 : };
87 :
88 0 : using all_vars =
89 : tmpl::append<required_primitive_vars, optional_primitive_vars>;
90 :
91 : // Input-file options
92 0 : struct Variables {
93 0 : using type = PrimitiveVars;
94 0 : static constexpr Options::String help =
95 : "Set of initial data variables from which the Valencia evolution "
96 : "variables are computed.";
97 : };
98 :
99 0 : struct DensityCutoff {
100 0 : using type = double;
101 0 : static constexpr Options::String help =
102 : "Where the density is below this cutoff the fluid variables are set to "
103 : "vacuum (zero density, pressure, energy and velocity, and unit Lorentz "
104 : "factor). "
105 : "During the evolution, atmosphere treatment will typically kick in and "
106 : "fix the value of the fluid variables in these regions. Therefore, "
107 : "it makes sense to set this density cutoff to the same value as the "
108 : "atmosphere density cutoff.";
109 0 : static constexpr double lower_bound() { return 0.; }
110 : };
111 :
112 0 : using options =
113 : tmpl::list<importers::OptionTags::VolumeData, Variables, DensityCutoff>;
114 :
115 0 : static constexpr Options::String help =
116 : "Numeric initial data loaded from volume data files";
117 :
118 0 : NumericInitialData() = default;
119 0 : NumericInitialData(const NumericInitialData& rhs) = default;
120 0 : NumericInitialData& operator=(const NumericInitialData& rhs) = default;
121 0 : NumericInitialData(NumericInitialData&& /*rhs*/) = default;
122 0 : NumericInitialData& operator=(NumericInitialData&& /*rhs*/) = default;
123 0 : ~NumericInitialData() = default;
124 :
125 : /// \cond
126 : explicit NumericInitialData(CkMigrateMessage* msg);
127 : using PUP::able::register_constructor;
128 : WRAPPED_PUPable_decl_template(NumericInitialData);
129 : /// \endcond
130 :
131 0 : std::unique_ptr<evolution::initial_data::InitialData> get_clone()
132 : const override {
133 : return std::make_unique<NumericInitialData>(*this);
134 : }
135 :
136 0 : NumericInitialData(importers::ImporterOptions importer_options,
137 : PrimitiveVars selected_variables, double density_cutoff);
138 :
139 0 : const importers::ImporterOptions& importer_options() const {
140 : return importer_options_;
141 : }
142 :
143 0 : const PrimitiveVars& selected_variables() const {
144 : return selected_variables_;
145 : }
146 :
147 0 : double density_cutoff() const { return density_cutoff_; }
148 :
149 0 : size_t volume_data_id() const;
150 :
151 : template <typename... AllTags>
152 0 : void select_for_import(
153 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> all_fields) const {
154 : // Select the subset of the available variables that we want to read from
155 : // the volume data file
156 : tmpl::for_each<primitive_vars_option_tags>([&all_fields,
157 : this](const auto option_tag_v) {
158 : using option_tag = tmpl::type_from<std::decay_t<decltype(option_tag_v)>>;
159 : using tag = typename option_tag::tag;
160 : static constexpr bool is_required = option_tag::is_required;
161 : const auto& selected_dataset_name = get<option_tag>(selected_variables_);
162 : if constexpr (is_required) {
163 : // Always select required tags for import
164 : get<importers::Tags::Selected<tag>>(*all_fields) =
165 : selected_dataset_name;
166 : } else {
167 : // Only select optional tags for import if a dataset name was
168 : // specified
169 : if (std::holds_alternative<std::string>(selected_dataset_name)) {
170 : get<importers::Tags::Selected<tag>>(*all_fields) =
171 : std::get<std::string>(selected_dataset_name);
172 : }
173 : }
174 : });
175 : }
176 :
177 : template <typename... AllTags>
178 0 : void set_initial_data(
179 : const gsl::not_null<Scalar<DataVector>*> rest_mass_density,
180 : const gsl::not_null<Scalar<DataVector>*> electron_fraction,
181 : const gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
182 : const gsl::not_null<tnsr::I<DataVector, 3>*> spatial_velocity,
183 : const gsl::not_null<tnsr::I<DataVector, 3>*> magnetic_field,
184 : const gsl::not_null<Scalar<DataVector>*> div_cleaning_field,
185 : const gsl::not_null<Scalar<DataVector>*> lorentz_factor,
186 : const gsl::not_null<Scalar<DataVector>*> pressure,
187 : const gsl::not_null<Scalar<DataVector>*> temperature,
188 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> numeric_data,
189 : const tnsr::II<DataVector, 3>& inv_spatial_metric,
190 : const EquationsOfState::EquationOfState<true, 3>& equation_of_state)
191 : const {
192 : // Rest mass density from dataset
193 : *rest_mass_density =
194 : std::move(get<hydro::Tags::RestMassDensity<DataVector>>(*numeric_data));
195 : const size_t num_points = get(*rest_mass_density).size();
196 : // Electron fraction from dataset or constant value
197 : if (const std::variant<double, std::string>& electron_fraction_selection =
198 : get<VarName<hydro::Tags::ElectronFraction<DataVector>,
199 : std::bool_constant<false>>>(selected_variables_);
200 : std::holds_alternative<std::string>(electron_fraction_selection)) {
201 : *electron_fraction = std::move(
202 : get<hydro::Tags::ElectronFraction<DataVector>>(*numeric_data));
203 : } else {
204 : const double constant_electron_fraction =
205 : std::get<double>(electron_fraction_selection);
206 : set_number_of_grid_points(electron_fraction, num_points);
207 : get(*electron_fraction) = constant_electron_fraction;
208 : }
209 : // Temperature from dataset or constant value
210 : if (const std::variant<double, std::string>& temperature_selection =
211 : get<VarName<hydro::Tags::Temperature<DataVector>,
212 : std::bool_constant<false>>>(selected_variables_);
213 : std::holds_alternative<std::string>(temperature_selection)) {
214 : *temperature =
215 : std::move(get<hydro::Tags::Temperature<DataVector>>(*numeric_data));
216 : } else {
217 : const double constant_temperature =
218 : std::get<double>(temperature_selection);
219 : set_number_of_grid_points(temperature, num_points);
220 : get(*temperature) = constant_temperature;
221 : }
222 : // Velocity and Lorentz factor from u_i dataset
223 : // W = 1 + W^2 v_i v^i
224 : // where W v_i = u_i, so we first raise the index on W v_i with the
225 : // spatial metric and then compute its magnitude. We use
226 : // `spatial_velocity` as intermediate memory buffer for W v^i.
227 : const auto& u_i = get<hydro::Tags::LowerSpatialFourVelocity<DataVector, 3>>(
228 : *numeric_data);
229 : raise_or_lower_index(spatial_velocity, u_i, inv_spatial_metric);
230 : dot_product(lorentz_factor, u_i, *spatial_velocity);
231 : get(*lorentz_factor) += 1.;
232 : for (size_t d = 0; d < 3; ++d) {
233 : spatial_velocity->get(d) /= get(*lorentz_factor);
234 : }
235 : // Specific internal energy, and pressure from EOS
236 : set_number_of_grid_points(specific_internal_energy, num_points);
237 : set_number_of_grid_points(pressure, num_points);
238 : for (size_t i = 0; i < num_points; ++i) {
239 : double& local_rest_mass_density = get(*rest_mass_density)[i];
240 : // Apply the EOS only where the density is above the cutoff, because the
241 : // fluid model breaks down in the zero-density limit
242 : if (local_rest_mass_density <= density_cutoff_) {
243 : local_rest_mass_density = 0.;
244 : get(*specific_internal_energy)[i] = 0.;
245 : get(*pressure)[i] = 0.;
246 : get(*temperature)[i] = 0.;
247 : // Also reset velocity and Lorentz factor below cutoff to be safe
248 : for (size_t d = 0; d < 3; ++d) {
249 : spatial_velocity->get(d)[i] = 0.;
250 : }
251 : get(*lorentz_factor)[i] = 1.;
252 : } else {
253 : get(*specific_internal_energy)[i] =
254 : get(equation_of_state
255 : .specific_internal_energy_from_density_and_temperature(
256 : Scalar<double>(local_rest_mass_density),
257 : Scalar<double>(get(*temperature)[i]),
258 : Scalar<double>(get(*electron_fraction)[i])));
259 : get(*pressure)[i] =
260 : get(equation_of_state.pressure_from_density_and_temperature(
261 : Scalar<double>(local_rest_mass_density),
262 : Scalar<double>(get(*temperature)[i]),
263 : Scalar<double>(get(*electron_fraction)[i])));
264 : }
265 : }
266 : // Magnetic field from dataset or constant value
267 : const std::variant<double, std::string>& magnetic_field_selection =
268 : get<VarName<hydro::Tags::MagneticField<DataVector, 3>,
269 : std::bool_constant<false>>>(selected_variables_);
270 : if (std::holds_alternative<std::string>(magnetic_field_selection)) {
271 : *magnetic_field = std::move(
272 : get<hydro::Tags::MagneticField<DataVector, 3>>(*numeric_data));
273 : } else {
274 : const double constant_magnetic_field =
275 : std::get<double>(magnetic_field_selection);
276 : if (constant_magnetic_field != 0.) {
277 : ERROR(
278 : "Choose a magnetic field dataset or set it to zero. "
279 : "Nonzero uniform magnetic fields cannot currently be chosen "
280 : "in the input file. Generate a dataset for the nonzero "
281 : "uniform magnetic field if you need to.");
282 : }
283 : set_number_of_grid_points(magnetic_field, num_points);
284 : std::fill(magnetic_field->begin(), magnetic_field->end(),
285 : constant_magnetic_field);
286 : }
287 : // Divergence cleaning field
288 : set_number_of_grid_points(div_cleaning_field, num_points);
289 : get(*div_cleaning_field) = 0.;
290 : }
291 :
292 0 : void pup(PUP::er& p) override;
293 :
294 0 : friend bool operator==(const NumericInitialData& lhs,
295 : const NumericInitialData& rhs);
296 :
297 : private:
298 0 : importers::ImporterOptions importer_options_{};
299 0 : PrimitiveVars selected_variables_{};
300 0 : double density_cutoff_{};
301 : };
302 :
303 0 : namespace Actions {
304 :
305 : /*!
306 : * \brief Dispatch loading numeric initial data from files.
307 : *
308 : * Place this action before
309 : * grmhd::ValenciaDivClean::Actions::SetNumericInitialData in the action list.
310 : * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
311 : * invoked by this action.
312 : */
313 1 : struct ReadNumericInitialData {
314 0 : using const_global_cache_tags =
315 : tmpl::list<evolution::initial_data::Tags::InitialData>;
316 :
317 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
318 : typename ArrayIndex, typename ActionList,
319 : typename ParallelComponent>
320 0 : static Parallel::iterable_action_return_t apply(
321 : db::DataBox<DbTagsList>& box,
322 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
323 : Parallel::GlobalCache<Metavariables>& cache,
324 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
325 : const ParallelComponent* const /*meta*/) {
326 : // Select the subset of the available variables that we want to read from
327 : // the volume data file
328 : const auto& initial_data = dynamic_cast<const NumericInitialData&>(
329 : db::get<evolution::initial_data::Tags::InitialData>(box));
330 : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
331 : importers::Tags::Selected, NumericInitialData::all_vars>>
332 : selected_fields{};
333 : initial_data.select_for_import(make_not_null(&selected_fields));
334 : // Dispatch loading the variables from the volume data file
335 : // - Not using `ckLocalBranch` here to make sure the simple action
336 : // invocation is asynchronous.
337 : auto& reader_component = Parallel::get_parallel_component<
338 : importers::ElementDataReader<Metavariables>>(cache);
339 : Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
340 : 3, NumericInitialData::all_vars, ParallelComponent>>(
341 : reader_component, initial_data.importer_options(),
342 : initial_data.volume_data_id(), std::move(selected_fields));
343 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
344 : }
345 : };
346 :
347 : /*!
348 : * \brief Receive numeric initial data loaded by
349 : * grmhd::ValenciaDivClean::Actions::ReadNumericInitialData.
350 : *
351 : * Place this action in the action list after
352 : * grmhd::ValenciaDivClean::Actions::ReadNumericInitialData to wait until the
353 : * data for this element has arrived, and then compute the remaining primitive
354 : * variables and store them in the DataBox to be used as initial data.
355 : *
356 : * This action modifies the tags listed in `hydro::grmhd_tags` in the DataBox
357 : * (i.e., the hydro primitives). It does not modify conservative variables, so
358 : * it relies on a primitive-to-conservative update in the action list before
359 : * the evolution can start.
360 : *
361 : * \requires This action requires that the (inverse) spatial metric is available
362 : * through the DataBox, so it should run after GR initial data has been loaded.
363 : *
364 : * \requires This action also requires an equation of state, which is retrieved
365 : * from the DataBox as `hydro::Tags::GrmhdEquationOfState`.
366 : */
367 1 : struct SetNumericInitialData {
368 0 : static constexpr size_t Dim = 3;
369 0 : using inbox_tags =
370 : tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
371 :
372 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
373 : typename ActionList, typename ParallelComponent>
374 0 : static Parallel::iterable_action_return_t apply(
375 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
376 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
377 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
378 : const ParallelComponent* const /*meta*/) {
379 : auto& inbox =
380 : tuples::get<importers::Tags::VolumeData<NumericInitialData::all_vars>>(
381 : inboxes);
382 : const auto& initial_data = dynamic_cast<const NumericInitialData&>(
383 : db::get<evolution::initial_data::Tags::InitialData>(box));
384 : const size_t volume_data_id = initial_data.volume_data_id();
385 : if (inbox.find(volume_data_id) == inbox.end()) {
386 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
387 : }
388 : auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
389 :
390 : const auto& inv_spatial_metric =
391 : db::get<gr::Tags::InverseSpatialMetric<DataVector, Dim>>(box);
392 : const auto& equation_of_state =
393 : db::get<hydro::Tags::GrmhdEquationOfState>(box);
394 :
395 : db::mutate<hydro::Tags::RestMassDensity<DataVector>,
396 : hydro::Tags::ElectronFraction<DataVector>,
397 : hydro::Tags::SpecificInternalEnergy<DataVector>,
398 : hydro::Tags::SpatialVelocity<DataVector, 3>,
399 : hydro::Tags::MagneticField<DataVector, 3>,
400 : hydro::Tags::DivergenceCleaningField<DataVector>,
401 : hydro::Tags::LorentzFactor<DataVector>,
402 : hydro::Tags::Pressure<DataVector>,
403 : hydro::Tags::Temperature<DataVector>>(
404 : [&initial_data, &numeric_data, &inv_spatial_metric, &equation_of_state](
405 : const gsl::not_null<Scalar<DataVector>*> rest_mass_density,
406 : const gsl::not_null<Scalar<DataVector>*> electron_fraction,
407 : const gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
408 : const gsl::not_null<tnsr::I<DataVector, 3>*> spatial_velocity,
409 : const gsl::not_null<tnsr::I<DataVector, 3>*> magnetic_field,
410 : const gsl::not_null<Scalar<DataVector>*> div_cleaning_field,
411 : const gsl::not_null<Scalar<DataVector>*> lorentz_factor,
412 : const gsl::not_null<Scalar<DataVector>*> pressure,
413 : const gsl::not_null<Scalar<DataVector>*> temperature) {
414 : initial_data.set_initial_data(
415 : rest_mass_density, electron_fraction, specific_internal_energy,
416 : spatial_velocity, magnetic_field, div_cleaning_field,
417 : lorentz_factor, pressure, temperature,
418 : make_not_null(&numeric_data), inv_spatial_metric,
419 : equation_of_state);
420 : },
421 : make_not_null(&box));
422 :
423 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
424 : }
425 : };
426 :
427 : } // namespace Actions
428 :
429 : } // namespace grmhd::ValenciaDivClean
|
