11#include <deal.II/base/config.h>
18#include <deal.II/base/parameter_acceptor.h>
19#include <deal.II/base/tensor.h>
25 namespace ScalarConservation
27 template <
int dim,
typename Number>
28 class HyperbolicSystemView;
42 static inline std::string
problem_name =
"Scalar conservation equation";
55 template <
int dim,
typename Number>
70 std::shared_ptr<Flux> selected_flux_;
72 template <
int dim,
typename Number>
83 template <
int dim,
typename Number>
92 : hyperbolic_system_(hyperbolic_system)
99 template <
int dim2,
typename Number2>
115 DEAL_II_ALWAYS_INLINE
inline const std::string &
flux()
const
117 return hyperbolic_system_.flux_;
123 const auto &
flux = hyperbolic_system_.selected_flux_;
136 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
142 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
170 using state_type = dealii::Tensor<1, problem_dimension, Number>;
182 std::array<std::string, problem_dimension>{
"u"};
189 std::array<std::string, problem_dimension>{
"u"};
195 dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
217 std::array<Number, n_precomputed_initial_values>;
230 std::array<std::string, n_precomputed_initial_values>{};
252 []() -> std::array<std::string, n_precomputed_values> {
253 if constexpr (dim == 1)
255 else if constexpr (dim == 2)
256 return {
"f_1",
"f_2",
"df_1",
"df_2"};
257 else if constexpr (dim == 3)
258 return {
"f_1",
"f_2",
"f_3",
"df_1",
"df_2",
"df_3"};
271 template <
typename DISPATCH,
typename SPARSITY>
273 const DISPATCH &dispatch_check,
278 unsigned int )
const;
347 template <
typename Lambda>
349 const dealii::types::boundary_id ,
351 const dealii::Tensor<1, dim, Number> & ,
352 const Lambda & )
const;
397 const unsigned int i,
403 const unsigned int *js,
413 const dealii::Tensor<1, dim, Number> &c_ij)
const;
421 const dealii::Tensor<1, dim, Number> &c_ij)
const =
delete;
433 const unsigned int i,
438 const unsigned int *js,
458 template <
typename ST>
470 return primitive_state;
487 template <
typename Lambda>
489 const Lambda & )
const
505 : ParameterAcceptor(subsection)
508 add_parameter(
"flux",
510 "The scalar flux. Valid names are given by any of the "
511 "subsections defined below");
519 const auto populate_functions = [
this]() {
520 bool initialized =
false;
521 for (
auto &it : flux_list_)
524 if (it->name() == flux_) {
526 it->parse_parameters_call_back();
527 problem_name =
"Scalar conservation equation (" + it->name() +
528 ": " + it->flux_formula() +
")";
533 AssertThrow(initialized,
535 "Could not find a flux description with name \"" +
539 ParameterAcceptor::parse_parameters_call_back.connect(populate_functions);
540 populate_functions();
544 template <
int dim,
typename Number>
545 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
548 const auto &flux = hyperbolic_system_.selected_flux_;
549 dealii::Tensor<1, dim, Number> result;
552 for (
unsigned int k = 0; k < dim; ++k) {
553 if constexpr (std::is_same_v<ScalarNumber, Number>) {
554 result[k] = flux->value(u, k);
556 for (
unsigned int s = 0; s < Number::size(); ++s) {
557 result[k][s] = flux->value(u[s], k);
566 template <
int dim,
typename Number>
567 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
569 const Number &u)
const
571 const auto &flux = hyperbolic_system_.selected_flux_;
572 dealii::Tensor<1, dim, Number> result;
575 for (
unsigned int k = 0; k < dim; ++k) {
576 if constexpr (std::is_same_v<ScalarNumber, Number>) {
577 result[k] = flux->gradient(u, k);
579 for (
unsigned int s = 0; s < Number::size(); ++s) {
580 result[k][s] = flux->gradient(u[s], k);
589 template <
int dim,
typename Number>
590 template <
typename DISPATCH,
typename SPARSITY>
591 DEAL_II_ALWAYS_INLINE
inline void
593 unsigned int cycle [[maybe_unused]],
594 const DISPATCH &dispatch_check,
596 const SPARSITY &sparsity_simd,
599 unsigned int right)
const
601 Assert(cycle == 0, dealii::ExcInternalError());
605 unsigned int stride_size = get_stride_size<Number>;
608 for (
unsigned int i = left; i < right; i += stride_size) {
611 const unsigned int row_length = sparsity_simd.row_length(i);
617 const auto U_i = U.template get_tensor<Number>(i);
618 const auto u_i = state(U_i);
620 const auto f_i = flux_function(u_i);
621 const auto df_i = flux_gradient_function(u_i);
625 for (
unsigned int k = 0; k < n_precomputed_values / 2; ++k) {
627 prec_i[dim + k] = df_i[k];
630 precomputed_values.template write_tensor<Number>(prec_i, i);
635 template <
int dim,
typename Number>
636 DEAL_II_ALWAYS_INLINE
inline Number
643 template <
int dim,
typename Number>
644 DEAL_II_ALWAYS_INLINE
inline Number
651 template <
int dim,
typename Number>
652 DEAL_II_ALWAYS_INLINE
inline Number
654 const Number &u)
const
660 template <
int dim,
typename Number>
661 DEAL_II_ALWAYS_INLINE
inline Number
663 const Number &u)
const
665 return std::abs(k - u);
669 template <
int dim,
typename Number>
670 DEAL_II_ALWAYS_INLINE
inline Number
672 const Number &k,
const Number &u)
const
674 constexpr auto gte = dealii::SIMDComparison::greater_than_or_equal;
676 return dealii::compare_and_apply_mask<gte>(u, k, Number(1.), Number(-1.));
680 template <
int dim,
typename Number>
681 template <
typename Lambda>
682 DEAL_II_ALWAYS_INLINE
inline auto
684 dealii::types::boundary_id
id,
686 const dealii::Tensor<1, dim, Number> & ,
687 const Lambda &get_dirichlet_data)
const ->
state_type
692 result = get_dirichlet_data();
697 dealii::ExcMessage(
"Invalid boundary ID »Boundary::slip«, slip "
698 "boundary conditions are unavailable for scalar "
699 "conservation equations."));
705 dealii::ExcMessage(
"Invalid boundary ID »Boundary::no_slip«, "
706 "no-slip boundary conditions are unavailable "
707 "for scalar conservation equations."));
713 dealii::ExcMessage(
"Invalid boundary ID »Boundary::dynamic«, "
714 "dynamic boundary conditions are unavailable "
715 "for scalar conservation equations."));
723 template <
int dim,
typename Number>
724 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
728 dealii::Tensor<1, dim, Number> result;
730 if constexpr (dim == 1) {
731 const auto &[f, df] = precomputed_state;
734 }
else if constexpr (dim == 2) {
735 const auto &[f_1, f_2, df_1, df_2] = precomputed_state;
739 }
else if constexpr (dim == 3) {
740 const auto &[f_1, f_2, f_3, df_1, df_2, df_3] = precomputed_state;
750 template <
int dim,
typename Number>
751 DEAL_II_ALWAYS_INLINE
inline dealii::Tensor<1, dim, Number>
755 dealii::Tensor<1, dim, Number> result;
757 if constexpr (dim == 1) {
758 const auto &[f, df] = precomputed_state;
761 }
else if constexpr (dim == 2) {
762 const auto &[f_1, f_2, df_1, df_2] = precomputed_state;
766 }
else if constexpr (dim == 3) {
767 const auto &[f_1, f_2, f_3, df_1, df_2, df_3] = precomputed_state;
777 template <
int dim,
typename Number>
778 DEAL_II_ALWAYS_INLINE
inline auto
782 const unsigned int i,
787 result[0] = construct_flux_tensor(
788 pv.template get_tensor<Number, precomputed_state_type>(i));
793 template <
int dim,
typename Number>
794 DEAL_II_ALWAYS_INLINE
inline auto
798 const unsigned int *js,
803 result[0] = construct_flux_tensor(
804 pv.template get_tensor<Number, precomputed_state_type>(js));
809 template <
int dim,
typename Number>
810 DEAL_II_ALWAYS_INLINE
inline auto
814 const dealii::Tensor<1, dim, Number> &c_ij)
const ->
state_type
std::array< Number, n_precomputed_values > precomputed_state_type
dealii::Tensor< 1, problem_dimension, Number > state_type
state_type apply_galilei_transform(const state_type &state, const Lambda &) const
static constexpr unsigned int n_precomputation_cycles
Number kruzkov_entropy_derivative(const Number &k, const Number &u) const
dealii::Tensor< 1, dim, Number > construct_flux_gradient_tensor(const precomputed_state_type &precomputed_state) const
DEAL_II_ALWAYS_INLINE const std::string & flux() const
static constexpr unsigned int n_precomputed_initial_values
static constexpr unsigned int problem_dimension
HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
DEAL_II_ALWAYS_INLINE dealii::Tensor< 1, dim, Number > flux_function(const Number &u) const
flux_contribution_type flux_contribution(const precomputed_vector_type &pv, const precomputed_initial_vector_type &, const unsigned int i, const state_type &) const
dealii::Tensor< 1, problem_dimension, dealii::Tensor< 1, dim, Number > > flux_type
state_type nodal_source(const precomputed_vector_type &pv, const unsigned int *js, const state_type &U_j, const ScalarNumber tau) const =delete
state_type from_primitive_state(const state_type &primitive_state) const
static const auto component_names
state_type apply_boundary_conditions(const dealii::types::boundary_id, const state_type &U, const dealii::Tensor< 1, dim, Number > &, const Lambda &) const
Number square_entropy_derivative(const Number &u) const
state_type nodal_source(const precomputed_vector_type &pv, const unsigned int i, const state_type &U_i, const ScalarNumber tau) const =delete
state_type to_primitive_state(const state_type &state) const
DEAL_II_ALWAYS_INLINE ScalarNumber derivative_approximation_delta() const
Number square_entropy(const Number &u) const
static const auto primitive_component_names
dealii::Tensor< 1, dim, Number > construct_flux_tensor(const precomputed_state_type &precomputed_state) const
static constexpr bool have_high_order_flux
typename get_value_type< Number >::type ScalarNumber
static Number state(const state_type &U)
Number kruzkov_entropy(const Number &k, const Number &u) const
DEAL_II_ALWAYS_INLINE dealii::Tensor< 1, dim, Number > flux_gradient_function(const Number &u) const
static constexpr unsigned int n_precomputed_values
dealii::Tensor< 1, problem_dimension, Number > state_type
std::array< Number, n_precomputed_values > precomputed_state_type
std::array< Number, n_precomputed_initial_values > precomputed_initial_state_type
void precomputation_loop(unsigned int, const DISPATCH &dispatch_check, precomputed_vector_type &, const SPARSITY &, const vector_type &, unsigned int, unsigned int) const
bool is_admissible(const state_type &) const
state_type expand_state(const ST &state) const
static const auto precomputed_names
state_type high_order_flux_divergence(const flux_contribution_type &, const flux_contribution_type &, const dealii::Tensor< 1, dim, Number > &c_ij) const =delete
flux_type flux_contribution_type
static const auto precomputed_initial_names
state_type flux_divergence(const flux_contribution_type &flux_i, const flux_contribution_type &flux_j, const dealii::Tensor< 1, dim, Number > &c_ij) const
static constexpr bool have_source_terms
static std::string problem_name
HyperbolicSystem(const std::string &subsection="/HyperbolicSystem")
std::set< std::shared_ptr< Flux > > flux_list_type
void populate_flux_list(flux_list_type &flux_list, const std::string &subsection)
DEAL_II_ALWAYS_INLINE FT add(const FT &flux_left_ij, const FT &flux_right_ij)
DEAL_II_ALWAYS_INLINE dealii::Tensor< 1, problem_dim, T > contract(const FT &flux_ij, const TT &c_ij)