SpECTRE Documentation Coverage Report
Current view: top level - ParallelAlgorithms/Amr/Criteria - Loehner.hpp Hit Total Coverage
Commit: 3c2e9d3ed337bca2146eee9de07432e292a38c3a Lines: 3 33 9.1 %
Date: 2024-06-11 22:56:19
Legend: Lines: hit not hit

          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 magnitude of
      37             :  * second derivatives
      38             :  *
      39             :  * This smoothness indicator is simply the L2 norm of the logical second
      40             :  * derivative of the tensor component in the given `dimension`:
      41             :  *
      42             :  * \begin{equation}
      43             :  * \epsilon_k =
      44             :  *   \sqrt{\frac{1}{N_\mathrm{points}} \sum_{p=1}^N_\mathrm{points}
      45             :  *     \left(\partial^2 u / \partial \xi_k^2\right)^2}
      46             :  * \end{equation}
      47             :  *
      48             :  * If the smoothness indicator is large in a direction, meaning the tensor
      49             :  * component has a large second derivative in that direction, the element should
      50             :  * be h-refined. If the smoothness indicator is small, the element should be
      51             :  * h-coarsened. A coarsing threshold of about a third of the refinement
      52             :  * threshold seems to work well, but this will need more testing.
      53             :  *
      54             :  * Note that it is not at all clear that a smoothness indicator based on the
      55             :  * magnitude of second derivatives is useful in a DG context. Smooth functions
      56             :  * with higher-order derivatives can be approximated just fine by higher-order
      57             :  * DG elements without the need for h-refinement. The use of second derivatives
      58             :  * to indicate the need for refinement originated in the finite element context
      59             :  * with linear elements. Other smoothness indicators might prove more useful for
      60             :  * DG elements, e.g. based on jumps or oscillations of the solution. We can also
      61             :  * explore applying the troubled-cell indicators (TCIs) used in hydrodynamic
      62             :  * simulations as h-refinement indicators.
      63             :  *
      64             :  * Specifically, this smoothness indicator is based on \cite Loehner1987 (hence
      65             :  * the name of the function), which is popular in the finite element community
      66             :  * and also used in a DG context by \cite Dumbser2013, Eq. (34), and by
      67             :  * \cite Renkhoff2023, Eq. (15). We make several modifications:
      68             :  *
      69             :  * - The original smoothness indicator is isotropic, i.e. it computes the norm
      70             :  *   over all (mixed) second derivatives. Here we compute an anisotropic
      71             :  *   indicator by computing second derivatives in each dimension separately
      72             :  *   and ignoring mixed derivatives.
      73             :  * - The original smoothness indicator is normalized by measures of the first
      74             :  *   derivative which don't generalize well to spectral elements. Therefore, we
      75             :  *   simplify the normalization to a standard relative and absolute tolerance.
      76             :  *   An alternative approach is proposed in \cite Renkhoff2023, Eq.(15), where
      77             :  *   the authors take the absolute value of the differentiation matrix and apply
      78             :  *   the resulting matrix to the absolute value of the data on the grid to
      79             :  *   compute the normalization. However, this quantity can produce quite large
      80             :  *   numbers and hence overestimates the smoothness by suppressing the second
      81             :  *   derivative.
      82             :  * - We compute the second derivative in logical coordinates. This seems
      83             :  *   easiest for spectral elements, but note that \cite Renkhoff2023 seem to
      84             :  *   use inertial coordinates.
      85             :  *
      86             :  * In addition to the above modifications, we can consider approximating the
      87             :  * second derivative using finite differences, as explored in the prototype
      88             :  * https://github.com/sxs-collaboration/dg-charm/blob/HpAmr/Evolution/HpAmr/LohnerRefiner.hpp.
      89             :  * This would allow falling back to the normalization used by Löhner and might
      90             :  * be cheaper to compute, but it requires an interpolation to the center and
      91             :  * maybe also to the faces, depending on the desired stencil.
      92             :  */
      93             : template <size_t Dim>
      94           1 : double loehner_smoothness_indicator(
      95             :     gsl::not_null<DataVector*> first_deriv_buffer,
      96             :     gsl::not_null<DataVector*> second_deriv_buffer,
      97             :     const DataVector& tensor_component, const Mesh<Dim>& mesh,
      98             :     size_t dimension);
      99             : template <size_t Dim>
     100           1 : std::array<double, Dim> loehner_smoothness_indicator(
     101             :     const DataVector& tensor_component, const Mesh<Dim>& mesh);
     102             : /// @}
     103             : 
     104             : namespace Loehner_detail {
     105             : template <size_t Dim>
     106             : void max_over_components(
     107             :     gsl::not_null<std::array<Flag, Dim>*> result,
     108             :     gsl::not_null<std::array<DataVector, 2>*> deriv_buffers,
     109             :     const DataVector& tensor_component, const Mesh<Dim>& mesh,
     110             :     double relative_tolerance, double absolute_tolerance,
     111             :     double coarsening_factor);
     112             : }
     113             : 
     114             : /*!
     115             :  * \brief h-refine the grid based on a smoothness indicator
     116             :  *
     117             :  * The smoothness indicator used here is based on the magnitude of second
     118             :  * derivatives. See `amr::Criteria::loehner_smoothness_indicator` for details
     119             :  * and caveats.
     120             :  *
     121             :  * \see amr::Criteria::loehner_smoothness_indicator
     122             :  */
     123             : template <size_t Dim, typename TensorTags>
     124           1 : class Loehner : public Criterion {
     125             :  public:
     126           0 :   struct VariablesToMonitor {
     127           0 :     using type = std::vector<std::string>;
     128           0 :     static constexpr Options::String help = {
     129             :         "The tensors to monitor for h-refinement."};
     130           0 :     static size_t lower_bound_on_size() { return 1; }
     131             :   };
     132           0 :   struct RelativeTolerance {
     133           0 :     using type = double;
     134           0 :     static constexpr Options::String help = {
     135             :         "If any tensor component has a second derivative magnitude above this "
     136             :         "value times the max of the absolute tensor component over the "
     137             :         "element, the element will be h-refined in that direction. "
     138             :         "Set to 0 to disable."};
     139           0 :     static double lower_bound() { return 0.; }
     140             :   };
     141           0 :   struct AbsoluteTolerance {
     142           0 :     using type = double;
     143           0 :     static constexpr Options::String help = {
     144             :         "If any tensor component has a second derivative magnitude above this "
     145             :         "value, the element will be h-refined in that direction. "
     146             :         "Set to 0 to disable."};
     147           0 :     static double lower_bound() { return 0.; }
     148             :   };
     149           0 :   struct CoarseningFactor {
     150           0 :     using type = double;
     151           0 :     static constexpr Options::String help = {
     152             :         "Factor applied to both relative and absolute tolerance to trigger "
     153             :         "h-coarsening. Set to 0 to disable h-coarsening altogether. "
     154             :         "Set closer to 1 to trigger h-coarsening more aggressively. "
     155             :         "Values too close to 1 risk that coarsened elements will immediately "
     156             :         "trigger h-refinement again. A reasonable value is 1/3."};
     157           0 :     static double lower_bound() { return 0.; }
     158           0 :     static double upper_bound() { return 1.; }
     159             :   };
     160             : 
     161           0 :   using options = tmpl::list<VariablesToMonitor, RelativeTolerance,
     162             :                              AbsoluteTolerance, CoarseningFactor>;
     163             : 
     164           0 :   static constexpr Options::String help = {
     165             :       "Refine the grid towards resolving an estimated error in the second "
     166             :       "derivative"};
     167             : 
     168           0 :   Loehner() = default;
     169             : 
     170           0 :   Loehner(std::vector<std::string> vars_to_monitor, double relative_tolerance,
     171             :           double absolute_tolerance, double coarsening_factor,
     172             :           const Options::Context& context = {});
     173             : 
     174             :   /// \cond
     175             :   explicit Loehner(CkMigrateMessage* msg);
     176             :   using PUP::able::register_constructor;
     177             :   WRAPPED_PUPable_decl_template(Loehner);  // NOLINT
     178             :   /// \endcond
     179             : 
     180           0 :   using compute_tags_for_observation_box = tmpl::list<>;
     181             : 
     182           0 :   using argument_tags = tmpl::list<::Tags::DataBox>;
     183             : 
     184             :   template <typename DbTagsList, typename Metavariables>
     185           0 :   std::array<Flag, Dim> operator()(const db::DataBox<DbTagsList>& box,
     186             :                                    Parallel::GlobalCache<Metavariables>& cache,
     187             :                                    const ElementId<Dim>& element_id) const;
     188             : 
     189           0 :   void pup(PUP::er& p) override;
     190             : 
     191             :  private:
     192           0 :   std::vector<std::string> vars_to_monitor_{};
     193           0 :   double relative_tolerance_ = std::numeric_limits<double>::signaling_NaN();
     194           0 :   double absolute_tolerance_ = std::numeric_limits<double>::signaling_NaN();
     195           0 :   double coarsening_factor_ = std::numeric_limits<double>::signaling_NaN();
     196             : };
     197             : 
     198             : // Out-of-line definitions
     199             : /// \cond
     200             : 
     201             : template <size_t Dim, typename TensorTags>
     202             : Loehner<Dim, TensorTags>::Loehner(std::vector<std::string> vars_to_monitor,
     203             :                                   const double relative_tolerance,
     204             :                                   const double absolute_tolerance,
     205             :                                   const double coarsening_factor,
     206             :                                   const Options::Context& context)
     207             :     : vars_to_monitor_(std::move(vars_to_monitor)),
     208             :       relative_tolerance_(relative_tolerance),
     209             :       absolute_tolerance_(absolute_tolerance),
     210             :       coarsening_factor_(coarsening_factor) {
     211             :   db::validate_selection<TensorTags>(vars_to_monitor_, context);
     212             :   if (relative_tolerance == 0. and absolute_tolerance == 0.) {
     213             :     PARSE_ERROR(
     214             :         context,
     215             :         "Must specify non-zero RelativeTolerance, AbsoluteTolerance, or both.");
     216             :   }
     217             : }
     218             : 
     219             : template <size_t Dim, typename TensorTags>
     220             : Loehner<Dim, TensorTags>::Loehner(CkMigrateMessage* msg) : Criterion(msg) {}
     221             : 
     222             : template <size_t Dim, typename TensorTags>
     223             : template <typename DbTagsList, typename Metavariables>
     224             : std::array<Flag, Dim> Loehner<Dim, TensorTags>::operator()(
     225             :     const db::DataBox<DbTagsList>& box,
     226             :     Parallel::GlobalCache<Metavariables>& /*cache*/,
     227             :     const ElementId<Dim>& /*element_id*/) const {
     228             :   auto result = make_array<Dim>(Flag::Undefined);
     229             :   const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
     230             :   // Check all tensors and all tensor components in turn. We take the
     231             :   // highest-priority refinement flag in each dimension, so if any tensor
     232             :   // component is non-smooth, the element will split in that dimension. And only
     233             :   // if all tensor components are smooth enough will elements join in that
     234             :   // dimension.
     235             :   std::array<DataVector, 2> deriv_buffers{};
     236             :   tmpl::for_each<TensorTags>(
     237             :       [&result, &box, &mesh, &deriv_buffers, this](const auto tag_v) {
     238             :         // Stop if we have already decided to refine every dimension
     239             :         if (result == make_array<Dim>(Flag::Split)) {
     240             :           return;
     241             :         }
     242             :         using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
     243             :         const std::string tag_name = db::tag_name<tag>();
     244             :         // Skip if this tensor is not being monitored
     245             :         if (not alg::found(vars_to_monitor_, tag_name)) {
     246             :           return;
     247             :         }
     248             :         const auto& tensor = db::get<tag>(box);
     249             :         for (const DataVector& tensor_component : tensor) {
     250             :           Loehner_detail::max_over_components(
     251             :               make_not_null(&result), make_not_null(&deriv_buffers),
     252             :               tensor_component, mesh, relative_tolerance_, absolute_tolerance_,
     253             :               coarsening_factor_);
     254             :         }
     255             :       });
     256             :   return result;
     257             : }
     258             : 
     259             : template <size_t Dim, typename TensorTags>
     260             : void Loehner<Dim, TensorTags>::pup(PUP::er& p) {
     261             :   p | vars_to_monitor_;
     262             :   p | relative_tolerance_;
     263             :   p | absolute_tolerance_;
     264             :   p | coarsening_factor_;
     265             : }
     266             : 
     267             : template <size_t Dim, typename TensorTags>
     268             : PUP::able::PUP_ID Loehner<Dim, TensorTags>::my_PUP_ID = 0;  // NOLINT
     269             : /// \endcond
     270             : 
     271             : }  // namespace amr::Criteria

Generated by: LCOV version 1.14