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             : 
       8             : #include "Domain/Structure/DirectionMap.hpp"
       9             : #include "Utilities/Gsl.hpp"
      10             : 
      11             : /// \cond
      12             : template <size_t Dim>
      13             : class Mesh;
      14             : /// \endcond
      15             : 
      16             : namespace fd {
      17             : /*!
      18             :  * \brief Apply a low-pass filter to the data.
      19             :  *
      20             :  * The filter fits a Legendre polynomial of degree equal to `fd_order` to each
      21             :  * grid point and its neighboring points, then subtracts out the highest mode
      22             :  * contribution at the grid point. This is inspired by a Heaviside filter for
      23             :  * Legendre polynomial-based spectral methods.
      24             :  *
      25             :  * The filter at different orders is given by:
      26             :  *
      27             :  * \f{align*}{
      28             :  * F^{(9)} u_i&= -\frac{35 \times 531441}{128} \left(
      29             :  *                \frac{1}{6406400} (u_{i-4} + u_{i+4}) -
      30             :  *                \frac{1}{800800} (u_{i-3} + u_{i+3}) +
      31             :  *                \frac{1}{228800} (u_{i-2} + u_{i+2}) -
      32             :  *                \frac{1}{114400} (u_{i-1} + u_{i+1}) +
      33             :  *                \frac{1}{91520} u_i\right) \\
      34             :  * F^{(7)} u_i&= \frac{5 \times 16807}{16} \left(
      35             :  *                \frac{1}{95040} (u_{i-3} + u_{i+3}) -
      36             :  *                \frac{1}{15840} (u_{i-2} + u_{i+2}) +
      37             :  *                \frac{1}{6336} (u_{i-1} + u_{i+1}) -
      38             :  *                \frac{1}{4572} u_i\right) \\
      39             :  * F^{(5)} u_i&= -\frac{3 \times 125}{8} \left(
      40             :  *                \frac{1}{336} (u_{i-2} + u_{i+2}) -
      41             :  *                \frac{1}{84} (u_{i-1} + u_{i+1}) +
      42             :  *                \frac{1}{56} u_i\right) \\
      43             :  * F^{(3)} u_i&= \frac{1 \times 3}{2} \left(
      44             :  *                \frac{1}{4} (u_{i-1} + u_{i+1}) -
      45             :  *                \frac{1}{2} u_i\right)
      46             :  * \f}
      47             :  *
      48             :  * \note The \f$F^{(11)}\f$ filter isn't implemented yet.
      49             :  *
      50             :  * \note The argument \f$\epsilon\f$ controls how much of the mode is filter
      51             :  * out. \f$\epsilon=1\f$ means the highest mode is completely filtered while
      52             :  * \f$\epsilon=0\f$ means it's not at all filtered.
      53             :  */
      54             : template <size_t Dim>
      55           1 : void low_pass_filter(
      56             :     gsl::not_null<gsl::span<double>*> filtered_data,
      57             :     const gsl::span<const double>& volume_vars,
      58             :     const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
      59             :     const Mesh<Dim>& volume_mesh, size_t number_of_variables, size_t fd_order,
      60             :     double epsilon);
      61             : 
      62             : /*!
      63             :  * \brief Apply Kreiss-Oliger dissipation \cite Kreiss1973 of order `fd_order`
      64             :  * to the variables.
      65             :  *
      66             :  * Define the operators \f$D_+\f$ and \f$D_-\f$ be defined as:
      67             :  *
      68             :  * \f{align*}{
      69             :  * D_+f_i&=\frac{(f_{i+1}-f_i)}{\Delta x} \\
      70             :  * D_-f_i&=\frac{(f_i-f_{i-1})}{\Delta x}
      71             :  * \f}
      72             :  *
      73             :  * where the subscript \f$i\f$ refers to the grid index. The dissipation
      74             :  * operators are generally applied dimension-by-dimension, and so we have:
      75             :  *
      76             :  * \f{align*}{
      77             :  * \mathcal{D}^{(2m)}=-\frac{(-1)^m}{2^{2m}}\Delta x^{2m-1}
      78             :  * \epsilon(D_{+})^m(D_{-})^m
      79             :  * \f}
      80             :  *
      81             :  * where \f$\epsilon\f$ controls the amount of dissipation and is restricted to
      82             :  * \f$0\leq\epsilon\leq1\f$, and \f$m\f$ is the order of the finite difference
      83             :  * derivative so as not to spoil the accuracy of the scheme. That is, for
      84             :  * second order FD, \f$m=2\f$ and one should use \f$\mathcal{D}^{(4)}\f$.
      85             :  * However, this choice requires a larger stencil and whether or not this is
      86             :  * necessary also depends on when and how in the algorithm the operators are
      87             :  * applied.
      88             :  *
      89             :  * We arrive at the following operators:
      90             :  *
      91             :  * \f{align*}{
      92             :  *  \mathcal{D}^{(2)}f_i&=-\frac{\epsilon}{\Delta x}
      93             :  *   (f_{i+1} - 2f_i+f_{i-1}) \\
      94             :  *  \mathcal{D}^{(4)}f_i&=-\frac{\epsilon}{16\Delta x}
      95             :  *   (f_{i+2}-4f_{i+1}+6f_i-4f_{i-1}+f_{i-2}) \\
      96             :  *  \mathcal{D}^{(6)}f_i&=\frac{\epsilon}{64\Delta x}(f_{i+3}-6f_{i+2}+15f_{i+1}
      97             :  *    -20f_i+15f_{i-1}-6f_{i-2}+f_{i-3}) \\
      98             :  *  \mathcal{D}^{(8)}f_i&=-\frac{\epsilon}{256\Delta x}
      99             :  *    (f_{i+4}-8f_{i+3}+28f_{i+2}-56f_{i+1}+70f_{i}-56f_{i-1}+28f_{i-2}
     100             :  *    -8f_{i-3}+f_{i-4}) \\
     101             :  *  \mathcal{D}^{(10)}f_i&=\frac{\epsilon}{1024\Delta x}
     102             :  *    (f_{i+5}-10f_{i+4}+45f_{i+3}
     103             :  *    -120f_{i+2}+210f_{i+1}-252f_{i}+210f_{i-1}-120f_{i-2}+45f_{i-3}
     104             :  *    -10f_{i-4}+f_{i-5})
     105             :  * \f}
     106             :  *
     107             :  * \note This function applies \f$\Delta x \mathcal{D}^{(2m)}\f$.
     108             :  */
     109             : template <size_t Dim>
     110           1 : void kreiss_oliger_filter(
     111             :     gsl::not_null<gsl::span<double>*> filtered_data,
     112             :     const gsl::span<const double>& volume_vars,
     113             :     const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
     114             :     const Mesh<Dim>& volume_mesh, size_t number_of_variables, size_t fd_order,
     115             :     double epsilon);
     116             : }  // namespace fd