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 =
137 : tmpl::push_back<importers::ImporterOptions::tags_list, 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 : std::optional<double> observation_value_epsilon,
164 : bool enable_interpolation,
165 : std::variant<AdmVars, GhVars> selected_variables);
166 :
167 0 : const importers::ImporterOptions& importer_options() const {
168 : return importer_options_;
169 : }
170 :
171 0 : const std::variant<AdmVars, GhVars>& selected_variables() const {
172 : return selected_variables_;
173 : }
174 :
175 : /*!
176 : * \brief Unique identifier for loading this volume data
177 : *
178 : * Involves a hash of the type name and the volume data file names.
179 : */
180 1 : size_t volume_data_id() const;
181 :
182 : /*!
183 : * \brief Selects which of the `fields` to import based on the choices in the
184 : * input-file options
185 : *
186 : * The `fields` are all datasets that are available to import, represented by
187 : * `importers::Tags::Selected<Tag>` tags. We select only those that we need by
188 : * setting their dataset name.
189 : */
190 : template <typename... AllTags>
191 1 : void select_for_import(
192 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> fields) const {
193 : // Select the subset of the available variables that we want to read from
194 : // the volume data file
195 : std::visit(
196 : [&fields](const auto& vars) {
197 : // This lambda is invoked with the set of vars selected in the input
198 : // file, which map to the tensor names that should be read from the H5
199 : // file
200 : using selected_vars = std::decay_t<decltype(vars)>;
201 : // Get the mapped tensor name from the input file and select it in the
202 : // set of all possible vars.
203 : tmpl::for_each<typename selected_vars::tags_list>(
204 : [&fields, &vars](const auto tag_v) {
205 : using tag = typename std::decay_t<decltype(tag_v)>::type::tag;
206 : get<importers::Tags::Selected<tag>>(*fields) =
207 : get<VarName<tag>>(vars);
208 : });
209 : },
210 : selected_variables_);
211 : }
212 :
213 : /*!
214 : * \brief Set GH initial data given numeric data loaded from files
215 : *
216 : * The `numeric_data` contains the datasets selected above (and possibly
217 : * more). We either set the GH variables directly, or compute them from the
218 : * ADM variables.
219 : */
220 : template <typename... AllTags>
221 1 : void set_initial_data(
222 : const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
223 : const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
224 : const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
225 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> numeric_data,
226 : const Mesh<3>& mesh,
227 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
228 : Frame::Inertial>& inv_jacobian,
229 : const tnsr::I<DataVector, 3, Frame::Inertial> inertial_coords) const {
230 : if (std::holds_alternative<NumericInitialData::GhVars>(
231 : selected_variables_)) {
232 : // We have loaded the GH system variables from the file, so just move the
233 : // data for spacetime_metric and Pi into the DataBox directly, with no
234 : // conversion needed. Set Phi to the spatial derivative of the spacetime
235 : // metric to enforce the 3-index constraint.
236 : *spacetime_metric = std::move(
237 : get<gr::Tags::SpacetimeMetric<DataVector, 3>>(*numeric_data));
238 : *pi = std::move(get<Tags::Pi<DataVector, 3>>(*numeric_data));
239 : *phi = get<Tags::Phi<DataVector, 3>>(*numeric_data);
240 : } else if (std::holds_alternative<NumericInitialData::AdmVars>(
241 : selected_variables_)) {
242 : // We have loaded ADM variables from the file. Convert to GH variables.
243 : const auto& spatial_metric =
244 : get<gr::Tags::SpatialMetric<DataVector, 3>>(*numeric_data);
245 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(*numeric_data);
246 : const auto& shift = get<gr::Tags::Shift<DataVector, 3>>(*numeric_data);
247 : const auto& extrinsic_curvature =
248 : get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(*numeric_data);
249 :
250 : initial_gh_variables_from_adm(spacetime_metric, pi, phi, spatial_metric,
251 : lapse, shift, extrinsic_curvature, mesh,
252 : inv_jacobian, inertial_coords);
253 : } else {
254 : ERROR(
255 : "These initial data variables are not implemented yet. Please add "
256 : "an implementation to gh::NumericInitialData.");
257 : }
258 : }
259 :
260 0 : void pup(PUP::er& p) override;
261 :
262 0 : friend bool operator==(const NumericInitialData& lhs,
263 : const NumericInitialData& rhs);
264 :
265 : private:
266 0 : importers::ImporterOptions importer_options_;
267 0 : std::variant<AdmVars, GhVars> selected_variables_{};
268 : };
269 :
270 0 : namespace Actions {
271 :
272 : /*!
273 : * \brief Dispatch loading numeric initial data from files or set analytic
274 : * initial data.
275 : *
276 : * Place this action before
277 : * gh::Actions::ReceiveNumericInitialData in the action list.
278 : * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
279 : * invoked by this action.
280 : * Analytic initial data is set directly by this action and terminates the
281 : * phase.
282 : */
283 1 : struct SetInitialData {
284 0 : using const_global_cache_tags =
285 : tmpl::list<evolution::initial_data::Tags::InitialData>;
286 :
287 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
288 : typename ArrayIndex, typename ActionList,
289 : typename ParallelComponent>
290 0 : static Parallel::iterable_action_return_t apply(
291 : db::DataBox<DbTagsList>& box,
292 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
293 : Parallel::GlobalCache<Metavariables>& cache,
294 : const ArrayIndex& array_index, const ActionList /*meta*/,
295 : const ParallelComponent* const parallel_component) {
296 : // Dispatch to the correct `apply` overload based on type of initial data
297 : using initial_data_classes =
298 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
299 : evolution::initial_data::InitialData>;
300 : return call_with_dynamic_type<Parallel::iterable_action_return_t,
301 : initial_data_classes>(
302 : &db::get<evolution::initial_data::Tags::InitialData>(box),
303 : [&box, &cache, &array_index,
304 : ¶llel_component](const auto* const initial_data) {
305 : return apply(make_not_null(&box), *initial_data, cache, array_index,
306 : parallel_component);
307 : });
308 : }
309 :
310 : private:
311 : // Numeric initial data
312 : template <typename DbTagsList, typename Metavariables, typename ArrayIndex,
313 : typename ParallelComponent>
314 0 : static Parallel::iterable_action_return_t apply(
315 : const gsl::not_null<db::DataBox<DbTagsList>*> /*box*/,
316 : const NumericInitialData& initial_data,
317 : Parallel::GlobalCache<Metavariables>& cache,
318 : const ArrayIndex& array_index, const ParallelComponent* const /*meta*/) {
319 : // If we are using GH Numeric ID, then we don't have to set Pi and Phi since
320 : // we are reading them in. Also we only need to mutate this tag once so do
321 : // it on the first element.
322 : if (is_zeroth_element(array_index) and
323 : std::holds_alternative<NumericInitialData::GhVars>(
324 : initial_data.selected_variables())) {
325 : Parallel::mutate<Tags::SetPiAndPhiFromConstraints,
326 : gh::gauges::SetPiAndPhiFromConstraintsCacheMutator>(
327 : cache, false);
328 : }
329 :
330 : // Select the subset of the available variables that we want to read from
331 : // the volume data file
332 : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
333 : importers::Tags::Selected, NumericInitialData::all_vars>>
334 : selected_fields{};
335 : initial_data.select_for_import(make_not_null(&selected_fields));
336 : // Dispatch loading the variables from the volume data file
337 : // - Not using `ckLocalBranch` here to make sure the simple action
338 : // invocation is asynchronous.
339 : auto& reader_component = Parallel::get_parallel_component<
340 : importers::ElementDataReader<Metavariables>>(cache);
341 : Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
342 : 3, NumericInitialData::all_vars, ParallelComponent>>(
343 : reader_component, initial_data.importer_options(),
344 : initial_data.volume_data_id(), std::move(selected_fields));
345 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
346 : }
347 :
348 : // "AnalyticData"-type initial data
349 : template <typename DbTagsList, typename InitialData, typename Metavariables,
350 : typename ArrayIndex, typename ParallelComponent>
351 0 : static Parallel::iterable_action_return_t apply(
352 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
353 : const InitialData& initial_data,
354 : Parallel::GlobalCache<Metavariables>& /*cache*/,
355 : const ArrayIndex& /*array_index*/,
356 : const ParallelComponent* const /*meta*/) {
357 : static constexpr size_t Dim = Metavariables::volume_dim;
358 :
359 : // Get ADM variables from analytic data / solution
360 : const auto& x =
361 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box);
362 : const auto adm_vars = evolution::Initialization::initial_data(
363 : initial_data, x, db::get<::Tags::Time>(*box),
364 : tmpl::list<gr::Tags::SpatialMetric<DataVector, Dim>,
365 : gr::Tags::Lapse<DataVector>,
366 : gr::Tags::Shift<DataVector, Dim>,
367 : gr::Tags::ExtrinsicCurvature<DataVector, Dim>>{});
368 : const auto& spatial_metric =
369 : get<gr::Tags::SpatialMetric<DataVector, Dim>>(adm_vars);
370 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(adm_vars);
371 : const auto& shift = get<gr::Tags::Shift<DataVector, Dim>>(adm_vars);
372 : const auto& extrinsic_curvature =
373 : get<gr::Tags::ExtrinsicCurvature<DataVector, Dim>>(adm_vars);
374 :
375 : // Compute GH vars from ADM vars
376 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
377 : const auto& inv_jacobian =
378 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
379 : Frame::Inertial>>(*box);
380 : const auto& inertial_coords =
381 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box);
382 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, Dim>,
383 : Tags::Pi<DataVector, Dim>, Tags::Phi<DataVector, Dim>>(
384 : &gh::initial_gh_variables_from_adm<Dim>, box, spatial_metric, lapse,
385 : shift, extrinsic_curvature, mesh, inv_jacobian, inertial_coords);
386 :
387 : // No need to import numeric initial data, so we terminate the phase by
388 : // pausing the algorithm on this element
389 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
390 : }
391 : };
392 :
393 : /*!
394 : * \brief Receive numeric initial data loaded by gh::Actions::SetInitialData.
395 : *
396 : * Place this action in the action list after
397 : * gh::Actions::SetInitialData to wait until the data
398 : * for this element has arrived, and then transform the data to GH variables and
399 : * store it in the DataBox to be used as initial data.
400 : *
401 : * This action modifies the following tags in the DataBox:
402 : * - gr::Tags::SpacetimeMetric<DataVector, 3>
403 : * - gh::Tags::Pi<DataVector, 3>
404 : * - gh::Tags::Phi<DataVector, 3>
405 : */
406 1 : struct ReceiveNumericInitialData {
407 0 : static constexpr size_t Dim = 3;
408 0 : using inbox_tags =
409 : tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
410 :
411 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
412 : typename ActionList, typename ParallelComponent>
413 0 : static Parallel::iterable_action_return_t apply(
414 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
415 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
416 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
417 : const ParallelComponent* const /*meta*/) {
418 : auto& inbox =
419 : tuples::get<importers::Tags::VolumeData<NumericInitialData::all_vars>>(
420 : inboxes);
421 : const auto& initial_data = dynamic_cast<const NumericInitialData&>(
422 : db::get<evolution::initial_data::Tags::InitialData>(box));
423 : const auto& volume_data_id = initial_data.volume_data_id();
424 : if (inbox.find(volume_data_id) == inbox.end()) {
425 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
426 : }
427 : auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
428 :
429 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
430 : const auto& inv_jacobian =
431 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
432 : Frame::Inertial>>(box);
433 : const auto& inertial_coords =
434 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box);
435 :
436 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
437 : Tags::Pi<DataVector, 3>, Tags::Phi<DataVector, 3>>(
438 : [&initial_data, &numeric_data, &mesh, &inv_jacobian, &inertial_coords](
439 : const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
440 : const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
441 : const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi) {
442 : initial_data.set_initial_data(spacetime_metric, pi, phi,
443 : make_not_null(&numeric_data), mesh,
444 : inv_jacobian, inertial_coords);
445 : },
446 : make_not_null(&box));
447 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
448 : }
449 : };
450 :
451 : } // namespace Actions
452 : } // namespace gh
|
