Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <limits>
9 : #include <pup.h>
10 : #include <string>
11 : #include <vector>
12 :
13 : #include "DataStructures/DataBox/DataBox.hpp"
14 : #include "DataStructures/DataBox/DataBoxTag.hpp"
15 : #include "DataStructures/DataBox/ValidateSelection.hpp"
16 : #include "DataStructures/DataVector.hpp"
17 : #include "DataStructures/Tensor/Tensor.hpp"
18 : #include "Domain/Amr/Flag.hpp"
19 : #include "Domain/Tags.hpp"
20 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
21 : #include "Options/Context.hpp"
22 : #include "Options/ParseError.hpp"
23 : #include "Options/String.hpp"
24 : #include "ParallelAlgorithms/Amr/Criteria/Criterion.hpp"
25 : #include "Utilities/TMPL.hpp"
26 :
27 : /// \cond
28 : template <size_t>
29 : class ElementId;
30 : /// \endcond
31 :
32 : namespace amr::Criteria {
33 :
34 : /// @{
35 : /*!
36 : * \brief Computes an anisotropic smoothness indicator based on the power in the
37 : * highest modes
38 : *
39 : * This smoothness indicator is the L2 norm of the tensor component with the
40 : * lowest N - M modes filtered out, where N is the number of grid points in the
41 : * given dimension and M is `num_highest_modes`.
42 : *
43 : * This strategy is similar to the Persson troubled-cell indicator (see
44 : * `evolution::dg::subcell::persson_tci`) with a few modifications:
45 : *
46 : * - The indicator is computed in each dimension separately for an anisotropic
47 : * measure.
48 : * - The number of highest modes to keep can be chosen as a parameter.
49 : * - We don't normalize by the L2 norm of the unfiltered data $u$ here. This
50 : * function just returns the L2 norm of the filtered data.
51 : */
52 : template <size_t Dim>
53 1 : double persson_smoothness_indicator(
54 : gsl::not_null<DataVector*> filtered_component_buffer,
55 : const DataVector& tensor_component, const Mesh<Dim>& mesh, size_t dimension,
56 : size_t num_highest_modes);
57 : template <size_t Dim>
58 1 : std::array<double, Dim> persson_smoothness_indicator(
59 : const DataVector& tensor_component, const Mesh<Dim>& mesh,
60 : size_t num_highest_modes);
61 : /// @}
62 :
63 : namespace Persson_detail {
64 : template <size_t Dim>
65 : void max_over_components(gsl::not_null<std::array<Flag, Dim>*> result,
66 : gsl::not_null<DataVector*> buffer,
67 : const DataVector& tensor_component,
68 : const Mesh<Dim>& mesh, size_t num_highest_modes,
69 : double alpha, double absolute_tolerance,
70 : double coarsening_factor);
71 : }
72 :
73 : /*!
74 : * \brief h-refine the grid based on power in the highest modes
75 : *
76 : * \see persson_smoothness_indicator
77 : */
78 : template <size_t Dim, typename TensorTags>
79 1 : class Persson : public Criterion {
80 : public:
81 0 : struct VariablesToMonitor {
82 0 : using type = std::vector<std::string>;
83 0 : static constexpr Options::String help = {
84 : "The tensors to monitor for h-refinement."};
85 0 : static size_t lower_bound_on_size() { return 1; }
86 : };
87 0 : struct NumHighestModes {
88 0 : using type = size_t;
89 0 : static constexpr Options::String help = {
90 : "Number of highest modes to monitor the power of."};
91 0 : static size_t lower_bound() { return 1; }
92 : };
93 0 : struct Exponent {
94 0 : using type = double;
95 0 : static constexpr Options::String help = {
96 : "The exponent at which the modes should decrease. "
97 : "Corresponds to a \"relative tolerance\" of N^(-alpha), where N is the "
98 : "number of grid points minus 'NumHighestModes'. "
99 : "If any tensor component has power in the highest modes above this "
100 : "value times the max of the absolute tensor component over the "
101 : "element, the element will be h-refined in that direction."};
102 0 : static double lower_bound() { return 0.; }
103 : };
104 0 : struct AbsoluteTolerance {
105 0 : using type = double;
106 0 : static constexpr Options::String help = {
107 : "If any tensor component has a power in the highest modes above this "
108 : "value, the element will be h-refined in that direction. "
109 : "Set to 0 to disable."};
110 0 : static double lower_bound() { return 0.; }
111 : };
112 0 : struct CoarseningFactor {
113 0 : using type = double;
114 0 : static constexpr Options::String help = {
115 : "Factor applied to both relative and absolute tolerance to trigger "
116 : "h-coarsening. Set to 0 to disable h-coarsening altogether. "
117 : "Set closer to 1 to trigger h-coarsening more aggressively. "
118 : "Values too close to 1 risk that coarsened elements will immediately "
119 : "trigger h-refinement again. A reasonable value is 0.1."};
120 0 : static double lower_bound() { return 0.; }
121 0 : static double upper_bound() { return 1.; }
122 : };
123 :
124 0 : using options = tmpl::list<VariablesToMonitor, NumHighestModes, Exponent,
125 : AbsoluteTolerance, CoarseningFactor>;
126 :
127 0 : static constexpr Options::String help = {
128 : "Refine the grid so the power in the highest modes stays below the "
129 : "tolerance"};
130 :
131 0 : Persson() = default;
132 :
133 0 : Persson(std::vector<std::string> vars_to_monitor,
134 : const size_t num_highest_modes, double alpha,
135 : double absolute_tolerance, double coarsening_factor,
136 : const Options::Context& context = {});
137 :
138 : /// \cond
139 : explicit Persson(CkMigrateMessage* msg);
140 : using PUP::able::register_constructor;
141 : WRAPPED_PUPable_decl_template(Persson); // NOLINT
142 : /// \endcond
143 :
144 0 : std::string observation_name() override { return "Persson"; }
145 :
146 0 : using compute_tags_for_observation_box = tmpl::list<>;
147 :
148 0 : using argument_tags = tmpl::list<::Tags::DataBox>;
149 :
150 : template <typename DbTagsList, typename Metavariables>
151 0 : std::array<Flag, Dim> operator()(const db::DataBox<DbTagsList>& box,
152 : Parallel::GlobalCache<Metavariables>& cache,
153 : const ElementId<Dim>& element_id) const;
154 :
155 0 : void pup(PUP::er& p) override;
156 :
157 : private:
158 0 : std::vector<std::string> vars_to_monitor_{};
159 0 : size_t num_highest_modes_{};
160 0 : double alpha_ = std::numeric_limits<double>::signaling_NaN();
161 0 : double absolute_tolerance_ = std::numeric_limits<double>::signaling_NaN();
162 0 : double coarsening_factor_ = std::numeric_limits<double>::signaling_NaN();
163 : };
164 :
165 : // Out-of-line definitions
166 : /// \cond
167 :
168 : template <size_t Dim, typename TensorTags>
169 : Persson<Dim, TensorTags>::Persson(std::vector<std::string> vars_to_monitor,
170 : const size_t num_highest_modes,
171 : const double alpha,
172 : const double absolute_tolerance,
173 : const double coarsening_factor,
174 : const Options::Context& context)
175 : : vars_to_monitor_(std::move(vars_to_monitor)),
176 : num_highest_modes_(num_highest_modes),
177 : alpha_(alpha),
178 : absolute_tolerance_(absolute_tolerance),
179 : coarsening_factor_(coarsening_factor) {
180 : db::validate_selection<TensorTags>(vars_to_monitor_, context);
181 : }
182 :
183 : template <size_t Dim, typename TensorTags>
184 : Persson<Dim, TensorTags>::Persson(CkMigrateMessage* msg) : Criterion(msg) {}
185 :
186 : template <size_t Dim, typename TensorTags>
187 : template <typename DbTagsList, typename Metavariables>
188 : std::array<Flag, Dim> Persson<Dim, TensorTags>::operator()(
189 : const db::DataBox<DbTagsList>& box,
190 : Parallel::GlobalCache<Metavariables>& /*cache*/,
191 : const ElementId<Dim>& /*element_id*/) const {
192 : auto result = make_array<Dim>(Flag::Undefined);
193 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
194 : // Check all tensors and all tensor components in turn. We take the
195 : // highest-priority refinement flag in each dimension, so if any tensor
196 : // component is non-smooth, the element will split in that dimension. And only
197 : // if all tensor components are smooth enough will elements join in that
198 : // dimension.
199 : DataVector buffer(mesh.number_of_grid_points());
200 : tmpl::for_each<TensorTags>(
201 : [&result, &box, &mesh, &buffer, this](const auto tag_v) {
202 : // Stop if we have already decided to refine every dimension
203 : if (result == make_array<Dim>(Flag::Split)) {
204 : return;
205 : }
206 : using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
207 : const std::string tag_name = db::tag_name<tag>();
208 : // Skip if this tensor is not being monitored
209 : if (not alg::found(vars_to_monitor_, tag_name)) {
210 : return;
211 : }
212 : const auto& tensor = db::get<tag>(box);
213 : for (const DataVector& tensor_component : tensor) {
214 : Persson_detail::max_over_components(
215 : make_not_null(&result), make_not_null(&buffer), tensor_component,
216 : mesh, num_highest_modes_, alpha_, absolute_tolerance_,
217 : coarsening_factor_);
218 : }
219 : });
220 : return result;
221 : }
222 :
223 : template <size_t Dim, typename TensorTags>
224 : void Persson<Dim, TensorTags>::pup(PUP::er& p) {
225 : p | vars_to_monitor_;
226 : p | num_highest_modes_;
227 : p | alpha_;
228 : p | absolute_tolerance_;
229 : p | coarsening_factor_;
230 : }
231 :
232 : template <size_t Dim, typename TensorTags>
233 : PUP::able::PUP_ID Persson<Dim, TensorTags>::my_PUP_ID = 0; // NOLINT
234 : /// \endcond
235 :
236 : } // namespace amr::Criteria
|