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