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