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