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