12#include <compile_time_options.h>
18 namespace ShallowWater
20 template <
typename ScalarNumber =
double>
25 : ParameterAcceptor(subsection)
29 "iterations", iterations_,
"Number of limiter iterations");
31 if constexpr (std::is_same<ScalarNumber, double>::value)
32 newton_tolerance_ = 1.e-10;
34 newton_tolerance_ = 1.e-4;
35 add_parameter(
"newton tolerance",
37 "Tolerance for the quadratic newton stopping criterion");
39 newton_max_iterations_ = 2;
40 add_parameter(
"newton max iterations",
41 newton_max_iterations_,
42 "Maximal number of quadratic newton iterations performed "
45 relaxation_factor_ = ScalarNumber(1.);
46 add_parameter(
"relaxation factor",
48 "Factor for scaling the relaxation window with r_i = "
49 "factor * (m_i/|Omega|)^(1.5/d).");
51 limit_on_kinetic_energy_ =
false;
52 add_parameter(
"limit on kinetic energy",
53 limit_on_kinetic_energy_,
54 "Limit on kinetic energy");
56 limit_on_square_velocity_ =
true;
57 add_parameter(
"limit on square velocity",
58 limit_on_square_velocity_,
59 "Limit on square velocity");
71 unsigned int iterations_;
72 ScalarNumber newton_tolerance_;
73 unsigned int newton_max_iterations_;
74 ScalarNumber relaxation_factor_;
76 bool limit_on_kinetic_energy_;
77 bool limit_on_square_velocity_;
86 template <
int dim,
typename Number =
double>
149 using Bounds = std::array<Number, n_bounds>;
158 : hyperbolic_system(hyperbolic_system)
159 , parameters(parameters)
160 , precomputed_values(precomputed_values)
167 void reset(
const unsigned int i,
178 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
179 const Number &beta_ij,
200 const Number t_min = Number(0.),
201 const Number t_max = Number(1.));
237 Number h_relaxation_numerator;
238 Number kin_relaxation_numerator;
239 Number v2_relaxation_numerator;
240 Number relaxation_denominator;
253 template <
int dim,
typename Number>
254 DEAL_II_ALWAYS_INLINE
inline void
261 auto &[h_min, h_max, h_small, kin_max, v2_max] = bounds_;
263 h_min = Number(std::numeric_limits<ScalarNumber>::max());
265 h_small = Number(0.);
266 kin_max = Number(0.);
269 h_relaxation_numerator = Number(0.);
270 kin_relaxation_numerator = Number(0.);
271 v2_relaxation_numerator = Number(0.);
272 relaxation_denominator = Number(0.);
276 template <
int dim,
typename Number>
281 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
282 const Number &beta_ij,
285 const auto view = hyperbolic_system.view<dim, Number>();
289 const auto f_star_ij = view.f(U_star_ij);
290 const auto f_star_ji = view.f(U_star_ji);
293 (U_star_ij + U_star_ji +
294 contract(
add(f_star_ij, -f_star_ji), scaled_c_ij));
296 U_ij_bar += affine_shift;
300 auto &[h_min, h_max, h_small, kin_max, v2_max] = bounds_;
302 const auto h_bar_ij = view.water_depth(U_ij_bar);
303 h_min = std::min(h_min, h_bar_ij);
304 h_max = std::max(h_max, h_bar_ij);
306 const auto kin_bar_ij = view.kinetic_energy(U_ij_bar);
307 kin_max = std::max(kin_max, kin_bar_ij);
309 const auto v_bar_ij = view.momentum(U_ij_bar) *
310 view.inverse_water_depth_mollified(U_ij_bar);
311 const auto v2_bar_ij = v_bar_ij.norm_square();
312 v2_max = std::max(v2_max, v2_bar_ij);
316 relaxation_denominator += std::abs(beta_ij);
318 const auto h_i = view.water_depth(U_i);
319 const auto h_j = view.water_depth(U_j);
320 h_relaxation_numerator += beta_ij * (h_i + h_j);
322 const auto kin_i = view.kinetic_energy(U_i);
323 const auto kin_j = view.kinetic_energy(U_j);
324 kin_relaxation_numerator += beta_ij * (kin_i + kin_j);
327 view.momentum(U_i) * view.inverse_water_depth_mollified(U_i);
329 view.momentum(U_j) * view.inverse_water_depth_mollified(U_j);
330 v2_relaxation_numerator +=
331 beta_ij * (-vel_i.norm_square() + vel_j.norm_square());
335 template <
int dim,
typename Number>
336 DEAL_II_ALWAYS_INLINE
inline auto
339 const auto view = hyperbolic_system.view<dim, Number>();
341 auto relaxed_bounds = bounds_;
342 auto &[h_min, h_max, h_small, kin_max, v2_max] = relaxed_bounds;
346 Number r_i = std::sqrt(hd_i);
347 if constexpr (dim == 2)
348 r_i = dealii::Utilities::fixed_power<3>(std::sqrt(r_i));
349 else if constexpr (dim == 1)
350 r_i = dealii::Utilities::fixed_power<3>(r_i);
351 r_i *= parameters.relaxation_factor();
353 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
356 std::abs(h_relaxation_numerator) /
357 (relaxation_denominator + Number(eps));
359 h_min = std::max((Number(1.) - r_i) * h_min, h_min - h_relaxed);
360 h_max = std::min((Number(1.) + r_i) * h_max, h_max + h_relaxed);
363 std::abs(kin_relaxation_numerator) /
364 (relaxation_denominator + Number(eps));
366 kin_max = std::min((Number(1.) + r_i) * kin_max, kin_max + kin_relaxed);
369 std::abs(v2_relaxation_numerator) /
370 (relaxation_denominator + Number(eps));
372 v2_max = std::min((Number(1.) + r_i) * v2_max, v2_max + v2_relaxed);
377 if constexpr (dim == 2)
378 r_i = std::sqrt(hd_i);
379 r_i *= view.dry_state_relaxation_factor();
381 h_small = view.reference_water_depth() * r_i;
383 return relaxed_bounds;
387 template <
int dim,
typename Number>
388 DEAL_II_ALWAYS_INLINE
inline bool
394 AssertThrow(
false, dealii::ExcNotImplemented());
static constexpr unsigned int n_precomputed_values
dealii::Tensor< 1, problem_dimension, Number > state_type
std::tuple< state_type, Number > flux_contribution_type
ACCESSOR_READ_ONLY(limit_on_square_velocity)
ACCESSOR_READ_ONLY(iterations)
ACCESSOR_READ_ONLY(relaxation_factor)
ACCESSOR_READ_ONLY(newton_tolerance)
ACCESSOR_READ_ONLY(newton_max_iterations)
LimiterParameters(const std::string &subsection="/Limiter")
ACCESSOR_READ_ONLY(limit_on_kinetic_energy)
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.))
typename View::flux_contribution_type flux_contribution_type
std::array< Number, n_bounds > Bounds
static constexpr unsigned int n_precomputed_values
void accumulate(const state_type &U_j, const state_type &U_star_ij, const state_type &U_star_ji, const dealii::Tensor< 1, dim, Number > &scaled_c_ij, const Number &beta_ij, const state_type &affine_shift)
typename get_value_type< Number >::type ScalarNumber
Bounds bounds(const Number hd_i) const
void reset(const unsigned int i, const state_type &U_i, const flux_contribution_type &flux_i)
typename View::state_type state_type
Limiter(const HyperbolicSystem &hyperbolic_system, const Parameters ¶meters, const MultiComponentVector< ScalarNumber, n_precomputed_values > &precomputed_values)
static bool is_in_invariant_domain(const HyperbolicSystem &, const Bounds &, const state_type &)
static constexpr unsigned int n_bounds
LimiterParameters< ScalarNumber > Parameters
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)