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 : * Dimensions with Cartoon bases are not filtered as there is nothing to filter.
49 : *
50 : * \note The \f$F^{(11)}\f$ filter isn't implemented yet.
51 : *
52 : * \note The argument \f$\epsilon\f$ controls how much of the mode is filter
53 : * out. \f$\epsilon=1\f$ means the highest mode is completely filtered while
54 : * \f$\epsilon=0\f$ means it's not at all filtered.
55 : */
56 : template <size_t Dim>
57 1 : void low_pass_filter(
58 : gsl::not_null<gsl::span<double>*> filtered_data,
59 : const gsl::span<const double>& volume_vars,
60 : const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
61 : const Mesh<Dim>& volume_mesh, size_t number_of_variables, size_t fd_order,
62 : double epsilon);
63 :
64 : /*!
65 : * \brief Apply Kreiss-Oliger dissipation \cite Kreiss1973 of order `fd_order`
66 : * to the variables.
67 : *
68 : * Define the operators \f$D_+\f$ and \f$D_-\f$ be defined as:
69 : *
70 : * \f{align*}{
71 : * D_+f_i&=\frac{(f_{i+1}-f_i)}{\Delta x} \\
72 : * D_-f_i&=\frac{(f_i-f_{i-1})}{\Delta x}
73 : * \f}
74 : *
75 : * where the subscript \f$i\f$ refers to the grid index. The dissipation
76 : * operators are generally applied dimension-by-dimension, and so we have:
77 : *
78 : * \f{align*}{
79 : * \mathcal{D}^{(2m)}=-\frac{(-1)^m}{2^{2m}}\Delta x^{2m-1}
80 : * \epsilon(D_{+})^m(D_{-})^m
81 : * \f}
82 : *
83 : * where \f$\epsilon\f$ controls the amount of dissipation and is restricted to
84 : * \f$0\leq\epsilon\leq1\f$, and \f$m\f$ is the order of the finite difference
85 : * derivative so as not to spoil the accuracy of the scheme. That is, for
86 : * second order FD, \f$m=2\f$ and one should use \f$\mathcal{D}^{(4)}\f$.
87 : * However, this choice requires a larger stencil and whether or not this is
88 : * necessary also depends on when and how in the algorithm the operators are
89 : * applied.
90 : *
91 : * We arrive at the following operators:
92 : *
93 : * \f{align*}{
94 : * \mathcal{D}^{(2)}f_i&=-\frac{\epsilon}{\Delta x}
95 : * (f_{i+1} - 2f_i+f_{i-1}) \\
96 : * \mathcal{D}^{(4)}f_i&=-\frac{\epsilon}{16\Delta x}
97 : * (f_{i+2}-4f_{i+1}+6f_i-4f_{i-1}+f_{i-2}) \\
98 : * \mathcal{D}^{(6)}f_i&=\frac{\epsilon}{64\Delta x}(f_{i+3}-6f_{i+2}+15f_{i+1}
99 : * -20f_i+15f_{i-1}-6f_{i-2}+f_{i-3}) \\
100 : * \mathcal{D}^{(8)}f_i&=-\frac{\epsilon}{256\Delta x}
101 : * (f_{i+4}-8f_{i+3}+28f_{i+2}-56f_{i+1}+70f_{i}-56f_{i-1}+28f_{i-2}
102 : * -8f_{i-3}+f_{i-4}) \\
103 : * \mathcal{D}^{(10)}f_i&=\frac{\epsilon}{1024\Delta x}
104 : * (f_{i+5}-10f_{i+4}+45f_{i+3}
105 : * -120f_{i+2}+210f_{i+1}-252f_{i}+210f_{i-1}-120f_{i-2}+45f_{i-3}
106 : * -10f_{i-4}+f_{i-5})
107 : * \f}
108 : *
109 : * Dimensions with Cartoon bases are not filtered as there is nothing to filter.
110 : *
111 : * \note This function applies \f$\Delta x \mathcal{D}^{(2m)}\f$.
112 : */
113 : template <size_t Dim>
114 1 : void kreiss_oliger_filter(
115 : gsl::not_null<gsl::span<double>*> filtered_data,
116 : const gsl::span<const double>& volume_vars,
117 : const DirectionMap<Dim, gsl::span<const double>>& ghost_cell_vars,
118 : const Mesh<Dim>& volume_mesh, size_t number_of_variables, size_t fd_order,
119 : double epsilon);
120 : } // namespace fd
|
