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/TensorYlm/Filter.hpp"
20 : #include "NumericalAlgorithms/TensorYlm/TensorYlm.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`. `SphericalShell` is the `Filters::Filter`-based implementation
75 : * that plugs into the filtering action and supports per-block selection
76 : * together with independent volume- and boundary-filtering cadences. It is
77 : * intended for spherical-shell blocks, which store Spherepack-normalized
78 : * 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 : ///
162 : /// \note Currently the check for whether to filter on every `N` steps is done
163 : /// relative to the start of the current Slab. This means that for GTS,
164 : /// independent of the value of `N` for every `N` steps, every step has a
165 : /// filter applied since GTS has one step per slab.
166 1 : struct VolumeFilterEveryNSteps {
167 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
168 0 : static constexpr Options::String help = {
169 : "Enable the volume filter on every N steps. 'None' to disable."};
170 : };
171 :
172 : /// \brief Apply the boundary correction filter once every `N` steps. `None`
173 : /// (`std::nullopt`) disables the every-N-steps trigger.
174 : ///
175 : /// \note Currently the check for whether to filter on every `N` steps is done
176 : /// relative to the start of the current Slab. This means that for GTS,
177 : /// independent of the value of `N` for every `N` steps, every step has a
178 : /// filter applied since GTS has one step per slab.
179 1 : struct BoundaryCorrectionFilterEveryNSteps {
180 0 : using type = Options::Auto<size_t, Options::AutoLabel::None>;
181 0 : static constexpr Options::String help = {
182 : "Enable the boundary filter on every N steps. 'None' to disable."};
183 : };
184 :
185 0 : using options =
186 : tmpl::list<NumModesToKill, AngularHalfPower, RadialHalfPower, Enable,
187 : BlocksToFilter, VolumeFilterOnSubstep,
188 : BoundaryCorrectionFilterOnSubstep, VolumeFilterEveryNSteps,
189 : BoundaryCorrectionFilterEveryNSteps>;
190 :
191 0 : static constexpr Options::String help = {
192 : "A spherical-shell filter applying a top-ell Heaviside cutoff in the "
193 : "angular direction with optional smooth exponential roll-offs in both "
194 : "the angular ell direction and the radial direction."};
195 :
196 0 : SphericalShell() = default;
197 :
198 0 : SphericalShell(
199 : size_t num_modes_to_kill, std::optional<size_t> angular_half_power,
200 : std::optional<size_t> radial_half_power, bool enable,
201 : const std::optional<std::vector<std::string>>& blocks_to_filter,
202 : bool volume_filter_on_substep, bool boundary_filter_on_substep,
203 : std::optional<size_t> volume_filter_every_n_steps,
204 : std::optional<size_t> boundary_filter_every_n_steps,
205 : const Options::Context& context = {});
206 :
207 0 : WRAPPED_PUPable_decl_base_template( // NOLINT
208 : SINGLE_ARG(Filter<3, TagList>), SphericalShell);
209 0 : explicit SphericalShell(CkMigrateMessage* msg) : Filter<3, TagList>(msg) {}
210 :
211 : // NOLINTNEXTLINE(google-runtime-references)
212 0 : void pup(PUP::er& p) override;
213 :
214 1 : std::unique_ptr<Filter<3, TagList>> get_clone() const override;
215 :
216 1 : bool apply_volume_filter_on_substep() const override;
217 1 : bool apply_volume_filter_on_this_step(size_t step_number) const override;
218 :
219 1 : bool apply_boundary_filter_on_substep() const override;
220 1 : bool apply_boundary_filter_on_this_step(size_t step_number) const override;
221 :
222 1 : bool need_jacobians() const override { return true; }
223 :
224 0 : bool supports_mesh(const Mesh<3>& mesh) const override;
225 :
226 1 : std::string name() const override { return "SphericalShell"; }
227 :
228 1 : const std::optional<std::vector<size_t>>& blocks_to_filter() const override;
229 :
230 1 : void set_blocks_to_filter(
231 : const std::vector<std::string>& all_block_names,
232 : const std::unordered_map<std::string, std::unordered_set<std::string>>&
233 : block_groups) override;
234 :
235 0 : void apply_in_volume(
236 : gsl::not_null<Variables<TagList>*> vars, const Mesh<3>& mesh,
237 : const std::optional<
238 : InverseJacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
239 : inv_jac_grid_to_inertial,
240 : const std::optional<
241 : Jacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
242 : jac_grid_to_inertial) const override;
243 :
244 0 : void apply_on_boundary(
245 : gsl::not_null<Variables<TagList>*> vars, const Mesh<2>& mesh,
246 : const std::optional<
247 : InverseJacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
248 : inv_jac_grid_to_inertial,
249 : const std::optional<
250 : Jacobian<DataVector, 3, Frame::Grid, Frame::Inertial>>&
251 : jac_grid_to_inertial) const override;
252 :
253 0 : bool is_equal(const Filter<3, TagList>& other) const override;
254 :
255 : private:
256 : template <typename LocalTagList>
257 : // NOLINTNEXTLINE(readability-redundant-declaration)
258 0 : friend bool operator==(const SphericalShell<LocalTagList>& lhs,
259 : const SphericalShell<LocalTagList>& rhs);
260 :
261 0 : size_t num_modes_to_kill_{0};
262 0 : std::optional<size_t> angular_half_power_{std::nullopt};
263 0 : std::optional<size_t> radial_half_power_{std::nullopt};
264 0 : bool enable_{true};
265 0 : std::optional<std::vector<std::string>> blocks_and_groups_to_filter_{};
266 0 : std::optional<std::vector<size_t>> blocks_to_filter_{};
267 0 : bool volume_filter_on_substep_{false};
268 0 : bool boundary_filter_on_substep_{false};
269 0 : std::optional<size_t> volume_filter_every_n_steps_{std::nullopt};
270 0 : std::optional<size_t> boundary_filter_every_n_steps_{std::nullopt};
271 :
272 : // Use Spherepack normalization because the variables are stored as Spherepack
273 : // modes
274 0 : static constexpr ylm::TensorYlm::CoefficientNormalization normalization_ =
275 : ylm::TensorYlm::CoefficientNormalization::Spherepack;
276 : // Caches and memory buffers
277 : // NOLINTNEXTLINE(spectre-mutable)
278 0 : mutable size_t cached_l_max_{0};
279 : // NOLINTNEXTLINE(spectre-mutable)
280 0 : mutable ylm::TensorYlm::FilterMatrixHolder filter_matrices_{};
281 : // NOLINTNEXTLINE(spectre-mutable)
282 0 : mutable size_t cached_radial_extents_{0};
283 : // NOLINTNEXTLINE(spectre-mutable)
284 0 : mutable Matrix cached_radial_filter_matrix_{};
285 : // NOLINTNEXTLINE(spectre-mutable)
286 0 : mutable Variables<TagList> temp_storage_{};
287 : };
288 :
289 : template <typename TagList>
290 0 : bool operator==(const SphericalShell<TagList>& lhs,
291 : const SphericalShell<TagList>& rhs);
292 :
293 : template <typename TagList>
294 0 : bool operator!=(const SphericalShell<TagList>& lhs,
295 : const SphericalShell<TagList>& rhs);
296 : } // namespace Filters
|
