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/Matrix.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "DataStructures/Tensor/TypeAliases.hpp"
17 : #include "DataStructures/Variables.hpp"
18 : #include "NumericalAlgorithms/LinearOperators/Filters/Filter.hpp"
19 : #include "NumericalAlgorithms/SphericalHarmonics/TensorYlm.hpp"
20 : #include "NumericalAlgorithms/SphericalHarmonics/TensorYlmFilter.hpp"
21 : #include "Options/Auto.hpp"
22 : #include "Options/Context.hpp"
23 : #include "Options/String.hpp"
24 : #include "Utilities/Gsl.hpp"
25 :
26 : /// \cond
27 : template <size_t Dim>
28 : class Mesh;
29 : /// \endcond
30 :
31 : namespace Filters {
32 : /*!
33 : * \ingroup DiscontinuousGalerkinGroup
34 : * \brief A modal filter for spherical-shell elements: a top-$\ell$ Heaviside
35 : * cutoff in the angular direction plus optional smooth exponential roll-offs
36 : * in both the angular $\ell$ direction and the radial direction.
37 : *
38 : * Concrete implementation of `Filters::Filter` for spherical-shell elements,
39 : * driven by the DG filtering action. See `Filters::Filter` for the framing of
40 : * volume vs. boundary application, the substep / every-N-steps cadence
41 : * controls, and the `blocks_to_filter` semantics.
42 : *
43 : * For each component of the tensors in `TagList`, the filter rescales the
44 : * Spherepack-normalized angular modal coefficients $c_{\ell'}$ as
45 : *
46 : * \f{align*}{
47 : * c_{\ell'} \to c_{\ell'} \exp\!\left[-36 \left(\frac{\ell'}
48 : * {\ell^+_{\mathrm{cut}}+1}\right)^{2\sigma_a}\right],
49 : * \f}
50 : *
51 : * where $\ell^+_{\mathrm{cut}} = \ell_{\mathrm{max}} -$ `NumModesToKill` is
52 : * the largest angular mode that is retained and $\sigma_a$ is the
53 : * `AngularHalfPower` option. With the fixed coefficient 36 and $\sigma_a$ in
54 : * the typical range 28-32 the angular filter is smooth below
55 : * $\ell^+_{\mathrm{cut}}$ and reduces to a sharp Heaviside cutoff as
56 : * $\sigma_a \to \infty$. When `AngularHalfPower` is `None`, only the
57 : * Heaviside cutoff is applied. See `ylm::TensorYlm` for the derivation of
58 : * the underlying angular filter.
59 : *
60 : * When `RadialHalfPower` has a value $\sigma_r$, an additional 1D
61 : * exponential filter is applied independently along each radial column. The
62 : * radial nodal coefficients $c_i$ are rescaled in modal space as
63 : *
64 : * \f{align*}{
65 : * c_i \to c_i \exp\!\left[-36 \left(\frac{i}{N_r}\right)^{2\sigma_r}\right],
66 : * \f}
67 : *
68 : * where $N_r$ is the radial basis degree (radial extent minus one). When
69 : * `RadialHalfPower` is `None`, the radial direction is left untouched.
70 : *
71 : * #### Design decision:
72 : *
73 : * - The exponential coefficient is hardcoded to 36, matching the choice in
74 : * `Hypercube` in `Cube.hpp` and `Filters::Exponential`. `SphericalShell` is
75 : * the `Filters::Filter`-based implementation that plugs into the filtering
76 : * action and supports per-block selection together with independent volume-
77 : * and boundary-filtering cadences. It is intended for spherical-shell
78 : * blocks, which store Spherepack-normalized spherical-harmonic modes.
79 : */
80 : template <typename TagList>
81 1 : class SphericalShell : public Filter<3, TagList> {
82 : public:
83 : /// \brief The number of top $\ell$ modes to set to zero.
84 1 : struct NumModesToKill {
85 0 : using type = size_t;
86 0 : static constexpr Options::String help =
87 : "The number of top ell modes to set to zero.";
88 : };
89 :
90 : /*!
91 : * \brief Half of the exponent $\sigma_a$ in the smooth exponential roll-off
92 : * applied to the angular $\ell$ modes below the top-$\ell$ cutoff.
93 : *
94 : * \f{align*}{
95 : * c_{\ell'} \to c_{\ell'} \exp\left[-36 \left(\frac{\ell'}
96 : * {\ell^+_{\mathrm{cut}}+1}\right)^{2\sigma_a}\right]
97 : * \f}
98 : *
99 : * If `None`, only the Heaviside top-$\ell$ cutoff is applied to the
100 : * angular modes.
101 : */
102 1 : struct AngularHalfPower {
103 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
104 0 : static constexpr Options::String help =
105 : "The half-power sigma for the angular ell-mode exponential roll-off. "
106 : "If None, only the top-ell Heaviside cutoff is applied.";
107 : };
108 :
109 : /*!
110 : * \brief Half of the exponent $\sigma_r$ in the smooth exponential
111 : * roll-off applied to the radial modal coefficients.
112 : *
113 : * \f{align*}{
114 : * c_i \to c_i \exp\left[-36 \left(\frac{i}{N_r}\right)^{2\sigma_r}\right]
115 : * \f}
116 : *
117 : * where $N_r$ is the radial basis degree. If `None`, the radial direction
118 : * is not filtered.
119 : */
120 1 : struct RadialHalfPower {
121 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
122 0 : static constexpr Options::String help =
123 : "The half-power sigma for the radial exponential filter. "
124 : "If None, no radial filtering is applied.";
125 : };
126 :
127 : /// \brief Enable (true) or disable (false) the filter
128 1 : struct Enable {
129 0 : using type = bool;
130 0 : static constexpr Options::String help = {"Enable the filter"};
131 : };
132 :
133 : /// \brief Which blocks the filter should be applied to.
134 1 : struct BlocksToFilter {
135 0 : using type =
136 : Options::Auto<std::vector<std::string>, Options::AutoLabel::All>;
137 0 : static constexpr Options::String help = {
138 : "List of blocks or block groups to apply filtering to. All other "
139 : "blocks will have no filtering. You can also specify 'All' to do "
140 : "filtering in all blocks of the domain that are spherical shells."};
141 : };
142 :
143 : /// \brief Apply the volume filter inside every Runge-Kutta substep
144 : /// instead of only at whole-step boundaries.
145 1 : struct VolumeFilterOnSubstep {
146 0 : using type = bool;
147 0 : static constexpr Options::String help = {
148 : "Enable the volume filter on every substep."};
149 : };
150 :
151 : /// \brief Apply the boundary correction filter inside every Runge-Kutta
152 : /// substep instead of only at whole-step boundaries.
153 1 : struct BoundaryCorrectionFilterOnSubstep {
154 0 : using type = bool;
155 0 : static constexpr Options::String help = {
156 : "Enable the boundary filter on every substep."};
157 : };
158 :
159 : /// \brief Apply the volume filter once every `N` steps. `None`
160 : /// (`std::nullopt`) disables the every-N-steps trigger.
161 1 : struct VolumeFilterEveryNSteps {
162 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
163 0 : static constexpr Options::String help = {
164 : "Enable the volume filter on every N steps. 'None' to disable."};
165 : };
166 :
167 : /// \brief Apply the boundary correction filter once every `N` steps. `None`
168 : /// (`std::nullopt`) disables the every-N-steps trigger.
169 1 : struct BoundaryCorrectionFilterEveryNSteps {
170 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
171 0 : static constexpr Options::String help = {
172 : "Enable the boundary filter on every N steps. 'None' to disable."};
173 : };
174 :
175 0 : using options =
176 : tmpl::list<NumModesToKill, AngularHalfPower, RadialHalfPower, Enable,
177 : BlocksToFilter, VolumeFilterOnSubstep,
178 : BoundaryCorrectionFilterOnSubstep, VolumeFilterEveryNSteps,
179 : BoundaryCorrectionFilterEveryNSteps>;
180 :
181 0 : static constexpr Options::String help = {
182 : "A spherical-shell filter applying a top-ell Heaviside cutoff in the "
183 : "angular direction with optional smooth exponential roll-offs in both "
184 : "the angular ell direction and the radial direction."};
185 :
186 0 : SphericalShell() = default;
187 :
188 0 : SphericalShell(
189 : size_t num_modes_to_kill, std::optional<size_t> angular_half_power,
190 : std::optional<size_t> radial_half_power, bool enable,
191 : const std::optional<std::vector<std::string>>& blocks_to_filter,
192 : bool volume_filter_on_substep, bool boundary_filter_on_substep,
193 : std::optional<size_t> volume_filter_every_n_steps,
194 : std::optional<size_t> boundary_filter_every_n_steps,
195 : const Options::Context& context = {});
196 :
197 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
198 : SINGLE_ARG(Filter<3, TagList>), SphericalShell);
199 0 : explicit SphericalShell(CkMigrateMessage* msg) : Filter<3, TagList>(msg) {}
200 :
201 : // NOLINTNEXTLINE(google-runtime-references)
202 0 : void pup(PUP::er& p) override;
203 :
204 1 : std::unique_ptr<Filter<3, TagList>> get_clone() const override;
205 :
206 1 : bool apply_volume_filter_on_substep() const override;
207 1 : bool apply_volume_filter_on_this_step(size_t step_number) const override;
208 :
209 1 : bool apply_boundary_filter_on_substep() const override;
210 1 : bool apply_boundary_filter_on_this_step(size_t step_number) const override;
211 :
212 1 : bool need_jacobians() const override { return true; }
213 :
214 0 : bool supports_mesh(const Mesh<3>& mesh) const override;
215 :
216 1 : const std::optional<std::vector<size_t>>& blocks_to_filter() const override;
217 :
218 1 : void set_blocks_to_filter(
219 : const std::vector<std::string>& all_block_names,
220 : const std::unordered_map<std::string, std::unordered_set<std::string>>&
221 : block_groups) override;
222 :
223 0 : void apply_in_volume(
224 : gsl::not_null<Variables<TagList>*> vars, const Mesh<3>& mesh,
225 : const std::optional<
226 : InverseJacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
227 : inv_jac_grid_to_inertial,
228 : const std::optional<
229 : Jacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
230 : jac_grid_to_inertial) const override;
231 :
232 0 : void apply_on_boundary(
233 : gsl::not_null<Variables<TagList>*> vars, const Mesh<2>& 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 : bool is_equal(const Filter<3, TagList>& other) const override;
242 :
243 : private:
244 : template <typename LocalTagList>
245 : // NOLINTNEXTLINE(readability-redundant-declaration)
246 0 : friend bool operator==(const SphericalShell<LocalTagList>& lhs,
247 : const SphericalShell<LocalTagList>& rhs);
248 :
249 0 : size_t num_modes_to_kill_{0};
250 0 : std::optional<size_t> angular_half_power_{std::nullopt};
251 0 : std::optional<size_t> radial_half_power_{std::nullopt};
252 0 : bool enable_{true};
253 0 : std::optional<std::vector<std::string>> blocks_and_groups_to_filter_{};
254 0 : std::optional<std::vector<size_t>> blocks_to_filter_{};
255 0 : bool volume_filter_on_substep_{false};
256 0 : bool boundary_filter_on_substep_{false};
257 0 : std::optional<size_t> volume_filter_every_n_steps_{std::nullopt};
258 0 : std::optional<size_t> boundary_filter_every_n_steps_{std::nullopt};
259 :
260 : // Use Spherepack normalization because the variables are stored as Spherepack
261 : // modes
262 0 : static constexpr ylm::TensorYlm::CoefficientNormalization normalization_ =
263 : ylm::TensorYlm::CoefficientNormalization::Spherepack;
264 : // Caches and memory buffers
265 : // NOLINTNEXTLINE(spectre-mutable)
266 0 : mutable size_t cached_l_max_{0};
267 : // NOLINTNEXTLINE(spectre-mutable)
268 0 : mutable ylm::TensorYlm::FilterMatrixHolder filter_matrices_{};
269 : // NOLINTNEXTLINE(spectre-mutable)
270 0 : mutable size_t cached_radial_extents_{0};
271 : // NOLINTNEXTLINE(spectre-mutable)
272 0 : mutable Matrix cached_radial_filter_matrix_{};
273 : // NOLINTNEXTLINE(spectre-mutable)
274 0 : mutable Variables<TagList> temp_storage_{};
275 : };
276 :
277 : template <typename TagList>
278 0 : bool operator==(const SphericalShell<TagList>& lhs,
279 : const SphericalShell<TagList>& rhs);
280 :
281 : template <typename TagList>
282 0 : bool operator!=(const SphericalShell<TagList>& lhs,
283 : const SphericalShell<TagList>& rhs);
284 : } // namespace Filters
|
