8#include <compile_time_options.h>
12#include <compile_time_options.h>
13#include <deal.II/base/config.h>
22 template <
typename ScalarNumber =
double>
27 : ParameterAcceptor(subsection)
31 "iterations", iterations_,
"Number of limiter iterations");
33 if constexpr (std::is_same<ScalarNumber, double>::value)
34 newton_tolerance_ = 1.e-10;
36 newton_tolerance_ = 1.e-4;
37 add_parameter(
"newton tolerance",
39 "Tolerance for the quadratic newton stopping criterion");
41 newton_max_iterations_ = 2;
42 add_parameter(
"newton max iterations",
43 newton_max_iterations_,
44 "Maximal number of quadratic newton iterations performed "
47 relaxation_factor_ = ScalarNumber(1.);
48 add_parameter(
"relaxation factor",
50 "Factor for scaling the relaxation window with r_i = "
51 "factor * (m_i/|Omega|)^(1.5/d).");
60 unsigned int iterations_;
61 ScalarNumber newton_tolerance_;
62 unsigned int newton_max_iterations_;
63 ScalarNumber relaxation_factor_;
97 template <
int dim,
typename Number =
double>
136 using Bounds = std::array<Number, n_bounds>;
144 : hyperbolic_system(hyperbolic_system)
145 , parameters(parameters)
146 , precomputed_values(precomputed_values)
163 const Bounds &bounds_right)
const;
198 void reset(
const unsigned int i,
209 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
239 const Number t_min = Number(0.),
240 const Number t_max = Number(1.))
const;
255 Number rho_relaxation_numerator;
256 Number rho_relaxation_denominator;
270 template <
int dim,
typename Number>
271 DEAL_II_ALWAYS_INLINE
inline auto
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);
280 return { rho_i, rho_i, s_i};
284 template <
int dim,
typename Number>
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;
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)};
297 template <
int dim,
typename Number>
298 DEAL_II_ALWAYS_INLINE
inline auto
300 const Number &hd)
const ->
Bounds
302 auto relaxed_bounds = bounds;
303 auto &[rho_min, rho_max, s_min] = relaxed_bounds;
307 Number r = std::sqrt(hd);
308 if constexpr (dim == 2)
309 r = dealii::Utilities::fixed_power<3>(std::sqrt(r));
310 else if constexpr (dim == 1)
311 r = dealii::Utilities::fixed_power<3>(r);
312 r *= parameters.relaxation_factor();
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));
319 return relaxed_bounds;
323 template <
int dim,
typename Number>
324 DEAL_II_ALWAYS_INLINE
inline void
333 auto &[rho_min, rho_max, s_min] = bounds_;
335 rho_min = Number(std::numeric_limits<ScalarNumber>::max());
336 rho_max = Number(0.);
337 s_min = Number(std::numeric_limits<ScalarNumber>::max());
341 rho_relaxation_numerator = Number(0.);
342 rho_relaxation_denominator = Number(0.);
343 s_interp_max = Number(0.);
347 template <
int dim,
typename Number>
349 const unsigned int *js,
352 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
359 Assert(std::max(affine_shift.norm(), Number(0.)) == Number(0.),
360 dealii::ExcNotImplemented());
362 const auto view = hyperbolic_system.view<dim, Number>();
365 auto &[rho_min, rho_max, s_min] = bounds_;
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);
374 const auto rho_ij_bar =
375 ScalarNumber(0.5) * (rho_i + rho_j + (m_i - m_j) * scaled_c_ij) +
378 rho_min = std::min(rho_min, rho_ij_bar);
379 rho_max = std::max(rho_max, rho_ij_bar);
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);
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);
392 const Number s_interp =
394 s_interp_max = std::max(s_interp_max, s_interp);
398 template <
int dim,
typename Number>
399 DEAL_II_ALWAYS_INLINE
inline auto
402 const auto &[rho_min, rho_max, s_min] = bounds_;
404 auto relaxed_bounds = fully_relax_bounds(bounds_, hd_i);
405 auto &[rho_min_relaxed, rho_max_relaxed, s_min_relaxed] = relaxed_bounds;
409 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
411 const auto rho_relaxation =
413 std::abs(rho_relaxation_numerator) /
414 (std::abs(rho_relaxation_denominator) + Number(eps));
416 const auto entropy_relaxation =
417 parameters.relaxation_factor() * (s_interp_max - s_min);
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);
423 return relaxed_bounds;
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
flux_type flux_contribution_type
ACCESSOR_READ_ONLY(iterations)
ACCESSOR_READ_ONLY(newton_max_iterations)
ACCESSOR_READ_ONLY(newton_tolerance)
ACCESSOR_READ_ONLY(relaxation_factor)
LimiterParameters(const std::string &subsection="/Limiter")
Limiter(const HyperbolicSystem &hyperbolic_system, const Parameters ¶meters, const PrecomputedVector &precomputed_values)
void reset(const unsigned int i, const state_type &U_i, const flux_contribution_type &flux_i)
typename View::PrecomputedVector PrecomputedVector
static constexpr auto problem_dimension
Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const
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)
typename View::precomputed_type precomputed_type
typename View::ScalarNumber ScalarNumber
typename View::state_type state_type
LimiterParameters< ScalarNumber > Parameters
Bounds bounds(const Number hd_i) const
static constexpr unsigned int n_bounds
std::array< Number, n_bounds > Bounds
typename View::flux_contribution_type flux_contribution_type
Bounds combine_bounds(const Bounds &bounds_left, const Bounds &bounds_right) const
Bounds projection_bounds_from_state(const unsigned int i, const state_type &U_i) const