10#include <compile_time_options.h>
19 namespace ShallowWater
21 template <
typename ScalarNumber =
double>
26 : ParameterAcceptor(subsection)
30 "iterations", iterations_,
"Number of limiter iterations");
32 if constexpr (std::is_same_v<ScalarNumber, double>)
33 newton_tolerance_ = 1.e-10;
35 newton_tolerance_ = 1.e-4;
36 add_parameter(
"newton tolerance",
38 "Tolerance for the quadratic newton stopping criterion");
40 newton_max_iterations_ = 2;
41 add_parameter(
"newton max iterations",
42 newton_max_iterations_,
43 "Maximal number of quadratic newton iterations performed "
46 relaxation_factor_ = ScalarNumber(1.);
47 add_parameter(
"relaxation factor",
49 "Factor for scaling the relaxation window with r_i = "
50 "factor * (m_i/|Omega|)^(1.5/d).");
59 unsigned int iterations_;
60 ScalarNumber newton_tolerance_;
61 unsigned int newton_max_iterations_;
62 ScalarNumber relaxation_factor_;
71 template <
int dim,
typename Number =
double>
109 using Bounds = std::array<Number, n_bounds>;
117 : hyperbolic_system(hyperbolic_system)
118 , parameters(parameters)
119 , precomputed_values(precomputed_values)
136 const Bounds &bounds_right)
const;
170 void reset(
const unsigned int i,
181 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
202 const Number t_min = Number(0.),
203 const Number t_max = Number(1.))
const;
220 Number h_relaxation_numerator;
221 Number v2_relaxation_numerator;
222 Number relaxation_denominator;
235 template <
int dim,
typename Number>
236 DEAL_II_ALWAYS_INLINE
inline auto
240 const auto view = hyperbolic_system.view<dim, Number>();
241 const auto h_i = view.water_depth(U_i);
243 view.momentum(U_i) * view.inverse_water_depth_mollified(U_i);
244 const auto v2_i = v_i.norm_square();
246 return { h_i, h_i, v2_i};
250 template <
int dim,
typename Number>
251 DEAL_II_ALWAYS_INLINE
inline auto
255 const auto &[h_min_l, h_max_l, v2_max_l] = bounds_l;
256 const auto &[h_min_r, h_max_r, v2_max_r] = bounds_r;
258 return {std::min(h_min_l, h_min_r),
259 std::max(h_max_l, h_max_r),
260 std::max(v2_max_l, v2_max_r)};
264 template <
int dim,
typename Number>
265 DEAL_II_ALWAYS_INLINE
inline auto
267 const Number &hd)
const ->
Bounds
269 auto relaxed_bounds = bounds;
270 auto &[h_min, h_max, v2_max] = relaxed_bounds;
274 Number r = std::sqrt(hd);
275 if constexpr (dim == 2)
276 r = dealii::Utilities::fixed_power<3>(std::sqrt(r));
277 else if constexpr (dim == 1)
278 r = dealii::Utilities::fixed_power<3>(r);
279 r *= parameters.relaxation_factor();
281 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
282 h_min *= std::max((Number(1.) - r), Number(eps));
283 h_max *= (Number(1.) + r);
284 v2_max *= (Number(1.) + r);
286 return relaxed_bounds;
290 template <
int dim,
typename Number>
291 DEAL_II_ALWAYS_INLINE
inline void
298 auto &[h_min, h_max, v2_max] = bounds_;
300 h_min = Number(std::numeric_limits<ScalarNumber>::max());
304 h_relaxation_numerator = Number(0.);
305 v2_relaxation_numerator = Number(0.);
306 relaxation_denominator = Number(0.);
310 template <
int dim,
typename Number>
315 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
318 const auto view = hyperbolic_system.view<dim, Number>();
322 const auto f_star_ij = view.f(U_star_ij);
323 const auto f_star_ji = view.f(U_star_ji);
326 const auto U_ij_bar =
328 (U_star_ij + U_star_ji +
329 contract(
add(f_star_ij, -f_star_ji), scaled_c_ij)) +
334 auto &[h_min, h_max, v2_max] = bounds_;
336 const auto h_bar_ij = view.water_depth(U_ij_bar);
337 h_min = std::min(h_min, h_bar_ij);
338 h_max = std::max(h_max, h_bar_ij);
340 const auto v_bar_ij = view.momentum(U_ij_bar) *
341 view.inverse_water_depth_mollified(U_ij_bar);
342 const auto v2_bar_ij = v_bar_ij.norm_square();
343 v2_max = std::max(v2_max, v2_bar_ij);
348 const auto beta_ij = Number(1.);
350 relaxation_denominator += std::abs(beta_ij);
352 const auto h_i = view.water_depth(U_i);
353 const auto h_j = view.water_depth(U_j);
354 h_relaxation_numerator += beta_ij * (h_i + h_j);
357 view.momentum(U_i) * view.inverse_water_depth_mollified(U_i);
359 view.momentum(U_j) * view.inverse_water_depth_mollified(U_j);
360 v2_relaxation_numerator +=
361 beta_ij * (-vel_i.norm_square() + vel_j.norm_square());
365 template <
int dim,
typename Number>
366 DEAL_II_ALWAYS_INLINE
inline auto
369 const auto &[h_min, h_max, v2_max] = bounds_;
371 auto relaxed_bounds = fully_relax_bounds(bounds_, hd_i);
372 auto &[h_min_relaxed, h_max_relaxed, v2_max_relaxed] = relaxed_bounds;
376 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
378 const Number h_relaxed =
380 std::abs(h_relaxation_numerator) /
381 (relaxation_denominator + Number(eps));
383 const Number v2_relaxed =
385 std::abs(v2_relaxation_numerator) /
386 (relaxation_denominator + Number(eps));
388 h_min_relaxed = std::max(h_min_relaxed, h_min - h_relaxed);
389 h_max_relaxed = std::min(h_max_relaxed, h_max + h_relaxed);
390 v2_max_relaxed = std::min(v2_max_relaxed, v2_max + v2_relaxed);
392 return relaxed_bounds;
typename get_value_type< Number >::type ScalarNumber
Vectors::MultiComponentVector< ScalarNumber, n_precomputed_values > PrecomputedVector
dealii::Tensor< 1, problem_dimension, Number > state_type
std::array< Number, n_precomputed_values > precomputed_type
static constexpr unsigned int problem_dimension
std::tuple< state_type, Number > flux_contribution_type
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")
Bounds projection_bounds_from_state(const unsigned int i, const state_type &U_i) const
typename View::ScalarNumber ScalarNumber
typename View::flux_contribution_type flux_contribution_type
Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const
std::array< Number, n_bounds > Bounds
typename View::PrecomputedVector PrecomputedVector
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 state_type &affine_shift)
Bounds bounds(const Number hd_i) 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
Limiter(const HyperbolicSystem &hyperbolic_system, const Parameters ¶meters, const PrecomputedVector &precomputed_values)
typename View::precomputed_type precomputed_type
Bounds combine_bounds(const Bounds &bounds_left, const Bounds &bounds_right) const
void reset(const unsigned int i, const state_type &U_i, const flux_contribution_type &flux_i)
typename View::state_type state_type
static constexpr auto problem_dimension
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)