Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <limits>
       9             : #include <optional>
      10             : #include <pup.h>
      11             : #include <string>
      12             : #include <vector>
      13             : 
      14             : #include "DataStructures/DataBox/ObservationBox.hpp"
      15             : #include "DataStructures/DataBox/ValidateSelection.hpp"
      16             : #include "DataStructures/DataVector.hpp"
      17             : #include "DataStructures/Tags/TempTensor.hpp"
      18             : #include "DataStructures/TempBuffer.hpp"
      19             : #include "DataStructures/Tensor/Tensor.hpp"
      20             : #include "Domain/Amr/Flag.hpp"
      21             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      22             : #include "Options/Context.hpp"
      23             : #include "Options/ParseError.hpp"
      24             : #include "Options/String.hpp"
      25             : #include "ParallelAlgorithms/Amr/Criteria/Criterion.hpp"
      26             : #include "ParallelAlgorithms/Amr/Criteria/Type.hpp"
      27             : #include "ParallelAlgorithms/Events/Tags.hpp"
      28             : #include "Utilities/Algorithm.hpp"
      29             : #include "Utilities/TMPL.hpp"
      30             : 
      31             : /// \cond
      32             : template <size_t>
      33             : class ElementId;
      34             : /// \endcond
      35             : 
      36           1 : namespace amr::Criteria {
      37             : 
      38             : namespace Constraints_detail {
      39             : 
      40             : /*!
      41             :  * \brief Computes the (squared) normalization factor $N_\hat{k}^2$ (see
      42             :  * `logical_constraints` below)
      43             :  */
      44             : template <size_t Dim>
      45             : void normalization_factor_square(
      46             :     gsl::not_null<tnsr::i<DataVector, Dim, Frame::ElementLogical>*> result,
      47             :     const Jacobian<DataVector, Dim, Frame::ElementLogical, Frame::Inertial>&
      48             :         jacobian);
      49             : 
      50             : /*!
      51             :  * \brief Computes a measure of the constraints in each logical direction of the
      52             :  * grid
      53             :  *
      54             :  * \requires `constraints_tensor` to be a tensor of rank 2 or higher. The first
      55             :  * index must be a lower spatial index that originates from a derivative.
      56             :  *
      57             :  * We follow \cite Szilagyi2014fna, Eq. (62)-(64) to compute:
      58             :  *
      59             :  * \begin{align}
      60             :  * \Epsilon_\hat{k} &= \frac{1}{N_\hat{k}} \sqrt{\sum_{a\ldots} \left(\sum_{i}
      61             :  *   \frac{\partial x^i}{\partial x^\hat{k}} C_{ia\ldots}\right)^2} \\
      62             :  * N_\hat{k} &= \sqrt{\sum_{i} \left(\frac{\partial x^i}{\partial x^\hat{k}}
      63             :  *   \right)^2}
      64             :  * \end{align}
      65             :  *
      66             :  * This transform the first lower spatial index of the tensor to the
      67             :  * element-logical frame, then takes an L2 norm over all remaining indices.
      68             :  * The (squared) normalization factor $N_\hat{k}^2$ is computed by
      69             :  * `normalization_factor` and passed in as an argument.
      70             :  */
      71             : template <size_t Dim, typename TensorType>
      72             : void logical_constraints(
      73             :     gsl::not_null<tnsr::i<DataVector, Dim, Frame::ElementLogical>*> result,
      74             :     gsl::not_null<DataVector*> buffer, const TensorType& constraints_tensor,
      75             :     const Jacobian<DataVector, Dim, Frame::ElementLogical, Frame::Inertial>&
      76             :         jacobian,
      77             :     const tnsr::i<DataVector, Dim, Frame::ElementLogical>&
      78             :         normalization_factor_square);
      79             : 
      80             : /*!
      81             :  * \brief Apply the AMR criterion to one of the constraints
      82             :  *
      83             :  * The `result` is the current decision in each dimension based on the previous
      84             :  * constraints. This function will update the flags if necessary. It takes
      85             :  * the "max" of the current and new flags, where the "highest" flag is
      86             :  * `Flag::IncreaseResolution`, followed by `Flag::DoNothing`, and then
      87             :  * `Flag::DecreaseResolution`.
      88             :  */
      89             : template <size_t Dim>
      90             : void max_over_components(
      91             :     gsl::not_null<std::array<Flag, Dim>*> result,
      92             :     const tnsr::i<DataVector, Dim, Frame::ElementLogical>& logical_constraints,
      93             :     double abs_target, double coarsening_factor);
      94             : 
      95             : }  // namespace Constraints_detail
      96             : 
      97             : /*!
      98             :  * \brief Refine the grid towards the target constraint violation
      99             :  *
     100             :  * - If any constraint is above the target value, the element will be p-refined.
     101             :  * - If all constraints are below the target times the "coarsening factor" the
     102             :  *   element will be p-coarsened.
     103             :  *
     104             :  * This criterion is based on Sec. 6.1.4 in \cite Szilagyi2014fna .
     105             :  *
     106             :  * If the coarsening factor turns out to be hard to choose, then we can try to
     107             :  * eliminate it by projecting the variables to a lower polynomial order before
     108             :  * computing constraints, or something like that.
     109             :  *
     110             :  * \tparam Dim Spatial dimension of the grid
     111             :  * \tparam TensorTags List of tags of the constraints to be monitored. These
     112             :  * must be tensors of rank 2 or higher. The first index must be a lower spatial
     113             :  * index that originates from a derivative.
     114             :  */
     115             : template <size_t Dim, typename TensorTags>
     116           1 : class Constraints : public Criterion {
     117             :  public:
     118           0 :   struct ConstraintsToMonitor {
     119           0 :     using type = std::vector<std::string>;
     120           0 :     static constexpr Options::String help = {"The constraints to monitor."};
     121           0 :     static size_t lower_bound_on_size() { return 1; }
     122             :   };
     123           0 :   struct AbsoluteTarget {
     124           0 :     using type = double;
     125           0 :     static constexpr Options::String help = {
     126             :         "The absolute target constraint violation. If any constraint is above "
     127             :         "this value, the element will be p-refined."};
     128           0 :     static double lower_bound() { return 0.; }
     129             :   };
     130           0 :   struct CoarseningFactor {
     131           0 :     using type = double;
     132           0 :     static constexpr Options::String help = {
     133             :         "If all constraints are below the 'AbsoluteTarget' times this factor, "
     134             :         "the element will be p-coarsened. "
     135             :         "A reasonable value is 0.1."};
     136           0 :     static double lower_bound() { return 0.; }
     137           0 :     static double upper_bound() { return 1.; }
     138             :   };
     139             : 
     140           0 :   using options =
     141             :       tmpl::list<ConstraintsToMonitor, AbsoluteTarget, CoarseningFactor>;
     142             : 
     143           0 :   static constexpr Options::String help = {
     144             :       "Refine the grid towards the target constraint violation"};
     145             : 
     146           0 :   Constraints() = default;
     147             : 
     148           0 :   Constraints(std::vector<std::string> vars_to_monitor, double abs_target,
     149             :               double coarsening_factor, const Options::Context& context = {});
     150             : 
     151             :   /// \cond
     152             :   explicit Constraints(CkMigrateMessage* msg);
     153             :   using PUP::able::register_constructor;
     154             :   WRAPPED_PUPable_decl_template(Constraints);  // NOLINT
     155             :   /// \endcond
     156             : 
     157           0 :   using compute_tags_for_observation_box = tmpl::append<
     158             :       TensorTags, tmpl::list<domain::Tags::JacobianCompute<
     159             :                                  Dim, Frame::ElementLogical, Frame::Inertial>,
     160             :                              Events::Tags::ObserverJacobianCompute<
     161             :                                  Dim, Frame::ElementLogical, Frame::Inertial>>>;
     162             : 
     163           0 :   Type type() override { return Type::p; }
     164             : 
     165           0 :   std::string observation_name() override { return "Constraints"; }
     166             : 
     167           0 :   using argument_tags = tmpl::list<::Tags::ObservationBox>;
     168             : 
     169             :   template <typename ComputeTagsList, typename DataBoxType,
     170             :             typename Metavariables>
     171           0 :   std::array<Flag, Dim> operator()(
     172             :       const ObservationBox<ComputeTagsList, DataBoxType>& box,
     173             :       Parallel::GlobalCache<Metavariables>& cache,
     174             :       const ElementId<Dim>& element_id) const;
     175             : 
     176           0 :   void pup(PUP::er& p) override;
     177             : 
     178             :  private:
     179           0 :   std::vector<std::string> vars_to_monitor_{};
     180           0 :   double abs_target_ = std::numeric_limits<double>::signaling_NaN();
     181           0 :   double coarsening_factor_ = std::numeric_limits<double>::signaling_NaN();
     182             : };
     183             : 
     184             : // Out-of-line definitions
     185             : /// \cond
     186             : 
     187             : template <size_t Dim, typename TensorTags>
     188             : Constraints<Dim, TensorTags>::Constraints(
     189             :     std::vector<std::string> vars_to_monitor, const double abs_target,
     190             :     const double coarsening_factor, const Options::Context& context)
     191             :     : vars_to_monitor_(std::move(vars_to_monitor)),
     192             :       abs_target_(abs_target),
     193             :       coarsening_factor_(coarsening_factor) {
     194             :   db::validate_selection<TensorTags>(vars_to_monitor_, context);
     195             : }
     196             : 
     197             : template <size_t Dim, typename TensorTags>
     198             : Constraints<Dim, TensorTags>::Constraints(CkMigrateMessage* msg)
     199             :     : Criterion(msg) {}
     200             : 
     201             : template <size_t Dim, typename TensorTags>
     202             : template <typename ComputeTagsList, typename DataBoxType,
     203             :           typename Metavariables>
     204             : std::array<Flag, Dim> Constraints<Dim, TensorTags>::operator()(
     205             :     const ObservationBox<ComputeTagsList, DataBoxType>& box,
     206             :     Parallel::GlobalCache<Metavariables>& /*cache*/,
     207             :     const ElementId<Dim>& /*element_id*/) const {
     208             :   auto result = make_array<Dim>(Flag::Undefined);
     209             :   const auto& jacobian =
     210             :       get<Events::Tags::ObserverJacobian<Dim, Frame::ElementLogical,
     211             :                                          Frame::Inertial>>(box);
     212             :   // Set up memory buffers
     213             :   const size_t num_points = jacobian.begin()->size();
     214             :   TempBuffer<tmpl::list<::Tags::Tempi<0, Dim, Frame::ElementLogical>,
     215             :                         ::Tags::Tempi<1, Dim, Frame::ElementLogical>,
     216             :                         ::Tags::TempScalar<2>>>
     217             :       buffer{num_points};
     218             :   auto& normalization_factor_square =
     219             :       get<::Tags::Tempi<0, Dim, Frame::ElementLogical>>(buffer);
     220             :   auto& logical_constraints =
     221             :       get<::Tags::Tempi<1, Dim, Frame::ElementLogical>>(buffer);
     222             :   auto& scalar_buffer = get(get<::Tags::TempScalar<2>>(buffer));
     223             :   Constraints_detail::normalization_factor_square(
     224             :       make_not_null(&normalization_factor_square), jacobian);
     225             :   // Check all constraints in turn
     226             :   tmpl::for_each<TensorTags>([&result, &box, &jacobian,
     227             :                               &normalization_factor_square,
     228             :                               &logical_constraints, &scalar_buffer,
     229             :                               this](const auto tag_v) {
     230             :     // Stop if we have already decided to refine every dimension
     231             :     if (result == make_array<Dim>(Flag::IncreaseResolution)) {
     232             :       return;
     233             :     }
     234             :     using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
     235             :     const std::string tag_name = db::tag_name<tag>();
     236             :     // Skip if this tensor is not being monitored
     237             :     if (not alg::found(vars_to_monitor_, tag_name)) {
     238             :       return;
     239             :     }
     240             :     Constraints_detail::logical_constraints(
     241             :         make_not_null(&logical_constraints), make_not_null(&scalar_buffer),
     242             :         get<tag>(box), jacobian, normalization_factor_square);
     243             :     Constraints_detail::max_over_components(make_not_null(&result),
     244             :                                             logical_constraints, abs_target_,
     245             :                                             coarsening_factor_);
     246             :   });
     247             :   return result;
     248             : }
     249             : 
     250             : template <size_t Dim, typename TensorTags>
     251             : void Constraints<Dim, TensorTags>::pup(PUP::er& p) {
     252             :   p | vars_to_monitor_;
     253             :   p | abs_target_;
     254             :   p | coarsening_factor_;
     255             : }
     256             : 
     257             : template <size_t Dim, typename TensorTags>
     258             : PUP::able::PUP_ID Constraints<Dim, TensorTags>::my_PUP_ID = 0;  // NOLINT
     259             : /// \endcond
     260             : 
     261             : }  // namespace amr::Criteria