Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : #include "DataStructures/DataVector.hpp"
       9             : #include "DataStructures/Index.hpp"
      10             : #include "DataStructures/Matrix.hpp"
      11             : #include "DataStructures/Tensor/Tensor.hpp"
      12             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      13             : #include "Utilities/Gsl.hpp"
      14             : 
      15             : namespace evolution::dg::subcell {
      16             : namespace detail {
      17             : template <size_t Dim>
      18             : bool persson_tci_impl(gsl::not_null<DataVector*> filtered_component,
      19             :                       const DataVector& component, const Mesh<Dim>& dg_mesh,
      20             :                       double alpha, size_t num_highest_modes);
      21             : }  // namespace detail
      22             : 
      23             : /*!
      24             :  * \brief Troubled cell indicator using the idea of spectral falloff by
      25             :  *    \cite Persson2006sub
      26             :  *
      27             :  * Consider a discontinuity sensing quantity \f$U\f$, which is typically a
      28             :  * scalar but could be a tensor of any rank. Let \f$U\f$ have the 1d spectral
      29             :  * decomposition (generalization to higher-dimensional tensor product bases is
      30             :  * done dimension-by-dimension):
      31             :  *
      32             :  * \f{align*}{
      33             :  *   U(x)=\sum_{i=0}^{N}c_i P_i(x),
      34             :  * \f}
      35             :  *
      36             :  * where \f$P_i(x)\f$ are the basis functions, in our case the Legendre
      37             :  * polynomials, and \f$c_i\f$ are the spectral coefficients. We then define a
      38             :  * filtered solution \f$\hat{U}\f$ as
      39             :  *
      40             :  * \f{align*}{
      41             :  *   \hat{U}(x)=\sum_{i=N+1-M}^{N} c_i P_i(x).
      42             :  * \f}
      43             :  *
      44             :  * where $M$ is the number of highest modes to include in the filtered solution.
      45             :  *
      46             :  * Note that when an exponential filter is being used to deal with aliasing,
      47             :  * lower modes can be included in \f$\hat{U}\f$. The main goal of \f$\hat{U}\f$
      48             :  * is to measure how much power is in the highest modes, which are the modes
      49             :  * responsible for Gibbs phenomena.
      50             :  *
      51             :  * A cell is troubled if
      52             :  *
      53             :  * \f{align*}{
      54             :  *    \frac{(\hat{U}, \hat{U})}{(U, U)} > (N + 1 - M)^{-\alpha}
      55             :  * \f}
      56             :  *
      57             :  * where \f$(\cdot,\cdot)\f$ is an inner product, which we take to be the
      58             :  * Euclidean \f$L_2\f$ norm (i.e. we do not divide by the number of grid points
      59             :  * since that cancels out anyway) computed as
      60             :  *
      61             :  * \f{align*}{
      62             :  *    (U, U) \equiv \sqrt{\sum_{i=0}^N U_i^2}
      63             :  * \f}
      64             :  *
      65             :  * where $U_i$ are nodal values of the quantity $U$ at grid points.
      66             :  *
      67             :  * Typically, \f$\alpha=4.0\f$ and $M=1$ is a good choice.
      68             :  *
      69             :  */
      70             : template <size_t Dim, typename SymmList, typename IndexList>
      71           1 : bool persson_tci(const Tensor<DataVector, SymmList, IndexList>& tensor,
      72             :                  const Mesh<Dim>& dg_mesh, const double alpha,
      73             :                  const size_t num_highest_modes) {
      74             :   DataVector filtered_component(dg_mesh.number_of_grid_points());
      75             :   for (size_t component_index = 0; component_index < tensor.size();
      76             :        ++component_index) {
      77             :     if (detail::persson_tci_impl(make_not_null(&filtered_component),
      78             :                                  tensor[component_index], dg_mesh, alpha,
      79             :                                  num_highest_modes)) {
      80             :       return true;
      81             :     }
      82             :   }
      83             :   return false;
      84             : }
      85             : }  // namespace evolution::dg::subcell