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