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