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 deserialize<type>(serialize<type>(analytic_solution).data());
      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