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