PerssonTci.hpp
1 // 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"
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, double zero_cutoff) noexcept;
21 } // namespace detail
22 
23 /*!
24  * \brief Troubled cell indicator using spectral falloff of \cite Persson2006sub
25  *
26  * Consider a discontinuity sensing quantity \f$U\f$, which is typically a
27  * scalar but could be a tensor of any rank. Let \f$U\f$ have the 1d spectral
28  * decomposition (generalization to higher-dimensional tensor product bases is
29  * done dimension-by-dimension):
30  *
31  * \f{align*}{
32  * U(x)=\sum_{i=0}^{N}c_i P_i(x),
33  * \f}
34  *
35  * where \f$P_i(x)\f$ are the basis functions, in our case the Legendre
36  * polynomials, and \f$c_i\f$ are the spectral coefficients. We then define a
37  * filtered solution \f$\hat{U}\f$ as
38  *
39  * \f{align*}{
40  * \hat{U}(x)=c_N P_N(x).
41  * \f}
42  *
43  * Note that when an exponential filter is being used to deal with aliasing,
44  * lower modes can be included in \f$\hat{U}\f$. The main goal of \f$\hat{U}\f$
45  * is to measure how much power is in the highest modes, which are the modes
46  * responsible for Gibbs phenomena. We define the discontinuity indicator
47  * \f$s^\Omega\f$ as
48  *
49  * \f{align*}{
50  * s^\Omega=\log_{10}\left(\frac{(\hat{U}, \hat{U})}{(U, U)}\right),
51  * \f}
52  *
53  * where \f$(\cdot,\cdot)\f$ is an inner product, which we take to be the
54  * Euclidean \f$L_2\f$ norm (i.e. we do not divide by the number of grid points
55  * since that cancels out anyway). A cell is troubled if
56  * \f$s^\Omega > -\alpha \log_{10}(N)\f$. Typically, \f$\alpha=4\f$ is a good
57  * choice.
58  *
59  * The parameter `zero_cutoff` is used to avoid division and logarithms of small
60  * numbers, which can be wildly fluctuating because of roundoff errors.
61  * We do not check the TCI for tensor components when \f$L_2(\hat{U}) \leq
62  * \epsilon L_2(U)\f$, where \f$\epsilon\f$ is the `zero_cutoff`. If all
63  * components are skipped the TCI returns `false`, i.e. the cell is not
64  * troubled.
65  */
66 template <size_t Dim, typename SymmList, typename IndexList>
67 bool persson_tci(const Tensor<DataVector, SymmList, IndexList>& tensor,
68  const Mesh<Dim>& dg_mesh, const double alpha,
69  const double zero_cutoff) noexcept {
70  DataVector filtered_component(dg_mesh.number_of_grid_points());
71  for (size_t component_index = 0; component_index < tensor.size();
72  ++component_index) {
73  if (detail::persson_tci_impl(make_not_null(&filtered_component),
74  tensor[component_index], dg_mesh, alpha,
75  zero_cutoff)) {
76  return true;
77  }
78  }
79  return false;
80 }
81 } // namespace evolution::dg::subcell
cstddef
Index.hpp
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
evolution::dg::subcell
Implementation of a generic finite volume/conservative finite difference subcell limiter.
Definition: Actions.hpp:6
evolution::dg::subcell::persson_tci
bool persson_tci(const Tensor< DataVector, SymmList, IndexList > &tensor, const Mesh< Dim > &dg_mesh, const double alpha, const double zero_cutoff) noexcept
Troubled cell indicator using spectral falloff of .
Definition: PerssonTci.hpp:67
Gsl.hpp
Tensor.hpp
Matrix.hpp
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13