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 <pup.h>
8 : #include <string>
9 : #include <variant>
10 :
11 : #include "DataStructures/DataBox/DataBox.hpp"
12 : #include "DataStructures/DataBox/Tag.hpp"
13 : #include "DataStructures/DataVector.hpp"
14 : #include "DataStructures/TaggedTuple.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "Domain/Structure/ElementId.hpp"
17 : #include "Domain/Tags.hpp"
18 : #include "Evolution/Initialization/InitialData.hpp"
19 : #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp"
20 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
21 : #include "IO/Importers/Actions/ReadVolumeData.hpp"
22 : #include "IO/Importers/ElementDataReader.hpp"
23 : #include "IO/Importers/Tags.hpp"
24 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
25 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
26 : #include "Parallel/AlgorithmExecution.hpp"
27 : #include "Parallel/GlobalCache.hpp"
28 : #include "Parallel/Invoke.hpp"
29 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
30 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
31 : #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
32 : #include "Utilities/CallWithDynamicType.hpp"
33 : #include "Utilities/ErrorHandling/Error.hpp"
34 : #include "Utilities/Gsl.hpp"
35 : #include "Utilities/Serialization/CharmPupable.hpp"
36 : #include "Utilities/TMPL.hpp"
37 :
38 : /// \cond
39 : namespace Tags {
40 : struct Time;
41 : } // namespace Tags
42 : /// \endcond
43 :
44 : namespace gh {
45 :
46 : /*!
47 : * \brief Compute initial GH variables from ADM variables
48 : *
49 : * - The spacetime metric is assembled from the spatial metric, lapse, and
50 : * shift. See `gr::spacetime_metric` for details.
51 : * - Phi is set to the numerical derivative of the spacetime metric. This
52 : * ensures that the 3-index constraint is initially satisfied.
53 : * - Pi is computed by choosing the time derivatives of lapse and shift to be
54 : * zero. The `gh::gauges::SetPiAndPhiFromConstraints` mutator exists to
55 : * override Pi later in the algorithm (it should be combined with this
56 : * function).
57 : */
58 : template <size_t Dim>
59 1 : void initial_gh_variables_from_adm(
60 : gsl::not_null<tnsr::aa<DataVector, Dim>*> spacetime_metric,
61 : gsl::not_null<tnsr::aa<DataVector, Dim>*> pi,
62 : gsl::not_null<tnsr::iaa<DataVector, Dim>*> phi,
63 : const tnsr::ii<DataVector, Dim>& spatial_metric,
64 : const Scalar<DataVector>& lapse, const tnsr::I<DataVector, Dim>& shift,
65 : const tnsr::ii<DataVector, Dim>& extrinsic_curvature, const Mesh<Dim>& mesh,
66 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
67 : Frame::Inertial>& inv_jacobian,
68 : const tnsr::I<DataVector, Dim, Frame::Inertial>& inertial_coords);
69 :
70 : /*!
71 : * \brief Numeric initial data loaded from volume data files
72 : *
73 : * This class can be factory-created in the input file to start an evolution
74 : * from numeric initial data. It selects the set of variables to load from
75 : * the volume data file (ADM or GH variables).
76 : */
77 1 : class NumericInitialData : public evolution::initial_data::InitialData {
78 : public:
79 : /// Name of a variable in the volume data file
80 : template <typename Tag>
81 1 : struct VarName {
82 0 : using tag = Tag;
83 0 : static std::string name() { return db::tag_name<Tag>(); }
84 0 : using type = std::string;
85 0 : static constexpr Options::String help =
86 : "Name of the variable in the volume data file";
87 : };
88 :
89 : // These are the sets of variables that we support loading from volume data
90 : // files:
91 : // - ADM variables
92 0 : using adm_vars =
93 : tmpl::list<gr::Tags::SpatialMetric<DataVector, 3>,
94 : gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>,
95 : gr::Tags::ExtrinsicCurvature<DataVector, 3>>;
96 0 : struct AdmVars : tuples::tagged_tuple_from_typelist<
97 : db::wrap_tags_in<VarName, adm_vars>> {
98 0 : static constexpr Options::String help =
99 : "ADM variables: 'Lapse', 'Shift', 'SpatialMetric' and "
100 : "'ExtrinsicCurvature'. The initial GH variables will be computed "
101 : "from these numeric fields, as well as their numeric spatial "
102 : "derivatives on the computational grid. The GH variable Pi will be set "
103 : "to satisfy the gauge constraint using the evolution gauge. The GH "
104 : "variable Phi will be set to satisfy the 3-index constraint.";
105 0 : using options = tags_list;
106 : using TaggedTuple::TaggedTuple;
107 : };
108 :
109 : // - Generalized harmonic variables
110 0 : using gh_vars = tmpl::list<gr::Tags::SpacetimeMetric<DataVector, 3>,
111 : Tags::Pi<DataVector, 3>, Tags::Phi<DataVector, 3>>;
112 0 : struct GhVars
113 : : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<VarName, gh_vars>> {
114 0 : static constexpr Options::String help =
115 : "GH variables: 'SpacetimeMetric', 'Pi', and 'Phi'. These variables are "
116 : "used to set the initial data directly.";
117 0 : using options = tags_list;
118 : using TaggedTuple::TaggedTuple;
119 : };
120 :
121 : // Collect all variables that we support loading from volume data files.
122 : // Remember to `tmpl::remove_duplicates` when adding overlapping sets of
123 : // vars.
124 0 : using all_vars = tmpl::append<adm_vars, gh_vars>;
125 :
126 : // Input-file options
127 0 : struct Variables {
128 : // The user can supply any of these choices of variables in the input
129 : // file
130 0 : using type = std::variant<AdmVars, GhVars>;
131 0 : static constexpr Options::String help =
132 : "Set of initial data variables from which the generalized harmonic "
133 : "system variables are computed.";
134 : };
135 :
136 0 : using options = tmpl::list<importers::OptionTags::VolumeData, Variables>;
137 :
138 0 : static constexpr Options::String help =
139 : "Numeric initial data loaded from volume data files";
140 :
141 0 : NumericInitialData() = default;
142 0 : NumericInitialData(const NumericInitialData& rhs) = default;
143 0 : NumericInitialData& operator=(const NumericInitialData& rhs) = default;
144 0 : NumericInitialData(NumericInitialData&& /*rhs*/) = default;
145 0 : NumericInitialData& operator=(NumericInitialData&& /*rhs*/) = default;
146 0 : ~NumericInitialData() = default;
147 :
148 : /// \cond
149 : explicit NumericInitialData(CkMigrateMessage* msg);
150 : using PUP::able::register_constructor;
151 : WRAPPED_PUPable_decl_template(NumericInitialData);
152 : /// \endcond
153 :
154 0 : std::unique_ptr<evolution::initial_data::InitialData> get_clone()
155 : const override {
156 : return std::make_unique<NumericInitialData>(*this);
157 : }
158 :
159 0 : NumericInitialData(importers::ImporterOptions importer_options,
160 : std::variant<AdmVars, GhVars> selected_variables);
161 :
162 0 : const importers::ImporterOptions& importer_options() const {
163 : return importer_options_;
164 : }
165 :
166 0 : const std::variant<AdmVars, GhVars>& selected_variables() const {
167 : return selected_variables_;
168 : }
169 :
170 : /*!
171 : * \brief Unique identifier for loading this volume data
172 : *
173 : * Involves a hash of the type name and the volume data file names.
174 : */
175 1 : size_t volume_data_id() const;
176 :
177 : /*!
178 : * \brief Selects which of the `fields` to import based on the choices in the
179 : * input-file options
180 : *
181 : * The `fields` are all datasets that are available to import, represented by
182 : * `importers::Tags::Selected<Tag>` tags. We select only those that we need by
183 : * setting their dataset name.
184 : */
185 : template <typename... AllTags>
186 1 : void select_for_import(
187 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> fields) const {
188 : // Select the subset of the available variables that we want to read from
189 : // the volume data file
190 : std::visit(
191 : [&fields](const auto& vars) {
192 : // This lambda is invoked with the set of vars selected in the input
193 : // file, which map to the tensor names that should be read from the H5
194 : // file
195 : using selected_vars = std::decay_t<decltype(vars)>;
196 : // Get the mapped tensor name from the input file and select it in the
197 : // set of all possible vars.
198 : tmpl::for_each<typename selected_vars::tags_list>(
199 : [&fields, &vars](const auto tag_v) {
200 : using tag = typename std::decay_t<decltype(tag_v)>::type::tag;
201 : get<importers::Tags::Selected<tag>>(*fields) =
202 : get<VarName<tag>>(vars);
203 : });
204 : },
205 : selected_variables_);
206 : }
207 :
208 : /*!
209 : * \brief Set GH initial data given numeric data loaded from files
210 : *
211 : * The `numeric_data` contains the datasets selected above (and possibly
212 : * more). We either set the GH variables directly, or compute them from the
213 : * ADM variables.
214 : */
215 : template <typename... AllTags>
216 1 : void set_initial_data(
217 : const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
218 : const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
219 : const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
220 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> numeric_data,
221 : const Mesh<3>& mesh,
222 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
223 : Frame::Inertial>& inv_jacobian,
224 : const tnsr::I<DataVector, 3, Frame::Inertial> inertial_coords) const {
225 : if (std::holds_alternative<NumericInitialData::GhVars>(
226 : selected_variables_)) {
227 : // We have loaded the GH system variables from the file, so just move the
228 : // data for spacetime_metric and Pi into the DataBox directly, with no
229 : // conversion needed. Set Phi to the spatial derivative of the spacetime
230 : // metric to enforce the 3-index constraint.
231 : *spacetime_metric = std::move(
232 : get<gr::Tags::SpacetimeMetric<DataVector, 3>>(*numeric_data));
233 : *pi = std::move(get<Tags::Pi<DataVector, 3>>(*numeric_data));
234 : *phi = get<Tags::Phi<DataVector, 3>>(*numeric_data);
235 : } else if (std::holds_alternative<NumericInitialData::AdmVars>(
236 : selected_variables_)) {
237 : // We have loaded ADM variables from the file. Convert to GH variables.
238 : const auto& spatial_metric =
239 : get<gr::Tags::SpatialMetric<DataVector, 3>>(*numeric_data);
240 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(*numeric_data);
241 : const auto& shift = get<gr::Tags::Shift<DataVector, 3>>(*numeric_data);
242 : const auto& extrinsic_curvature =
243 : get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(*numeric_data);
244 :
245 : initial_gh_variables_from_adm(spacetime_metric, pi, phi, spatial_metric,
246 : lapse, shift, extrinsic_curvature, mesh,
247 : inv_jacobian, inertial_coords);
248 : } else {
249 : ERROR(
250 : "These initial data variables are not implemented yet. Please add "
251 : "an implementation to gh::NumericInitialData.");
252 : }
253 : }
254 :
255 0 : void pup(PUP::er& p) override;
256 :
257 0 : friend bool operator==(const NumericInitialData& lhs,
258 : const NumericInitialData& rhs);
259 :
260 : private:
261 0 : importers::ImporterOptions importer_options_;
262 0 : std::variant<AdmVars, GhVars> selected_variables_{};
263 : };
264 :
265 0 : namespace Actions {
266 :
267 : /*!
268 : * \brief Dispatch loading numeric initial data from files or set analytic
269 : * initial data.
270 : *
271 : * Place this action before
272 : * gh::Actions::ReceiveNumericInitialData in the action list.
273 : * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
274 : * invoked by this action.
275 : * Analytic initial data is set directly by this action and terminates the
276 : * phase.
277 : */
278 1 : struct SetInitialData {
279 0 : using const_global_cache_tags =
280 : tmpl::list<evolution::initial_data::Tags::InitialData>;
281 :
282 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
283 : typename ArrayIndex, typename ActionList,
284 : typename ParallelComponent>
285 0 : static Parallel::iterable_action_return_t apply(
286 : db::DataBox<DbTagsList>& box,
287 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
288 : Parallel::GlobalCache<Metavariables>& cache,
289 : const ArrayIndex& array_index, const ActionList /*meta*/,
290 : const ParallelComponent* const parallel_component) {
291 : // Dispatch to the correct `apply` overload based on type of initial data
292 : using initial_data_classes =
293 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
294 : evolution::initial_data::InitialData>;
295 : return call_with_dynamic_type<Parallel::iterable_action_return_t,
296 : initial_data_classes>(
297 : &db::get<evolution::initial_data::Tags::InitialData>(box),
298 : [&box, &cache, &array_index,
299 : ¶llel_component](const auto* const initial_data) {
300 : return apply(make_not_null(&box), *initial_data, cache, array_index,
301 : parallel_component);
302 : });
303 : }
304 :
305 : private:
306 : // Numeric initial data
307 : template <typename DbTagsList, typename Metavariables, typename ArrayIndex,
308 : typename ParallelComponent>
309 0 : static Parallel::iterable_action_return_t apply(
310 : const gsl::not_null<db::DataBox<DbTagsList>*> /*box*/,
311 : const NumericInitialData& initial_data,
312 : Parallel::GlobalCache<Metavariables>& cache,
313 : const ArrayIndex& array_index, const ParallelComponent* const /*meta*/) {
314 : // If we are using GH Numeric ID, then we don't have to set Pi and Phi since
315 : // we are reading them in. Also we only need to mutate this tag once so do
316 : // it on the first element.
317 : if (is_zeroth_element(array_index) and
318 : std::holds_alternative<NumericInitialData::GhVars>(
319 : initial_data.selected_variables())) {
320 : Parallel::mutate<Tags::SetPiAndPhiFromConstraints,
321 : gh::gauges::SetPiAndPhiFromConstraintsCacheMutator>(
322 : cache, false);
323 : }
324 :
325 : // Select the subset of the available variables that we want to read from
326 : // the volume data file
327 : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
328 : importers::Tags::Selected, NumericInitialData::all_vars>>
329 : selected_fields{};
330 : initial_data.select_for_import(make_not_null(&selected_fields));
331 : // Dispatch loading the variables from the volume data file
332 : // - Not using `ckLocalBranch` here to make sure the simple action
333 : // invocation is asynchronous.
334 : auto& reader_component = Parallel::get_parallel_component<
335 : importers::ElementDataReader<Metavariables>>(cache);
336 : Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
337 : 3, NumericInitialData::all_vars, ParallelComponent>>(
338 : reader_component, initial_data.importer_options(),
339 : initial_data.volume_data_id(), std::move(selected_fields));
340 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
341 : }
342 :
343 : // "AnalyticData"-type initial data
344 : template <typename DbTagsList, typename InitialData, typename Metavariables,
345 : typename ArrayIndex, typename ParallelComponent>
346 0 : static Parallel::iterable_action_return_t apply(
347 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
348 : const InitialData& initial_data,
349 : Parallel::GlobalCache<Metavariables>& /*cache*/,
350 : const ArrayIndex& /*array_index*/,
351 : const ParallelComponent* const /*meta*/) {
352 : static constexpr size_t Dim = Metavariables::volume_dim;
353 :
354 : // Get ADM variables from analytic data / solution
355 : const auto& x =
356 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box);
357 : const auto adm_vars = evolution::Initialization::initial_data(
358 : initial_data, x, db::get<::Tags::Time>(*box),
359 : tmpl::list<gr::Tags::SpatialMetric<DataVector, Dim>,
360 : gr::Tags::Lapse<DataVector>,
361 : gr::Tags::Shift<DataVector, Dim>,
362 : gr::Tags::ExtrinsicCurvature<DataVector, Dim>>{});
363 : const auto& spatial_metric =
364 : get<gr::Tags::SpatialMetric<DataVector, Dim>>(adm_vars);
365 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(adm_vars);
366 : const auto& shift = get<gr::Tags::Shift<DataVector, Dim>>(adm_vars);
367 : const auto& extrinsic_curvature =
368 : get<gr::Tags::ExtrinsicCurvature<DataVector, Dim>>(adm_vars);
369 :
370 : // Compute GH vars from ADM vars
371 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
372 : const auto& inv_jacobian =
373 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
374 : Frame::Inertial>>(*box);
375 : const auto& inertial_coords =
376 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box);
377 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, Dim>,
378 : Tags::Pi<DataVector, Dim>, Tags::Phi<DataVector, Dim>>(
379 : &gh::initial_gh_variables_from_adm<Dim>, box, spatial_metric, lapse,
380 : shift, extrinsic_curvature, mesh, inv_jacobian, inertial_coords);
381 :
382 : // No need to import numeric initial data, so we terminate the phase by
383 : // pausing the algorithm on this element
384 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
385 : }
386 : };
387 :
388 : /*!
389 : * \brief Receive numeric initial data loaded by gh::Actions::SetInitialData.
390 : *
391 : * Place this action in the action list after
392 : * gh::Actions::SetInitialData to wait until the data
393 : * for this element has arrived, and then transform the data to GH variables and
394 : * store it in the DataBox to be used as initial data.
395 : *
396 : * This action modifies the following tags in the DataBox:
397 : * - gr::Tags::SpacetimeMetric<DataVector, 3>
398 : * - gh::Tags::Pi<DataVector, 3>
399 : * - gh::Tags::Phi<DataVector, 3>
400 : */
401 1 : struct ReceiveNumericInitialData {
402 0 : static constexpr size_t Dim = 3;
403 0 : using inbox_tags =
404 : tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
405 :
406 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
407 : typename ActionList, typename ParallelComponent>
408 0 : static Parallel::iterable_action_return_t apply(
409 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
410 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
411 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
412 : const ParallelComponent* const /*meta*/) {
413 : auto& inbox =
414 : tuples::get<importers::Tags::VolumeData<NumericInitialData::all_vars>>(
415 : inboxes);
416 : const auto& initial_data = dynamic_cast<const NumericInitialData&>(
417 : db::get<evolution::initial_data::Tags::InitialData>(box));
418 : const auto& volume_data_id = initial_data.volume_data_id();
419 : if (inbox.find(volume_data_id) == inbox.end()) {
420 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
421 : }
422 : auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
423 :
424 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
425 : const auto& inv_jacobian =
426 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
427 : Frame::Inertial>>(box);
428 : const auto& inertial_coords =
429 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box);
430 :
431 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
432 : Tags::Pi<DataVector, 3>, Tags::Phi<DataVector, 3>>(
433 : [&initial_data, &numeric_data, &mesh, &inv_jacobian, &inertial_coords](
434 : const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
435 : const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
436 : const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi) {
437 : initial_data.set_initial_data(spacetime_metric, pi, phi,
438 : make_not_null(&numeric_data), mesh,
439 : inv_jacobian, inertial_coords);
440 : },
441 : make_not_null(&box));
442 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
443 : }
444 : };
445 :
446 : } // namespace Actions
447 : } // namespace gh
|
