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 <pup.h>
      10             : #include <string>
      11             : #include <vector>
      12             : 
      13             : #include "DataStructures/DataBox/DataBox.hpp"
      14             : #include "DataStructures/DataBox/DataBoxTag.hpp"
      15             : #include "DataStructures/DataBox/ValidateSelection.hpp"
      16             : #include "DataStructures/DataVector.hpp"
      17             : #include "DataStructures/Tensor/Tensor.hpp"
      18             : #include "Domain/Amr/Flag.hpp"
      19             : #include "Domain/Tags.hpp"
      20             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      21             : #include "Options/Context.hpp"
      22             : #include "Options/ParseError.hpp"
      23             : #include "Options/String.hpp"
      24             : #include "ParallelAlgorithms/Amr/Criteria/Criterion.hpp"
      25             : #include "ParallelAlgorithms/Amr/Criteria/Type.hpp"
      26             : #include "Utilities/TMPL.hpp"
      27             : 
      28             : /// \cond
      29             : template <size_t>
      30             : class ElementId;
      31             : /// \endcond
      32             : 
      33             : namespace amr::Criteria {
      34             : 
      35             : /// @{
      36             : /*!
      37             :  * \brief Computes an anisotropic smoothness indicator based on the power in the
      38             :  * highest modes
      39             :  *
      40             :  * This smoothness indicator is the L2 norm of the tensor component with the
      41             :  * lowest N - M modes filtered out, where N is the number of grid points in the
      42             :  * given dimension and M is `num_highest_modes`.
      43             :  *
      44             :  * This strategy is similar to the Persson troubled-cell indicator (see
      45             :  * `evolution::dg::subcell::persson_tci`) with a few modifications:
      46             :  *
      47             :  * - The indicator is computed in each dimension separately for an anisotropic
      48             :  *   measure.
      49             :  * - The number of highest modes to keep can be chosen as a parameter.
      50             :  * - We don't normalize by the L2 norm of the unfiltered data $u$ here. This
      51             :  *   function just returns the L2 norm of the filtered data.
      52             :  */
      53             : template <typename VectorType, size_t Dim>
      54           1 : double persson_smoothness_indicator(
      55             :     gsl::not_null<VectorType*> filtered_component_buffer,
      56             :     const VectorType& tensor_component, const Mesh<Dim>& mesh, size_t dimension,
      57             :     size_t num_highest_modes);
      58             : template <typename VectorType, size_t Dim>
      59           1 : std::array<double, Dim> persson_smoothness_indicator(
      60             :     const VectorType& tensor_component, const Mesh<Dim>& mesh,
      61             :     size_t num_highest_modes);
      62             : /// @}
      63             : 
      64             : namespace Persson_detail {
      65             : template <typename VectorType, size_t Dim>
      66             : void max_over_components(gsl::not_null<std::array<Flag, Dim>*> result,
      67             :                          gsl::not_null<VectorType*> buffer,
      68             :                          const VectorType& tensor_component,
      69             :                          const Mesh<Dim>& mesh, size_t num_highest_modes,
      70             :                          double alpha, double absolute_tolerance,
      71             :                          double coarsening_factor);
      72             : }
      73             : 
      74             : /*!
      75             :  * \brief h-refine the grid based on power in the highest modes
      76             :  *
      77             :  * \see persson_smoothness_indicator
      78             :  */
      79             : template <size_t Dim, typename TensorTags>
      80           1 : class Persson : public Criterion {
      81             :  public:
      82           0 :   struct VariablesToMonitor {
      83           0 :     using type = std::vector<std::string>;
      84           0 :     static constexpr Options::String help = {
      85             :         "The tensors to monitor for h-refinement."};
      86           0 :     static size_t lower_bound_on_size() { return 1; }
      87             :   };
      88           0 :   struct NumHighestModes {
      89           0 :     using type = size_t;
      90           0 :     static constexpr Options::String help = {
      91             :         "Number of highest modes to monitor the power of."};
      92           0 :     static size_t lower_bound() { return 1; }
      93             :   };
      94           0 :   struct Exponent {
      95           0 :     using type = double;
      96           0 :     static constexpr Options::String help = {
      97             :         "The exponent at which the modes should decrease. "
      98             :         "Corresponds to a \"relative tolerance\" of N^(-alpha), where N is the "
      99             :         "number of grid points minus 'NumHighestModes'. "
     100             :         "If any tensor component has power in the highest modes above this "
     101             :         "value times the max of the absolute tensor component over the "
     102             :         "element, the element will be h-refined in that direction."};
     103           0 :     static double lower_bound() { return 0.; }
     104             :   };
     105           0 :   struct AbsoluteTolerance {
     106           0 :     using type = double;
     107           0 :     static constexpr Options::String help = {
     108             :         "If any tensor component has a power in the highest modes above this "
     109             :         "value, the element will be h-refined in that direction. "
     110             :         "Set to 0 to disable."};
     111           0 :     static double lower_bound() { return 0.; }
     112             :   };
     113           0 :   struct CoarseningFactor {
     114           0 :     using type = double;
     115           0 :     static constexpr Options::String help = {
     116             :         "Factor applied to both relative and absolute tolerance to trigger "
     117             :         "h-coarsening. Set to 0 to disable h-coarsening altogether. "
     118             :         "Set closer to 1 to trigger h-coarsening more aggressively. "
     119             :         "Values too close to 1 risk that coarsened elements will immediately "
     120             :         "trigger h-refinement again. A reasonable value is 0.1."};
     121           0 :     static double lower_bound() { return 0.; }
     122           0 :     static double upper_bound() { return 1.; }
     123             :   };
     124             : 
     125           0 :   using options = tmpl::list<VariablesToMonitor, NumHighestModes, Exponent,
     126             :                              AbsoluteTolerance, CoarseningFactor>;
     127             : 
     128           0 :   static constexpr Options::String help = {
     129             :       "Refine the grid so the power in the highest modes stays below the "
     130             :       "tolerance"};
     131             : 
     132           0 :   Persson() = default;
     133             : 
     134           0 :   Persson(std::vector<std::string> vars_to_monitor,
     135             :           const size_t num_highest_modes, double alpha,
     136             :           double absolute_tolerance, double coarsening_factor,
     137             :           const Options::Context& context = {});
     138             : 
     139             :   /// \cond
     140             :   explicit Persson(CkMigrateMessage* msg);
     141             :   using PUP::able::register_constructor;
     142             :   WRAPPED_PUPable_decl_template(Persson);  // NOLINT
     143             :   /// \endcond
     144             : 
     145           0 :   Type type() override { return Type::h; }
     146             : 
     147           0 :   std::string observation_name() override { return "Persson"; }
     148             : 
     149           0 :   using compute_tags_for_observation_box = tmpl::list<>;
     150             : 
     151           0 :   using argument_tags = tmpl::list<::Tags::DataBox>;
     152             : 
     153             :   template <typename DbTagsList, typename Metavariables>
     154           0 :   std::array<Flag, Dim> operator()(const db::DataBox<DbTagsList>& box,
     155             :                                    Parallel::GlobalCache<Metavariables>& cache,
     156             :                                    const ElementId<Dim>& element_id) const;
     157             : 
     158           0 :   void pup(PUP::er& p) override;
     159             : 
     160             :  private:
     161           0 :   std::vector<std::string> vars_to_monitor_{};
     162           0 :   size_t num_highest_modes_{};
     163           0 :   double alpha_ = std::numeric_limits<double>::signaling_NaN();
     164           0 :   double absolute_tolerance_ = std::numeric_limits<double>::signaling_NaN();
     165           0 :   double coarsening_factor_ = std::numeric_limits<double>::signaling_NaN();
     166             : };
     167             : 
     168             : // Out-of-line definitions
     169             : /// \cond
     170             : 
     171             : template <size_t Dim, typename TensorTags>
     172             : Persson<Dim, TensorTags>::Persson(std::vector<std::string> vars_to_monitor,
     173             :                                   const size_t num_highest_modes,
     174             :                                   const double alpha,
     175             :                                   const double absolute_tolerance,
     176             :                                   const double coarsening_factor,
     177             :                                   const Options::Context& context)
     178             :     : vars_to_monitor_(std::move(vars_to_monitor)),
     179             :       num_highest_modes_(num_highest_modes),
     180             :       alpha_(alpha),
     181             :       absolute_tolerance_(absolute_tolerance),
     182             :       coarsening_factor_(coarsening_factor) {
     183             :   db::validate_selection<TensorTags>(vars_to_monitor_, context);
     184             : }
     185             : 
     186             : template <size_t Dim, typename TensorTags>
     187             : Persson<Dim, TensorTags>::Persson(CkMigrateMessage* msg) : Criterion(msg) {}
     188             : 
     189             : template <size_t Dim, typename TensorTags>
     190             : template <typename DbTagsList, typename Metavariables>
     191             : std::array<Flag, Dim> Persson<Dim, TensorTags>::operator()(
     192             :     const db::DataBox<DbTagsList>& box,
     193             :     Parallel::GlobalCache<Metavariables>& /*cache*/,
     194             :     const ElementId<Dim>& /*element_id*/) const {
     195             :   auto result = make_array<Dim>(Flag::Undefined);
     196             :   const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
     197             :   using VectorType = typename Variables<TensorTags>::vector_type;
     198             :   // Check all tensors and all tensor components in turn. We take the
     199             :   // highest-priority refinement flag in each dimension, so if any tensor
     200             :   // component is non-smooth, the element will split in that dimension. And only
     201             :   // if all tensor components are smooth enough will elements join in that
     202             :   // dimension.
     203             :   VectorType buffer(mesh.number_of_grid_points());
     204             :   tmpl::for_each<TensorTags>(
     205             :       [&result, &box, &mesh, &buffer, this](const auto tag_v) {
     206             :         // Stop if we have already decided to refine every dimension
     207             :         if (result == make_array<Dim>(Flag::Split)) {
     208             :           return;
     209             :         }
     210             :         using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
     211             :         const std::string tag_name = db::tag_name<tag>();
     212             :         // Skip if this tensor is not being monitored
     213             :         if (not alg::found(vars_to_monitor_, tag_name)) {
     214             :           return;
     215             :         }
     216             :         const auto& tensor = db::get<tag>(box);
     217             :         for (const VectorType& tensor_component : tensor) {
     218             :           Persson_detail::max_over_components(
     219             :               make_not_null(&result), make_not_null(&buffer), tensor_component,
     220             :               mesh, num_highest_modes_, alpha_, absolute_tolerance_,
     221             :               coarsening_factor_);
     222             :         }
     223             :       });
     224             :   return result;
     225             : }
     226             : 
     227             : template <size_t Dim, typename TensorTags>
     228             : void Persson<Dim, TensorTags>::pup(PUP::er& p) {
     229             :   p | vars_to_monitor_;
     230             :   p | num_highest_modes_;
     231             :   p | alpha_;
     232             :   p | absolute_tolerance_;
     233             :   p | coarsening_factor_;
     234             : }
     235             : 
     236             : template <size_t Dim, typename TensorTags>
     237             : PUP::able::PUP_ID Persson<Dim, TensorTags>::my_PUP_ID = 0;  // NOLINT
     238             : /// \endcond
     239             : 
     240             : }  // namespace amr::Criteria