ExponentialFilter.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <pup.h>
8 #include <string>
9 
10 #include "Options/Options.hpp"
11 #include "Utilities/TMPL.hpp"
12 
13 /// \cond
14 class Matrix;
15 template <size_t Dim>
16 class Mesh;
17 
18 /// \endcond
19 
20 namespace Filters {
21 
22 /*!
23  * \ingroup DiscontinuousGalerkinGroup
24  * \brief A cached exponential filter.
25  *
26  * Applies an exponential filter in each logical direction to each component of
27  * the tensors `TagsToFilter`. The exponential filter rescales the 1d modal
28  * coefficients \f$c_i\f$ as:
29  *
30  * \f{align*}{
31  * c_i\to c_i \exp\left[-\alpha_{\mathrm{ef}}
32  * \left(\frac{i}{N}\right)^{2\beta_{\mathrm{ef}}}\right]
33  * \f}
34  *
35  * where \f$N\f$ is the basis order, \f$\alpha_{\mathrm{ef}}\f$ determines how
36  * much the coefficients are rescaled, and \f$\beta_{\mathrm{ef}}\f$ (given by
37  * the `HalfPower` option) determines how aggressive/broad the filter is (lower
38  * values means filtering more coefficients). Setting
39  * \f$\alpha_{\mathrm{ef}}=36\f$ results in effectively zeroing the highest
40  * coefficient (in practice it gets rescaled by machine epsilon). The same
41  * \f$\alpha_{\mathrm{ef}}\f$ and \f$\beta_{\mathrm{ef}}\f$ are used in each
42  * logical direction. For a discussion of filtering see section 5.3 of
43  * \cite HesthavenWarburton.
44  *
45  * #### Design decision:
46  *
47  * - The reason for the `size_t` template parameter is to allow for different
48  * `Alpha` and `HalfPower` parameters for different tensors while still being
49  * able to cache the matrices. If different `Alpha` or `HalfPower` parameters
50  * are desired for filtering different tensors, then multiple filters must be
51  * inserted into the ConstGlobalCache with different `FilterIndex` values. In
52  * the input file these will be specified as `ExpFilterFILTER_INDEX`, e.g.
53  * - Filtering:
54  * - ExpFilter0:
55  * - Alpha: 12
56  * - HalfPower: 32
57  * - ExpFilter1:
58  * - Alpha: 36
59  * - HalfPower: 32
60  */
61 template <size_t FilterIndex>
62 class Exponential {
63  public:
64  /// \brief The value of `exp(-alpha)` is what the highest modal coefficient is
65  /// rescaled by.
66  struct Alpha {
67  using type = double;
68  static constexpr OptionString help =
69  "exp(-alpha) is rescaling of highest coefficient";
70  static type lower_bound() noexcept { return 0.0; }
71  };
72 
73  /*!
74  * \brief Half of the exponent in the exponential.
75  *
76  * \f{align*}{
77  * c_i\to c_i \exp\left[-\alpha \left(\frac{i}{N}\right)^{2m}\right]
78  * \f}
79  */
80  struct HalfPower {
81  using type = unsigned;
82  static constexpr OptionString help =
83  "Half of the exponent in the generalized Gaussian";
84  static type lower_bound() noexcept { return 1; }
85  };
86 
87  /// \brief Turn the filter off
88  ///
89  /// This option exists to temporarily disable the filter for debugging
90  /// purposes. For problems where filtering is not needed, the preferred
91  /// approach is to not compile the filter into the executable.
93  using type = bool;
94  static type default_value() noexcept { return false; }
95  static constexpr OptionString help = {"Disable the filter"};
96  };
97 
98  using options = tmpl::list<Alpha, HalfPower, DisableForDebugging>;
99  static constexpr OptionString help = {"An exponential filter."};
100  static std::string name() noexcept {
101  return "ExpFilter" + std::to_string(FilterIndex);
102  }
103 
104  Exponential() = default;
105 
106  Exponential(double alpha, unsigned half_power,
107  bool disable_for_debugging) noexcept;
108 
109  /// A cached matrix used to apply the filter to the given mesh
110  const Matrix& filter_matrix(const Mesh<1>& mesh) const noexcept;
111 
112  bool disable_for_debugging() const noexcept { return disable_for_debugging_; }
113 
114  // NOLINTNEXTLINE(google-runtime-references)
115  void pup(PUP::er& p) noexcept;
116 
117  private:
118  template <size_t LocalFilterIndex>
119  // NOLINTNEXTLINE(google-runtime-references)
120  friend bool operator==(const Exponential<LocalFilterIndex>& lhs,
121  const Exponential<LocalFilterIndex>& rhs) noexcept;
122 
123  double alpha_{36.0};
124  unsigned half_power_{16};
125  bool disable_for_debugging_{false};
126 };
127 
128 template <size_t LocalFilterIndex>
129 bool operator==(const Exponential<LocalFilterIndex>& lhs,
130  const Exponential<LocalFilterIndex>& rhs) noexcept;
131 
132 template <size_t FilterIndex>
133 bool operator!=(const Exponential<FilterIndex>& lhs,
134  const Exponential<FilterIndex>& rhs) noexcept;
135 } // namespace Filters
const Matrix & filter_matrix(const Mesh< 1 > &mesh) const noexcept
A cached matrix used to apply the filter to the given mesh.
Half of the exponent in the exponential.
Definition: ExponentialFilter.hpp:80
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
Turn the filter off.
Definition: ExponentialFilter.hpp:92
Defines classes and functions for making classes creatable from input files.
A cached exponential filter.
Definition: ExponentialFilter.hpp:62
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:29
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
Wraps the template metaprogramming library used (brigand)
Definition: ExponentialFilter.hpp:20
The value of exp(-alpha) is what the highest modal coefficient is rescaled by.
Definition: ExponentialFilter.hpp:66