ryujin 2.1.1 revision ef0fcd4010d109b860652ece4a7b8963fb7d46b1
indicator.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2024 by the ryujin authors
4//
5
6#pragma once
7
8#include <compile_time_options.h>
9
10#include "hyperbolic_system.h"
11
13#include <simd.h>
14
15#include <deal.II/base/parameter_acceptor.h>
16#include <deal.II/base/vectorization.h>
17
18namespace ryujin
19{
20 namespace Euler
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
81 template <int dim, typename Number = double>
83 {
84 public:
89
91
93
95
96 using state_type = typename View::state_type;
97
98 using flux_type = typename View::flux_type;
99
101
103
105
107
125
129 Indicator(const HyperbolicSystem &hyperbolic_system,
130 const Parameters &parameters,
131 const PrecomputedVector &precomputed_values)
132 : hyperbolic_system(hyperbolic_system)
133 , parameters(parameters)
134 , precomputed_values(precomputed_values)
135 {
136 }
137
142 void reset(const unsigned int i, const state_type &U_i);
143
148 void accumulate(const unsigned int *js,
149 const state_type &U_j,
150 const dealii::Tensor<1, dim, Number> &c_ij);
151
155 Number alpha(const Number h_i) const;
156
158
159 private:
164
165 const HyperbolicSystem &hyperbolic_system;
166 const Parameters &parameters;
167 const PrecomputedVector &precomputed_values;
168
169 Number rho_i_inverse = 0.;
170 Number eta_i = 0.;
171 flux_type f_i;
172 state_type d_eta_i;
173
174 Number left = 0.;
175 state_type right;
176
178 };
179
180
181 /*
182 * -------------------------------------------------------------------------
183 * Inline definitions
184 * -------------------------------------------------------------------------
185 */
186
187
188 template <int dim, typename Number>
189 DEAL_II_ALWAYS_INLINE inline void
190 Indicator<dim, Number>::reset(const unsigned int i, const state_type &U_i)
191 {
192 /* Entropy viscosity commutator: */
193
194 const auto view = hyperbolic_system.view<dim, Number>();
195
196 const auto &[new_s_i, new_eta_i] =
197 precomputed_values.template get_tensor<Number, precomputed_type>(i);
198
199 const auto rho_i = view.density(U_i);
200 rho_i_inverse = Number(1.) / rho_i;
201 eta_i = new_eta_i;
202
203 d_eta_i = view.harten_entropy_derivative(U_i);
204 d_eta_i[0] -= eta_i * rho_i_inverse;
205 f_i = view.f(U_i);
206
207 left = 0.;
208 right = 0.;
209 }
210
211
212 template <int dim, typename Number>
213 DEAL_II_ALWAYS_INLINE inline void Indicator<dim, Number>::accumulate(
214 const unsigned int *js,
215 const state_type &U_j,
216 const dealii::Tensor<1, dim, Number> &c_ij)
217 {
218 /* Entropy viscosity commutator: */
219
220 const auto view = hyperbolic_system.view<dim, Number>();
221
222 const auto &[s_j, eta_j] =
223 precomputed_values.template get_tensor<Number, precomputed_type>(js);
224
225 const auto rho_j = view.density(U_j);
226 const auto rho_j_inverse = Number(1.) / rho_j;
227
228 const auto m_j = view.momentum(U_j);
229 const auto f_j = view.f(U_j);
230
231 const auto entropy_flux =
232 (eta_j * rho_j_inverse - eta_i * rho_i_inverse) * (m_j * c_ij);
233
234 left += entropy_flux;
235 for (unsigned int k = 0; k < problem_dimension; ++k) {
236 const auto component = (f_j[k] - f_i[k]) * c_ij;
237 right[k] += component;
238 }
239 }
240
241
242 template <int dim, typename Number>
243 DEAL_II_ALWAYS_INLINE inline Number
244 Indicator<dim, Number>::alpha(const Number hd_i) const
245 {
246 /* Entropy viscosity commutator: */
247
248 Number numerator = left;
249 Number denominator = std::abs(left);
250 for (unsigned int k = 0; k < problem_dimension; ++k) {
251 numerator -= d_eta_i[k] * right[k];
252 denominator += std::abs(d_eta_i[k] * right[k]);
253 }
254
255 const auto quotient =
256 std::abs(numerator) / (denominator + hd_i * std::abs(eta_i));
257
258 return std::min(Number(1.), parameters.evc_factor() * quotient);
259 }
260 } // namespace Euler
261} // namespace ryujin
typename get_value_type< Number >::type ScalarNumber
dealii::Tensor< 1, problem_dimension, dealii::Tensor< 1, dim, Number > > flux_type
Vectors::MultiComponentVector< ScalarNumber, n_precomputed_values > PrecomputedVector
static constexpr unsigned int problem_dimension
dealii::Tensor< 1, problem_dimension, Number > state_type
std::array< Number, n_precomputed_values > precomputed_type
IndicatorParameters(const std::string &subsection="/Indicator")
Definition: indicator.h:26
void accumulate(const unsigned int *js, const state_type &U_j, const dealii::Tensor< 1, dim, Number > &c_ij)
Definition: indicator.h:213
typename View::state_type state_type
Definition: indicator.h:96
IndicatorParameters< ScalarNumber > Parameters
Definition: indicator.h:104
Number alpha(const Number h_i) const
Definition: indicator.h:244
typename View::flux_type flux_type
Definition: indicator.h:98
typename View::precomputed_type precomputed_type
Definition: indicator.h:100
typename View::ScalarNumber ScalarNumber
Definition: indicator.h:92
typename View::PrecomputedVector PrecomputedVector
Definition: indicator.h:102
Indicator(const HyperbolicSystem &hyperbolic_system, const Parameters &parameters, const PrecomputedVector &precomputed_values)
Definition: indicator.h:129
void reset(const unsigned int i, const state_type &U_i)
Definition: indicator.h:190
static constexpr auto problem_dimension
Definition: indicator.h:94