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