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 magnitude of
38 : * second derivatives
39 : *
40 : * This smoothness indicator is simply the L2 norm of the logical second
41 : * derivative of the tensor component in the given `dimension`:
42 : *
43 : * \begin{equation}
44 : * \epsilon_k =
45 : * \sqrt{\frac{1}{N_{\mathrm{points}}} \sum_{p=1}^{N_\mathrm{points}}
46 : * \left(\partial^2 u / \partial \xi_k^2\right)^2}
47 : * \end{equation}
48 : *
49 : * If the smoothness indicator is large in a direction, meaning the tensor
50 : * component has a large second derivative in that direction, the element should
51 : * be h-refined. If the smoothness indicator is small, the element should be
52 : * h-coarsened. A coarsing threshold of about a third of the refinement
53 : * threshold seems to work well, but this will need more testing.
54 : *
55 : * Note that it is not at all clear that a smoothness indicator based on the
56 : * magnitude of second derivatives is useful in a DG context. Smooth functions
57 : * with higher-order derivatives can be approximated just fine by higher-order
58 : * DG elements without the need for h-refinement. The use of second derivatives
59 : * to indicate the need for refinement originated in the finite element context
60 : * with linear elements. Other smoothness indicators might prove more useful for
61 : * DG elements, e.g. based on jumps or oscillations of the solution. We can also
62 : * explore applying the troubled-cell indicators (TCIs) used in hydrodynamic
63 : * simulations as h-refinement indicators.
64 : *
65 : * Specifically, this smoothness indicator is based on \cite Loehner1987 (hence
66 : * the name of the function), which is popular in the finite element community
67 : * and also used in a DG context by \cite Dumbser2013, Eq. (34), and by
68 : * \cite Renkhoff2023, Eq. (15). We make several modifications:
69 : *
70 : * - The original smoothness indicator is isotropic, i.e. it computes the norm
71 : * over all (mixed) second derivatives. Here we compute an anisotropic
72 : * indicator by computing second derivatives in each dimension separately
73 : * and ignoring mixed derivatives.
74 : * - The original smoothness indicator is normalized by measures of the first
75 : * derivative which don't generalize well to spectral elements. Therefore, we
76 : * simplify the normalization to a standard relative and absolute tolerance.
77 : * An alternative approach is proposed in \cite Renkhoff2023, Eq.(15), where
78 : * the authors take the absolute value of the differentiation matrix and apply
79 : * the resulting matrix to the absolute value of the data on the grid to
80 : * compute the normalization. However, this quantity can produce quite large
81 : * numbers and hence overestimates the smoothness by suppressing the second
82 : * derivative.
83 : * - We compute the second derivative in logical coordinates. This seems
84 : * easiest for spectral elements, but note that \cite Renkhoff2023 seem to
85 : * use inertial coordinates.
86 : *
87 : * In addition to the above modifications, we can consider approximating the
88 : * second derivative using finite differences, as explored in the prototype
89 : * https://github.com/sxs-collaboration/dg-charm/blob/HpAmr/Evolution/HpAmr/LohnerRefiner.hpp.
90 : * This would allow falling back to the normalization used by Löhner and might
91 : * be cheaper to compute, but it requires an interpolation to the center and
92 : * maybe also to the faces, depending on the desired stencil.
93 : */
94 : template <typename VectorType, size_t Dim>
95 1 : double loehner_smoothness_indicator(
96 : gsl::not_null<VectorType*> first_deriv_buffer,
97 : gsl::not_null<VectorType*> second_deriv_buffer,
98 : const VectorType& tensor_component, const Mesh<Dim>& mesh,
99 : size_t dimension);
100 : template <typename VectorType, size_t Dim>
101 1 : std::array<double, Dim> loehner_smoothness_indicator(
102 : const VectorType& tensor_component, const Mesh<Dim>& mesh);
103 : /// @}
104 :
105 : namespace Loehner_detail {
106 : template <typename VectorType, size_t Dim>
107 : void max_over_components(
108 : gsl::not_null<std::array<Flag, Dim>*> result,
109 : gsl::not_null<std::array<VectorType, 2>*> deriv_buffers,
110 : const VectorType& tensor_component, const Mesh<Dim>& mesh,
111 : double relative_tolerance, double absolute_tolerance,
112 : double coarsening_factor);
113 : }
114 :
115 : /*!
116 : * \brief h-refine the grid based on a smoothness indicator
117 : *
118 : * The smoothness indicator used here is based on the magnitude of second
119 : * derivatives. See `amr::Criteria::loehner_smoothness_indicator` for details
120 : * and caveats.
121 : *
122 : * \see amr::Criteria::loehner_smoothness_indicator
123 : */
124 : template <size_t Dim, typename TensorTags>
125 1 : class Loehner : public Criterion {
126 : public:
127 0 : struct VariablesToMonitor {
128 0 : using type = std::vector<std::string>;
129 0 : static constexpr Options::String help = {
130 : "The tensors to monitor for h-refinement."};
131 0 : static size_t lower_bound_on_size() { return 1; }
132 : };
133 0 : struct RelativeTolerance {
134 0 : using type = double;
135 0 : static constexpr Options::String help = {
136 : "If any tensor component has a second derivative magnitude above this "
137 : "value times the max of the absolute tensor component over the "
138 : "element, the element will be h-refined in that direction. "
139 : "Set to 0 to disable."};
140 0 : static double lower_bound() { return 0.; }
141 : };
142 0 : struct AbsoluteTolerance {
143 0 : using type = double;
144 0 : static constexpr Options::String help = {
145 : "If any tensor component has a second derivative magnitude above this "
146 : "value, the element will be h-refined in that direction. "
147 : "Set to 0 to disable."};
148 0 : static double lower_bound() { return 0.; }
149 : };
150 0 : struct CoarseningFactor {
151 0 : using type = double;
152 0 : static constexpr Options::String help = {
153 : "Factor applied to both relative and absolute tolerance to trigger "
154 : "h-coarsening. Set to 0 to disable h-coarsening altogether. "
155 : "Set closer to 1 to trigger h-coarsening more aggressively. "
156 : "Values too close to 1 risk that coarsened elements will immediately "
157 : "trigger h-refinement again. A reasonable value is 1/3."};
158 0 : static double lower_bound() { return 0.; }
159 0 : static double upper_bound() { return 1.; }
160 : };
161 :
162 0 : using options = tmpl::list<VariablesToMonitor, RelativeTolerance,
163 : AbsoluteTolerance, CoarseningFactor>;
164 :
165 0 : static constexpr Options::String help = {
166 : "Refine the grid towards resolving an estimated error in the second "
167 : "derivative"};
168 :
169 0 : Loehner() = default;
170 :
171 0 : Loehner(std::vector<std::string> vars_to_monitor, double relative_tolerance,
172 : double absolute_tolerance, double coarsening_factor,
173 : const Options::Context& context = {});
174 :
175 : /// \cond
176 : explicit Loehner(CkMigrateMessage* msg);
177 : using PUP::able::register_constructor;
178 : WRAPPED_PUPable_decl_template(Loehner); // NOLINT
179 : /// \endcond
180 :
181 0 : Type type() override { return Type::h; }
182 :
183 0 : std::string observation_name() override { return "Loehner"; }
184 :
185 0 : using compute_tags_for_observation_box = tmpl::list<>;
186 :
187 0 : using argument_tags = tmpl::list<::Tags::DataBox>;
188 :
189 : template <typename DbTagsList, typename Metavariables>
190 0 : std::array<Flag, Dim> operator()(const db::DataBox<DbTagsList>& box,
191 : Parallel::GlobalCache<Metavariables>& cache,
192 : const ElementId<Dim>& element_id) const;
193 :
194 0 : void pup(PUP::er& p) override;
195 :
196 : private:
197 0 : std::vector<std::string> vars_to_monitor_{};
198 0 : double relative_tolerance_ = std::numeric_limits<double>::signaling_NaN();
199 0 : double absolute_tolerance_ = std::numeric_limits<double>::signaling_NaN();
200 0 : double coarsening_factor_ = std::numeric_limits<double>::signaling_NaN();
201 : };
202 :
203 : // Out-of-line definitions
204 : /// \cond
205 :
206 : template <size_t Dim, typename TensorTags>
207 : Loehner<Dim, TensorTags>::Loehner(std::vector<std::string> vars_to_monitor,
208 : const double relative_tolerance,
209 : const double absolute_tolerance,
210 : const double coarsening_factor,
211 : const Options::Context& context)
212 : : vars_to_monitor_(std::move(vars_to_monitor)),
213 : relative_tolerance_(relative_tolerance),
214 : absolute_tolerance_(absolute_tolerance),
215 : coarsening_factor_(coarsening_factor) {
216 : db::validate_selection<TensorTags>(vars_to_monitor_, context);
217 : if (relative_tolerance == 0. and absolute_tolerance == 0.) {
218 : PARSE_ERROR(
219 : context,
220 : "Must specify non-zero RelativeTolerance, AbsoluteTolerance, or both.");
221 : }
222 : }
223 :
224 : template <size_t Dim, typename TensorTags>
225 : Loehner<Dim, TensorTags>::Loehner(CkMigrateMessage* msg) : Criterion(msg) {}
226 :
227 : template <size_t Dim, typename TensorTags>
228 : template <typename DbTagsList, typename Metavariables>
229 : std::array<Flag, Dim> Loehner<Dim, TensorTags>::operator()(
230 : const db::DataBox<DbTagsList>& box,
231 : Parallel::GlobalCache<Metavariables>& /*cache*/,
232 : const ElementId<Dim>& /*element_id*/) const {
233 : auto result = make_array<Dim>(Flag::Undefined);
234 : const auto& mesh = db::get<domain::Tags::Mesh<Dim>>(box);
235 : using VectorType = typename Variables<TensorTags>::vector_type;
236 : // Check all tensors and all tensor components in turn. We take the
237 : // highest-priority refinement flag in each dimension, so if any tensor
238 : // component is non-smooth, the element will split in that dimension. And only
239 : // if all tensor components are smooth enough will elements join in that
240 : // dimension.
241 : std::array<VectorType, 2> deriv_buffers{};
242 : tmpl::for_each<TensorTags>(
243 : [&result, &box, &mesh, &deriv_buffers, this](const auto tag_v) {
244 : // Stop if we have already decided to refine every dimension
245 : if (result == make_array<Dim>(Flag::Split)) {
246 : return;
247 : }
248 : using tag = tmpl::type_from<std::decay_t<decltype(tag_v)>>;
249 : const std::string tag_name = db::tag_name<tag>();
250 : // Skip if this tensor is not being monitored
251 : if (not alg::found(vars_to_monitor_, tag_name)) {
252 : return;
253 : }
254 : const auto& tensor = db::get<tag>(box);
255 : for (const VectorType& tensor_component : tensor) {
256 : Loehner_detail::max_over_components(
257 : make_not_null(&result), make_not_null(&deriv_buffers),
258 : tensor_component, mesh, relative_tolerance_, absolute_tolerance_,
259 : coarsening_factor_);
260 : }
261 : });
262 : return result;
263 : }
264 :
265 : template <size_t Dim, typename TensorTags>
266 : void Loehner<Dim, TensorTags>::pup(PUP::er& p) {
267 : p | vars_to_monitor_;
268 : p | relative_tolerance_;
269 : p | absolute_tolerance_;
270 : p | coarsening_factor_;
271 : }
272 :
273 : template <size_t Dim, typename TensorTags>
274 : PUP::able::PUP_ID Loehner<Dim, TensorTags>::my_PUP_ID = 0; // NOLINT
275 : /// \endcond
276 :
277 : } // namespace amr::Criteria
|
