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 <optional>
8 : #include <string>
9 : #include <variant>
10 :
11 : #include "DataStructures/DataBox/DataBox.hpp"
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/TaggedTuple.hpp"
14 : #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "Domain/Structure/ElementId.hpp"
17 : #include "Evolution/DgSubcell/ActiveGrid.hpp"
18 : #include "Evolution/DgSubcell/Tags/ActiveGrid.hpp"
19 : #include "Evolution/DgSubcell/Tags/Coordinates.hpp"
20 : #include "Evolution/DgSubcell/Tags/Jacobians.hpp"
21 : #include "Evolution/DgSubcell/Tags/Mesh.hpp"
22 : #include "Evolution/NumericInitialData.hpp"
23 : #include "Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp"
24 : #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp"
25 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
26 : #include "Evolution/Systems/GrMhd/ValenciaDivClean/Actions/NumericInitialData.hpp"
27 : #include "IO/Importers/Actions/ReadVolumeData.hpp"
28 : #include "IO/Importers/ElementDataReader.hpp"
29 : #include "IO/Importers/Tags.hpp"
30 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
31 : #include "Options/String.hpp"
32 : #include "Parallel/AlgorithmExecution.hpp"
33 : #include "Parallel/GlobalCache.hpp"
34 : #include "Parallel/Invoke.hpp"
35 : #include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp"
36 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
37 : #include "PointwiseFunctions/Hydro/EquationsOfState/EquationOfState.hpp"
38 : #include "PointwiseFunctions/Hydro/Tags.hpp"
39 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
40 : #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
41 : #include "Utilities/ErrorHandling/Error.hpp"
42 : #include "Utilities/Gsl.hpp"
43 : #include "Utilities/TMPL.hpp"
44 :
45 : /// \cond
46 : namespace Tags {
47 : struct Time;
48 : } // namespace Tags
49 : /// \endcond
50 :
51 1 : namespace grmhd::GhValenciaDivClean {
52 :
53 : /*!
54 : * \brief Numeric initial data loaded from volume data files
55 : *
56 : * This class can be factory-created in the input file to start an evolution
57 : * from numeric initial data. It selects the set of GH variables to load from
58 : * the volume data file (ADM or GH variables), the set of hydro variables, and
59 : * allows to set constant values for some of the hydro variables. This class
60 : * mostly combines the `gh::NumericInitialData` and
61 : * `grmhd::ValenciaDivClean::NumericInitialData` classes.
62 : */
63 1 : class NumericInitialData : public evolution::initial_data::InitialData,
64 : public evolution::NumericInitialData {
65 : private:
66 0 : using GhNumericId = gh::NumericInitialData;
67 0 : using HydroNumericId = grmhd::ValenciaDivClean::NumericInitialData;
68 :
69 : public:
70 0 : using all_vars =
71 : tmpl::append<GhNumericId::all_vars, HydroNumericId::all_vars>;
72 :
73 0 : struct GhVariables : GhNumericId::Variables {};
74 0 : struct HydroVariables : HydroNumericId::Variables {};
75 :
76 0 : using options =
77 : tmpl::push_back<importers::ImporterOptions::tags_list, GhVariables,
78 : HydroVariables, HydroNumericId::DensityCutoff>;
79 :
80 0 : static constexpr Options::String help =
81 : "Numeric initial data loaded from volume data files";
82 :
83 0 : NumericInitialData() = default;
84 0 : NumericInitialData(const NumericInitialData& rhs) = default;
85 0 : NumericInitialData& operator=(const NumericInitialData& rhs) = default;
86 0 : NumericInitialData(NumericInitialData&& /*rhs*/) = default;
87 0 : NumericInitialData& operator=(NumericInitialData&& /*rhs*/) = default;
88 0 : ~NumericInitialData() = default;
89 :
90 : /// \cond
91 : explicit NumericInitialData(CkMigrateMessage* msg);
92 : using PUP::able::register_constructor;
93 : WRAPPED_PUPable_decl_template(NumericInitialData);
94 : /// \endcond
95 :
96 0 : std::unique_ptr<evolution::initial_data::InitialData> get_clone()
97 : const override {
98 : return std::make_unique<NumericInitialData>(*this);
99 : }
100 :
101 0 : NumericInitialData(
102 : std::string file_glob, std::string subfile_name,
103 : std::variant<double, importers::ObservationSelector> observation_value,
104 : std::optional<double> observation_value_epsilon,
105 : bool enable_interpolation,
106 : typename GhNumericId::Variables::type gh_selected_variables,
107 : typename HydroNumericId::Variables::type hydro_selected_variables,
108 : double density_cutoff);
109 :
110 0 : const importers::ImporterOptions& importer_options() const {
111 : return gh_numeric_id_.importer_options();
112 : }
113 :
114 0 : const GhNumericId& gh_numeric_id() const { return gh_numeric_id_; }
115 :
116 0 : const HydroNumericId& hydro_numeric_id() const { return hydro_numeric_id_; }
117 :
118 0 : size_t volume_data_id() const;
119 :
120 : template <typename... AllTags>
121 0 : void select_for_import(
122 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> fields) const {
123 : gh_numeric_id_.select_for_import(fields);
124 : hydro_numeric_id_.select_for_import(fields);
125 : }
126 :
127 : template <typename... AllTags, size_t ThermodynamicDim>
128 0 : void set_initial_data(
129 : const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
130 : const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
131 : const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
132 : const gsl::not_null<Scalar<DataVector>*> rest_mass_density,
133 : const gsl::not_null<Scalar<DataVector>*> electron_fraction,
134 : const gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
135 : const gsl::not_null<tnsr::I<DataVector, 3>*> spatial_velocity,
136 : const gsl::not_null<tnsr::I<DataVector, 3>*> magnetic_field,
137 : const gsl::not_null<Scalar<DataVector>*> div_cleaning_field,
138 : const gsl::not_null<Scalar<DataVector>*> lorentz_factor,
139 : const gsl::not_null<Scalar<DataVector>*> pressure,
140 : const gsl::not_null<Scalar<DataVector>*> temperature,
141 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> numeric_data,
142 : const Mesh<3>& mesh,
143 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
144 : Frame::Inertial>& inv_jacobian,
145 : const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords,
146 : const EquationsOfState::EquationOfState<true, ThermodynamicDim>&
147 : equation_of_state) const {
148 : gh_numeric_id_.set_initial_data(spacetime_metric, pi, phi, numeric_data,
149 : mesh, inv_jacobian, inertial_coords);
150 : const auto spatial_metric = gr::spatial_metric(*spacetime_metric);
151 : const auto inv_spatial_metric =
152 : determinant_and_inverse(spatial_metric).second;
153 : hydro_numeric_id_.set_initial_data(
154 : rest_mass_density, electron_fraction, specific_internal_energy,
155 : spatial_velocity, magnetic_field, div_cleaning_field, lorentz_factor,
156 : pressure, temperature, numeric_data, inv_spatial_metric,
157 : equation_of_state);
158 : }
159 :
160 0 : void pup(PUP::er& p) override;
161 :
162 0 : friend bool operator==(const NumericInitialData& lhs,
163 : const NumericInitialData& rhs);
164 :
165 : private:
166 0 : GhNumericId gh_numeric_id_{};
167 0 : HydroNumericId hydro_numeric_id_{};
168 : };
169 :
170 0 : namespace Actions {
171 :
172 : /*!
173 : * \brief Dispatch loading numeric initial data from files.
174 : *
175 : * Place this action before
176 : * grmhd::GhValenciaDivClean::Actions::SetNumericInitialData in the action list.
177 : * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
178 : * invoked by this action.
179 : */
180 1 : struct SetInitialData {
181 0 : using const_global_cache_tags =
182 : tmpl::list<evolution::initial_data::Tags::InitialData>;
183 :
184 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
185 : typename ArrayIndex, typename ActionList,
186 : typename ParallelComponent>
187 0 : static Parallel::iterable_action_return_t apply(
188 : db::DataBox<DbTagsList>& box,
189 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
190 : Parallel::GlobalCache<Metavariables>& cache,
191 : const ArrayIndex& array_index, const ActionList /*meta*/,
192 : const ParallelComponent* const parallel_component) {
193 : // Dispatch to the correct `apply` overload based on type of initial data
194 : using initial_data_classes =
195 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
196 : evolution::initial_data::InitialData>;
197 : return call_with_dynamic_type<Parallel::iterable_action_return_t,
198 : initial_data_classes>(
199 : &db::get<evolution::initial_data::Tags::InitialData>(box),
200 : [&box, &cache, &array_index,
201 : ¶llel_component](const auto* const initial_data) {
202 : return apply(make_not_null(&box), *initial_data, cache, array_index,
203 : parallel_component);
204 : });
205 : }
206 :
207 : private:
208 0 : static constexpr size_t Dim = 3;
209 :
210 : // Numeric initial data
211 : template <typename DbTagsList, typename Metavariables, typename ArrayIndex,
212 : typename ParallelComponent>
213 0 : static Parallel::iterable_action_return_t apply(
214 : const gsl::not_null<db::DataBox<DbTagsList>*> /*box*/,
215 : const NumericInitialData& initial_data,
216 : Parallel::GlobalCache<Metavariables>& cache,
217 : const ArrayIndex& array_index, const ParallelComponent* const /*meta*/) {
218 : // If we are using GH Numeric ID, then we don't have to set Pi and Phi since
219 : // we are reading them in. Also we only need to mutate this tag once so do
220 : // it on the first element.
221 : if (is_zeroth_element(array_index) and
222 : std::holds_alternative<gh::NumericInitialData::GhVars>(
223 : initial_data.gh_numeric_id().selected_variables())) {
224 : Parallel::mutate<gh::Tags::SetPiAndPhiFromConstraints,
225 : gh::gauges::SetPiAndPhiFromConstraintsCacheMutator>(
226 : cache, false);
227 : }
228 : // Select the subset of the available variables that we want to read from
229 : // the volume data file
230 : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
231 : importers::Tags::Selected, NumericInitialData::all_vars>>
232 : selected_fields{};
233 : initial_data.select_for_import(make_not_null(&selected_fields));
234 : // Dispatch loading the variables from the volume data file
235 : // - Not using `ckLocalBranch` here to make sure the simple action
236 : // invocation is asynchronous.
237 : auto& reader_component = Parallel::get_parallel_component<
238 : importers::ElementDataReader<Metavariables>>(cache);
239 : Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
240 : 3, NumericInitialData::all_vars, ParallelComponent>>(
241 : reader_component, initial_data.importer_options(),
242 : initial_data.volume_data_id(), std::move(selected_fields));
243 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
244 : }
245 :
246 : // "AnalyticData"-type initial data
247 : template <typename DbTagsList, typename InitialData, typename Metavariables,
248 : typename ArrayIndex, typename ParallelComponent>
249 0 : static Parallel::iterable_action_return_t apply(
250 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
251 : const InitialData& initial_data,
252 : Parallel::GlobalCache<Metavariables>& /*cache*/,
253 : const ArrayIndex& /*array_index*/,
254 : const ParallelComponent* const /*meta*/) {
255 : // Get ADM + hydro variables from analytic data / solution
256 : const auto& [coords, mesh, inv_jacobian] = [&box]() {
257 : if constexpr (db::tag_is_retrievable_v<
258 : evolution::dg::subcell::Tags::ActiveGrid,
259 : db::DataBox<DbTagsList>>) {
260 : const bool on_subcell =
261 : db::get<evolution::dg::subcell::Tags::ActiveGrid>(*box) ==
262 : evolution::dg::subcell::ActiveGrid::Subcell;
263 : if (on_subcell) {
264 : return std::forward_as_tuple(
265 : db::get<evolution::dg::subcell::Tags::Coordinates<
266 : Dim, Frame::Inertial>>(*box),
267 : db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(*box),
268 : db::get<evolution::dg::subcell::fd::Tags::
269 : InverseJacobianLogicalToInertial<Dim>>(*box));
270 : }
271 : }
272 : return std::forward_as_tuple(
273 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box),
274 : db::get<domain::Tags::Mesh<Dim>>(*box),
275 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
276 : Frame::Inertial>>(*box));
277 : }();
278 : auto vars = evolution::Initialization::initial_data(
279 : initial_data, coords, db::get<::Tags::Time>(*box),
280 : tmpl::append<tmpl::list<gr::Tags::SpatialMetric<DataVector, 3>,
281 : gr::Tags::Lapse<DataVector>,
282 : gr::Tags::Shift<DataVector, 3>,
283 : gr::Tags::ExtrinsicCurvature<DataVector, 3>>,
284 : hydro::grmhd_tags<DataVector>>{});
285 : const auto& spatial_metric =
286 : get<gr::Tags::SpatialMetric<DataVector, 3>>(vars);
287 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(vars);
288 : const auto& shift = get<gr::Tags::Shift<DataVector, 3>>(vars);
289 : const auto& extrinsic_curvature =
290 : get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(vars);
291 :
292 : // Compute GH vars from ADM vars
293 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
294 : gh::Tags::Pi<DataVector, 3>, gh::Tags::Phi<DataVector, 3>>(
295 : &gh::initial_gh_variables_from_adm<3>, box, spatial_metric, lapse,
296 : shift, extrinsic_curvature, mesh, inv_jacobian, coords);
297 :
298 : // Move hydro vars directly into the DataBox
299 : tmpl::for_each<hydro::grmhd_tags<DataVector>>(
300 : [&box, &vars](const auto tag_v) {
301 : using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
302 : Initialization::mutate_assign<tmpl::list<tag>>(
303 : box, std::move(get<tag>(vars)));
304 : });
305 :
306 : // No need to import numeric initial data, so we terminate the phase by
307 : // pausing the algorithm on this element
308 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
309 : }
310 : };
311 :
312 : /*!
313 : * \brief Receive numeric initial data loaded by
314 : * grmhd::GhValenciaDivClean::Actions::SetInitialData.
315 : *
316 : * Place this action in the action list after
317 : * grmhd::GhValenciaDivClean::Actions::SetInitialData to wait until the
318 : * data for this element has arrived, and then compute the GH variables and the
319 : * remaining primitive variables and store them in the DataBox to be used as
320 : * initial data.
321 : *
322 : * This action modifies the GH system tags (spacetime metric, pi, phi) and the
323 : * tags listed in `hydro::grmhd_tags` in the DataBox (i.e., the hydro
324 : * primitives). It does not modify conservative variables, so it relies on a
325 : * primitive-to-conservative update in the action list before the evolution can
326 : * start.
327 : *
328 : * \requires This action requires an equation of state, which is retrieved from
329 : * the DataBox as `hydro::Tags::GrmhdEquationOfState`.
330 : */
331 1 : struct ReceiveNumericInitialData {
332 0 : static constexpr size_t Dim = 3;
333 0 : using inbox_tags =
334 : tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
335 :
336 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
337 : typename ActionList, typename ParallelComponent>
338 0 : static Parallel::iterable_action_return_t apply(
339 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
340 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
341 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
342 : const ParallelComponent* const /*meta*/) {
343 : auto& inbox =
344 : tuples::get<importers::Tags::VolumeData<NumericInitialData::all_vars>>(
345 : inboxes);
346 : const auto& initial_data = dynamic_cast<const NumericInitialData&>(
347 : db::get<evolution::initial_data::Tags::InitialData>(box));
348 : const size_t volume_data_id = initial_data.volume_data_id();
349 : if (inbox.find(volume_data_id) == inbox.end()) {
350 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
351 : }
352 : auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
353 :
354 : const auto& [coords, mesh, inv_jacobian] = [&box]() {
355 : if constexpr (db::tag_is_retrievable_v<
356 : evolution::dg::subcell::Tags::ActiveGrid,
357 : db::DataBox<DbTagsList>>) {
358 : const bool on_subcell =
359 : db::get<evolution::dg::subcell::Tags::ActiveGrid>(box) ==
360 : evolution::dg::subcell::ActiveGrid::Subcell;
361 : if (on_subcell) {
362 : return std::forward_as_tuple(
363 : db::get<evolution::dg::subcell::Tags::Coordinates<
364 : Dim, Frame::Inertial>>(box),
365 : db::get<evolution::dg::subcell::Tags::Mesh<Dim>>(box),
366 : db::get<evolution::dg::subcell::fd::Tags::
367 : InverseJacobianLogicalToInertial<Dim>>(box));
368 : }
369 : }
370 : return std::forward_as_tuple(
371 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box),
372 : db::get<domain::Tags::Mesh<Dim>>(box),
373 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
374 : Frame::Inertial>>(box));
375 : }();
376 : const auto& equation_of_state =
377 : db::get<hydro::Tags::GrmhdEquationOfState>(box);
378 :
379 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
380 : gh::Tags::Pi<DataVector, 3>, gh::Tags::Phi<DataVector, 3>,
381 : hydro::Tags::RestMassDensity<DataVector>,
382 : hydro::Tags::ElectronFraction<DataVector>,
383 : hydro::Tags::SpecificInternalEnergy<DataVector>,
384 : hydro::Tags::SpatialVelocity<DataVector, 3>,
385 : hydro::Tags::MagneticField<DataVector, 3>,
386 : hydro::Tags::DivergenceCleaningField<DataVector>,
387 : hydro::Tags::LorentzFactor<DataVector>,
388 : hydro::Tags::Pressure<DataVector>,
389 : hydro::Tags::Temperature<DataVector>>(
390 : [&initial_data, &numeric_data, &equation_of_state](
391 : const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
392 : const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
393 : const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
394 : const gsl::not_null<Scalar<DataVector>*> rest_mass_density,
395 : const gsl::not_null<Scalar<DataVector>*> electron_fraction,
396 : const gsl::not_null<Scalar<DataVector>*> specific_internal_energy,
397 : const gsl::not_null<tnsr::I<DataVector, 3>*> spatial_velocity,
398 : const gsl::not_null<tnsr::I<DataVector, 3>*> magnetic_field,
399 : const gsl::not_null<Scalar<DataVector>*> div_cleaning_field,
400 : const gsl::not_null<Scalar<DataVector>*> lorentz_factor,
401 : const gsl::not_null<Scalar<DataVector>*> pressure,
402 : const gsl::not_null<Scalar<DataVector>*> temperature,
403 : const auto& local_mesh, const auto& local_inv_jacobian,
404 : const auto& local_inertial_coords) {
405 : initial_data.set_initial_data(
406 : spacetime_metric, pi, phi, rest_mass_density, electron_fraction,
407 : specific_internal_energy, spatial_velocity, magnetic_field,
408 : div_cleaning_field, lorentz_factor, pressure, temperature,
409 : make_not_null(&numeric_data), local_mesh, local_inv_jacobian,
410 : local_inertial_coords, equation_of_state);
411 : },
412 : make_not_null(&box), mesh, inv_jacobian, coords);
413 :
414 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
415 : }
416 : };
417 :
418 : } // namespace Actions
419 :
420 : } // namespace grmhd::GhValenciaDivClean
|
