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