10#include <compile_time_options.h>
14#include <deal.II/base/parameter_acceptor.h>
15#include <deal.II/base/vectorization.h>
22 template <
typename ScalarNumber =
double>
27 : ParameterAcceptor(subsection)
29 evc_factor_ = ScalarNumber(1.);
30 add_parameter(
"evc factor",
32 "Factor for scaling the entropy viscocity commuator");
38 ScalarNumber evc_factor_;
81 template <
int dim,
typename Number =
double>
157 : hyperbolic_system(hyperbolic_system)
158 , parameters(parameters)
159 , precomputed_values(precomputed_values)
175 const dealii::Tensor<1, dim, Number> &c_ij);
180 Number
alpha(
const Number h_i)
const;
197 Number rho_i_inverse = 0.;
218 template <
int dim,
typename Number>
219 DEAL_II_ALWAYS_INLINE
inline void
222 const auto view = hyperbolic_system.view<dim, Number>();
224 if (!view.compute_strict_bounds())
229 const auto &[p_i, gamma_min_i, s_i, new_eta_i] =
231 .template get_tensor<Number, precomputed_state_type>(i);
233 gamma_min = gamma_min_i;
235 const auto rho_i = view.density(U_i);
236 rho_i_inverse = Number(1.) / rho_i;
239 d_eta_i = view.surrogate_harten_entropy_derivative(U_i, eta_i, gamma_min);
240 d_eta_i[0] -= eta_i * rho_i_inverse;
242 const auto surrogate_p_i = view.surrogate_pressure(U_i, gamma_min);
243 f_i = view.f(U_i, surrogate_p_i);
250 template <
int dim,
typename Number>
252 const unsigned int * ,
254 const dealii::Tensor<1, dim, Number> &c_ij)
256 const auto view = hyperbolic_system.view<dim, Number>();
258 if (!view.compute_strict_bounds())
263 const auto eta_j = view.surrogate_harten_entropy(U_j, gamma_min);
265 const auto rho_j = view.density(U_j);
266 const auto rho_j_inverse = Number(1.) / rho_j;
268 const auto m_j = view.momentum(U_j);
270 const auto surrogate_p_j = view.surrogate_pressure(U_j, gamma_min);
271 const auto f_j = view.f(U_j, surrogate_p_j);
273 left += (eta_j * rho_j_inverse - eta_i * rho_i_inverse) * (m_j * c_ij);
274 for (
unsigned int k = 0; k < problem_dimension; ++k)
275 right[k] += (f_j[k] - f_i[k]) * c_ij;
279 template <
int dim,
typename Number>
280 DEAL_II_ALWAYS_INLINE
inline Number
283 const auto view = hyperbolic_system.view<dim, Number>();
285 if (!view.compute_strict_bounds())
288 Number numerator = left;
289 Number denominator = std::abs(left);
290 for (
unsigned int k = 0; k < problem_dimension; ++k) {
291 numerator -= d_eta_i[k] * right[k];
292 denominator += std::abs(d_eta_i[k] * right[k]);
295 const auto quotient =
296 std::abs(numerator) / (denominator + hd_i * std::abs(eta_i));
298 return std::min(Number(1.), parameters.evc_factor() * quotient);
dealii::Tensor< 1, problem_dimension, Number > state_type
std::array< Number, n_precomputed_values > precomputed_state_type
typename get_value_type< Number >::type ScalarNumber
static constexpr unsigned int n_precomputed_values
static constexpr unsigned int problem_dimension
flux_type flux_contribution_type
dealii::Tensor< 1, problem_dimension, dealii::Tensor< 1, dim, Number > > flux_type
ACCESSOR_READ_ONLY(evc_factor)
IndicatorParameters(const std::string &subsection="/Indicator")
typename View::flux_type flux_type
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
typename View::flux_contribution_type flux_contribution_type
Number alpha(const Number h_i) const
void reset(const unsigned int i, const state_type &U_i)
IndicatorParameters< ScalarNumber > Parameters
Indicator(const HyperbolicSystem &hyperbolic_system, const Parameters ¶meters, const MultiComponentVector< ScalarNumber, n_precomputed_values > &precomputed_values)
static constexpr unsigned int n_precomputed_values
static constexpr unsigned int problem_dimension
typename View::precomputed_state_type precomputed_state_type
typename View::ScalarNumber ScalarNumber