10#include <compile_time_options.h>
19 template <
typename ScalarNumber =
double>
24 : ParameterAcceptor(subsection)
28 "iterations", iterations_,
"Number of limiter iterations");
30 if constexpr (std::is_same<ScalarNumber, double>::value)
31 newton_tolerance_ = 1.e-10;
33 newton_tolerance_ = 1.e-4;
34 add_parameter(
"newton tolerance",
36 "Tolerance for the quadratic newton stopping criterion");
38 newton_max_iterations_ = 2;
39 add_parameter(
"newton max iterations",
40 newton_max_iterations_,
41 "Maximal number of quadratic newton iterations performed "
44 relaxation_factor_ = ScalarNumber(1.);
45 add_parameter(
"relaxation factor",
47 "Factor for scaling the relaxation window with r_i = "
48 "factor * (m_i/|Omega|)^(1.5/d).");
57 unsigned int iterations_;
58 ScalarNumber newton_tolerance_;
59 unsigned int newton_max_iterations_;
60 ScalarNumber relaxation_factor_;
96 template <
int dim,
typename Number =
double>
168 using Bounds = std::array<Number, n_bounds>;
177 : hyperbolic_system(hyperbolic_system)
178 , parameters(parameters)
179 , precomputed_values(precomputed_values)
186 void reset(
const unsigned int i,
197 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
198 const Number beta_ij);
226 const Number t_min = Number(0.),
227 const Number t_max = Number(1.));
261 Number rho_relaxation_numerator;
262 Number rho_relaxation_denominator;
276 template <
int dim,
typename Number>
277 DEAL_II_ALWAYS_INLINE
inline void
287 auto &[rho_min, rho_max, s_min, gamma_min] = bounds_;
289 rho_min = Number(std::numeric_limits<ScalarNumber>::max());
290 rho_max = Number(0.);
291 s_min = Number(std::numeric_limits<ScalarNumber>::max());
293 const auto &[p_i, gamma_min_i, s_i, eta_i] =
295 .template get_tensor<Number, precomputed_state_type>(i);
297 gamma_min = gamma_min_i;
301 rho_relaxation_numerator = Number(0.);
302 rho_relaxation_denominator = Number(0.);
303 s_interp_max = Number(0.);
307 template <
int dim,
typename Number>
309 const unsigned int *js,
312 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
313 const Number beta_ij)
315 const auto view = hyperbolic_system.view<dim, Number>();
318 auto &[rho_min, rho_max, s_min, gamma_min] = bounds_;
320 const auto rho_i = view.density(U_i);
321 const auto rho_j = view.density(U_j);
323 const auto U_ij_bar =
327 const auto rho_ij_bar = view.density(U_ij_bar);
331 rho_min = std::min(rho_min, rho_ij_bar);
332 rho_max = std::max(rho_max, rho_ij_bar);
336 rho_relaxation_numerator += beta_ij * (rho_i + rho_j);
337 rho_relaxation_denominator += std::abs(beta_ij);
341 if (view.compute_strict_bounds()) {
352 const auto s_j = view.surrogate_specific_entropy(U_j, gamma_min);
354 const auto s_ij_bar =
355 view.surrogate_specific_entropy(U_ij_bar, gamma_min);
357 const Number s_interp = view.surrogate_specific_entropy(
360 s_min = std::min(s_min, s_j);
361 s_min = std::min(s_min, s_ij_bar);
362 s_interp_max = std::max(s_interp_max, s_interp);
371 const auto [p_j, gamma_min_j, s_j, eta_j] =
373 .template get_tensor<Number, precomputed_state_type>(js);
375 const auto s_ij_bar =
376 view.surrogate_specific_entropy(U_ij_bar, gamma_min);
378 s_min = std::min(s_min, s_j);
379 s_min = std::min(s_min, s_ij_bar);
380 s_interp_max = std::max(s_interp_max, s_ij_bar);
385 template <
int dim,
typename Number>
386 DEAL_II_ALWAYS_INLINE
inline auto
389 const auto view = hyperbolic_system.view<dim, Number>();
391 auto relaxed_bounds = bounds_;
392 auto &[rho_min, rho_max, s_min, gamma_min] = relaxed_bounds;
396 Number r_i = std::sqrt(hd_i);
397 if constexpr (dim == 2)
398 r_i = dealii::Utilities::fixed_power<3>(std::sqrt(r_i));
399 else if constexpr (dim == 1)
400 r_i = dealii::Utilities::fixed_power<3>(r_i);
401 r_i *= parameters.relaxation_factor();
403 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
404 const Number rho_relaxation =
405 std::abs(rho_relaxation_numerator) /
406 (std::abs(rho_relaxation_denominator) + Number(eps));
408 rho_min = std::max((Number(1.) - r_i) * rho_min,
411 rho_max = std::min((Number(1.) + r_i) * rho_max,
414 s_min = std::max((Number(1.) - r_i) * s_min,
415 Number(2.) * s_min - s_interp_max);
423 const auto numerator = (gamma_min + Number(1.)) * rho_max;
424 const auto interpolation_b = view.eos_interpolation_b();
425 const auto denominator =
426 gamma_min - Number(1.) +
ScalarNumber(2.) * interpolation_b * rho_max;
427 const auto upper_bound = numerator / denominator;
429 rho_max = std::min(upper_bound, rho_max);
431 return relaxed_bounds;
435 template <
int dim,
typename Number>
436 DEAL_II_ALWAYS_INLINE
inline bool
442 AssertThrow(
false, dealii::ExcNotImplemented());
dealii::Tensor< 1, problem_dimension, Number > state_type
std::array< Number, n_precomputed_values > precomputed_state_type
typename get_value_type< Number >::type ScalarNumber
static constexpr unsigned int n_precomputed_values
static constexpr unsigned int problem_dimension
flux_type flux_contribution_type
LimiterParameters(const std::string &subsection="/Limiter")
ACCESSOR_READ_ONLY(iterations)
ACCESSOR_READ_ONLY(newton_max_iterations)
ACCESSOR_READ_ONLY(relaxation_factor)
ACCESSOR_READ_ONLY(newton_tolerance)
typename View::state_type state_type
static constexpr unsigned int n_bounds
typename View::ScalarNumber ScalarNumber
static constexpr unsigned int n_precomputed_values
static constexpr unsigned int problem_dimension
typename View::precomputed_state_type precomputed_state_type
typename View::flux_contribution_type flux_contribution_type
void reset(const unsigned int i, const state_type &U_i, const flux_contribution_type &flux_i)
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 Number beta_ij)
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.))
Limiter(const HyperbolicSystem &hyperbolic_system, const Parameters ¶meters, const MultiComponentVector< ScalarNumber, n_precomputed_values > &precomputed_values)
std::array< Number, n_bounds > Bounds
Bounds bounds(const Number hd_i) const
LimiterParameters< ScalarNumber > Parameters
static bool is_in_invariant_domain(const HyperbolicSystem &hyperbolic_system, const Bounds &bounds, const state_type &U)
DEAL_II_ALWAYS_INLINE FT add(const FT &flux_left_ij, const FT &flux_right_ij)
DEAL_II_ALWAYS_INLINE dealii::Tensor< 1, problem_dim, T > contract(const FT &flux_ij, const TT &c_ij)