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