Filtering.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <array>
7 #include <cstddef>
8 #include <string>
9 #include <tuple>
10 
12 #include "DataStructures/DataVector.hpp"
14 #include "Domain/Mesh.hpp"
15 #include "Domain/Tags.hpp"
16 #include "ErrorHandling/Error.hpp"
17 #include "NumericalAlgorithms/LinearOperators/ApplyMatrices.hpp"
18 #include "NumericalAlgorithms/Spectral/Filtering.hpp"
20 #include "Options/Options.hpp"
22 #include "Utilities/Gsl.hpp"
23 #include "Utilities/StdHelpers.hpp"
24 #include "Utilities/TMPL.hpp"
25 #include "Utilities/TypeTraits.hpp"
26 
27 namespace OptionTags {
28 /*!
29  * \ingroup OptionGroupsGroup
30  * \brief Groups the filtering configurations in the input file.
31  */
33  static std::string name() noexcept { return "Filtering"; }
34  static constexpr OptionString help = "Options for filtering";
35 };
36 } // namespace OptionTags
37 
38 namespace dg {
39 namespace Actions {
40 namespace ExponentialFilter_detail {
41 template <bool SameSize, bool SameList>
42 struct FilterAllEvolvedVars {
43  template <typename EvolvedVarsTagList, typename... FilterTags>
45 };
46 
47 template <>
48 struct FilterAllEvolvedVars<true, false> {
49  template <typename EvolvedVarsTagList, typename... FilterTags>
50  using f =
51  std::integral_constant<bool, tmpl2::flat_all_v<tmpl::list_contains_v<
52  EvolvedVarsTagList, FilterTags>...>>;
53 };
54 
55 template <>
56 struct FilterAllEvolvedVars<true, true> {
57  template <typename EvolvedVarsTagList, typename... FilterTags>
59 };
60 } // namespace ExponentialFilter_detail
61 
62 /// \cond
63 template <size_t FilterIndex, typename TagsToFilterList>
64 struct ExponentialFilter;
65 /// \endcond
66 
67 /*!
68  * \ingroup DiscontinuousGalerkinGroup
69  * \brief Applies an exponential filter to the specified tags.
70  *
71  * Applies an exponential filter in each logical direction to each component of
72  * the tensors `TagsToFilter`. The exponential filter rescales the 1d modal
73  * coefficients \f$c_i\f$ as:
74  *
75  * \f{align*}{
76  * c_i\to c_i \exp\left[-\alpha_{\mathrm{ef}}
77  * \left(\frac{i}{N}\right)^{2\beta_{\mathrm{ef}}}\right]
78  * \f}
79  *
80  * where \f$N\f$ is the basis order, \f$\alpha_{\mathrm{ef}}\f$ determines how
81  * much the coefficients are rescaled, and \f$\beta_{\mathrm{ef}}\f$ (given by
82  * the `HalfPower` option) determines how aggressive/broad the filter is (lower
83  * values means filtering more coefficients). Setting
84  * \f$\alpha_{\mathrm{ef}}=36\f$ results in effectively zeroing the highest
85  * coefficient (in practice it gets rescaled by machine epsilon). The same
86  * \f$\alpha_{\mathrm{ef}}\f$ and \f$\beta_{\mathrm{ef}}\f$ are used in each
87  * logical direction. For a discussion of filtering see section 5.3 of
88  * \cite HesthavenWarburton.
89  *
90  * If different `Alpha` or `HalfPower` parameters are desired for different tags
91  * then multiple `ExponentialFilter` actions must be inserted into the action
92  * list with different `FilterIndex` values. In the input file these will be
93  * specified as `ExpFilterFILTER_INDEX`, e.g. `ExpFilter0` and `ExpFilter1`.
94  * Below is an example of an action list with two different exponential filters:
95  *
96  * \snippet DiscontinuousGalerkin/Test_Filtering.cpp action_list_example
97  *
98  * #### Action properties
99  *
100  * Uses:
101  * - ConstGlobalCache:
102  * - `ExponentialFilter`
103  * - DataBox:
104  * - `Tags::Mesh`
105  * - DataBox changes:
106  * - Adds: nothing
107  * - Removes: nothing
108  * - Modifies:
109  * - `TagsToFilter`
110  * - System:
111  * - `volume_dim`
112  * - `variables_tag`
113  *
114  * #### Design decision:
115  *
116  * - The reason for the `size_t` template parameter is to allow for different
117  * `Alpha` and `HalfPower` parameters for different tensors while still being
118  * able to cache the matrices.
119  */
120 template <size_t FilterIndex, typename... TagsToFilter>
121 class ExponentialFilter<FilterIndex, tmpl::list<TagsToFilter...>> {
122  public:
123  using const_global_cache_tags = tmpl::list<ExponentialFilter>;
124 
125  /// \brief The value of `exp(-alpha)` is what the highest modal coefficient is
126  /// rescaled by.
127  struct Alpha {
128  using type = double;
129  static constexpr OptionString help =
130  "exp(-alpha) is rescaling of highest coefficient";
131  static type lower_bound() noexcept { return 0.0; }
132  };
133  /*!
134  * \brief Half of the exponent in the exponential.
135  *
136  * \f{align*}{
137  * c_i\to c_i \exp\left[-\alpha \left(\frac{i}{N}\right)^{2m}\right]
138  * \f}
139  */
140  struct HalfPower {
141  using type = unsigned;
142  static constexpr OptionString help =
143  "Half of the exponent in the generalized Gaussian";
144  static type lower_bound() noexcept { return 1; }
145  };
146  /// \brief Turn the filter off
147  ///
148  /// This option exists to temporarily disable the filter for debugging
149  /// purposes. For problems where filtering is not needed, the preferred
150  /// approach is to not compile the filter into the executable.
151  struct DisableForDebugging {
152  using type = bool;
153  static type default_value() noexcept { return false; }
154  static constexpr OptionString help = {"Disable the filter"};
155  };
156 
157  using type = ExponentialFilter;
158  static std::string name() noexcept {
159  return "ExpFilter" + std::to_string(FilterIndex);
160  }
162 
163  using options = tmpl::list<Alpha, HalfPower, DisableForDebugging>;
164  static constexpr OptionString help = {"An exponential filter."};
165 
166  ExponentialFilter() = default;
167 
168  ExponentialFilter(double alpha, unsigned half_power,
169  bool disable_for_debugging) noexcept;
170 
171  // Action part of the class
172  template <typename DbTags, typename... InboxTags, typename ArrayIndex,
173  typename ActionList, typename ParallelComponent,
174  typename Metavariables>
176  db::DataBox<DbTags>& box,
177  const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
179  const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
180  const ParallelComponent* const /*meta*/) noexcept {
181  constexpr size_t volume_dim = Metavariables::system::volume_dim;
182  using evolved_vars_tag = typename Metavariables::system::variables_tag;
183  using evolved_vars_tags_list = typename evolved_vars_tag::tags_list;
184  const ExponentialFilter& filter_helper =
185  Parallel::get<ExponentialFilter>(cache);
186  if (UNLIKELY(filter_helper.disable_for_debugging_)) {
187  return {std::move(box)};
188  }
189  const Mesh<volume_dim> mesh = db::get<::Tags::Mesh<volume_dim>>(box);
190  const Matrix empty{};
191  auto filter = make_array<volume_dim>(std::cref(empty));
192  for (size_t d = 0; d < volume_dim; d++) {
193  gsl::at(filter, d) =
194  std::cref(filter_helper.filter_matrix(mesh.slice_through(d)));
195  }
196 
197  // In the case that the tags we are filtering are all the evolved variables
198  // we filter the entire Variables at once to be more efficient. This case is
199  // the first branch of the `if-else`.
200  if (ExponentialFilter_detail::FilterAllEvolvedVars<
201  sizeof...(TagsToFilter) ==
202  tmpl::size<evolved_vars_tags_list>::value,
203  cpp17::is_same_v<evolved_vars_tags_list,
204  tmpl::list<TagsToFilter...>>>::
205  template f<evolved_vars_tags_list, TagsToFilter...>::value) {
206  db::mutate<typename Metavariables::system::variables_tag>(
207  make_not_null(&box),
208  [&filter](
209  const gsl::not_null<
211  vars,
212  const auto& local_mesh) noexcept {
213  *vars = apply_matrices(filter, *vars, local_mesh.extents());
214  },
215  mesh);
216  } else {
217  db::mutate<TagsToFilter...>(
218  make_not_null(&box),
219  [&filter](const gsl::not_null<
220  db::item_type<TagsToFilter>*>... tensors_to_filter,
221  const auto& local_mesh) noexcept {
222  DataVector temp(local_mesh.number_of_grid_points(), 0.0);
223  const auto helper =
224  [&local_mesh, &filter, &temp ](const auto tensor) noexcept {
225  for (auto& component : *tensor) {
226  temp = 0.0;
227  apply_matrices(make_not_null(&temp), filter, component,
228  local_mesh.extents());
229  component = temp;
230  }
231  };
232  EXPAND_PACK_LEFT_TO_RIGHT(helper(tensors_to_filter));
233  },
234  mesh);
235  }
236  return {std::move(box)};
237  }
238 
239  // NOLINTNEXTLINE(google-runtime-references)
240  void pup(PUP::er& p) noexcept;
241 
242  bool operator==(const ExponentialFilter& rhs) const noexcept;
243 
244  private:
245  const Matrix& filter_matrix(const Mesh<1>& mesh) const noexcept {
246  // All these switch gymnastics are to translate the runtime mesh values into
247  // compile time template parameters. This is used so we have a sparse lazy
248  // cache where we only cache matrices that we actually need for the meshes
249  // used in the simulation.
250  switch (mesh.basis(0)) {
251  case Spectral::Basis::Legendre:
252  switch (mesh.quadrature(0)) {
253  case Spectral::Quadrature::Gauss:
254  return filter_matrix_impl<Spectral::Basis::Legendre,
255  Spectral::Quadrature::Gauss>(
256  mesh,
258  Spectral::Basis::Legendre>>{});
259  case Spectral::Quadrature::GaussLobatto:
260  return filter_matrix_impl<Spectral::Basis::Legendre,
261  Spectral::Quadrature::GaussLobatto>(
262  mesh,
264  Spectral::Basis::Legendre>>{});
265  default:
266  ERROR("Missing quadrature in exponential filter matrix action.");
267  }
268  case Spectral::Basis::Chebyshev:
269  switch (mesh.quadrature(0)) {
270  case Spectral::Quadrature::Gauss:
271  return filter_matrix_impl<Spectral::Basis::Chebyshev,
272  Spectral::Quadrature::Gauss>(
273  mesh,
275  Spectral::Basis::Chebyshev>>{});
276  case Spectral::Quadrature::GaussLobatto:
277  return filter_matrix_impl<Spectral::Basis::Chebyshev,
278  Spectral::Quadrature::GaussLobatto>(
279  mesh.slice_through(0),
281  Spectral::Basis::Chebyshev>>{});
282  default:
283  ERROR("Missing quadrature in exponential filter matrix action.");
284  }
285  default:
286  ERROR("Missing basis in exponential filter matrix action.");
287  };
288  }
289 
290  template <Spectral::Basis BasisType, Spectral::Quadrature QuadratureType,
291  size_t I>
292  static const Matrix& filter_matrix_cache(const Mesh<1>& mesh,
293  const double alpha,
294  const unsigned half_power) noexcept {
295  static Matrix matrix =
296  Spectral::filtering::exponential_filter(mesh, alpha, half_power);
297  return matrix;
298  }
299 
300  template <Spectral::Basis BasisType, Spectral::Quadrature QuadratureType,
301  size_t... Is>
302  const Matrix& filter_matrix_impl(const Mesh<1>& mesh,
303  std::index_sequence<Is...> /*meta*/) const
304  noexcept {
305  static const std::array<const Matrix& (*)(const Mesh<1>&, double, unsigned),
306  sizeof...(Is)>
307  cache{{&filter_matrix_cache<BasisType, QuadratureType, Is>...}};
308  return gsl::at(cache, mesh.extents(0))(mesh, alpha_, half_power_);
309  }
310 
311  double alpha_{36.0};
312  unsigned half_power_{16};
313  bool disable_for_debugging_{false};
314 };
315 
316 template <size_t FilterIndex, typename... TagsToFilter>
317 ExponentialFilter<FilterIndex, tmpl::list<TagsToFilter...>>::ExponentialFilter(
318  const double alpha, const unsigned half_power,
319  const bool disable_for_debugging) noexcept
320  : alpha_(alpha),
321  half_power_(half_power),
322  disable_for_debugging_(disable_for_debugging) {}
323 
324 template <size_t FilterIndex, typename... TagsToFilter>
325 void ExponentialFilter<FilterIndex, tmpl::list<TagsToFilter...>>::pup(
326  PUP::er& p) noexcept {
327  p | alpha_;
328  p | half_power_;
329  p | disable_for_debugging_;
330 }
331 
332 template <size_t FilterIndex, typename... TagsToFilter>
333 bool ExponentialFilter<FilterIndex, tmpl::list<TagsToFilter...>>::operator==(
334  const ExponentialFilter& rhs) const noexcept {
335  return alpha_ == rhs.alpha_ and half_power_ == rhs.half_power_ and
336  disable_for_debugging_ == rhs.disable_for_debugging_;
337 }
338 
339 template <size_t FilterIndex, typename TagsToFilterList>
340 bool operator!=(
341  const ExponentialFilter<FilterIndex, TagsToFilterList>& lhs,
342  const ExponentialFilter<FilterIndex, TagsToFilterList>& rhs) noexcept {
343  return not(lhs == rhs);
344 }
345 } // namespace Actions
346 } // namespace dg
Quadrature
The choice of quadrature method to compute integration weights.
Definition: Spectral.hpp:68
Defines class Matrix.
Defines helper functions for the standard library.
const std::array< Spectral::Quadrature, Dim > & quadrature() const noexcept
The quadrature chosen in each dimension of the grid.
Definition: Mesh.hpp:155
#define ERROR(m)
prints an error message to the standard error stream and aborts the program.
Definition: Error.hpp:35
void mutate(const gsl::not_null< DataBox< TagList > *> box, Invokable &&invokable, Args &&... args) noexcept
Allows changing the state of one or more non-computed elements in the DataBox.
Definition: DataBox.hpp:1099
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:561
constexpr bool flat_all_v
A non-short-circuiting logical AND between bools &#39;B"".
Definition: TMPL.hpp:504
Definition: Digraph.hpp:11
#define UNLIKELY(x)
Definition: Gsl.hpp:72
Definition: Filtering.hpp:38
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
void apply_matrices(const gsl::not_null< Variables< VariableTags > *> result, const std::array< MatrixType, Dim > &matrices, const Variables< VariableTags > &u, const Index< Dim > &extents) noexcept
Multiply by matrices in each dimension.
Definition: ApplyMatrices.hpp:62
Defines classes and functions for making classes creatable from input files.
constexpr size_t maximum_number_of_points
Maximum number of allowed collocation points.
Definition: Spectral.hpp:96
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1609
Mesh< sizeof...(D)> slice_through(D... d) const noexcept
Returns a Mesh with the dimensions d, ... present (zero-indexed).
Definition: Mesh.hpp:186
Groups the filtering configurations in the input file.
Definition: Filtering.hpp:32
Matrix exponential_filter(const Mesh< 1 > &mesh, const double alpha, const unsigned half_power) noexcept
Returns a Matrix by which to multiply the nodal coefficients to apply a stable exponential filter...
Definition: Filtering.cpp:15
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:27
Basis
The choice of basis functions for computing collocation points and weights.
Definition: Spectral.hpp:53
An associative container that is indexed by structs.
Definition: TaggedTuple.hpp:272
Defines classes and functions used for manipulating DataBox&#39;s.
constexpr bool is_same_v
Variable template for is_same.
Definition: TypeTraits.hpp:221
Definition: Strahlkorper.hpp:167
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
Definition: InterpolationTargetWedgeSectionTorus.hpp:24
A Charm++ chare that caches constant data once per Charm++ node.
Definition: ConstGlobalCache.hpp:76
Defines the class template Mesh.
const Index< Dim > & extents() const noexcept
The number of grid points in each dimension of the grid.
Definition: Mesh.hpp:107
Stores a collection of function values.
Definition: DataVector.hpp:42
Wraps the template metaprogramming library used (brigand)
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that is returned by the Tag. If it is a base tag then a TagList must be passed as a seco...
Definition: DataBoxTag.hpp:410
Defines functions and classes from the GSL.
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, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:863
Defines tags related to domain quantities.
Definition: SolvePoissonProblem.hpp:37
Defines class template ConstGlobalCache.
Defines macro ERROR.
Defines type traits, some of which are future STL type_traits header.
const std::array< Spectral::Basis, Dim > & basis() const noexcept
The basis chosen in each dimension of the grid.
Definition: Mesh.hpp:141
Definition: ComputeTimeDerivative.hpp:28
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid...
Definition: Gsl.hpp:124