ryujin 2.1.1 revision c38afcfdec63b34961924f05408760247496a1f0
indicator.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
13
14#include <deal.II/base/parameter_acceptor.h>
15#include <deal.II/base/vectorization.h>
16
17
18namespace ryujin
19{
20 namespace ShallowWater
21 {
22 template <typename ScalarNumber = double>
23 class IndicatorParameters : public dealii::ParameterAcceptor
24 {
25 public:
26 IndicatorParameters(const std::string &subsection = "/Indicator")
27 : ParameterAcceptor(subsection)
28 {
29 evc_factor_ = ScalarNumber(1.);
30 add_parameter("evc factor",
31 evc_factor_,
32 "Factor for scaling the entropy viscocity commuator");
33 }
34
35 ACCESSOR_READ_ONLY(evc_factor);
36
37 private:
38 ScalarNumber evc_factor_;
39 };
40
41
48 template <int dim, typename Number = double>
50 {
51 public:
56
58
60
62
63 using state_type = typename View::state_type;
64
65 using flux_type = typename View::flux_type;
66
68
70
72
74
92
96 Indicator(const HyperbolicSystem &hyperbolic_system,
97 const Parameters &parameters,
98 const PrecomputedVector &precomputed_values)
99 : hyperbolic_system(hyperbolic_system)
100 , parameters(parameters)
101 , precomputed_values(precomputed_values)
102 {
103 }
104
109 void reset(const unsigned int /*i*/, const state_type &U_i);
110
115 void accumulate(const unsigned int *js,
116 const state_type &U_j,
117 const dealii::Tensor<1, dim, Number> &c_ij);
118
122 Number alpha(const Number h_i);
123
125
126 private:
131
132 const HyperbolicSystem &hyperbolic_system;
133 const Parameters &parameters;
134 const PrecomputedVector &precomputed_values;
135
136 Number h_i = 0.;
137 Number eta_i = 0.;
138 flux_type f_i;
139 state_type d_eta_i;
140 Number pressure_i = 0.;
141
142 Number left = 0.;
143 state_type right;
145 };
146
147
148 /*
149 * -------------------------------------------------------------------------
150 * Inline definitions
151 * -------------------------------------------------------------------------
152 */
153
154
155 template <int dim, typename Number>
156 DEAL_II_ALWAYS_INLINE inline void
157 Indicator<dim, Number>::reset(const unsigned int i, const state_type &U_i)
158 {
159 /* entropy viscosity commutator: */
160
161 const auto view = hyperbolic_system.view<dim, Number>();
162
163 const auto &[eta_m, h_star] =
164 precomputed_values.template get_tensor<Number, precomputed_type>(i);
165
166 h_i = view.water_depth(U_i);
167 eta_i = eta_m;
168 d_eta_i = view.mathematical_entropy_derivative(U_i);
169 f_i = view.f(U_i);
170 pressure_i = view.pressure(U_i);
171
172 left = 0.;
173 right = 0.;
174 }
175
176
177 template <int dim, typename Number>
178 DEAL_II_ALWAYS_INLINE inline void Indicator<dim, Number>::accumulate(
179 const unsigned int *js,
180 const state_type &U_j,
181 const dealii::Tensor<1, dim, Number> &c_ij)
182 {
183 /* entropy viscosity commutator: */
184
185 const auto view = hyperbolic_system.view<dim, Number>();
186
187 const auto &[eta_j, h_star_j] =
188 precomputed_values.template get_tensor<Number, precomputed_type>(js);
189
190 const auto velocity_j =
191 view.momentum(U_j) * view.inverse_water_depth_sharp(U_j);
192 const auto f_j = view.f(U_j);
193 const auto pressure_j = view.pressure(U_j);
194
195 left += (eta_j + pressure_j) * (velocity_j * c_ij);
196
197 for (unsigned int k = 0; k < problem_dimension; ++k)
198 right[k] += (f_j[k] - f_i[k]) * c_ij;
199 }
200
201
202 template <int dim, typename Number>
203 DEAL_II_ALWAYS_INLINE inline Number
205 {
206 Number my_sum = 0.;
207 for (unsigned int k = 0; k < problem_dimension; ++k) {
208 my_sum += d_eta_i[k] * right[k];
209 }
210
211 Number numerator = std::abs(left - my_sum);
212 Number denominator = std::abs(left) + std::abs(my_sum);
213
214 const auto regularization =
215 Number(100. * std::numeric_limits<ScalarNumber>::min());
216
217 const auto quotient =
218 std::abs(numerator) /
219 (denominator + std::max(hd_i * std::abs(eta_i), regularization));
220
221 return std::min(Number(1.), parameters.evc_factor() * quotient);
222 }
223
224
225 } // namespace ShallowWater
226} // namespace ryujin
dealii::Tensor< 1, problem_dimension, dealii::Tensor< 1, dim, Number > > flux_type
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
IndicatorParameters(const std::string &subsection="/Indicator")
Definition: indicator.h:26
typename View::PrecomputedVector PrecomputedVector
Definition: indicator.h:69
typename View::state_type state_type
Definition: indicator.h:63
typename View::flux_type flux_type
Definition: indicator.h:65
typename View::precomputed_type precomputed_type
Definition: indicator.h:67
void reset(const unsigned int, const state_type &U_i)
Definition: indicator.h:157
void accumulate(const unsigned int *js, const state_type &U_j, const dealii::Tensor< 1, dim, Number > &c_ij)
Definition: indicator.h:178
static constexpr auto problem_dimension
Definition: indicator.h:61
typename View::ScalarNumber ScalarNumber
Definition: indicator.h:59
Indicator(const HyperbolicSystem &hyperbolic_system, const Parameters &parameters, const PrecomputedVector &precomputed_values)
Definition: indicator.h:96
Number alpha(const Number h_i)
Definition: indicator.h:204
IndicatorParameters< ScalarNumber > Parameters
Definition: indicator.h:71