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 <pup.h>
9 : #include <string>
10 : #include <variant>
11 :
12 : #include "DataStructures/DataBox/DataBox.hpp"
13 : #include "DataStructures/DataVector.hpp"
14 : #include "DataStructures/TaggedTuple.hpp"
15 : #include "DataStructures/Tensor/EagerMath/DeterminantAndInverse.hpp"
16 : #include "DataStructures/Tensor/Tensor.hpp"
17 : #include "Domain/Structure/ElementId.hpp"
18 : #include "Domain/Tags.hpp"
19 : #include "Evolution/Initialization/InitialData.hpp"
20 : #include "Evolution/NumericInitialData.hpp"
21 : #include "Evolution/Systems/CurvedScalarWave/Actions/SetInitialData.hpp"
22 : #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
23 : #include "Evolution/Systems/GeneralizedHarmonic/Actions/SetInitialData.hpp"
24 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
25 : #include "IO/Importers/Actions/ReadVolumeData.hpp"
26 : #include "IO/Importers/ElementDataReader.hpp"
27 : #include "IO/Importers/Tags.hpp"
28 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
29 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
30 : #include "Options/String.hpp"
31 : #include "Parallel/AlgorithmExecution.hpp"
32 : #include "Parallel/GlobalCache.hpp"
33 : #include "Parallel/Invoke.hpp"
34 : #include "PointwiseFunctions/GeneralRelativity/SpatialMetric.hpp"
35 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
36 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
37 : #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
38 : #include "Utilities/CallWithDynamicType.hpp"
39 : #include "Utilities/ErrorHandling/Error.hpp"
40 : #include "Utilities/Gsl.hpp"
41 : #include "Utilities/Serialization/CharmPupable.hpp"
42 : #include "Utilities/TMPL.hpp"
43 :
44 1 : namespace ScalarTensor {
45 :
46 : /*!
47 : * \brief Numeric initial data loaded from volume data files
48 : */
49 1 : class NumericInitialData : public evolution::initial_data::InitialData,
50 : public evolution::NumericInitialData {
51 : private:
52 0 : using GhNumericId = gh::NumericInitialData;
53 0 : using ScalarNumericId = CurvedScalarWave::NumericInitialData;
54 :
55 : public:
56 0 : using all_vars =
57 : tmpl::append<GhNumericId::all_vars, ScalarNumericId::all_vars>;
58 :
59 0 : struct GhVariables : GhNumericId::Variables {};
60 0 : struct ScalarVariables : ScalarNumericId::Variables {};
61 :
62 0 : using options = tmpl::push_back<importers::ImporterOptions::tags_list,
63 : GhVariables, ScalarVariables>;
64 :
65 0 : static constexpr Options::String help =
66 : "Numeric initial data for the Scalar Tensor system loaded from volume "
67 : "data files";
68 :
69 0 : NumericInitialData() = default;
70 0 : NumericInitialData(const NumericInitialData& rhs) = default;
71 0 : NumericInitialData& operator=(const NumericInitialData& rhs) = default;
72 0 : NumericInitialData(NumericInitialData&& /*rhs*/) = default;
73 0 : NumericInitialData& operator=(NumericInitialData&& /*rhs*/) = default;
74 0 : ~NumericInitialData() override = default;
75 :
76 : /// \cond
77 : explicit NumericInitialData(CkMigrateMessage* msg);
78 : using PUP::able::register_constructor;
79 : WRAPPED_PUPable_decl_template(NumericInitialData);
80 : /// \endcond
81 :
82 0 : std::unique_ptr<evolution::initial_data::InitialData> get_clone()
83 : const override {
84 : return std::make_unique<NumericInitialData>(*this);
85 : }
86 :
87 0 : NumericInitialData(
88 : std::string file_glob, std::string subfile_name,
89 : std::variant<double, importers::ObservationSelector> observation_value,
90 : std::optional<double> observation_value_epsilon,
91 : bool enable_interpolation,
92 : typename GhNumericId::Variables::type gh_selected_variables,
93 : typename ScalarNumericId::Variables::type hydro_selected_variables);
94 :
95 0 : const importers::ImporterOptions& importer_options() const {
96 : return gh_numeric_id_.importer_options();
97 : }
98 :
99 0 : const GhNumericId& gh_numeric_id() const { return gh_numeric_id_; }
100 :
101 0 : const ScalarNumericId& scalar_numeric_id() const {
102 : return scalar_numeric_id_;
103 : }
104 :
105 0 : size_t volume_data_id() const;
106 :
107 : template <typename... AllTags>
108 0 : void select_for_import(
109 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> fields) const {
110 : gh_numeric_id_.select_for_import(fields);
111 : scalar_numeric_id_.select_for_import(fields);
112 : }
113 :
114 : template <typename... AllTags>
115 0 : void set_initial_data(
116 : const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
117 : const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
118 : const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
119 : const gsl::not_null<Scalar<DataVector>*> psi_scalar,
120 : const gsl::not_null<Scalar<DataVector>*> pi_scalar,
121 : const gsl::not_null<tnsr::i<DataVector, 3>*> phi_scalar,
122 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> numeric_data,
123 : const Mesh<3>& mesh,
124 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
125 : Frame::Inertial>& inv_jacobian,
126 : const tnsr::I<DataVector, 3, Frame::Inertial>& inertial_coords) const {
127 : gh_numeric_id_.set_initial_data(spacetime_metric, pi, phi, numeric_data,
128 : mesh, inv_jacobian, inertial_coords);
129 : scalar_numeric_id_.set_initial_data(psi_scalar, pi_scalar, phi_scalar,
130 : numeric_data);
131 : }
132 :
133 0 : void pup(PUP::er& p) override;
134 :
135 0 : friend bool operator==(const NumericInitialData& lhs,
136 : const NumericInitialData& rhs);
137 :
138 : private:
139 0 : GhNumericId gh_numeric_id_{};
140 0 : ScalarNumericId scalar_numeric_id_{};
141 : };
142 :
143 0 : namespace Actions {
144 :
145 : /*!
146 : * \brief Dispatch loading numeric initial data from files.
147 : *
148 : * Place this action before
149 : * ScalarTensor::Actions::SetNumericInitialData in the action list.
150 : * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
151 : * invoked by this action.
152 : */
153 1 : struct SetInitialData {
154 0 : using const_global_cache_tags =
155 : tmpl::list<evolution::initial_data::Tags::InitialData>;
156 :
157 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
158 : typename ArrayIndex, typename ActionList,
159 : typename ParallelComponent>
160 0 : static Parallel::iterable_action_return_t apply(
161 : db::DataBox<DbTagsList>& box,
162 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
163 : Parallel::GlobalCache<Metavariables>& cache,
164 : const ArrayIndex& array_index, const ActionList /*meta*/,
165 : const ParallelComponent* const parallel_component) {
166 : // Dispatch to the correct `apply` overload based on type of initial data
167 : using initial_data_classes =
168 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
169 : evolution::initial_data::InitialData>;
170 : return call_with_dynamic_type<Parallel::iterable_action_return_t,
171 : initial_data_classes>(
172 : &db::get<evolution::initial_data::Tags::InitialData>(box),
173 : [&box, &cache, &array_index,
174 : ¶llel_component](const auto* const initial_data) {
175 : return apply(make_not_null(&box), *initial_data, cache, array_index,
176 : parallel_component);
177 : });
178 : }
179 :
180 : private:
181 0 : static constexpr size_t Dim = 3;
182 :
183 : // Numeric initial data
184 : template <typename DbTagsList, typename Metavariables, typename ArrayIndex,
185 : typename ParallelComponent>
186 0 : static Parallel::iterable_action_return_t apply(
187 : const gsl::not_null<db::DataBox<DbTagsList>*> /*box*/,
188 : const NumericInitialData& initial_data,
189 : Parallel::GlobalCache<Metavariables>& cache,
190 : const ArrayIndex& array_index, const ParallelComponent* const /*meta*/) {
191 : // If we are using GH Numeric ID, then we don't have to set Pi and Phi since
192 : // we are reading them in. Also we only need to mutate this tag once so do
193 : // it on the first element.
194 : if (is_zeroth_element(array_index) and
195 : std::holds_alternative<gh::NumericInitialData::GhVars>(
196 : initial_data.gh_numeric_id().selected_variables())) {
197 : Parallel::mutate<gh::Tags::SetPiAndPhiFromConstraints,
198 : gh::gauges::SetPiAndPhiFromConstraintsCacheMutator>(
199 : cache, false);
200 : }
201 : // Select the subset of the available variables that we want to read from
202 : // the volume data file
203 : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
204 : importers::Tags::Selected, NumericInitialData::all_vars>>
205 : selected_fields{};
206 : initial_data.select_for_import(make_not_null(&selected_fields));
207 : // Dispatch loading the variables from the volume data file
208 : // - Not using `ckLocalBranch` here to make sure the simple action
209 : // invocation is asynchronous.
210 : auto& reader_component = Parallel::get_parallel_component<
211 : importers::ElementDataReader<Metavariables>>(cache);
212 : Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
213 : 3, NumericInitialData::all_vars, ParallelComponent>>(
214 : reader_component, initial_data.importer_options(),
215 : initial_data.volume_data_id(), std::move(selected_fields));
216 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
217 : }
218 :
219 : // "AnalyticData"-type initial data
220 : template <typename DbTagsList, typename InitialData, typename Metavariables,
221 : typename ArrayIndex, typename ParallelComponent>
222 0 : static Parallel::iterable_action_return_t apply(
223 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
224 : const InitialData& initial_data,
225 : Parallel::GlobalCache<Metavariables>& /*cache*/,
226 : const ArrayIndex& /*array_index*/,
227 : const ParallelComponent* const /*meta*/) {
228 : // Get ADM + scalar variables from analytic data / solution
229 : const auto& [coords, mesh, inv_jacobian] = [&box]() {
230 : return std::forward_as_tuple(
231 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box),
232 : db::get<domain::Tags::Mesh<Dim>>(*box),
233 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
234 : Frame::Inertial>>(*box));
235 : }();
236 : auto vars = evolution::Initialization::initial_data(
237 : initial_data, coords, db::get<::Tags::Time>(*box),
238 : tmpl::append<tmpl::list<gr::Tags::SpatialMetric<DataVector, 3>,
239 : gr::Tags::Lapse<DataVector>,
240 : gr::Tags::Shift<DataVector, 3>,
241 : gr::Tags::ExtrinsicCurvature<DataVector, 3>>,
242 : // Don't use the scalar gradient
243 : tmpl::list<CurvedScalarWave::Tags::Psi,
244 : CurvedScalarWave::Tags::Pi>>{});
245 : const auto& spatial_metric =
246 : get<gr::Tags::SpatialMetric<DataVector, 3>>(vars);
247 : const auto& lapse = get<gr::Tags::Lapse<DataVector>>(vars);
248 :
249 : const auto& shift = get<gr::Tags::Shift<DataVector, 3>>(vars);
250 :
251 : const auto& extrinsic_curvature =
252 : get<gr::Tags::ExtrinsicCurvature<DataVector, 3>>(vars);
253 :
254 : // Compute GH vars from ADM vars
255 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
256 : gh::Tags::Pi<DataVector, 3>, gh::Tags::Phi<DataVector, 3>>(
257 : &gh::initial_gh_variables_from_adm<3>, box, spatial_metric, lapse,
258 : shift, extrinsic_curvature, mesh, inv_jacobian, coords);
259 :
260 : // Move scalar variables and compute gradient
261 : db::mutate<CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi,
262 : CurvedScalarWave::Tags::Phi<3>>(
263 : [&vars](const gsl::not_null<Scalar<DataVector>*> psi_scalar,
264 : const gsl::not_null<Scalar<DataVector>*> pi_scalar,
265 : const gsl::not_null<tnsr::i<DataVector, 3>*> phi_scalar,
266 : const Mesh<3>& local_mesh,
267 : const InverseJacobian<DataVector, 3_st, Frame::ElementLogical,
268 : Frame::Inertial>& local_inv_jacobian) {
269 : *psi_scalar = std::move(get<CurvedScalarWave::Tags::Psi>(vars));
270 : *pi_scalar = std::move(get<CurvedScalarWave::Tags::Pi>(vars));
271 : // Set Phi to the numerical spatial derivative of the scalar
272 : partial_derivative(phi_scalar, *psi_scalar, local_mesh,
273 : local_inv_jacobian);
274 : },
275 : box, mesh, inv_jacobian);
276 :
277 : // No need to import numeric initial data, so we terminate the phase by
278 : // pausing the algorithm on this element
279 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
280 : }
281 : };
282 :
283 : /*!
284 : * \brief Receive numeric initial data loaded by
285 : * ScalarTensor::Actions::SetInitialData.
286 : */
287 1 : struct ReceiveNumericInitialData {
288 0 : static constexpr size_t Dim = 3;
289 0 : using inbox_tags =
290 : tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
291 :
292 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
293 : typename ActionList, typename ParallelComponent>
294 0 : static Parallel::iterable_action_return_t apply(
295 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
296 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
297 : const ElementId<Dim>& /*element_id*/, const ActionList /*meta*/,
298 : const ParallelComponent* const /*meta*/) {
299 : auto& inbox =
300 : tuples::get<importers::Tags::VolumeData<NumericInitialData::all_vars>>(
301 : inboxes);
302 : const auto& initial_data = dynamic_cast<const NumericInitialData&>(
303 : db::get<evolution::initial_data::Tags::InitialData>(box));
304 : const size_t volume_data_id = initial_data.volume_data_id();
305 : if (inbox.find(volume_data_id) == inbox.end()) {
306 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
307 : }
308 : auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
309 :
310 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
311 : const auto& inv_jacobian =
312 : db::get<domain::Tags::InverseJacobian<Dim, Frame::ElementLogical,
313 : Frame::Inertial>>(box);
314 : const auto& inertial_coords =
315 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(box);
316 :
317 : db::mutate<gr::Tags::SpacetimeMetric<DataVector, 3>,
318 : gh::Tags::Pi<DataVector, 3>, gh::Tags::Phi<DataVector, 3>,
319 : CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi,
320 : CurvedScalarWave::Tags::Phi<3>>(
321 : [&initial_data, &numeric_data, &mesh, &inv_jacobian, &inertial_coords](
322 : const gsl::not_null<tnsr::aa<DataVector, 3>*> spacetime_metric,
323 : const gsl::not_null<tnsr::aa<DataVector, 3>*> pi,
324 : const gsl::not_null<tnsr::iaa<DataVector, 3>*> phi,
325 : const gsl::not_null<Scalar<DataVector>*> psi_scalar,
326 : const gsl::not_null<Scalar<DataVector>*> pi_scalar,
327 : const gsl::not_null<tnsr::i<DataVector, 3>*> phi_scalar) {
328 : initial_data.set_initial_data(spacetime_metric, pi, phi,
329 :
330 : psi_scalar, pi_scalar, phi_scalar,
331 :
332 : make_not_null(&numeric_data), mesh,
333 : inv_jacobian, inertial_coords);
334 : },
335 : make_not_null(&box));
336 :
337 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
338 : }
339 : };
340 :
341 : } // namespace Actions
342 :
343 : } // namespace ScalarTensor
|
