10#include <compile_time_options.h>
14#include <deal.II/base/parameter_acceptor.h>
15#include <deal.II/base/vectorization.h>
21 template <
typename ScalarNumber =
double>
26 : ParameterAcceptor(subsection)
28 evc_factor_ = ScalarNumber(1.);
29 add_parameter(
"evc factor",
31 "Factor for scaling the entropy viscocity commuator");
37 ScalarNumber evc_factor_;
80 template <
int dim,
typename Number =
double>
131 : hyperbolic_system(hyperbolic_system)
132 , parameters(parameters)
133 , precomputed_values(precomputed_values)
149 const dealii::Tensor<1, dim, Number> &c_ij);
154 Number
alpha(
const Number h_i)
const;
168 Number rho_i_inverse = 0.;
187 template <
int dim,
typename Number>
188 DEAL_II_ALWAYS_INLINE
inline void
193 const auto view = hyperbolic_system.view<dim, Number>();
195 const auto &[new_s_i, new_eta_i] =
196 precomputed_values.template get_tensor<Number, precomputed_type>(i);
198 const auto rho_i = view.density(U_i);
199 rho_i_inverse = Number(1.) / rho_i;
202 d_eta_i = view.harten_entropy_derivative(U_i);
203 d_eta_i[0] -= eta_i * rho_i_inverse;
211 template <
int dim,
typename Number>
213 const unsigned int *js,
215 const dealii::Tensor<1, dim, Number> &c_ij)
219 const auto view = hyperbolic_system.view<dim, Number>();
221 const auto &[s_j, eta_j] =
222 precomputed_values.template get_tensor<Number, precomputed_type>(js);
224 const auto rho_j = view.density(U_j);
225 const auto rho_j_inverse = Number(1.) / rho_j;
227 const auto m_j = view.momentum(U_j);
228 const auto f_j = view.f(U_j);
230 const auto entropy_flux =
231 (eta_j * rho_j_inverse - eta_i * rho_i_inverse) * (m_j * c_ij);
233 left += entropy_flux;
234 for (
unsigned int k = 0; k < problem_dimension; ++k) {
235 const auto component = (f_j[k] - f_i[k]) * c_ij;
236 right[k] += component;
241 template <
int dim,
typename Number>
242 DEAL_II_ALWAYS_INLINE
inline Number
247 Number numerator = left;
248 Number denominator = std::abs(left);
249 for (
unsigned int k = 0; k < problem_dimension; ++k) {
250 numerator -= d_eta_i[k] * right[k];
251 denominator += std::abs(d_eta_i[k] * right[k]);
254 const auto quotient =
255 std::abs(numerator) / (denominator + hd_i * std::abs(eta_i));
257 return std::min(Number(1.), parameters.evc_factor() * quotient);
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")
ACCESSOR_READ_ONLY(evc_factor)
void accumulate(const unsigned int *js, const state_type &U_j, const dealii::Tensor< 1, dim, Number > &c_ij)
typename View::state_type state_type
IndicatorParameters< ScalarNumber > Parameters
Number alpha(const Number h_i) const
typename View::flux_type flux_type
typename View::precomputed_type precomputed_type
typename View::ScalarNumber ScalarNumber
typename View::PrecomputedVector PrecomputedVector
Indicator(const HyperbolicSystem &hyperbolic_system, const Parameters ¶meters, const PrecomputedVector &precomputed_values)
void reset(const unsigned int i, const state_type &U_i)
static constexpr auto problem_dimension