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/Tensor/Tensor.hpp"
15 : #include "Domain/Structure/ElementId.hpp"
16 : #include "Domain/Tags.hpp"
17 : #include "Evolution/Initialization/InitialData.hpp"
18 : #include "Evolution/Systems/GeneralizedHarmonic/GaugeSourceFunctions/SetPiAndPhiFromConstraints.hpp"
19 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
20 : #include "IO/Importers/Actions/ReadVolumeData.hpp"
21 : #include "IO/Importers/ElementDataReader.hpp"
22 : #include "IO/Importers/Tags.hpp"
23 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
24 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
25 : #include "Parallel/AlgorithmExecution.hpp"
26 : #include "Parallel/GlobalCache.hpp"
27 : #include "Parallel/Invoke.hpp"
28 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
29 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
30 : #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
31 : #include "Utilities/CallWithDynamicType.hpp"
32 : #include "Utilities/ErrorHandling/Error.hpp"
33 : #include "Utilities/Gsl.hpp"
34 : #include "Utilities/Serialization/CharmPupable.hpp"
35 : #include "Utilities/TMPL.hpp"
36 : #include "Utilities/TaggedTuple.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 :
69 : /*!
70 : * \brief Numeric initial data loaded from volume data files
71 : *
72 : * This class can be factory-created in the input file to start an evolution
73 : * from numeric initial data. It selects the set of variables to load from
74 : * the volume data file (ADM or GH variables).
75 : */
76 1 : class NumericInitialData : public evolution::initial_data::InitialData {
77 : public:
78 : /// Name of a variable in the volume data file
79 : template <typename Tag>
80 1 : struct VarName {
81 0 : using tag = Tag;
82 0 : static std::string name() { return db::tag_name<Tag>(); }
83 0 : using type = std::string;
84 0 : static constexpr Options::String help =
85 : "Name of the variable in the volume data file";
86 : };
87 :
88 : // These are the sets of variables that we support loading from volume data
89 : // files:
90 : // - ADM variables
91 0 : using adm_vars =
92 : tmpl::list<gr::Tags::SpatialMetric<DataVector, 3>,
93 : gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, 3>,
94 : gr::Tags::ExtrinsicCurvature<DataVector, 3>>;
95 0 : struct AdmVars : tuples::tagged_tuple_from_typelist<
96 : db::wrap_tags_in<VarName, adm_vars>> {
97 0 : static constexpr Options::String help =
98 : "ADM variables: 'Lapse', 'Shift', 'SpatialMetric' and "
99 : "'ExtrinsicCurvature'. The initial GH variables will be computed "
100 : "from these numeric fields, as well as their numeric spatial "
101 : "derivatives on the computational grid. The GH variable Pi will be set "
102 : "to satisfy the gauge constraint using the evolution gauge. The GH "
103 : "variable Phi will be set to satisfy the 3-index constraint.";
104 0 : using options = tags_list;
105 : using TaggedTuple::TaggedTuple;
106 : };
107 :
108 : // - Generalized harmonic variables
109 0 : using gh_vars = tmpl::list<gr::Tags::SpacetimeMetric<DataVector, 3>,
110 : Tags::Pi<DataVector, 3>, Tags::Phi<DataVector, 3>>;
111 0 : struct GhVars
112 : : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<VarName, gh_vars>> {
113 0 : static constexpr Options::String help =
114 : "GH variables: 'SpacetimeMetric', 'Pi', and 'Phi'. These variables are "
115 : "used to set the initial data directly.";
116 0 : using options = tags_list;
117 : using TaggedTuple::TaggedTuple;
118 : };
119 :
120 : // Collect all variables that we support loading from volume data files.
121 : // Remember to `tmpl::remove_duplicates` when adding overlapping sets of
122 : // vars.
123 0 : using all_vars = tmpl::append<adm_vars, gh_vars>;
124 :
125 : // Input-file options
126 0 : struct Variables {
127 : // The user can supply any of these choices of variables in the input
128 : // file
129 0 : using type = std::variant<AdmVars, GhVars>;
130 0 : static constexpr Options::String help =
131 : "Set of initial data variables from which the generalized harmonic "
132 : "system variables are computed.";
133 : };
134 :
135 0 : using options =
136 : tmpl::push_back<importers::ImporterOptions::tags_list, 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(
160 : std::string file_glob, std::string subfile_name,
161 : std::variant<double, importers::ObservationSelector> observation_value,
162 : std::optional<double> observation_value_epsilon,
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 : *phi = get<Tags::Phi<DataVector, 3>>(*numeric_data);
238 : } else if (std::holds_alternative<NumericInitialData::AdmVars>(
239 : selected_variables_)) {
240 : // We have loaded ADM variables from the file. Convert to GH variables.
241 : const auto& spatial_metric =
242 : get<gr::Tags::SpatialMetric<DataVector, 3>>(*numeric_data);
243 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(*numeric_data);
244 : const auto& shift = get<gr::Tags::Shift<DataVector, 3>>(*numeric_data);
245 : const auto& extrinsic_curvature =
246 : get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(*numeric_data);
247 :
248 : initial_gh_variables_from_adm(spacetime_metric, pi, phi, spatial_metric,
249 : lapse, shift, extrinsic_curvature, mesh,
250 : inv_jacobian);
251 : } else {
252 : ERROR(
253 : "These initial data variables are not implemented yet. Please add "
254 : "an implementation to gh::NumericInitialData.");
255 : }
256 : }
257 :
258 0 : void pup(PUP::er& p) override;
259 :
260 0 : friend bool operator==(const NumericInitialData& lhs,
261 : const NumericInitialData& rhs);
262 :
263 : private:
264 0 : importers::ImporterOptions importer_options_;
265 0 : std::variant<AdmVars, GhVars> selected_variables_{};
266 : };
267 :
268 0 : namespace Actions {
269 :
270 : /*!
271 : * \brief Dispatch loading numeric initial data from files or set analytic
272 : * initial data.
273 : *
274 : * Place this action before
275 : * gh::Actions::ReceiveNumericInitialData in the action list.
276 : * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
277 : * invoked by this action.
278 : * Analytic initial data is set directly by this action and terminates the
279 : * phase.
280 : */
281 1 : struct SetInitialData {
282 0 : using const_global_cache_tags =
283 : tmpl::list<evolution::initial_data::Tags::InitialData>;
284 :
285 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
286 : typename ArrayIndex, typename ActionList,
287 : typename ParallelComponent>
288 0 : static Parallel::iterable_action_return_t apply(
289 : db::DataBox<DbTagsList>& box,
290 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
291 : Parallel::GlobalCache<Metavariables>& cache,
292 : const ArrayIndex& array_index, const ActionList /*meta*/,
293 : const ParallelComponent* const parallel_component) {
294 : // Dispatch to the correct `apply` overload based on type of initial data
295 : using initial_data_classes =
296 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
297 : evolution::initial_data::InitialData>;
298 : return call_with_dynamic_type<Parallel::iterable_action_return_t,
299 : initial_data_classes>(
300 : &db::get<evolution::initial_data::Tags::InitialData>(box),
301 : [&box, &cache, &array_index,
302 : ¶llel_component](const auto* const initial_data) {
303 : return apply(make_not_null(&box), *initial_data, cache, array_index,
304 : parallel_component);
305 : });
306 : }
307 :
308 : private:
309 : // Numeric initial data
310 : template <typename DbTagsList, typename Metavariables, typename ArrayIndex,
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 ArrayIndex& array_index, const ParallelComponent* const /*meta*/) {
317 : // If we are using GH Numeric ID, then we don't have to set Pi and Phi since
318 : // we are reading them in. Also we only need to mutate this tag once so do
319 : // it on the first element.
320 : if (is_zeroth_element(array_index) and
321 : std::holds_alternative<NumericInitialData::GhVars>(
322 : initial_data.selected_variables())) {
323 : Parallel::mutate<Tags::SetPiAndPhiFromConstraints,
324 : gh::gauges::SetPiAndPhiFromConstraintsCacheMutator>(
325 : cache, false);
326 : }
327 :
328 : // Select the subset of the available variables that we want to read from
329 : // the volume data file
330 : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
331 : importers::Tags::Selected, NumericInitialData::all_vars>>
332 : selected_fields{};
333 : initial_data.select_for_import(make_not_null(&selected_fields));
334 : // Dispatch loading the variables from the volume data file
335 : // - Not using `ckLocalBranch` here to make sure the simple action
336 : // invocation is asynchronous.
337 : auto& reader_component = Parallel::get_parallel_component<
338 : importers::ElementDataReader<Metavariables>>(cache);
339 : Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
340 : 3, NumericInitialData::all_vars, ParallelComponent>>(
341 : reader_component, initial_data.importer_options(),
342 : initial_data.volume_data_id(), std::move(selected_fields));
343 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
344 : }
345 :
346 : // "AnalyticData"-type initial data
347 : template <typename DbTagsList, typename InitialData, typename Metavariables,
348 : typename ArrayIndex, typename ParallelComponent>
349 0 : static Parallel::iterable_action_return_t apply(
350 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
351 : const InitialData& initial_data,
352 : Parallel::GlobalCache<Metavariables>& /*cache*/,
353 : const ArrayIndex& /*array_index*/,
354 : const ParallelComponent* const /*meta*/) {
355 : static constexpr size_t Dim = Metavariables::volume_dim;
356 :
357 : // Get ADM variables from analytic data / solution
358 : const auto& x =
359 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box);
360 : const auto adm_vars = evolution::Initialization::initial_data(
361 : initial_data, x, db::get<::Tags::Time>(*box),
362 : tmpl::list<gr::Tags::SpatialMetric<DataVector, Dim>,
363 : gr::Tags::Lapse<DataVector>,
364 : gr::Tags::Shift<DataVector, Dim>,
365 : gr::Tags::ExtrinsicCurvature<DataVector, Dim>>{});
366 : const auto& spatial_metric =
367 : get<gr::Tags::SpatialMetric<DataVector, Dim>>(adm_vars);
368 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(adm_vars);
369 : const auto& shift = get<gr::Tags::Shift<DataVector, Dim>>(adm_vars);
370 : const auto& extrinsic_curvature =
371 : get<gr::Tags::ExtrinsicCurvature<DataVector, Dim>>(adm_vars);
372 :
373 : // Compute GH vars from ADM vars
374 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(*box);
375 : const auto& inv_jacobian =
376 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
377 : Frame::Inertial>>(*box);
378 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, Dim>,
379 : Tags::Pi<DataVector, Dim>, Tags::Phi<DataVector, Dim>>(
380 : &gh::initial_gh_variables_from_adm<Dim>, box, spatial_metric, lapse,
381 : shift, extrinsic_curvature, mesh, inv_jacobian);
382 :
383 : // No need to import numeric initial data, so we terminate the phase by
384 : // pausing the algorithm on this element
385 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
386 : }
387 : };
388 :
389 : /*!
390 : * \brief Receive numeric initial data loaded by gh::Actions::SetInitialData.
391 : *
392 : * Place this action in the action list after
393 : * gh::Actions::SetInitialData to wait until the data
394 : * for this element has arrived, and then transform the data to GH variables and
395 : * store it in the DataBox to be used as initial data.
396 : *
397 : * This action modifies the following tags in the DataBox:
398 : * - gr::Tags::SpacetimeMetric<DataVector, 3>
399 : * - gh::Tags::Pi<DataVector, 3>
400 : * - gh::Tags::Phi<DataVector, 3>
401 : */
402 1 : struct ReceiveNumericInitialData {
403 0 : static constexpr size_t Dim = 3;
404 0 : using inbox_tags =
405 : tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
406 :
407 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
408 : typename ActionList, typename ParallelComponent>
409 0 : static Parallel::iterable_action_return_t apply(
410 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
411 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
412 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
413 : const ParallelComponent* const /*meta*/) {
414 : auto& inbox =
415 : tuples::get<importers::Tags::VolumeData<NumericInitialData::all_vars>>(
416 : inboxes);
417 : const auto& initial_data = dynamic_cast<const NumericInitialData&>(
418 : db::get<evolution::initial_data::Tags::InitialData>(box));
419 : const auto& volume_data_id = initial_data.volume_data_id();
420 : if (inbox.find(volume_data_id) == inbox.end()) {
421 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
422 : }
423 : auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
424 :
425 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
426 : const auto& inv_jacobian =
427 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
428 : Frame::Inertial>>(box);
429 :
430 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
431 : Tags::Pi<DataVector, 3>, Tags::Phi<DataVector, 3>>(
432 : [&initial_data, &numeric_data, &mesh, &inv_jacobian](
433 : const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
434 : const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
435 : const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi) {
436 : initial_data.set_initial_data(spacetime_metric, pi, phi,
437 : make_not_null(&numeric_data), mesh,
438 : inv_jacobian);
439 : },
440 : make_not_null(&box));
441 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
442 : }
443 : };
444 :
445 : } // namespace Actions
446 : } // namespace gh
|
