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