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/TaggedTuple.hpp"
14 : #include "DataStructures/Tensor/EagerMath/DotProduct.hpp"
15 : #include "DataStructures/Tensor/EagerMath/RaiseOrLowerIndex.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/Systems/CurvedScalarWave/System.hpp"
21 : #include "Evolution/Systems/CurvedScalarWave/Tags.hpp"
22 : #include "Evolution/Systems/ScalarWave/System.hpp"
23 : #include "IO/Importers/Actions/ReadVolumeData.hpp"
24 : #include "IO/Importers/ElementDataReader.hpp"
25 : #include "IO/Importers/Tags.hpp"
26 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
27 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
28 : #include "Options/String.hpp"
29 : #include "Parallel/AlgorithmExecution.hpp"
30 : #include "Parallel/GlobalCache.hpp"
31 : #include "Parallel/Invoke.hpp"
32 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
33 : #include "PointwiseFunctions/InitialDataUtilities/InitialData.hpp"
34 : #include "PointwiseFunctions/InitialDataUtilities/Tags/InitialData.hpp"
35 : #include "Utilities/CallWithDynamicType.hpp"
36 : #include "Utilities/ErrorHandling/Error.hpp"
37 : #include "Utilities/Gsl.hpp"
38 : #include "Utilities/Serialization/CharmPupable.hpp"
39 : #include "Utilities/SetNumberOfGridPoints.hpp"
40 : #include "Utilities/TMPL.hpp"
41 :
42 : namespace CurvedScalarWave {
43 :
44 : /*!
45 : * \brief Numeric initial data loaded from volume data files
46 : */
47 1 : class NumericInitialData : public evolution::initial_data::InitialData {
48 : public:
49 : template <typename Tag>
50 0 : struct VarName {
51 0 : using tag = Tag;
52 0 : static std::string name() { return db::tag_name<Tag>(); }
53 0 : using type = std::string;
54 0 : static constexpr Options::String help =
55 : "Name of the variable in the volume data file";
56 : };
57 :
58 : // These are the scalar variables that we support loading from volume
59 : // data files
60 0 : using all_vars =
61 : tmpl::list<CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi,
62 : CurvedScalarWave::Tags::Phi<3>>;
63 0 : using optional_primitive_vars = tmpl::list<>;
64 :
65 0 : struct ScalarVars : tuples::tagged_tuple_from_typelist<
66 : db::wrap_tags_in<VarName, all_vars>> {
67 0 : static constexpr Options::String help =
68 : "Scalar variables: 'Psi', 'Pi' and 'Phi'.";
69 0 : using options = tags_list;
70 : using TaggedTuple::TaggedTuple;
71 : };
72 :
73 : // Input-file options
74 0 : struct Variables {
75 0 : using type = ScalarVars;
76 0 : static constexpr Options::String help =
77 : "Set of initial data variables for the CurvedScalarWave system.";
78 : };
79 :
80 0 : using options = tmpl::list<importers::OptionTags::VolumeData, Variables>;
81 :
82 0 : static constexpr Options::String help =
83 : "Numeric initial data loaded from volume data files";
84 :
85 0 : NumericInitialData() = default;
86 0 : NumericInitialData(const NumericInitialData& rhs) = default;
87 0 : NumericInitialData& operator=(const NumericInitialData& rhs) = default;
88 0 : NumericInitialData(NumericInitialData&& /*rhs*/) = default;
89 0 : NumericInitialData& operator=(NumericInitialData&& /*rhs*/) = default;
90 0 : ~NumericInitialData() override = default;
91 :
92 : /// \cond
93 : explicit NumericInitialData(CkMigrateMessage* msg);
94 : using PUP::able::register_constructor;
95 : WRAPPED_PUPable_decl_template(NumericInitialData);
96 : /// \endcond
97 :
98 0 : std::unique_ptr<evolution::initial_data::InitialData> get_clone()
99 : const override {
100 : return std::make_unique<NumericInitialData>(*this);
101 : }
102 :
103 0 : NumericInitialData(importers::ImporterOptions importer_options,
104 : ScalarVars selected_variables);
105 :
106 0 : const importers::ImporterOptions& importer_options() const;
107 :
108 0 : const ScalarVars& selected_variables() const { return selected_variables_; }
109 :
110 0 : size_t volume_data_id() const;
111 :
112 : template <typename... AllTags>
113 0 : void select_for_import(
114 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*> fields) const {
115 : // Select the subset of the available variables that we want to read from
116 : // the volume data file
117 : using selected_vars = std::decay_t<decltype(selected_variables_)>;
118 : tmpl::for_each<typename selected_vars::tags_list>(
119 : [&fields, this](const auto tag_v) {
120 : using tag = typename std::decay_t<decltype(tag_v)>::type::tag;
121 : get<importers::Tags::Selected<tag>>(*fields) =
122 : get<VarName<tag>>(selected_variables_);
123 : });
124 : }
125 :
126 : template <typename... AllTags>
127 0 : void set_initial_data(const gsl::not_null<Scalar<DataVector>*> psi_scalar,
128 : const gsl::not_null<Scalar<DataVector>*> pi_scalar,
129 : const gsl::not_null<tnsr::i<DataVector, 3>*> phi_scalar,
130 : const gsl::not_null<tuples::TaggedTuple<AllTags...>*>
131 : numeric_data) const {
132 : *psi_scalar = std::move(get<CurvedScalarWave::Tags::Psi>(*numeric_data));
133 : *pi_scalar = std::move(get<CurvedScalarWave::Tags::Pi>(*numeric_data));
134 : *phi_scalar = std::move(get<CurvedScalarWave::Tags::Phi<3>>(*numeric_data));
135 : }
136 :
137 0 : void pup(PUP::er& p) override;
138 :
139 0 : friend bool operator==(const NumericInitialData& lhs,
140 : const NumericInitialData& rhs);
141 :
142 : private:
143 0 : importers::ImporterOptions importer_options_{};
144 0 : ScalarVars selected_variables_{};
145 : };
146 :
147 : namespace Actions {
148 :
149 : /*!
150 : * \brief Dispatch loading numeric initial data from files.
151 : *
152 : * Place this action before
153 : * CurvedScalarWave::Actions::SetNumericInitialData in the action list.
154 : * See importers::Actions::ReadAllVolumeDataAndDistribute for details, which is
155 : * invoked by this action.
156 : */
157 1 : struct SetInitialData {
158 0 : using const_global_cache_tags =
159 : tmpl::list<evolution::initial_data::Tags::InitialData>;
160 :
161 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
162 : typename ArrayIndex, typename ActionList,
163 : typename ParallelComponent>
164 0 : static Parallel::iterable_action_return_t apply(
165 : db::DataBox<DbTagsList>& box,
166 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
167 : Parallel::GlobalCache<Metavariables>& cache,
168 : const ArrayIndex& array_index, const ActionList /*meta*/,
169 : const ParallelComponent* const parallel_component) {
170 : // Dispatch to the correct `apply` overload based on type of initial data
171 : using initial_data_classes =
172 : tmpl::at<typename Metavariables::factory_creation::factory_classes,
173 : evolution::initial_data::InitialData>;
174 : return call_with_dynamic_type<Parallel::iterable_action_return_t,
175 : initial_data_classes>(
176 : &db::get<evolution::initial_data::Tags::InitialData>(box),
177 : [&box, &cache, &array_index,
178 : ¶llel_component](const auto* const initial_data) {
179 : return apply(make_not_null(&box), *initial_data, cache, array_index,
180 : parallel_component);
181 : });
182 : }
183 :
184 : private:
185 : // Numeric initial data
186 : template <typename DbTagsList, typename Metavariables, typename ArrayIndex,
187 : typename ParallelComponent>
188 0 : static Parallel::iterable_action_return_t apply(
189 : const gsl::not_null<db::DataBox<DbTagsList>*> /*box*/,
190 : const NumericInitialData& initial_data,
191 : Parallel::GlobalCache<Metavariables>& cache,
192 : const ArrayIndex& /*array_index*/,
193 : const ParallelComponent* const /*meta*/) {
194 : // Select the subset of the available variables that we want to read from
195 : // the volume data file
196 : tuples::tagged_tuple_from_typelist<db::wrap_tags_in<
197 : importers::Tags::Selected, NumericInitialData::all_vars>>
198 : selected_fields{};
199 : initial_data.select_for_import(make_not_null(&selected_fields));
200 : auto& reader_component = Parallel::get_parallel_component<
201 : importers::ElementDataReader<Metavariables>>(cache);
202 : Parallel::simple_action<importers::Actions::ReadAllVolumeDataAndDistribute<
203 : 3, NumericInitialData::all_vars, ParallelComponent>>(
204 : reader_component, initial_data.importer_options(),
205 : initial_data.volume_data_id(), std::move(selected_fields));
206 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
207 : }
208 :
209 : // "AnalyticData"-type initial data
210 : template <typename DbTagsList, typename InitialData, typename Metavariables,
211 : typename ArrayIndex, typename ParallelComponent>
212 0 : static Parallel::iterable_action_return_t apply(
213 : const gsl::not_null<db::DataBox<DbTagsList>*> box,
214 : const InitialData& initial_data,
215 : Parallel::GlobalCache<Metavariables>& /*cache*/,
216 : const ArrayIndex& /*array_index*/,
217 : const ParallelComponent* const /*meta*/) {
218 : static constexpr size_t Dim = Metavariables::volume_dim;
219 : using flat_variables_tag = typename ScalarWave::System<Dim>::variables_tag;
220 : using curved_variables_tag =
221 : typename CurvedScalarWave::System<Dim>::variables_tag;
222 : const auto inertial_coords =
223 : db::get<domain::Tags::Coordinates<Dim, Frame::Inertial>>(*box);
224 : const double initial_time = db::get<::Tags::Time>(*box);
225 : if constexpr (is_analytic_data_v<InitialData> or
226 : is_analytic_solution_v<InitialData>) {
227 : if constexpr (tmpl::list_contains_v<typename InitialData::tags,
228 : CurvedScalarWave::Tags::Psi>) {
229 : const auto curved_initial_data =
230 : evolution::Initialization::initial_data(
231 : initial_data, inertial_coords, initial_time,
232 : typename curved_variables_tag::tags_list{});
233 :
234 : db::mutate<typename CurvedScalarWave::System<Dim>::variables_tag>(
235 : [&curved_initial_data](
236 : const gsl::not_null<typename curved_variables_tag::type*>
237 : evolved_vars) {
238 : evolved_vars->assign_subset(curved_initial_data);
239 : },
240 : box);
241 : } else {
242 : static_assert(tmpl::list_contains_v<typename InitialData::tags,
243 : ScalarWave::Tags::Psi>,
244 : "The initial data class must either calculate ScalarWave "
245 : "or CurvedScalarWave variables.");
246 : const auto flat_initial_data = evolution::Initialization::initial_data(
247 : initial_data, inertial_coords, initial_time,
248 : typename flat_variables_tag::tags_list{});
249 : const auto& shift = db::get<gr::Tags::Shift<DataVector, Dim>>(*box);
250 : const auto& lapse = db::get<gr::Tags::Lapse<DataVector>>(*box);
251 : const auto shift_dot_dpsi = dot_product(
252 : shift, get<ScalarWave::Tags::Phi<Dim>>(flat_initial_data));
253 : db::mutate<typename CurvedScalarWave::System<Dim>::variables_tag>(
254 : [&flat_initial_data, &shift_dot_dpsi,
255 : &lapse](const gsl::not_null<typename curved_variables_tag::type*>
256 : evolved_vars) {
257 : get<CurvedScalarWave::Tags::Psi>(*evolved_vars) =
258 : get<ScalarWave::Tags::Psi>(flat_initial_data);
259 : get<CurvedScalarWave::Tags::Phi<Dim>>(*evolved_vars) =
260 : get<ScalarWave::Tags::Phi<Dim>>(flat_initial_data);
261 : get(get<CurvedScalarWave::Tags::Pi>(*evolved_vars)) =
262 : (get(shift_dot_dpsi) +
263 : get(get<ScalarWave::Tags::Pi>(flat_initial_data))) /
264 : get(lapse);
265 : },
266 : box);
267 : }
268 : } else {
269 : ERROR(
270 : "Trying to use "
271 : "'evolution::Initialization::Actions::SetVariables' with a "
272 : "class that's not marked as analytic solution or analytic "
273 : "data. To support numeric initial data, add a "
274 : "system-specific initialization routine to your executable.");
275 : }
276 :
277 : return {Parallel::AlgorithmExecution::Pause, std::nullopt};
278 : }
279 : };
280 :
281 : /*!
282 : * \brief Receive numeric initial data loaded by
283 : * CurvedScalarWave::Actions::ReadNumericInitialData.
284 : */
285 1 : struct ReceiveNumericInitialData {
286 0 : static constexpr size_t Dim = 3;
287 0 : using inbox_tags =
288 : tmpl::list<importers::Tags::VolumeData<NumericInitialData::all_vars>>;
289 :
290 : template <typename DbTagsList, typename... InboxTags, typename Metavariables,
291 : typename ArrayIndex, typename ActionList,
292 : typename ParallelComponent>
293 0 : static Parallel::iterable_action_return_t apply(
294 : db::DataBox<DbTagsList>& box, tuples::TaggedTuple<InboxTags...>& inboxes,
295 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
296 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
297 : const ParallelComponent* const /*meta*/) {
298 : if constexpr (Metavariables::volume_dim != Dim) {
299 : ERROR(
300 : "CurvedScalarWave numeric initial data currently requires a 3D "
301 : "domain.");
302 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
303 : } else {
304 : auto& inbox = tuples::get<
305 : importers::Tags::VolumeData<NumericInitialData::all_vars>>(inboxes);
306 : const auto& initial_data = dynamic_cast<const NumericInitialData&>(
307 : db::get<evolution::initial_data::Tags::InitialData>(box));
308 : const size_t volume_data_id = initial_data.volume_data_id();
309 : if (inbox.find(volume_data_id) == inbox.end()) {
310 : return {Parallel::AlgorithmExecution::Retry, std::nullopt};
311 : }
312 : auto numeric_data = std::move(inbox.extract(volume_data_id).mapped());
313 :
314 : db::mutate<CurvedScalarWave::Tags::Psi, CurvedScalarWave::Tags::Pi,
315 : CurvedScalarWave::Tags::Phi<Dim>>(
316 : [&initial_data, &numeric_data](
317 : const gsl::not_null<Scalar<DataVector>*> psi_scalar,
318 : const gsl::not_null<Scalar<DataVector>*> pi_scalar,
319 : const gsl::not_null<tnsr::i<DataVector, Dim>*> phi_scalar) {
320 : initial_data.set_initial_data(psi_scalar, pi_scalar, phi_scalar,
321 : make_not_null(&numeric_data));
322 : },
323 : make_not_null(&box));
324 :
325 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
326 : }
327 : }
328 : };
329 :
330 : } // namespace Actions
331 :
332 : } // namespace CurvedScalarWave
|
