Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <optional>
       8             : #include <pup.h>
       9             : #include <string>
      10             : #include <unordered_set>
      11             : 
      12             : #include "NumericalAlgorithms/Spectral/Parity.hpp"
      13             : #include "Options/Auto.hpp"
      14             : #include "Options/String.hpp"
      15             : #include "Utilities/TMPL.hpp"
      16             : 
      17             : /// \cond
      18             : class Matrix;
      19             : template <size_t Dim>
      20             : class Mesh;
      21             : 
      22             : /// \endcond
      23             : 
      24             : namespace Filters {
      25             : 
      26             : /*!
      27             :  * \ingroup DiscontinuousGalerkinGroup
      28             :  * \brief A cached exponential filter.
      29             :  *
      30             :  * Applies an exponential filter in each logical direction to each component of
      31             :  * the tensors `TagsToFilter`. The exponential filter rescales the 1d modal
      32             :  * coefficients \f$c_i\f$ as:
      33             :  *
      34             :  * \f{align*}{
      35             :  *  c_i\to c_i \exp\left[-\alpha_{\mathrm{ef}}
      36             :  *   \left(\frac{i}{N}\right)^{2\beta_{\mathrm{ef}}}\right]
      37             :  * \f}
      38             :  *
      39             :  * where \f$N\f$ is the basis degree (number of grid points per element per
      40             :  * dimension minus one), \f$\alpha_{\mathrm{ef}}\f$ determines how much the
      41             :  * coefficients are rescaled, and \f$\beta_{\mathrm{ef}}\f$ (given by the
      42             :  * `HalfPower` option) determines how aggressive/broad the filter is (lower
      43             :  * values means filtering more coefficients). Setting
      44             :  * \f$\alpha_{\mathrm{ef}}=36\f$ results in effectively zeroing the highest
      45             :  * coefficient (in practice it gets rescaled by machine epsilon). The same
      46             :  * \f$\alpha_{\mathrm{ef}}\f$ and \f$\beta_{\mathrm{ef}}\f$ are used in each
      47             :  * logical direction. For a discussion of filtering see section 5.3 of
      48             :  * \cite HesthavenWarburton.
      49             :  *
      50             :  * #### Design decision:
      51             :  *
      52             :  * - The reason for the `size_t` template parameter is to allow for different
      53             :  * `Alpha` and `HalfPower` parameters for different tensors while still being
      54             :  * able to cache the matrices.   If different `Alpha` or `HalfPower` parameters
      55             :  * are desired for filtering different tensors, then multiple filters must be
      56             :  * inserted into the GlobalCache with different `FilterIndex` values. In
      57             :  * the input file these will be specified as `ExpFilterFILTER_INDEX`, e.g.
      58             :  * \snippet LinearOperators/Test_Filtering.cpp multiple_exponential_filters
      59             :  */
      60             : template <size_t FilterIndex>
      61           1 : class Exponential {
      62             :  public:
      63             :   /// \brief The value of `exp(-alpha)` is what the highest modal coefficient is
      64             :   /// rescaled by.
      65           1 :   struct Alpha {
      66           0 :     using type = double;
      67           0 :     static constexpr Options::String help =
      68             :         "exp(-alpha) is rescaling of highest coefficient";
      69           0 :     static type lower_bound() { return 0.0; }
      70             :   };
      71             : 
      72             :   /*!
      73             :    * \brief Half of the exponent in the exponential.
      74             :    *
      75             :    * \f{align*}{
      76             :    *  c_i\to c_i \exp\left[-\alpha \left(\frac{i}{N}\right)^{2m}\right]
      77             :    * \f}
      78             :    */
      79           1 :   struct HalfPower {
      80           0 :     using type = unsigned;
      81           0 :     static constexpr Options::String help =
      82             :         "Half of the exponent in the generalized Gaussian";
      83           0 :     static type lower_bound() { return 1; }
      84             :   };
      85             : 
      86             :   /// \brief Turn the filter off
      87           1 :   struct Enable {
      88           0 :     using type = bool;
      89           0 :     static constexpr Options::String help = {"Enable the filter"};
      90             :   };
      91             : 
      92           0 :   struct BlocksToFilter {
      93           0 :     using type =
      94             :         Options::Auto<std::vector<std::string>, Options::AutoLabel::All>;
      95           0 :     static constexpr Options::String help = {
      96             :         "List of blocks or block groups to apply filtering to. All other "
      97             :         "blocks will have no filtering. You can also specify 'All' to do "
      98             :         "filtering in all blocks of the domain."};
      99             :   };
     100             : 
     101           0 :   using options = tmpl::list<Alpha, HalfPower, Enable, BlocksToFilter>;
     102           0 :   static constexpr Options::String help = {"An exponential filter."};
     103           0 :   static std::string name() {
     104             :     return "ExpFilter" + std::to_string(FilterIndex);
     105             :   }
     106             : 
     107           0 :   Exponential() = default;
     108             : 
     109           0 :   Exponential(double alpha, unsigned half_power, bool enable,
     110             :               const std::optional<std::vector<std::string>>& blocks_to_filter,
     111             :               const Options::Context& context = {});
     112             : 
     113             :   /// A cached matrix used to apply the filter to the given mesh
     114           1 :   const Matrix& filter_matrix(
     115             :       const Mesh<1>& mesh,
     116             :       Spectral::Parity parity = Spectral::Parity::Uninitialized) const;
     117             : 
     118           0 :   bool enable() const { return enable_; }
     119             : 
     120           0 :   const std::optional<std::unordered_set<std::string>>& blocks_to_filter()
     121             :       const {
     122             :     return blocks_to_filter_;
     123             :   }
     124             : 
     125             :   // NOLINTNEXTLINE(google-runtime-references)
     126           0 :   void pup(PUP::er& p);
     127             : 
     128             :  private:
     129             :   template <size_t LocalFilterIndex>
     130             :   // NOLINTNEXTLINE(readability-redundant-declaration)
     131           0 :   friend bool operator==(const Exponential<LocalFilterIndex>& lhs,
     132             :                          const Exponential<LocalFilterIndex>& rhs);
     133             : 
     134           0 :   double alpha_{36.0};
     135           0 :   unsigned half_power_{16};
     136           0 :   bool enable_{true};
     137           0 :   std::optional<std::unordered_set<std::string>> blocks_to_filter_{};
     138             : };
     139             : 
     140             : template <size_t LocalFilterIndex>
     141           0 : bool operator==(const Exponential<LocalFilterIndex>& lhs,
     142             :                 const Exponential<LocalFilterIndex>& rhs);
     143             : 
     144             : template <size_t FilterIndex>
     145           0 : bool operator!=(const Exponential<FilterIndex>& lhs,
     146             :                 const Exponential<FilterIndex>& rhs);
     147             : }  // namespace Filters