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