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_map>
11 : #include <unordered_set>
12 : #include <vector>
13 :
14 : #include "DataStructures/Tensor/TypeAliases.hpp"
15 : #include "NumericalAlgorithms/LinearOperators/Filters/Filter.hpp"
16 : #include "Options/Auto.hpp"
17 : #include "Options/Context.hpp"
18 : #include "Options/String.hpp"
19 : #include "Utilities/Gsl.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : /// \cond
23 : class DataVector;
24 : template <size_t Dim>
25 : class Mesh;
26 : class Matrix;
27 : template <typename TagsList>
28 : class Variables;
29 : /// \endcond
30 :
31 : namespace Filters {
32 : /*!
33 : * \ingroup DiscontinuousGalerkinGroup
34 : * \brief An exponential spectral filter applied in each logical direction of
35 : * a tensor-product (line, square, cube, ...) element.
36 : *
37 : * Concrete implementation of `Filters::Filter` for tensor-product
38 : * (hypercube) elements, driven by the DG filtering action. See
39 : * `Filters::Filter` for the framing of volume vs. boundary application, the
40 : * substep / every-N-steps cadence controls, and the `blocks_to_filter`
41 : * semantics.
42 : *
43 : * For each component of the tensors in `TagList`, the filter rescales the
44 : * 1-D modal coefficients \f$c_i\f$ in each logical direction as
45 : *
46 : * \f{align*}{
47 : * c_i \to c_i \exp\!\left[-36 \left(\frac{i}{N}\right)^{2m}\right],
48 : * \f}
49 : *
50 : * where \f$N\f$ is the basis degree (number of grid points per element per
51 : * dimension minus one) and \f$m\f$ is the `HalfPower` option. The same
52 : * coefficient and `HalfPower` are used in every logical direction. With the
53 : * fixed coefficient 36 the highest mode is rescaled by approximately machine
54 : * epsilon, i.e. effectively zeroed. For a discussion of filtering see
55 : * section 5.3 of \cite HesthavenWarburton.
56 : *
57 : * #### Design decision:
58 : *
59 : * The exponential coefficient is hardcoded to 36 since this is what has
60 : * worked well in practice for several decades in SpEC. If we ever use
61 : * quad or double-double types, we may want to try 72, but that is unlikely to
62 : * be necessary since 36 decreases the highest coefficient by 1e-16. I.e.,
63 : * this is not relative to the largest coefficient.
64 : */
65 : template <size_t Dim, typename TagList>
66 1 : class Hypercube : public Filter<Dim, TagList> {
67 : public:
68 : /*!
69 : * \brief Half of the exponent in the exponential.
70 : *
71 : * I.e., this is \f$m\f$ in
72 : *
73 : * \f{align*}{
74 : * c_i\to c_i \exp\left[-\alpha \left(\frac{i}{N}\right)^{2m}\right]
75 : * \f}
76 : */
77 1 : struct HalfPower {
78 0 : using type = unsigned;
79 0 : static constexpr Options::String help =
80 : "Half of the exponent in the generalized Gaussian";
81 0 : static type lower_bound() { return 1; }
82 : };
83 :
84 : /// \brief Enable the filter
85 1 : struct Enable {
86 0 : using type = bool;
87 0 : static constexpr Options::String help = {"Enable the filter"};
88 : };
89 :
90 : /// \brief Which blocks and block groups the filter should be applied to.
91 1 : 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 that are hypercubes."};
98 : };
99 :
100 : /// \brief Apply the volume filter inside every substep instead of only at
101 : /// step boundaries.
102 1 : struct VolumeFilterOnSubstep {
103 0 : using type = bool;
104 0 : static constexpr Options::String help = {
105 : "Enable the volume filter on every substep."};
106 : };
107 :
108 : /// \brief Apply the boundary correction filter inside every substep instead
109 : /// of only at step boundaries.
110 1 : struct BoundaryCorrectionFilterOnSubstep {
111 0 : using type = bool;
112 0 : static constexpr Options::String help = {
113 : "Enable the boundary filter on every substep."};
114 : };
115 :
116 : /// \brief Apply the volume filter once every `N` steps. `None`
117 : /// (`std::nullopt`) disables the every-N-steps trigger.
118 1 : struct VolumeFilterEveryNSteps {
119 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
120 0 : static constexpr Options::String help = {
121 : "Enable the volume filter on every N steps. 'None' to disable."};
122 : };
123 :
124 : /// \brief Apply the boundary correction filter once every `N` steps. `None`
125 : /// (`std::nullopt`) disables the every-N-steps trigger.
126 1 : struct BoundaryCorrectionFilterEveryNSteps {
127 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
128 0 : static constexpr Options::String help = {
129 : "Enable the boundary filter on every N steps. 'None' to disable."};
130 : };
131 :
132 0 : using options =
133 : tmpl::list<HalfPower, Enable, BlocksToFilter, VolumeFilterOnSubstep,
134 : BoundaryCorrectionFilterOnSubstep, VolumeFilterEveryNSteps,
135 : BoundaryCorrectionFilterEveryNSteps>;
136 :
137 0 : static constexpr Options::String help = {
138 : "An exponential filter applied in each direction of a line, square, or "
139 : "cube (hypercube)."};
140 :
141 0 : Hypercube();
142 :
143 0 : Hypercube(unsigned half_power, bool enable,
144 : const std::optional<std::vector<std::string>>& blocks_to_filter,
145 : bool volume_filter_on_substep, bool boundary_filter_on_substep,
146 : std::optional<size_t> volume_filter_every_n_steps,
147 : std::optional<size_t> boundary_filter_every_n_steps,
148 : const Options::Context& context = {});
149 :
150 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
151 : SINGLE_ARG(Filter<Dim, TagList>), Hypercube);
152 0 : explicit Hypercube(CkMigrateMessage* msg) : Filter<Dim, TagList>(msg) {}
153 :
154 : // NOLINTNEXTLINE(google-runtime-references)
155 0 : void pup(PUP::er& p) override;
156 :
157 1 : std::unique_ptr<Filter<Dim, TagList>> get_clone() const override;
158 :
159 1 : bool apply_volume_filter_on_substep() const override;
160 1 : bool apply_volume_filter_on_this_step(size_t step_number) const override;
161 :
162 1 : bool apply_boundary_filter_on_substep() const override;
163 1 : bool apply_boundary_filter_on_this_step(size_t step_number) const override;
164 :
165 1 : bool need_jacobians() const override { return false; }
166 :
167 1 : bool supports_mesh(const Mesh<Dim>& mesh) const override;
168 :
169 1 : const std::optional<std::vector<size_t>>& blocks_to_filter() const override;
170 :
171 1 : void set_blocks_to_filter(
172 : const std::vector<std::string>& all_block_names,
173 : const std::unordered_map<std::string, std::unordered_set<std::string>>&
174 : block_groups) override;
175 :
176 1 : void apply_in_volume(
177 : gsl::not_null<Variables<TagList>*> vars, const Mesh<Dim>& mesh,
178 : const std::optional<
179 : InverseJacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>>&
180 : inv_jac_grid_to_inertial,
181 : const std::optional<
182 : Jacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>>&
183 : jac_grid_to_inertial) const override;
184 :
185 1 : void apply_on_boundary(
186 : gsl::not_null<Variables<TagList>*> vars, const Mesh<Dim - 1>& mesh,
187 : const std::optional<
188 : InverseJacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>>&
189 : inv_jac_grid_to_inertial,
190 : const std::optional<
191 : Jacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>>&
192 : jac_grid_to_inertial) const override;
193 :
194 1 : bool is_equal(const Filter<Dim, TagList>& other) const override;
195 :
196 : private:
197 0 : const Matrix& filter_matrix(const Mesh<1>& mesh) const;
198 :
199 : template <size_t LocalDim, typename LocalTagList>
200 : // NOLINTNEXTLINE(readability-redundant-declaration)
201 0 : friend bool operator==(const Hypercube<LocalDim, LocalTagList>& lhs,
202 : const Hypercube<LocalDim, LocalTagList>& rhs);
203 :
204 0 : unsigned half_power_{0};
205 0 : bool enable_{true};
206 0 : std::optional<std::vector<std::string>> blocks_and_groups_to_filter_{};
207 0 : std::optional<std::vector<size_t>> blocks_to_filter_{};
208 0 : bool volume_filter_on_substep_{false};
209 0 : bool boundary_filter_on_substep_{false};
210 0 : std::optional<size_t> volume_filter_every_n_steps_{std::nullopt};
211 0 : std::optional<size_t> boundary_filter_every_n_steps_{std::nullopt};
212 : };
213 :
214 : template <size_t Dim, typename TagList>
215 0 : bool operator==(const Hypercube<Dim, TagList>& lhs,
216 : const Hypercube<Dim, TagList>& rhs);
217 :
218 : template <size_t Dim, typename TagList>
219 0 : bool operator!=(const Hypercube<Dim, TagList>& lhs,
220 : const Hypercube<Dim, TagList>& rhs);
221 : } // namespace Filters
|
