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