ryujin 2.1.1 revision c38afcfdec63b34961924f05408760247496a1f0
limiter.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0
3// [LANL Copyright Statement]
4// Copyright (C) 2023 - 2024 by the ryujin authors
5// Copyright (C) 2023 - 2024 by Triad National Security, LLC
6//
7
8#pragma once
9
10#include "hyperbolic_system.h"
11
12#include <compile_time_options.h>
14#include <newton.h>
15
16namespace ryujin
17{
18 namespace ShallowWater
19 {
20 template <typename ScalarNumber = double>
21 class LimiterParameters : public dealii::ParameterAcceptor
22 {
23 public:
24 LimiterParameters(const std::string &subsection = "/Limiter")
25 : ParameterAcceptor(subsection)
26 {
27 iterations_ = 2;
28 add_parameter(
29 "iterations", iterations_, "Number of limiter iterations");
30
31 if constexpr (std::is_same<ScalarNumber, double>::value)
32 newton_tolerance_ = 1.e-10;
33 else
34 newton_tolerance_ = 1.e-4;
35 add_parameter("newton tolerance",
36 newton_tolerance_,
37 "Tolerance for the quadratic newton stopping criterion");
38
39 newton_max_iterations_ = 2;
40 add_parameter("newton max iterations",
41 newton_max_iterations_,
42 "Maximal number of quadratic newton iterations performed "
43 "during limiting");
44
45 relaxation_factor_ = ScalarNumber(1.);
46 add_parameter("relaxation factor",
47 relaxation_factor_,
48 "Factor for scaling the relaxation window with r_i = "
49 "factor * (m_i/|Omega|)^(1.5/d).");
50
51 limit_on_kinetic_energy_ = false;
52 add_parameter("limit on kinetic energy",
53 limit_on_kinetic_energy_,
54 "Limit on kinetic energy");
55
56 limit_on_square_velocity_ = true;
57 add_parameter("limit on square velocity",
58 limit_on_square_velocity_,
59 "Limit on square velocity");
60 }
61
62 ACCESSOR_READ_ONLY(iterations);
63 ACCESSOR_READ_ONLY(newton_tolerance);
64 ACCESSOR_READ_ONLY(newton_max_iterations);
65 ACCESSOR_READ_ONLY(relaxation_factor);
66
67 ACCESSOR_READ_ONLY(limit_on_kinetic_energy);
68 ACCESSOR_READ_ONLY(limit_on_square_velocity);
69
70 private:
71 unsigned int iterations_;
72 ScalarNumber newton_tolerance_;
73 unsigned int newton_max_iterations_;
74 ScalarNumber relaxation_factor_;
75
76 bool limit_on_kinetic_energy_;
77 bool limit_on_square_velocity_;
78 };
79
80
86 template <int dim, typename Number = double>
87 class Limiter
88 {
89 public:
94
96
98
100
101 using state_type = typename View::state_type;
102
104
106
108
110
112
130
134 static constexpr unsigned int n_bounds = 5;
135
139 using Bounds = std::array<Number, n_bounds>;
140
144 Limiter(const HyperbolicSystem &hyperbolic_system,
145 const Parameters &parameters,
146 const PrecomputedVector &precomputed_values)
147 : hyperbolic_system(hyperbolic_system)
148 , parameters(parameters)
149 , precomputed_values(precomputed_values)
150 {
151 }
152
156 void reset(const unsigned int i,
157 const state_type &U_i,
158 const flux_contribution_type &flux_i);
159
164 void accumulate(const state_type &U_j,
165 const state_type &U_star_ij,
166 const state_type &U_star_ji,
167 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
168 const state_type &affine_shift);
169
173 Bounds bounds(const Number hd_i) const;
174
180 static Bounds combine_bounds(const Bounds &bounds_left,
181 const Bounds &bounds_right);
182
183 //*}
186
193 std::tuple<Number, bool> limit(const Bounds &bounds,
194 const state_type &U,
195 const state_type &P,
196 const Number t_min = Number(0.),
197 const Number t_max = Number(1.));
198
199 //*}
204
211 static bool
212 is_in_invariant_domain(const HyperbolicSystem & /*hyperbolic_system*/,
213 const Bounds & /*bounds*/,
214 const state_type & /*U*/);
215
216 private:
218
220
221 const HyperbolicSystem &hyperbolic_system;
222 const Parameters &parameters;
223 const PrecomputedVector &precomputed_values;
224
225 state_type U_i;
226
227 Bounds bounds_;
228
229 /* for relaxation */
230
231 Number h_relaxation_numerator;
232 Number kin_relaxation_numerator;
233 Number v2_relaxation_numerator;
234 Number relaxation_denominator;
235
237 };
238
239
240 /*
241 * -------------------------------------------------------------------------
242 * Inline definitions
243 * -------------------------------------------------------------------------
244 */
245
246
247 template <int dim, typename Number>
248 DEAL_II_ALWAYS_INLINE inline void
249 Limiter<dim, Number>::reset(unsigned int /*i*/,
250 const state_type &new_U_i,
251 const flux_contribution_type & /*new_flux_i*/)
252 {
253 U_i = new_U_i;
254
255 auto &[h_min, h_max, h_small, kin_max, v2_max] = bounds_;
256
257 h_min = Number(std::numeric_limits<ScalarNumber>::max());
258 h_max = Number(0.);
259 h_small = Number(0.);
260 kin_max = Number(0.);
261 v2_max = Number(0.);
262
263 h_relaxation_numerator = Number(0.);
264 kin_relaxation_numerator = Number(0.);
265 v2_relaxation_numerator = Number(0.);
266 relaxation_denominator = Number(0.);
267 }
268
269
270 template <int dim, typename Number>
271 DEAL_II_ALWAYS_INLINE inline void Limiter<dim, Number>::accumulate(
272 const state_type &U_j,
273 const state_type &U_star_ij,
274 const state_type &U_star_ji,
275 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
276 const state_type &affine_shift)
277 {
278 const auto view = hyperbolic_system.view<dim, Number>();
279
280 /* The bar states: */
281
282 const auto f_star_ij = view.f(U_star_ij);
283 const auto f_star_ji = view.f(U_star_ji);
284
285 /* bar state shifted by an affine shift: */
286 const auto U_ij_bar =
287 ScalarNumber(0.5) *
288 (U_star_ij + U_star_ji +
289 contract(add(f_star_ij, -f_star_ji), scaled_c_ij)) +
290 affine_shift;
291
292 /* Bounds: */
293
294 auto &[h_min, h_max, h_small, kin_max, v2_max] = bounds_;
295
296 const auto h_bar_ij = view.water_depth(U_ij_bar);
297 h_min = std::min(h_min, h_bar_ij);
298 h_max = std::max(h_max, h_bar_ij);
299
300 const auto kin_bar_ij = view.kinetic_energy(U_ij_bar);
301 kin_max = std::max(kin_max, kin_bar_ij);
302
303 const auto v_bar_ij = view.momentum(U_ij_bar) *
304 view.inverse_water_depth_mollified(U_ij_bar);
305 const auto v2_bar_ij = v_bar_ij.norm_square();
306 v2_max = std::max(v2_max, v2_bar_ij);
307
308 /* Relaxation: */
309
310 /* Use a uniform weight. */
311 const auto beta_ij = Number(1.);
312
313 relaxation_denominator += std::abs(beta_ij);
314
315 const auto h_i = view.water_depth(U_i);
316 const auto h_j = view.water_depth(U_j);
317 h_relaxation_numerator += beta_ij * (h_i + h_j);
318
319 const auto kin_i = view.kinetic_energy(U_i);
320 const auto kin_j = view.kinetic_energy(U_j);
321 kin_relaxation_numerator += beta_ij * (kin_i + kin_j);
322
323 const auto vel_i =
324 view.momentum(U_i) * view.inverse_water_depth_mollified(U_i);
325 const auto vel_j =
326 view.momentum(U_j) * view.inverse_water_depth_mollified(U_j);
327 v2_relaxation_numerator +=
328 beta_ij * (-vel_i.norm_square() + vel_j.norm_square());
329 }
330
331
332 template <int dim, typename Number>
333 DEAL_II_ALWAYS_INLINE inline auto
334 Limiter<dim, Number>::bounds(const Number hd_i) const -> Bounds
335 {
336 const auto view = hyperbolic_system.view<dim, Number>();
337
338 auto relaxed_bounds = bounds_;
339 auto &[h_min, h_max, h_small, kin_max, v2_max] = relaxed_bounds;
340
341 /* Use r_i = factor * (m_i / |Omega|) ^ (1.5 / d): */
342
343 Number r_i = std::sqrt(hd_i); // in 3D: ^ 3/6
344 if constexpr (dim == 2) //
345 r_i = dealii::Utilities::fixed_power<3>(std::sqrt(r_i)); // in 2D: ^ 3/4
346 else if constexpr (dim == 1) //
347 r_i = dealii::Utilities::fixed_power<3>(r_i); // in 1D: ^ 3/2
348 r_i *= parameters.relaxation_factor();
349
350 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
351
352 const Number h_relaxed = ScalarNumber(2.) *
353 std::abs(h_relaxation_numerator) /
354 (relaxation_denominator + Number(eps));
355
356 h_min = std::max((Number(1.) - r_i) * h_min, h_min - h_relaxed);
357 h_max = std::min((Number(1.) + r_i) * h_max, h_max + h_relaxed);
358
359 const Number kin_relaxed = ScalarNumber(2.) *
360 std::abs(kin_relaxation_numerator) /
361 (relaxation_denominator + Number(eps));
362
363 kin_max = std::min((Number(1.) + r_i) * kin_max, kin_max + kin_relaxed);
364
365 const Number v2_relaxed = ScalarNumber(2.) *
366 std::abs(v2_relaxation_numerator) /
367 (relaxation_denominator + Number(eps));
368
369 v2_max = std::min((Number(1.) + r_i) * v2_max, v2_max + v2_relaxed);
370
371 /* Use r_i = 0.2 * (m_i / |Omega|) ^ (1 / d): */
372
373 r_i = hd_i;
374 if constexpr (dim == 2)
375 r_i = std::sqrt(hd_i);
376 r_i *= view.dry_state_relaxation_factor();
377
378 h_small = view.reference_water_depth() * r_i;
379
380 return relaxed_bounds;
381 }
382
383
384 template <int dim, typename Number>
385 DEAL_II_ALWAYS_INLINE inline auto
387 const Bounds &bounds_r) -> Bounds
388 {
389 const auto &[h_min_l, h_max_l, h_small_l, k_max_l, v2_max_l] = bounds_l;
390 const auto &[h_min_r, h_max_r, h_small_r, k_max_r, v2_max_r] = bounds_r;
391
392 return {std::min(h_min_l, h_min_r),
393 std::max(h_max_l, h_max_r),
394 std::min(h_small_l, h_small_r),
395 std::max(k_max_l, h_max_r),
396 std::max(v2_max_l, v2_max_r)};
397 }
398
399
400 template <int dim, typename Number>
401 DEAL_II_ALWAYS_INLINE inline bool
403 const HyperbolicSystem & /*hyperbolic_system*/,
404 const Bounds & /*bounds*/,
405 const state_type & /*U*/)
406 {
407 AssertThrow(false, dealii::ExcNotImplemented());
408 __builtin_trap();
409 return true;
410 }
411
412
413 } // namespace ShallowWater
414} // namespace ryujin
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(limit_on_square_velocity)
ACCESSOR_READ_ONLY(newton_max_iterations)
LimiterParameters(const std::string &subsection="/Limiter")
Definition: limiter.h:24
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::ScalarNumber ScalarNumber
Definition: limiter.h:97
typename View::flux_contribution_type flux_contribution_type
Definition: limiter.h:103
std::array< Number, n_bounds > Bounds
Definition: limiter.h:139
typename View::PrecomputedVector PrecomputedVector
Definition: limiter.h:107
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)
Definition: limiter.h:271
Bounds bounds(const Number hd_i) const
Definition: limiter.h:334
Limiter(const HyperbolicSystem &hyperbolic_system, const Parameters &parameters, const PrecomputedVector &precomputed_values)
Definition: limiter.h:144
typename View::precomputed_type precomputed_type
Definition: limiter.h:105
void reset(const unsigned int i, const state_type &U_i, const flux_contribution_type &flux_i)
Definition: limiter.h:249
typename View::state_type state_type
Definition: limiter.h:101
static Bounds combine_bounds(const Bounds &bounds_left, const Bounds &bounds_right)
Definition: limiter.h:386
static bool is_in_invariant_domain(const HyperbolicSystem &, const Bounds &, const state_type &)
Definition: limiter.h:402
static constexpr auto problem_dimension
Definition: limiter.h:99
static constexpr unsigned int n_bounds
Definition: limiter.h:134
LimiterParameters< ScalarNumber > Parameters
Definition: limiter.h:109
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)