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 <memory>
8 : #include <optional>
9 : #include <string>
10 : #include <unordered_map>
11 : #include <unordered_set>
12 : #include <vector>
13 :
14 : #include "DataStructures/DataVector.hpp"
15 : #include "DataStructures/Matrix.hpp"
16 : #include "DataStructures/Tensor/TypeAliases.hpp"
17 : #include "NumericalAlgorithms/LinearOperators/Filters/Filter.hpp"
18 : #include "Options/Auto.hpp"
19 : #include "Options/Context.hpp"
20 : #include "Options/String.hpp"
21 : #include "Utilities/Serialization/CharmPupable.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 : /// \cond
25 : template <size_t Dim>
26 : class Mesh;
27 : template <typename TagsList>
28 : class Variables;
29 : namespace PUP {
30 : class er;
31 : } // namespace PUP
32 : namespace gsl {
33 : template <typename T>
34 : class not_null;
35 : } // namespace gsl
36 : /// \endcond
37 :
38 : namespace Filters {
39 : /*!
40 : * \ingroup DiscontinuousGalerkinGroup
41 : * \brief A modal filter for filled-cylinder elements: a smooth exponential
42 : * roll-off of the coupled ZernikeB2 radial-angular (disk) modes and an
43 : * independent exponential roll-off in the axial direction, plus an optional
44 : * top-mode Heaviside cutoff in the angular direction.
45 : *
46 : * Concrete implementation of `Filters::Filter` for the central filled-cylinder
47 : * block created by `domain::creators::AngularCylinder`, driven by the DG
48 : * filtering action. See `Filters::Filter` for the framing of volume vs.
49 : * boundary application, the substep / every-N-steps cadence controls, and the
50 : * `blocks_to_filter` semantics.
51 : *
52 : * A filled-cylinder element uses a ZernikeB2 basis in the radial direction
53 : * (logical dimension 0), a ZernikeB2 basis in the angular direction (logical
54 : * dimension 1), and a Legendre/Chebyshev basis in the axial \f$z\f$ direction
55 : * (logical dimension 2). Unlike the Fourier-angular `Filters::HollowCylinder`
56 : * shells, the ZernikeB2 radial and angular spectral spaces are intertwined (a
57 : * two-dimensional disk basis), so the radial and angular roll-offs cannot be
58 : * specified independently: a single `RadialAngularHalfPower` drives the
59 : * combined disk exponential filter, while `ZHalfPower` drives the axial filter.
60 : *
61 : * In the volume the disk and axial filters are applied by
62 : * `Spectral::filtering::zernike_b2_cylinder_filter`. The coefficient of the
63 : * exponential roll-off is hardcoded to 36, matching `Filters::Hypercube`,
64 : * `Filters::SphericalShell`, and `Filters::HollowCylinder`. When
65 : * `NumModesToKill` is nonzero, the top `NumModesToKill` angular (Fourier
66 : * \f$m\f$) modes are additionally set to zero by a Heaviside cutoff (the
67 : * \f$m=0\f$ mode is always retained).
68 : *
69 : * For boundary-correction filtering the face mesh is classified by basis and
70 : * quadrature:
71 : * - an axial face (both directions ZernikeB2) is a full disk and is filtered by
72 : * `Spectral::filtering::zernike_b2_disk_filter`;
73 : * - a mantle/radial face (ZernikeB2 with `Equiangular` quadrature in the
74 : * angular direction, Legendre/Chebyshev in \f$z\f$) is filtered by treating
75 : * the equiangular angular direction as a Fourier direction (an exponential
76 : * roll-off using `RadialAngularHalfPower` plus the `NumModesToKill` cutoff),
77 : * together with the axial filter;
78 : */
79 : template <typename TagList>
80 1 : class FilledCylinder : public Filter<3, TagList> {
81 : public:
82 : /// \brief The number of top angular (Fourier) \f$m\f$-modes to set to zero.
83 1 : struct NumModesToKill {
84 0 : using type = size_t;
85 0 : static constexpr Options::String help =
86 : "The number of top angular (Fourier) m-modes to set to zero.";
87 : };
88 :
89 : /*!
90 : * \brief Half of the exponent \f$m\f$ in the smooth exponential roll-off
91 : * applied to the coupled ZernikeB2 radial-angular (disk) modal coefficients
92 : * (logical dimensions 0 and 1).
93 : *
94 : * \f{align*}{
95 : * c \to c \exp\left[-36 \left(\frac{n}{N_r}\right)^{2m}\right]
96 : * \exp\left[-36 \left(\frac{k}{K}\right)^{2m}\right]
97 : * \f}
98 : *
99 : * where \f$n\f$ and \f$k\f$ are the radial and angular mode numbers and
100 : * \f$N_r\f$, \f$K\f$ are the corresponding maximum resolved modes. If `None`,
101 : * only the top-mode Heaviside cutoff (if any) is applied to the disk modes.
102 : */
103 1 : struct RadialAngularHalfPower {
104 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
105 0 : static constexpr Options::String help =
106 : "The half-power for the coupled ZernikeB2 radial-angular (disk) "
107 : "exponential roll-off. If None, only the top-mode cutoff is applied to "
108 : "the disk modes.";
109 : };
110 :
111 : /*!
112 : * \brief Half of the exponent \f$m\f$ in the smooth exponential roll-off
113 : * applied to the axial \f$z\f$ modal coefficients (logical dimension 2). If
114 : * `None`, the axial direction is not filtered.
115 : */
116 1 : struct ZHalfPower {
117 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
118 0 : static constexpr Options::String help =
119 : "The half-power for the axial (z) exponential filter. If None, no "
120 : "axial filtering is applied.";
121 : };
122 :
123 : /// \brief Enable (true) or disable (false) the filter
124 1 : struct Enable {
125 0 : using type = bool;
126 0 : static constexpr Options::String help = {"Enable the filter"};
127 : };
128 :
129 : /// \brief Which blocks the filter should be applied to.
130 1 : struct BlocksToFilter {
131 0 : using type =
132 : Options::Auto<std::vector<std::string>, Options::AutoLabel::All>;
133 0 : static constexpr Options::String help = {
134 : "List of blocks or block groups to apply filtering to. All other "
135 : "blocks will have no filtering. You can also specify 'All' to do "
136 : "filtering in all blocks of the domain that are filled cylinders."};
137 : };
138 :
139 : /// \brief Apply the volume filter inside every Runge-Kutta substep
140 : /// instead of only at whole-step boundaries.
141 1 : struct VolumeFilterOnSubstep {
142 0 : using type = bool;
143 0 : static constexpr Options::String help = {
144 : "Enable the volume filter on every substep."};
145 : };
146 :
147 : /// \brief Apply the boundary correction filter inside every Runge-Kutta
148 : /// substep instead of only at whole-step boundaries.
149 1 : struct BoundaryCorrectionFilterOnSubstep {
150 0 : using type = bool;
151 0 : static constexpr Options::String help = {
152 : "Enable the boundary filter on every substep."};
153 : };
154 :
155 : /// \brief Apply the volume filter once every `N` steps. `None`
156 : /// (`std::nullopt`) disables the every-N-steps trigger.
157 : ///
158 : /// \note Currently the check for whether to filter on every `N` steps is done
159 : /// relative to the start of the current Slab. This means that for GTS,
160 : /// independent of the value of `N` for every `N` steps, every step has a
161 : /// filter applied since GTS has one step per slab.
162 1 : struct VolumeFilterEveryNSteps {
163 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
164 0 : static constexpr Options::String help = {
165 : "Enable the volume filter on every N steps. 'None' to disable."};
166 : };
167 :
168 : /// \brief Apply the boundary correction filter once every `N` steps. `None`
169 : /// (`std::nullopt`) disables the every-N-steps trigger.
170 : ///
171 : /// \note Currently the check for whether to filter on every `N` steps is done
172 : /// relative to the start of the current Slab. This means that for GTS,
173 : /// independent of the value of `N` for every `N` steps, every step has a
174 : /// filter applied since GTS has one step per slab.
175 1 : struct BoundaryCorrectionFilterEveryNSteps {
176 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
177 0 : static constexpr Options::String help = {
178 : "Enable the boundary filter on every N steps. 'None' to disable."};
179 : };
180 :
181 0 : using options =
182 : tmpl::list<NumModesToKill, RadialAngularHalfPower, ZHalfPower, Enable,
183 : BlocksToFilter, VolumeFilterOnSubstep,
184 : BoundaryCorrectionFilterOnSubstep, VolumeFilterEveryNSteps,
185 : BoundaryCorrectionFilterEveryNSteps>;
186 :
187 0 : static constexpr Options::String help = {
188 : "A filled-cylinder filter applying a smooth exponential roll-off of the "
189 : "coupled ZernikeB2 radial-angular (disk) modes and an independent "
190 : "exponential roll-off in the axial direction, with an optional top-mode "
191 : "Heaviside cutoff in the angular direction."};
192 :
193 0 : FilledCylinder() = default;
194 :
195 0 : FilledCylinder(
196 : size_t num_modes_to_kill, std::optional<size_t> radial_angular_half_power,
197 : std::optional<size_t> z_half_power, bool enable,
198 : const std::optional<std::vector<std::string>>& blocks_to_filter,
199 : bool volume_filter_on_substep, bool boundary_filter_on_substep,
200 : std::optional<size_t> volume_filter_every_n_steps,
201 : std::optional<size_t> boundary_filter_every_n_steps,
202 : const Options::Context& context = {});
203 :
204 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
205 : SINGLE_ARG(Filter<3, TagList>), FilledCylinder);
206 0 : explicit FilledCylinder(CkMigrateMessage* msg) : Filter<3, TagList>(msg) {}
207 :
208 : // NOLINTNEXTLINE(google-runtime-references)
209 0 : void pup(PUP::er& p) override;
210 :
211 1 : std::unique_ptr<Filter<3, TagList>> get_clone() const override;
212 :
213 1 : bool apply_volume_filter_on_substep() const override;
214 1 : bool apply_volume_filter_on_this_step(size_t step_number) const override;
215 :
216 1 : bool apply_boundary_filter_on_substep() const override;
217 1 : bool apply_boundary_filter_on_this_step(size_t step_number) const override;
218 :
219 1 : bool need_jacobians() const override { return false; }
220 :
221 0 : bool supports_mesh(const Mesh<3>& mesh) const override;
222 :
223 1 : std::string name() const override { return "FilledCylinder"; }
224 :
225 1 : const std::optional<std::vector<size_t>>& blocks_to_filter() const override;
226 :
227 1 : void set_blocks_to_filter(
228 : const std::vector<std::string>& all_block_names,
229 : const std::unordered_map<std::string, std::unordered_set<std::string>>&
230 : block_groups) override;
231 :
232 0 : void apply_in_volume(
233 : gsl::not_null<Variables<TagList>*> vars, const Mesh<3>& mesh,
234 : const std::optional<
235 : InverseJacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
236 : inv_jac_grid_to_inertial,
237 : const std::optional<
238 : Jacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
239 : jac_grid_to_inertial) const override;
240 :
241 0 : void apply_on_boundary(
242 : gsl::not_null<Variables<TagList>*> vars, const Mesh<2>& mesh,
243 : const std::optional<
244 : InverseJacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
245 : inv_jac_grid_to_inertial,
246 : const std::optional<
247 : Jacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
248 : jac_grid_to_inertial) const override;
249 :
250 0 : bool is_equal(const Filter<3, TagList>& other) const override;
251 :
252 : private:
253 : template <typename LocalTagList>
254 : // NOLINTNEXTLINE(readability-redundant-declaration)
255 0 : friend bool operator==(const FilledCylinder<LocalTagList>& lhs,
256 : const FilledCylinder<LocalTagList>& rhs);
257 :
258 : // Caches a single filter matrix together with the parameters it was built
259 : // from: the 1-D extent, the half-power, and (for the angular direction) the
260 : // number of top modes killed. Rebuilt (replacing the previous entry) whenever
261 : // queried with any different parameter, keeping the memory footprint bounded
262 : // instead of growing with the number of distinct extents encountered.
263 0 : struct SingleExtentCache {
264 0 : std::optional<size_t> extent{std::nullopt};
265 0 : std::optional<size_t> half_power{std::nullopt};
266 0 : size_t num_modes_to_kill{0};
267 0 : Matrix matrix{};
268 : };
269 :
270 : // Returns the cached axial exponential filter matrix for a Legendre or
271 : // Chebyshev direction with the given half-power and 1-D mesh, stored in the
272 : // passed-in cache. Returns an empty (identity) matrix when the half-power is
273 : // None.
274 0 : const Matrix& exponential_filter_matrix(std::optional<size_t> half_power,
275 : const Mesh<1>& mesh_1d,
276 : SingleExtentCache& cache) const;
277 :
278 : // Returns the cached angular (Fourier) filter matrix combining the optional
279 : // exponential roll-off and the optional top-mode cutoff. On the mantle face
280 : // the equiangular ZernikeB2-angular direction is treated as Fourier, so the
281 : // caller passes a Fourier/Equiangular `Mesh<1>` with the matching extent.
282 : // Returns an empty (identity) matrix when neither is requested.
283 0 : const Matrix& angular_filter_matrix(const Mesh<1>& mesh_1d) const;
284 :
285 0 : size_t num_modes_to_kill_{0};
286 0 : std::optional<size_t> radial_angular_half_power_{std::nullopt};
287 0 : std::optional<size_t> z_half_power_{std::nullopt};
288 0 : bool enable_{true};
289 0 : std::optional<std::vector<std::string>> blocks_and_groups_to_filter_{};
290 0 : std::optional<std::vector<size_t>> blocks_to_filter_{};
291 0 : bool volume_filter_on_substep_{false};
292 0 : bool boundary_filter_on_substep_{false};
293 0 : std::optional<size_t> volume_filter_every_n_steps_{std::nullopt};
294 0 : std::optional<size_t> boundary_filter_every_n_steps_{std::nullopt};
295 :
296 : // Boundary-filter matrix caches, each bound to the single extent it was built
297 : // for and rebuilt when queried with a different extent. The angular filter
298 : // is used only by `apply_on_boundary` right now: the radial and angular
299 : // filter is applied uncached by
300 : // `Spectral::filtering::zernike_b2_cylinder_filter`, which intertwines the
301 : // radial and angular disk modes and so cannot reuse a per-direction 1-D
302 : // matrix. `cached_z_filter_` holds the axial (Legendre/Chebyshev) filter and
303 : // `cached_angular_filter_` the Fourier-treated mantle-angular filter. The
304 : // empty matrix is returned (and not cached) when a direction is not filtered.
305 : // NOLINTNEXTLINE(spectre-mutable)
306 0 : mutable SingleExtentCache cached_z_filter_{};
307 : // NOLINTNEXTLINE(spectre-mutable)
308 0 : mutable SingleExtentCache cached_angular_filter_{};
309 : // Returned by const reference to represent the identity for unfiltered
310 : // directions; never mutated, so it needs no `mutable`.
311 0 : Matrix empty_matrix_{};
312 : // Buffer used internally by `Spectral::filtering::zernike_b2_cylinder_filter`
313 : // NOLINTNEXTLINE(spectre-mutable)
314 0 : mutable DataVector temp_storage_{};
315 : };
316 :
317 : template <typename TagList>
318 0 : bool operator==(const FilledCylinder<TagList>& lhs,
319 : const FilledCylinder<TagList>& rhs);
320 :
321 : template <typename TagList>
322 0 : bool operator!=(const FilledCylinder<TagList>& lhs,
323 : const FilledCylinder<TagList>& rhs);
324 : } // namespace Filters
|
