ryujin 2.1.1 revision ef0fcd4010d109b860652ece4a7b8963fb7d46b1
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 <compile_time_options.h>
11
12#include "hyperbolic_system.h"
13
15#include <newton.h>
16
17namespace ryujin
18{
19 namespace ShallowWater
20 {
21 template <typename ScalarNumber = double>
22 class LimiterParameters : public dealii::ParameterAcceptor
23 {
24 public:
25 LimiterParameters(const std::string &subsection = "/Limiter")
26 : ParameterAcceptor(subsection)
27 {
28 iterations_ = 2;
29 add_parameter(
30 "iterations", iterations_, "Number of limiter iterations");
31
32 if constexpr (std::is_same_v<ScalarNumber, double>)
33 newton_tolerance_ = 1.e-10;
34 else
35 newton_tolerance_ = 1.e-4;
36 add_parameter("newton tolerance",
37 newton_tolerance_,
38 "Tolerance for the quadratic newton stopping criterion");
39
40 newton_max_iterations_ = 2;
41 add_parameter("newton max iterations",
42 newton_max_iterations_,
43 "Maximal number of quadratic newton iterations performed "
44 "during limiting");
45
46 relaxation_factor_ = ScalarNumber(1.);
47 add_parameter("relaxation factor",
48 relaxation_factor_,
49 "Factor for scaling the relaxation window with r_i = "
50 "factor * (m_i/|Omega|)^(1.5/d).");
51 }
52
53 ACCESSOR_READ_ONLY(iterations);
54 ACCESSOR_READ_ONLY(newton_tolerance);
55 ACCESSOR_READ_ONLY(newton_max_iterations);
56 ACCESSOR_READ_ONLY(relaxation_factor);
57
58 private:
59 unsigned int iterations_;
60 ScalarNumber newton_tolerance_;
61 unsigned int newton_max_iterations_;
62 ScalarNumber relaxation_factor_;
63 };
64
65
71 template <int dim, typename Number = double>
72 class Limiter
73 {
74 public:
79
81
83
85
86 using state_type = typename View::state_type;
87
89
91
93
95
97
104 static constexpr unsigned int n_bounds = 3;
105
109 using Bounds = std::array<Number, n_bounds>;
110
114 Limiter(const HyperbolicSystem &hyperbolic_system,
115 const Parameters &parameters,
116 const PrecomputedVector &precomputed_values)
117 : hyperbolic_system(hyperbolic_system)
118 , parameters(parameters)
119 , precomputed_values(precomputed_values)
120 {
121 }
122
127 Bounds projection_bounds_from_state(const unsigned int i,
128 const state_type &U_i) const;
129
135 Bounds combine_bounds(const Bounds &bounds_left,
136 const Bounds &bounds_right) const;
137
145 Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const;
146
148
166
170 void reset(const unsigned int i,
171 const state_type &U_i,
172 const flux_contribution_type &flux_i);
173
178 void accumulate(const state_type &U_j,
179 const state_type &U_star_ij,
180 const state_type &U_star_ji,
181 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
182 const state_type &affine_shift);
183
187 Bounds bounds(const Number hd_i) const;
188
189 //*}
192
199 std::tuple<Number, bool> limit(const Bounds &bounds,
200 const state_type &U,
201 const state_type &P,
202 const Number t_min = Number(0.),
203 const Number t_max = Number(1.)) const;
204
205 private:
207
209
210 const HyperbolicSystem &hyperbolic_system;
211 const Parameters &parameters;
212 const PrecomputedVector &precomputed_values;
213
214 state_type U_i;
215
216 Bounds bounds_;
217
218 /* for relaxation */
219
220 Number h_relaxation_numerator;
221 Number v2_relaxation_numerator;
222 Number relaxation_denominator;
223
225 };
226
227
228 /*
229 * -------------------------------------------------------------------------
230 * Inline definitions
231 * -------------------------------------------------------------------------
232 */
233
234
235 template <int dim, typename Number>
236 DEAL_II_ALWAYS_INLINE inline auto
238 const unsigned int /*i*/, const state_type &U_i) const -> Bounds
239 {
240 const auto view = hyperbolic_system.view<dim, Number>();
241 const auto h_i = view.water_depth(U_i);
242 const auto v_i =
243 view.momentum(U_i) * view.inverse_water_depth_mollified(U_i);
244 const auto v2_i = v_i.norm_square();
245
246 return {/*h_min*/ h_i, /*h_max*/ h_i, /*v2_max*/ v2_i};
247 }
248
249
250 template <int dim, typename Number>
251 DEAL_II_ALWAYS_INLINE inline auto
253 const Bounds &bounds_r) const -> Bounds
254 {
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;
257
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)};
261 }
262
263
264 template <int dim, typename Number>
265 DEAL_II_ALWAYS_INLINE inline auto
267 const Number &hd) const -> Bounds
268 {
269 auto relaxed_bounds = bounds;
270 auto &[h_min, h_max, v2_max] = relaxed_bounds;
271
272 /* Use r = factor * (m_i / |Omega|) ^ (1.5 / d): */
273
274 Number r = std::sqrt(hd); // in 3D: ^ 3/6
275 if constexpr (dim == 2) //
276 r = dealii::Utilities::fixed_power<3>(std::sqrt(r)); // in 2D: ^ 3/4
277 else if constexpr (dim == 1) //
278 r = dealii::Utilities::fixed_power<3>(r); // in 1D: ^ 3/2
279 r *= parameters.relaxation_factor();
280
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);
285
286 return relaxed_bounds;
287 }
288
289
290 template <int dim, typename Number>
291 DEAL_II_ALWAYS_INLINE inline void
292 Limiter<dim, Number>::reset(unsigned int /*i*/,
293 const state_type &new_U_i,
294 const flux_contribution_type & /*new_flux_i*/)
295 {
296 U_i = new_U_i;
297
298 auto &[h_min, h_max, v2_max] = bounds_;
299
300 h_min = Number(std::numeric_limits<ScalarNumber>::max());
301 h_max = Number(0.);
302 v2_max = Number(0.);
303
304 h_relaxation_numerator = Number(0.);
305 v2_relaxation_numerator = Number(0.);
306 relaxation_denominator = Number(0.);
307 }
308
309
310 template <int dim, typename Number>
311 DEAL_II_ALWAYS_INLINE inline void Limiter<dim, Number>::accumulate(
312 const state_type &U_j,
313 const state_type &U_star_ij,
314 const state_type &U_star_ji,
315 const dealii::Tensor<1, dim, Number> &scaled_c_ij,
316 const state_type &affine_shift)
317 {
318 const auto view = hyperbolic_system.view<dim, Number>();
319
320 /* The bar states: */
321
322 const auto f_star_ij = view.f(U_star_ij);
323 const auto f_star_ji = view.f(U_star_ji);
324
325 /* bar state shifted by an affine shift: */
326 const auto U_ij_bar =
327 ScalarNumber(0.5) *
328 (U_star_ij + U_star_ji +
329 contract(add(f_star_ij, -f_star_ji), scaled_c_ij)) +
330 affine_shift;
331
332 /* Bounds: */
333
334 auto &[h_min, h_max, v2_max] = bounds_;
335
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);
339
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);
344
345 /* Relaxation: */
346
347 /* Use a uniform weight. */
348 const auto beta_ij = Number(1.);
349
350 relaxation_denominator += std::abs(beta_ij);
351
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);
355
356 const auto vel_i =
357 view.momentum(U_i) * view.inverse_water_depth_mollified(U_i);
358 const auto vel_j =
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());
362 }
363
364
365 template <int dim, typename Number>
366 DEAL_II_ALWAYS_INLINE inline auto
367 Limiter<dim, Number>::bounds(const Number hd_i) const -> Bounds
368 {
369 const auto &[h_min, h_max, v2_max] = bounds_;
370
371 auto relaxed_bounds = fully_relax_bounds(bounds_, hd_i);
372 auto &[h_min_relaxed, h_max_relaxed, v2_max_relaxed] = relaxed_bounds;
373
374 /* Apply a stricter window: */
375
376 constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
377
378 const Number h_relaxed =
379 ScalarNumber(2. * parameters.relaxation_factor()) *
380 std::abs(h_relaxation_numerator) /
381 (relaxation_denominator + Number(eps));
382
383 const Number v2_relaxed =
384 ScalarNumber(2. * parameters.relaxation_factor()) *
385 std::abs(v2_relaxation_numerator) /
386 (relaxation_denominator + Number(eps));
387
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);
391
392 return relaxed_bounds;
393 }
394 } // namespace ShallowWater
395} // 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(newton_max_iterations)
LimiterParameters(const std::string &subsection="/Limiter")
Definition: limiter.h:25
Bounds projection_bounds_from_state(const unsigned int i, const state_type &U_i) const
Definition: limiter.h:237
typename View::ScalarNumber ScalarNumber
Definition: limiter.h:82
typename View::flux_contribution_type flux_contribution_type
Definition: limiter.h:88
Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const
Definition: limiter.h:266
std::array< Number, n_bounds > Bounds
Definition: limiter.h:109
typename View::PrecomputedVector PrecomputedVector
Definition: limiter.h:92
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:311
Bounds bounds(const Number hd_i) const
Definition: limiter.h:367
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 &parameters, const PrecomputedVector &precomputed_values)
Definition: limiter.h:114
typename View::precomputed_type precomputed_type
Definition: limiter.h:90
Bounds combine_bounds(const Bounds &bounds_left, const Bounds &bounds_right) const
Definition: limiter.h:252
void reset(const unsigned int i, const state_type &U_i, const flux_contribution_type &flux_i)
Definition: limiter.h:292
typename View::state_type state_type
Definition: limiter.h:86
static constexpr auto problem_dimension
Definition: limiter.h:84
static constexpr unsigned int n_bounds
Definition: limiter.h:104
LimiterParameters< ScalarNumber > Parameters
Definition: limiter.h:94
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)