Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <optional>
7 :
8 : #include "DataStructures/DataBox/PrefixHelpers.hpp"
9 : #include "DataStructures/DataBox/SubitemTag.hpp"
10 : #include "DataStructures/DataBox/Subitems.hpp"
11 : #include "DataStructures/DataBox/Tag.hpp"
12 : #include "DataStructures/Variables.hpp"
13 : #include "Options/String.hpp"
14 : #include "PointwiseFunctions/AnalyticData/Tags.hpp"
15 : #include "PointwiseFunctions/AnalyticSolutions/AnalyticSolution.hpp"
16 : #include "Utilities/PrettyType.hpp"
17 : #include "Utilities/Serialization/Serialize.hpp"
18 : #include "Utilities/TMPL.hpp"
19 :
20 : namespace OptionTags {
21 : /// \ingroup OptionGroupsGroup
22 : /// Holds the `OptionTags::AnalyticSolution` option in the input file
23 1 : struct AnalyticSolutionGroup {
24 0 : static std::string name() { return "AnalyticSolution"; }
25 0 : static constexpr Options::String help =
26 : "Analytic solution used for the initial data and errors";
27 : };
28 :
29 : /// \ingroup OptionTagsGroup
30 : /// The analytic solution, with the type of the analytic solution set as the
31 : /// template parameter
32 : template <typename SolutionType>
33 1 : struct AnalyticSolution {
34 0 : static std::string name() { return pretty_type::name<SolutionType>(); }
35 0 : static constexpr Options::String help = "Options for the analytic solution";
36 0 : using type = SolutionType;
37 0 : using group = AnalyticSolutionGroup;
38 : };
39 : } // namespace OptionTags
40 :
41 : namespace Tags {
42 : /// \ingroup OptionTagsGroup
43 : /// The analytic solution, with the type of the analytic solution set as the
44 : /// template parameter
45 : template <typename SolutionType>
46 1 : struct AnalyticSolution : db::SimpleTag {
47 0 : using type = SolutionType;
48 0 : using option_tags = tmpl::list<::OptionTags::AnalyticSolution<SolutionType>>;
49 :
50 0 : static constexpr bool pass_metavariables = false;
51 0 : static SolutionType create_from_options(
52 : const SolutionType& analytic_solution) {
53 : return serialize_and_deserialize<type>(analytic_solution);
54 : }
55 : };
56 :
57 : /// \ingroup DataBoxTagsGroup
58 : /// \brief Prefix indicating the analytic solution value for a quantity
59 : ///
60 : /// \snippet AnalyticSolutions/Test_Tags.cpp analytic_name
61 : template <typename Tag>
62 1 : struct Analytic : db::PrefixTag, db::SimpleTag {
63 0 : using type = std::optional<typename Tag::type>;
64 0 : using tag = Tag;
65 : };
66 :
67 : namespace detail {
68 : template <typename Tag>
69 : struct AnalyticImpl : db::PrefixTag, db::SimpleTag {
70 : using type = typename Tag::type;
71 : using tag = Tag;
72 : };
73 : } // namespace detail
74 :
75 : /*!
76 : * \brief The analytic solution of the `FieldTags`.
77 : *
78 : * The `std::optional` is a `nullopt` if there is no analytic solution.
79 : *
80 : * The individual `FieldTags` are added as `Tags::Analytic<Tag>` which holds a
81 : * `std::optional<Tensor>`.
82 : */
83 : template <typename FieldTags>
84 1 : struct AnalyticSolutions : db::SimpleTag {
85 0 : using field_tags = FieldTags;
86 0 : using type = std::optional<
87 : ::Variables<db::wrap_tags_in<detail::AnalyticImpl, FieldTags>>>;
88 : };
89 :
90 : namespace detail {
91 : // Check if the argument is a `::Tags::AnalyticSolutions` template, or derived
92 : // from it
93 : template <typename FieldTags>
94 : constexpr std::true_type is_analytic_solutions(
95 : ::Tags::AnalyticSolutions<FieldTags>&&) {
96 : return {};
97 : }
98 :
99 : constexpr std::false_type is_analytic_solutions(...) { return {}; }
100 :
101 : template <typename Tag>
102 : static constexpr bool is_analytic_solutions_v =
103 : decltype(is_analytic_solutions(std::declval<Tag>()))::value;
104 : } // namespace detail
105 :
106 : /*!
107 : * \brief The error of the `Tag` defined as `numerical - analytic`.
108 : *
109 : * The `std::optional` is a `nullopt` if no error was able to be computed, e.g.
110 : * if there is no analytic solution to compare to.
111 : */
112 : template <typename Tag>
113 1 : struct Error : db::PrefixTag, db::SimpleTag {
114 0 : using type = std::optional<typename Tag::type>;
115 0 : using tag = Tag;
116 : };
117 :
118 : namespace detail {
119 : template <typename Tag>
120 : struct ErrorImpl : db::PrefixTag, db::SimpleTag {
121 : using type = typename Tag::type;
122 : using tag = Tag;
123 : };
124 : } // namespace detail
125 :
126 : /*!
127 : * \brief The error of the `FieldTags`, defined as `numerical - analytic`.
128 : *
129 : * The `std::optional` is a `nullopt` if no error was able to be computed, e.g.
130 : * if there is no analytic solution to compare to.
131 : *
132 : * The individual `FieldTags` are added as `Tags::Error<Tag>` which holds a
133 : * `std::optional<Tensor>`.
134 : */
135 : template <typename FieldTags>
136 1 : struct Errors : db::SimpleTag {
137 0 : using field_tags = FieldTags;
138 0 : using type = std::optional<
139 : ::Variables<db::wrap_tags_in<detail::ErrorImpl, FieldTags>>>;
140 : };
141 :
142 : /// \cond
143 : template <typename FieldTagsList>
144 : struct ErrorsCompute;
145 : /// \endcond
146 :
147 : /*!
148 : * \brief Compute tag for computing the error from the `Tags::Analytic` of the
149 : * `FieldTags`.
150 : *
151 : * The error is defined as `numerical - analytic`.
152 : *
153 : * We use individual `Tensor`s rather than `Variables` of the `FieldTags`
154 : * because not all `FieldTags` are always stored in the same `Variables`.
155 : * For example, in the Generalized Harmonic-GRMHD combined system, the analytic
156 : * variables for the GH system are part of the `evolved_variables_tag` while the
157 : * GRMHD analytic variables are part of the `primitive_variables_tag`. A similar
158 : * issue arises in some elliptic systems. The main drawback of having to use the
159 : * tensor-by-tensor implementation instead of a `Variables` implementation is
160 : * the added loop complexity. However, this is offset by the reduced code
161 : * duplication and flexibility.
162 : */
163 : template <typename... FieldTags>
164 1 : struct ErrorsCompute<tmpl::list<FieldTags...>>
165 : : Errors<tmpl::list<FieldTags...>>, db::ComputeTag {
166 0 : using field_tags = tmpl::list<FieldTags...>;
167 0 : using base = Errors<tmpl::list<FieldTags...>>;
168 0 : using return_type = typename base::type;
169 0 : using argument_tags =
170 : tmpl::list<FieldTags..., ::Tags::Analytic<FieldTags>...>;
171 0 : static void function(
172 : const gsl::not_null<return_type*> errors,
173 : const typename FieldTags::type&... vars,
174 : const typename ::Tags::Analytic<FieldTags>::type&... analytic_vars) {
175 : if (const auto& first_analytic_var = get_first_argument(analytic_vars...);
176 : first_analytic_var.has_value()) {
177 : // Construct a Variables of the size of the DataVector of the first
178 : // analytic_vars tensor.
179 : *errors = typename return_type::value_type{
180 : first_analytic_var.value()[0].size()};
181 : const auto compute_error = [](const auto error_ptr, const auto& var,
182 : const auto& analytic_var) {
183 : for (size_t tensor_index = 0; tensor_index < var.size();
184 : ++tensor_index) {
185 : (*error_ptr)[tensor_index] =
186 : var[tensor_index] - analytic_var[tensor_index];
187 : }
188 : };
189 : EXPAND_PACK_LEFT_TO_RIGHT(compute_error(
190 : make_not_null(
191 : &get<::Tags::detail::ErrorImpl<FieldTags>>(errors->value())),
192 : vars, analytic_vars.value()));
193 : } else {
194 : *errors = std::nullopt;
195 : }
196 : }
197 : };
198 :
199 : namespace detail {
200 : // Check if the argument is a `::Tags::Errors` template, or derived
201 : // from it
202 : template <typename FieldTags>
203 : constexpr std::true_type is_errors(::Tags::Errors<FieldTags>&&) {
204 : return {};
205 : }
206 :
207 : constexpr std::false_type is_errors(...) { return {}; }
208 :
209 : template <typename Tag>
210 : static constexpr bool is_errors_v =
211 : decltype(is_errors(std::declval<Tag>()))::value;
212 : } // namespace detail
213 : } // namespace Tags
214 :
215 : /// \cond
216 : namespace Tags {
217 : template <typename Tag, typename ParentTag>
218 : struct Subitem<Tag, ParentTag,
219 : Requires<detail::is_analytic_solutions_v<ParentTag>>>
220 : : Tag, db::ComputeTag {
221 : using field_tags = typename ParentTag::field_tags;
222 : using base = Tag;
223 : using DataType = typename Tag::tag::type::type;
224 : using parent_tag = AnalyticSolutions<field_tags>;
225 : using argument_tags = tmpl::list<parent_tag>;
226 : static void function(const gsl::not_null<typename base::type*> result,
227 : const typename parent_tag::type& vars) {
228 : if (vars.has_value()) {
229 : result->reset();
230 : *result = typename base::type::value_type{};
231 : const auto& tensor =
232 : get<typename detail::AnalyticImpl<typename Tag::tag>>(*vars);
233 : const size_t num_components = tensor.size();
234 : for (size_t storage_index = 0; storage_index < num_components;
235 : ++storage_index) {
236 : (**result)[storage_index].set_data_ref(
237 : // NOLINTNEXTLINE
238 : make_not_null(&const_cast<DataType&>(tensor[storage_index])));
239 : }
240 : } else {
241 : *result = std::nullopt;
242 : }
243 : }
244 : };
245 :
246 : template <typename Tag, typename ParentTag>
247 : struct Subitem<Tag, ParentTag, Requires<detail::is_errors_v<ParentTag>>>
248 : : Tag, db::ComputeTag {
249 : using field_tags = typename ParentTag::field_tags;
250 : using base = Tag;
251 : using DataType = typename Tag::tag::type::type;
252 : using parent_tag = Errors<field_tags>;
253 : using argument_tags = tmpl::list<parent_tag>;
254 : static void function(const gsl::not_null<typename base::type*> result,
255 : const typename parent_tag::type& vars) {
256 : if (vars.has_value()) {
257 : result->reset();
258 : *result = typename base::type::value_type{};
259 : const auto& tensor =
260 : get<typename detail::ErrorImpl<typename Tag::tag>>(*vars);
261 : const size_t num_components = tensor.size();
262 : for (size_t storage_index = 0; storage_index < num_components;
263 : ++storage_index) {
264 : (**result)[storage_index].set_data_ref(
265 : // NOLINTNEXTLINE
266 : make_not_null(&const_cast<DataType&>(tensor[storage_index])));
267 : }
268 : } else {
269 : *result = std::nullopt;
270 : }
271 : }
272 : };
273 : } // namespace Tags
274 :
275 : namespace db {
276 : template <typename Tag>
277 : struct Subitems<Tag, Requires<Tags::detail::is_analytic_solutions_v<Tag>>> {
278 : using field_tags = typename Tag::field_tags;
279 : using type = db::wrap_tags_in<::Tags::Analytic, field_tags>;
280 :
281 : template <typename Subtag>
282 : static void create_item(
283 : const gsl::not_null<typename Tag::type*> parent_value,
284 : const gsl::not_null<typename Subtag::type*> sub_value) {
285 : if (parent_value->has_value()) {
286 : auto& tensor = get<::Tags::detail::AnalyticImpl<typename Subtag::tag>>(
287 : parent_value->value());
288 : sub_value->emplace();
289 : // Only update the Tensor if the Variables has changed its allocation
290 : if constexpr (not is_any_spin_weighted_v<
291 : typename Subtag::tag::type::type>) {
292 : if (tensor.begin()->data() != sub_value->value().begin()->data()) {
293 : for (auto tensor_component = tensor.begin(),
294 : sub_var_it = sub_value->value().begin();
295 : tensor_component != tensor.end();
296 : ++tensor_component, ++sub_var_it) {
297 : sub_var_it->set_data_ref(make_not_null(&*tensor_component));
298 : }
299 : }
300 : } else {
301 : if (tensor.begin()->data().data() !=
302 : sub_value->value().begin()->data().data()) {
303 : for (auto tensor_component = tensor.begin(),
304 : sub_var_it = sub_value->value().begin();
305 : tensor_component != tensor.end();
306 : ++tensor_component, ++sub_var_it) {
307 : sub_var_it->set_data_ref(make_not_null(&*tensor_component));
308 : }
309 : }
310 : }
311 : } else {
312 : *sub_value = std::nullopt;
313 : }
314 : }
315 : };
316 :
317 : template <typename Tag>
318 : struct Subitems<Tag, Requires<Tags::detail::is_errors_v<Tag>>> {
319 : using field_tags = typename Tag::field_tags;
320 : using type = db::wrap_tags_in<::Tags::Error, field_tags>;
321 :
322 : template <typename Subtag>
323 : static void create_item(
324 : const gsl::not_null<typename Tag::type*> parent_value,
325 : const gsl::not_null<typename Subtag::type*> sub_value) {
326 : if (parent_value->has_value()) {
327 : auto& tensor = get<::Tags::detail::ErrorImpl<typename Subtag::tag>>(
328 : parent_value->value());
329 : *sub_value = std::decay_t<decltype(tensor)>{};
330 : // Only update the Tensor if the Variables has changed its allocation
331 : if constexpr (not is_any_spin_weighted_v<
332 : typename Subtag::tag::type::type>) {
333 : if (tensor.begin()->data() != sub_value->value().begin()->data()) {
334 : for (auto tensor_component = tensor.begin(),
335 : sub_var_it = sub_value->value().begin();
336 : tensor_component != tensor.end();
337 : ++tensor_component, ++sub_var_it) {
338 : sub_var_it->set_data_ref(make_not_null(&*tensor_component));
339 : }
340 : }
341 : } else {
342 : if (tensor.begin()->data().data() !=
343 : sub_value->value().begin()->data().data()) {
344 : for (auto tensor_component = tensor.begin(),
345 : sub_var_it = sub_value->value().begin();
346 : tensor_component != tensor.end();
347 : ++tensor_component, ++sub_var_it) {
348 : sub_var_it->set_data_ref(make_not_null(&*tensor_component));
349 : }
350 : }
351 : }
352 : } else {
353 : *sub_value = std::nullopt;
354 : }
355 : }
356 : };
357 : } // namespace db
358 : /// \endcond
|
