ryujin 2.1.1 revision ed06d5bf0c7d511bfbe1eea34219071b2f01d008
limiter.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9
10#include "hyperbolic_system.h"
11
12#include <compile_time_options.h>
13#include <deal.II/base/config.h>
15#include <newton.h>
16#include <simd.h>
17
18namespace ryujin
19{
20 namespace Euler
21 {
22 template <typename ScalarNumber = double>
23 class LimiterParameters : public dealii::ParameterAcceptor
24 {
25 public:
26 LimiterParameters(const std::string &subsection = "/Limiter")
27 : ParameterAcceptor(subsection)
28 {
29 iterations_ = 2;
30 add_parameter(
31 "iterations", iterations_, "Number of limiter iterations");
32
33 if constexpr (std::is_same<ScalarNumber, double>::value)
34 newton_tolerance_ = 1.e-10;
35 else
36 newton_tolerance_ = 1.e-4;
37 add_parameter("newton tolerance",
38 newton_tolerance_,
39 "Tolerance for the quadratic newton stopping criterion");
40
41 newton_max_iterations_ = 2;
42 add_parameter("newton max iterations",
43 newton_max_iterations_,
44 "Maximal number of quadratic newton iterations performed "
45 "during limiting");
46
47 relaxation_factor_ = ScalarNumber(1.);
48 add_parameter("relaxation factor",
49 relaxation_factor_,
50 "Factor for scaling the relaxation window with r_i = "
51 "factor * (m_i/|Omega|)^(1.5/d).");
52 }
53
54 ACCESSOR_READ_ONLY(iterations);
55 ACCESSOR_READ_ONLY(newton_tolerance);
56 ACCESSOR_READ_ONLY(newton_max_iterations);
57 ACCESSOR_READ_ONLY(relaxation_factor);
58
59 private:
60 unsigned int iterations_;
61 ScalarNumber newton_tolerance_;
62 unsigned int newton_max_iterations_;
63 ScalarNumber relaxation_factor_;
64 };
65
66
97 template <int dim, typename Number = double>
98 class Limiter
99 {
100 public:
105
107
109
111
112 using state_type = typename View::state_type;
113
115
117
119
121
123
126 //
128
131 static constexpr unsigned int n_bounds = 3;
132
136 using Bounds = std::array<Number, n_bounds>;
137
141 Limiter(const HyperbolicSystem &hyperbolic_system,
142 const Parameters &parameters,
143 const PrecomputedVector &precomputed_values)
144 : hyperbolic_system(hyperbolic_system)
145 , parameters(parameters)
146 , precomputed_values(precomputed_values)
147 {
148 }
149
154 Bounds projection_bounds_from_state(const unsigned int i,
155 const state_type &U_i) const;
156
162 Bounds combine_bounds(const Bounds &bounds_left,
163 const Bounds &bounds_right) const;
164
173 Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const;
174
176
194
198 void reset(const unsigned int i,
199 const state_type &U_i,
200 const flux_contribution_type &flux_i);
201
206 void accumulate(const unsigned int *js,
207 const state_type &U_j,
208 const flux_contribution_type &flux_j,
209 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
210 const state_type &affine_shift);
211
215 Bounds bounds(const Number hd_i) const;
216
217 //*}
220
236 std::tuple<Number, bool> limit(const Bounds &bounds,
237 const state_type &U,
238 const state_type &P,
239 const Number t_min = Number(0.),
240 const Number t_max = Number(1.)) const;
241
242 private:
244
246
247 const HyperbolicSystem &hyperbolic_system;
248 const Parameters &parameters;
249 const PrecomputedVector &precomputed_values;
250
251 state_type U_i;
252
253 Bounds bounds_;
254
255 Number rho_relaxation_numerator;
256 Number rho_relaxation_denominator;
257 Number s_interp_max;
258
260 };
261
262
263 /*
264 * -------------------------------------------------------------------------
265 * Inline definitions
266 * -------------------------------------------------------------------------
267 */
268
269
270 template <int dim, typename Number>
271 DEAL_II_ALWAYS_INLINE inline auto
273 const unsigned int i, const state_type &U_i) const -> Bounds
274 {
275 const auto view = hyperbolic_system.view<dim, Number>();
276 const auto rho_i = view.density(U_i);
277 const auto &[s_i, eta_i] =
278 precomputed_values.template get_tensor<Number, precomputed_type>(i);
279
280 return {/*rho_min*/ rho_i, /*rho_max*/ rho_i, /*s_min*/ s_i};
281 }
282
283
284 template <int dim, typename Number>
285 DEAL_II_ALWAYS_INLINE inline auto Limiter<dim, Number>::combine_bounds(
286 const Bounds &bounds_left, const Bounds &bounds_right) const -> Bounds
287 {
288 const auto &[rho_min_l, rho_max_l, s_min_l] = bounds_left;
289 const auto &[rho_min_r, rho_max_r, s_min_r] = bounds_right;
290
291 return {std::min(rho_min_l, rho_min_r),
292 std::max(rho_max_l, rho_max_r),
293 std::min(s_min_l, s_min_r)};
294 }
295
296
297 template <int dim, typename Number>
298 DEAL_II_ALWAYS_INLINE inline auto
300 const Number &hd) const -> Bounds
301 {
302 auto relaxed_bounds = bounds;
303 auto &[rho_min, rho_max, s_min] = relaxed_bounds;
304
305 /* Use r = factor * (m_i / |Omega|) ^ (1.5 / d): */
306
307 Number r = std::sqrt(hd); // in 3D: ^ 3/6
308 if constexpr (dim == 2) //
309 r = dealii::Utilities::fixed_power<3>(std::sqrt(r)); // in 2D: ^ 3/4
310 else if constexpr (dim == 1) //
311 r = dealii::Utilities::fixed_power<3>(r); // in 1D: ^ 3/2
312 r *= parameters.relaxation_factor();
313
314 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
315 rho_min *= std::max(Number(1.) - r, Number(eps));
316 rho_max *= (Number(1.) + r);
317 s_min *= std::max(Number(1.) - r, Number(eps));
318
319 return relaxed_bounds;
320 }
321
322
323 template <int dim, typename Number>
324 DEAL_II_ALWAYS_INLINE inline void
325 Limiter<dim, Number>::reset(const unsigned int /*i*/,
326 const state_type &new_U_i,
327 const flux_contribution_type & /*new_flux_i*/)
328 {
329 U_i = new_U_i;
330
331 /* Bounds: */
332
333 auto &[rho_min, rho_max, s_min] = bounds_;
334
335 rho_min = Number(std::numeric_limits<ScalarNumber>::max());
336 rho_max = Number(0.);
337 s_min = Number(std::numeric_limits<ScalarNumber>::max());
338
339 /* Relaxation: */
340
341 rho_relaxation_numerator = Number(0.);
342 rho_relaxation_denominator = Number(0.);
343 s_interp_max = Number(0.);
344 }
345
346
347 template <int dim, typename Number>
348 DEAL_II_ALWAYS_INLINE inline void Limiter<dim, Number>::accumulate(
349 const unsigned int *js,
350 const state_type &U_j,
351 const flux_contribution_type & /*flux_j*/,
352 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
353 const state_type &affine_shift)
354 {
355 // TODO: Currently we only apply the affine_shift to U_ij_bar (which
356 // then enters all bounds), but we do not modify s_interp and
357 // rho_relaxation. When actually adding a source term to the Euler
358 // equations verify that this does the right thing.
359 Assert(std::max(affine_shift.norm(), Number(0.)) == Number(0.),
360 dealii::ExcNotImplemented());
361
362 const auto view = hyperbolic_system.view<dim, Number>();
363
364 /* Bounds: */
365 auto &[rho_min, rho_max, s_min] = bounds_;
366
367 const auto rho_i = view.density(U_i);
368 const auto m_i = view.momentum(U_i);
369 const auto rho_j = view.density(U_j);
370 const auto m_j = view.momentum(U_j);
371 const auto rho_affine_shift = view.density(affine_shift);
372
373 /* bar state shifted by an affine shift: */
374 const auto rho_ij_bar =
375 ScalarNumber(0.5) * (rho_i + rho_j + (m_i - m_j) * scaled_c_ij) +
376 rho_affine_shift;
377
378 rho_min = std::min(rho_min, rho_ij_bar);
379 rho_max = std::max(rho_max, rho_ij_bar);
380
381 const auto &[s_j, eta_j] =
382 precomputed_values.template get_tensor<Number, precomputed_type>(js);
383 s_min = std::min(s_min, s_j);
384
385 /* Relaxation: */
386
387 /* Use a uniform weight. */
388 const auto beta_ij = Number(1.);
389 rho_relaxation_numerator += beta_ij * (rho_i + rho_j);
390 rho_relaxation_denominator += std::abs(beta_ij);
391
392 const Number s_interp =
393 view.specific_entropy((U_i + U_j) * ScalarNumber(.5));
394 s_interp_max = std::max(s_interp_max, s_interp);
395 }
396
397
398 template <int dim, typename Number>
399 DEAL_II_ALWAYS_INLINE inline auto
400 Limiter<dim, Number>::bounds(const Number hd_i) const -> Bounds
401 {
402 const auto &[rho_min, rho_max, s_min] = bounds_;
403
404 auto relaxed_bounds = fully_relax_bounds(bounds_, hd_i);
405 auto &[rho_min_relaxed, rho_max_relaxed, s_min_relaxed] = relaxed_bounds;
406
407 /* Apply a stricter window: */
408
409 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
410
411 const auto rho_relaxation =
412 ScalarNumber(2. * parameters.relaxation_factor()) *
413 std::abs(rho_relaxation_numerator) /
414 (std::abs(rho_relaxation_denominator) + Number(eps));
415
416 const auto entropy_relaxation =
417 parameters.relaxation_factor() * (s_interp_max - s_min);
418
419 rho_min_relaxed = std::max(rho_min_relaxed, rho_min - rho_relaxation);
420 rho_max_relaxed = std::min(rho_max_relaxed, rho_max + rho_relaxation);
421 s_min_relaxed = std::max(s_min_relaxed, s_min - entropy_relaxation);
422
423 return relaxed_bounds;
424 }
425 } // namespace Euler
426} // namespace ryujin
typename get_value_type< Number >::type ScalarNumber
Vectors::MultiComponentVector< ScalarNumber, n_precomputed_values > PrecomputedVector
static constexpr unsigned int problem_dimension
dealii::Tensor< 1, problem_dimension, Number > state_type
std::array< Number, n_precomputed_values > precomputed_type
ACCESSOR_READ_ONLY(newton_max_iterations)
ACCESSOR_READ_ONLY(newton_tolerance)
ACCESSOR_READ_ONLY(relaxation_factor)
LimiterParameters(const std::string &subsection="/Limiter")
Definition: limiter.h:26
Limiter(const HyperbolicSystem &hyperbolic_system, const Parameters &parameters, const PrecomputedVector &precomputed_values)
Definition: limiter.h:141
void reset(const unsigned int i, const state_type &U_i, const flux_contribution_type &flux_i)
Definition: limiter.h:325
typename View::PrecomputedVector PrecomputedVector
Definition: limiter.h:118
static constexpr auto problem_dimension
Definition: limiter.h:110
Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const
Definition: limiter.h:299
std::tuple< Number, bool > limit(const Bounds &bounds, const state_type &U, const state_type &P, const Number t_min=Number(0.), const Number t_max=Number(1.)) const
void accumulate(const unsigned int *js, const state_type &U_j, const flux_contribution_type &flux_j, const dealii::Tensor< 1, dim, Number > &scaled_c_ij, const state_type &affine_shift)
Definition: limiter.h:348
typename View::precomputed_type precomputed_type
Definition: limiter.h:116
typename View::ScalarNumber ScalarNumber
Definition: limiter.h:108
typename View::state_type state_type
Definition: limiter.h:112
LimiterParameters< ScalarNumber > Parameters
Definition: limiter.h:120
Bounds bounds(const Number hd_i) const
Definition: limiter.h:400
static constexpr unsigned int n_bounds
Definition: limiter.h:131
std::array< Number, n_bounds > Bounds
Definition: limiter.h:136
typename View::flux_contribution_type flux_contribution_type
Definition: limiter.h:114
Bounds combine_bounds(const Bounds &bounds_left, const Bounds &bounds_right) const
Definition: limiter.h:285
Bounds projection_bounds_from_state(const unsigned int i, const state_type &U_i) const
Definition: limiter.h:272