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) 2023 - 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
18
19namespace ryujin
20{
21 namespace ScalarConservation
22 {
23 template <typename ScalarNumber = double>
24 class IndicatorParameters : public dealii::ParameterAcceptor
25 {
26 public:
27 IndicatorParameters(const std::string &subsection = "/Indicator")
28 : ParameterAcceptor(subsection)
29 {
30 evc_factor_ = ScalarNumber(1.);
31 add_parameter("evc factor",
32 evc_factor_,
33 "Factor for scaling the entropy viscocity commuator");
34 }
35
36 ACCESSOR_READ_ONLY(evc_factor);
37
38 private:
39 ScalarNumber evc_factor_;
40 };
41
42
49 template <int dim, typename Number = double>
51 {
52 public:
57
59
61
63
64 using state_type = typename View::state_type;
65
66 using flux_type = typename View::flux_type;
67
69
71
73
75
93
97 Indicator(const HyperbolicSystem &hyperbolic_system,
98 const Parameters &parameters,
99 const PrecomputedVector &precomputed_values)
100 : hyperbolic_system(hyperbolic_system)
101 , parameters(parameters)
102 , precomputed_values(precomputed_values)
103 {
104 }
105
110 void reset(const unsigned int i, const state_type &U_i);
111
116 void accumulate(const unsigned int *js,
117 const state_type &U_j,
118 const dealii::Tensor<1, dim, Number> &c_ij);
119
123 Number alpha(const Number h_i) const;
124
126
127 private:
132
133 const HyperbolicSystem &hyperbolic_system;
134 const Parameters &parameters;
135 const PrecomputedVector &precomputed_values;
136
137 Number u_i;
138 Number u_abs_max;
139 dealii::Tensor<1, dim, Number> f_i;
140 Number left;
141 Number right;
143 };
144
145
146 /*
147 * -------------------------------------------------------------------------
148 * Inline definitions
149 * -------------------------------------------------------------------------
150 */
151
152
153 template <int dim, typename Number>
154 DEAL_II_ALWAYS_INLINE inline void
155 Indicator<dim, Number>::reset(const unsigned int i, const state_type &U_i)
156 {
157 /* entropy viscosity commutator: */
158
159 const auto view = hyperbolic_system.view<dim, Number>();
160
161 const auto prec_i =
162 precomputed_values.template get_tensor<Number, precomputed_type>(i);
163
164 u_i = view.state(U_i);
165 u_abs_max = std::abs(u_i);
166 f_i = view.construct_flux_tensor(prec_i);
167 left = 0.;
168 right = 0.;
169 }
170
171
172 template <int dim, typename Number>
173 DEAL_II_ALWAYS_INLINE inline void Indicator<dim, Number>::accumulate(
174 const unsigned int *js,
175 const state_type &U_j,
176 const dealii::Tensor<1, dim, Number> &c_ij)
177 {
178 /* entropy viscosity commutator: */
179
180 const auto view = hyperbolic_system.view<dim, Number>();
181
182 const auto prec_j =
183 precomputed_values.template get_tensor<Number, precomputed_type>(js);
184
185 const auto u_j = view.state(U_j);
186 u_abs_max = std::max(u_abs_max, std::abs(u_j));
187 const auto d_eta_j = view.kruzkov_entropy_derivative(u_i, u_j);
188 const auto f_j = view.construct_flux_tensor(prec_j);
189
190 left += d_eta_j * (f_j * c_ij);
191 right += d_eta_j * (f_i * c_ij);
192 }
193
194
195 template <int dim, typename Number>
196 DEAL_II_ALWAYS_INLINE inline Number
197 Indicator<dim, Number>::alpha(const Number hd_i) const
198 {
199 Number numerator = left - right;
200 Number denominator = std::abs(left) + std::abs(right);
201
202 const auto regularization =
203 Number(100. * std::numeric_limits<ScalarNumber>::min());
204
205 const auto quotient =
206 std::abs(numerator) /
207 (denominator + std::max(hd_i * std::abs(u_abs_max), regularization));
208
209 return std::min(Number(1.), parameters.evc_factor() * quotient);
210 }
211
212 } // namespace ScalarConservation
213} // namespace ryujin
dealii::Tensor< 1, problem_dimension, dealii::Tensor< 1, dim, Number > > flux_type
std::array< Number, n_precomputed_values > precomputed_type
Vectors::MultiComponentVector< ScalarNumber, n_precomputed_values > PrecomputedVector
typename get_value_type< Number >::type ScalarNumber
dealii::Tensor< 1, problem_dimension, Number > state_type
IndicatorParameters(const std::string &subsection="/Indicator")
Definition: indicator.h:27
void reset(const unsigned int i, const state_type &U_i)
Definition: indicator.h:155
typename View::state_type state_type
Definition: indicator.h:64
typename View::PrecomputedVector PrecomputedVector
Definition: indicator.h:70
typename View::ScalarNumber ScalarNumber
Definition: indicator.h:60
void accumulate(const unsigned int *js, const state_type &U_j, const dealii::Tensor< 1, dim, Number > &c_ij)
Definition: indicator.h:173
static constexpr auto problem_dimension
Definition: indicator.h:62
typename View::flux_type flux_type
Definition: indicator.h:66
Indicator(const HyperbolicSystem &hyperbolic_system, const Parameters &parameters, const PrecomputedVector &precomputed_values)
Definition: indicator.h:97
Number alpha(const Number h_i) const
Definition: indicator.h:197
typename View::precomputed_type precomputed_type
Definition: indicator.h:68