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