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_;
94 template <
int dim,
typename Number =
double>
148 using Bounds = std::array<Number, n_bounds>;
156 : hyperbolic_system(hyperbolic_system)
157 , parameters(parameters)
158 , precomputed_values(precomputed_values)
165 void reset(
const unsigned int i,
176 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
177 const Number beta_ij,
206 const Number t_min = Number(0.),
207 const Number t_max = Number(1.));
238 Number rho_relaxation_numerator;
239 Number rho_relaxation_denominator;
249 template <
int dim,
typename Number>
250 DEAL_II_ALWAYS_INLINE
inline void
259 auto &[rho_min, rho_max, s_min] = bounds_;
261 rho_min = Number(std::numeric_limits<ScalarNumber>::max());
262 rho_max = Number(0.);
263 s_min = Number(std::numeric_limits<ScalarNumber>::max());
267 rho_relaxation_numerator = Number(0.);
268 rho_relaxation_denominator = Number(0.);
269 s_interp_max = Number(0.);
273 template <
int dim,
typename Number>
275 const unsigned int *js,
278 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
279 const Number beta_ij,
286 Assert(std::max(affine_shift.norm(), Number(0.)) == Number(0.),
287 dealii::ExcNotImplemented());
289 const auto view = hyperbolic_system.view<dim, Number>();
292 auto &[rho_min, rho_max, s_min] = bounds_;
294 const auto rho_i = view.density(U_i);
295 const auto m_i = view.momentum(U_i);
296 const auto rho_j = view.density(U_j);
297 const auto m_j = view.momentum(U_j);
298 const auto rho_affine_shift = view.density(affine_shift);
301 const auto rho_ij_bar =
302 ScalarNumber(0.5) * (rho_i + rho_j + (m_i - m_j) * scaled_c_ij) +
305 rho_min = std::min(rho_min, rho_ij_bar);
306 rho_max = std::max(rho_max, rho_ij_bar);
308 const auto &[s_j, eta_j] =
309 precomputed_values.template get_tensor<Number, precomputed_type>(js);
310 s_min = std::min(s_min, s_j);
314 rho_relaxation_numerator += beta_ij * (rho_i + rho_j);
315 rho_relaxation_denominator += std::abs(beta_ij);
317 const Number s_interp =
319 s_interp_max = std::max(s_interp_max, s_interp);
323 template <
int dim,
typename Number>
324 DEAL_II_ALWAYS_INLINE
inline auto
327 auto relaxed_bounds = bounds_;
328 auto &[rho_min, rho_max, s_min] = relaxed_bounds;
332 Number r_i = std::sqrt(hd_i);
333 if constexpr (dim == 2)
334 r_i = dealii::Utilities::fixed_power<3>(std::sqrt(r_i));
335 else if constexpr (dim == 1)
336 r_i = dealii::Utilities::fixed_power<3>(r_i);
337 r_i *= parameters.relaxation_factor();
339 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
340 const Number rho_relaxation =
341 std::abs(rho_relaxation_numerator) /
342 (std::abs(rho_relaxation_denominator) + Number(eps));
344 rho_min = std::max((Number(1.) - r_i) * rho_min,
346 rho_max = std::min((Number(1.) + r_i) * rho_max,
349 s_min = std::max((Number(1.) - r_i) * s_min,
350 Number(2.) * s_min - s_interp_max);
352 return relaxed_bounds;
356 template <
int dim,
typename Number>
357 DEAL_II_ALWAYS_INLINE
inline bool
363 AssertThrow(
false, dealii::ExcNotImplemented());
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
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, 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
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.))
static bool is_in_invariant_domain(const HyperbolicSystem &hyperbolic_system, const Bounds &bounds, const state_type &U)