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