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 <utility>
9 :
10 : #include "DataStructures/DataBox/DataBox.hpp"
11 : #include "DataStructures/DataVector.hpp"
12 : #include "DataStructures/Tensor/Tensor.hpp"
13 : #include "DataStructures/Variables.hpp"
14 : #include "DataStructures/VariablesTag.hpp"
15 : #include "Domain/Tags.hpp"
16 : #include "NumericalAlgorithms/LinearOperators/Filters/Filter.hpp"
17 : #include "NumericalAlgorithms/LinearOperators/Filters/None.hpp"
18 : #include "NumericalAlgorithms/LinearOperators/Filters/Tag.hpp"
19 : #include "Parallel/AlgorithmExecution.hpp"
20 : #include "Parallel/GlobalCache.hpp"
21 : #include "Time/Tags/StepNumberWithinSlab.hpp"
22 : #include "Utilities/Gsl.hpp"
23 :
24 : /// \cond
25 : template <size_t Dim>
26 : class Mesh;
27 : namespace tuples {
28 : template <typename... Tags>
29 : class TaggedTuple;
30 : } // namespace tuples
31 : /// \endcond
32 :
33 0 : namespace dg::Actions {
34 : /*!
35 : * \ingroup DiscontinuousGalerkinGroup
36 : * \brief Applies the element's spectral volume filter to the evolved variables.
37 : *
38 : * Retrieves the per-element filter from `Filters::Tags::SpectralFilter` from
39 : * the DataBox and applies the volume filter.
40 : *
41 : * The filter is skipped when:
42 : * - the filter is `Filters::None` (explicit no-op), or
43 : * - both `apply_volume_filter_on_substep()` and
44 : * `apply_volume_filter_on_this_step(step_number)` return `false`.
45 : *
46 : * The grid-to-inertial Jacobian and its inverse are only retrieved from the
47 : * DataBox when `Filters::Filter::need_jacobians()` returns `true`; otherwise
48 : * `std::nullopt` is passed for both arguments.
49 : *
50 : * \note Currently the check for whether to filter on every `N` steps is done
51 : * relative to the start of the current Slab. This means that for GTS,
52 : * independent of the value of `N` for every `N` steps, every step has a
53 : * filter applied since GTS has one step per slab.
54 : *
55 : * Uses:
56 : * - DataBox:
57 : * - `Filters::Tags::SpectralFilter<Dim, TagList>`
58 : * - `Tags::StepNumberWithinSlab`
59 : * - `domain::Tags::Mesh<Dim>`
60 : * - `domain::Tags::InverseJacobian<Dim, Frame::Grid, Frame::Inertial>`
61 : * (only when `need_jacobians()` is `true`)
62 : * - `domain::Tags::Jacobian<Dim, Frame::Grid, Frame::Inertial>`
63 : * (only when `need_jacobians()` is `true`)
64 : * - DataBox changes:
65 : * - Modifies: `Metavariables::system::variables_tag`
66 : * - System:
67 : * - `volume_dim`
68 : * - `variables_tag`
69 : */
70 : template <size_t Dim, typename TagList>
71 1 : struct SpectralFilter {
72 : template <typename DbTags, typename... InboxTags, typename ArrayIndex,
73 : typename ActionList, typename ParallelComponent,
74 : typename Metavariables>
75 0 : static Parallel::iterable_action_return_t apply(
76 : db::DataBox<DbTags>& box,
77 : const tuples::TaggedTuple<InboxTags...>& /*inboxes*/,
78 : const Parallel::GlobalCache<Metavariables>& /*cache*/,
79 : const ArrayIndex& /*array_index*/, const ActionList /*meta*/,
80 : const ParallelComponent* const /*meta*/) {
81 : const auto& filter =
82 : db::get<Filters::Tags::SpectralFilter<Dim, TagList>>(box);
83 :
84 : const Mesh<Dim>& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
85 : if (not filter.supports_mesh(mesh)) {
86 : ERROR("The filter '"
87 : << filter.name() << "' on element "
88 : << db::get<domain::Tags::Element<Dim>>(box).id() << " with mesh "
89 : << mesh
90 : << " does not support that mesh. This is a bug in initialization "
91 : "or if AMR projection during h-refinement incorrectly adjusted "
92 : "the filter.");
93 : }
94 :
95 : // None is explicitly a no-op; skip immediately.
96 : if (dynamic_cast<const Filters::None<Dim, TagList>*>(&filter) != nullptr) {
97 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
98 : }
99 :
100 : const auto step_number =
101 : static_cast<size_t>(db::get<::Tags::StepNumberWithinSlab>(box));
102 : if (not(filter.apply_volume_filter_on_substep() or
103 : filter.apply_volume_filter_on_this_step(step_number))) {
104 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
105 : }
106 :
107 : const size_t num_pts = mesh.number_of_grid_points();
108 :
109 : std::optional<
110 : InverseJacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>>
111 : inv_jac{std::nullopt};
112 : std::optional<Jacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>> jac{
113 : std::nullopt};
114 : if (filter.need_jacobians()) {
115 : if constexpr (db::tag_is_retrievable_v<
116 : domain::Tags::InverseJacobian<Dim, Frame::Grid,
117 : Frame::Inertial>,
118 : db::DataBox<DbTags>>) {
119 : const auto& src = db::get<
120 : domain::Tags::InverseJacobian<Dim, Frame::Grid, Frame::Inertial>>(
121 : box);
122 : inv_jac =
123 : InverseJacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>{};
124 : for (size_t i = 0; i < src.size(); ++i) {
125 : make_const_view(make_not_null(&std::as_const((*inv_jac)[i])), src[i],
126 : 0, num_pts);
127 : }
128 : } else {
129 : ERROR(
130 : "SpectralFilter: Missing "
131 : "domain::Tags::InverseJacobian<Dim, Grid, Inertial> in the "
132 : "DataBox. This should be present for both time-dependent and "
133 : "time-independent evolutions.");
134 : }
135 : if constexpr (db::tag_is_retrievable_v<
136 : domain::Tags::Jacobian<Dim, Frame::Grid,
137 : Frame::Inertial>,
138 : db::DataBox<DbTags>>) {
139 : const auto& src =
140 : db::get<domain::Tags::Jacobian<Dim, Frame::Grid, Frame::Inertial>>(
141 : box);
142 : jac = Jacobian<DataVector, Dim, Frame::Grid, Frame::Inertial>{};
143 : for (size_t i = 0; i < src.size(); ++i) {
144 : make_const_view(make_not_null(&std::as_const((*jac)[i])), src[i], 0,
145 : num_pts);
146 : }
147 : } else {
148 : ERROR(
149 : "SpectralFilter: Missing "
150 : "domain::Tags::Jacobian<Dim, Grid, Inertial> in the "
151 : "DataBox. This should be present for both time-dependent and "
152 : "time-independent evolutions.");
153 : }
154 : }
155 :
156 : db::mutate<::Tags::Variables<TagList>>(
157 : [&filter, &mesh, &inv_jac,
158 : &jac](const gsl::not_null<Variables<TagList>*> vars) {
159 : filter.apply_in_volume(vars, mesh, inv_jac, jac);
160 : },
161 : make_not_null(&box));
162 :
163 : return {Parallel::AlgorithmExecution::Continue, std::nullopt};
164 : }
165 : };
166 : } // namespace dg::Actions
|
