From 67fe0f5003008a114c5d15fe107d02b40c4f4726 Mon Sep 17 00:00:00 2001 From: Eric Tovar Date: Tue, 10 Feb 2026 20:10:09 -0700 Subject: [PATCH 1/5] MultiSpeciesEuler: add files for new module --- source/multi_euler_ideal/CMakeLists.txt | 17 + source/multi_euler_ideal/description.h | 84 + source/multi_euler_ideal/equation_dispatch.cc | 18 + source/multi_euler_ideal/hyperbolic_system.h | 1637 +++++++++++++++++ source/multi_euler_ideal/indicator.h | 276 +++ .../initial_state_contrast.h | 96 + .../initial_state_exact_riemann_solution.h | 536 ++++++ .../initial_state_function.h | 149 ++ .../initial_state_icf_like.h | 207 +++ .../initial_state_library.cc | 14 + .../initial_state_library.template.h | 41 + .../initial_state_library_multi_euler.h | 46 + .../initial_state_radial_contrast.h | 127 ++ .../initial_state_shock_bubble.h | 183 ++ .../initial_state_smooth_wave.h | 107 ++ .../initial_state_three_state_contrast.h | 131 ++ source/multi_euler_ideal/instantiate.h | 18 + source/multi_euler_ideal/limiter.cc | 25 + source/multi_euler_ideal/limiter.h | 549 ++++++ source/multi_euler_ideal/limiter.template.h | 502 +++++ source/multi_euler_ideal/riemann_solver.cc | 25 + source/multi_euler_ideal/riemann_solver.h | 273 +++ .../riemann_solver.template.h | 661 +++++++ 23 files changed, 5722 insertions(+) create mode 100644 source/multi_euler_ideal/CMakeLists.txt create mode 100644 source/multi_euler_ideal/description.h create mode 100644 source/multi_euler_ideal/equation_dispatch.cc create mode 100644 source/multi_euler_ideal/hyperbolic_system.h create mode 100644 source/multi_euler_ideal/indicator.h create mode 100644 source/multi_euler_ideal/initial_state_contrast.h create mode 100644 source/multi_euler_ideal/initial_state_exact_riemann_solution.h create mode 100644 source/multi_euler_ideal/initial_state_function.h create mode 100644 source/multi_euler_ideal/initial_state_icf_like.h create mode 100644 source/multi_euler_ideal/initial_state_library.cc create mode 100644 source/multi_euler_ideal/initial_state_library.template.h create mode 100644 source/multi_euler_ideal/initial_state_library_multi_euler.h create mode 100644 source/multi_euler_ideal/initial_state_radial_contrast.h create mode 100644 source/multi_euler_ideal/initial_state_shock_bubble.h create mode 100644 source/multi_euler_ideal/initial_state_smooth_wave.h create mode 100644 source/multi_euler_ideal/initial_state_three_state_contrast.h create mode 100644 source/multi_euler_ideal/instantiate.h create mode 100644 source/multi_euler_ideal/limiter.cc create mode 100644 source/multi_euler_ideal/limiter.h create mode 100644 source/multi_euler_ideal/limiter.template.h create mode 100644 source/multi_euler_ideal/riemann_solver.cc create mode 100644 source/multi_euler_ideal/riemann_solver.h create mode 100644 source/multi_euler_ideal/riemann_solver.template.h diff --git a/source/multi_euler_ideal/CMakeLists.txt b/source/multi_euler_ideal/CMakeLists.txt new file mode 100644 index 00000000..7b169e0d --- /dev/null +++ b/source/multi_euler_ideal/CMakeLists.txt @@ -0,0 +1,17 @@ +## +## SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +## Copyright (C) 2023 - 2025 by the ryujin authors +## Copyright (C) 2025 by Triad National Security, LLC +## + +add_library(obj_multi_euler_ideal OBJECT + equation_dispatch.cc + initial_state_library.cc + limiter.cc + riemann_solver.cc + ) +set_target_properties(obj_multi_euler_ideal PROPERTIES LINKER_LANGUAGE CXX) +deal_ii_setup_target(obj_multi_euler_ideal) +target_link_libraries(obj_multi_euler_ideal obj_common ${EXTERNAL_TARGETS}) +# Propagate the current source directory with PUBLIC visibility +target_include_directories(obj_multi_euler_ideal PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/source/multi_euler_ideal/description.h b/source/multi_euler_ideal/description.h new file mode 100644 index 00000000..9ee0ef15 --- /dev/null +++ b/source/multi_euler_ideal/description.h @@ -0,0 +1,84 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include + +#include "../stub_parabolic_module.h" +#include "../stub_parabolic_system.h" +#include "hyperbolic_system.h" +#include "indicator.h" +#include "limiter.h" +#include "riemann_solver.h" + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + /** + * A struct that contains all equation specific classes describing the + * chosen hyperbolic system, the indicator, the limiter and + * (approximate) Riemann solver. + * + * The compressible multi-species Euler equations of gas dynamics for + * an ideal gas mixture under the assumption of thermal-mechanical + * equilibrium. Each species \f$k\f$ is modeled by an ideal gas with + * specific heat capacities \f$c_{p,k}\f$ and \f$c_{v,k}\f$. + * + * The governing equations are: + * \f{align} + * \partial_t (\alpha_k \rho_k) + \nabla \cdot (\alpha_k \rho_k + * \mathbf{v}) &= 0, \quad k = 1, \ldots, n_s, \\ + * \partial_t \mathbf{m} + \nabla \cdot (\mathbf{v} \otimes \mathbf{m} + * + p \mathbb{I}_d) &= 0, \\ + * \partial_t E + \nabla \cdot ((E + p) \mathbf{v}) &= 0, + * \f} + * where \f$\alpha_k \rho_k\f$ are the partial densities, \f$\mathbf{m} + * = \rho \mathbf{v}\f$ is the mixture momentum, \f$E\f$ is the mixture + * total energy, and \f$p\f$ is the mixture pressure. The mixture + * pressure is given by \f$p = (\bar{\gamma} - 1) \varepsilon\f$ where + * \f$\bar{\gamma} = \bar{c}_p / \bar{c}_v\f$ is the mixture ratio of + * specific heats. + * + * The numerical scheme ensures preservation of an invariant domain + * \f$\mathcal{A} = \{ \mathbf{u} : \alpha_k \rho_k \geq 0,\; + * \varepsilon(\mathbf{u}) > 0,\; s(\mathbf{u}) \geq s_{\min} \}\f$, + * where \f$s(\mathbf{u})\f$ is the mixture specific entropy. + * + * The parabolic subsystem is chosen to be the identity (inviscid flow). + * + * @note See "Second-order invariant-domain preserving approximation to + * the multi-species Euler equations" by Clayton, Dzanic, and Tovar + * (arXiv:2505.09581) for full details of the mathematical formulation + * and numerical scheme. + * + * @ingroup MultiSpeciesEulerEquations + */ + struct Description { + using HyperbolicSystem = MultiSpeciesEuler::HyperbolicSystem; + + template + using HyperbolicSystemView = + MultiSpeciesEuler::HyperbolicSystemView; + + using ParabolicSystem = ryujin::StubParabolicSystem; + + template + using ParabolicModule = + ryujin::StubParabolicModule; + + template + using Indicator = MultiSpeciesEuler::Indicator; + + template + using Limiter = MultiSpeciesEuler::Limiter; + + template + using RiemannSolver = MultiSpeciesEuler::RiemannSolver; + }; + } // namespace MultiSpeciesEuler +} // namespace ryujin diff --git a/source/multi_euler_ideal/equation_dispatch.cc b/source/multi_euler_ideal/equation_dispatch.cc new file mode 100644 index 00000000..a17fc589 --- /dev/null +++ b/source/multi_euler_ideal/equation_dispatch.cc @@ -0,0 +1,18 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2024 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#include "description.h" + +#include +#include + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + Dispatch dispatch_instance("multi species euler"); + } // namespace MultiSpeciesEuler +} // namespace ryujin diff --git a/source/multi_euler_ideal/hyperbolic_system.h b/source/multi_euler_ideal/hyperbolic_system.h new file mode 100644 index 00000000..9da92bbc --- /dev/null +++ b/source/multi_euler_ideal/hyperbolic_system.h @@ -0,0 +1,1637 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2020 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + /* + * For various divisions in the multi-species module we have a + * mathematical guarantee that the numerator and denominator are + * nonnegative and the limit (of zero numerator and denominator) must + * converge to zero. The following function takes care of rounding + * issues when computing such quotients by (a) avoiding division by + * zero and (b) ensuring non-negativity of the result. + */ + template + DEAL_II_ALWAYS_INLINE inline Number safe_division(const Number &numerator, + const Number &denominator) + { + using ScalarNumber = typename get_value_type::type; + constexpr ScalarNumber min = std::numeric_limits::min(); + + return std::max(numerator, Number(0.)) / + std::max(denominator, Number(min)); + } + + + template + class HyperbolicSystemView; + + /** + * The compressible multi-species Euler equations for a gas mixture + * under the assumption of thermal-mechanical equilibrium. Each species + * \f$k\f$ is modeled by an ideal gas with specific heat capacities + * \f$c_{p,k}\f$ and \f$c_{v,k}\f$ at constant pressure and volume, + * respectively. + * + * We have a \f$(n_s + 1 + d)\f$ dimensional state space + * \f[ + * \mathbf{u} = [(\alpha_k \rho_k)_{k=1}^{n_s}, \mathbf{m}, E]^T + * \in \mathbb{R}^{n_s + 1 + d}, + * \f] + * where \f$\alpha_k \rho_k\f$ denotes the partial densities, + * \f$\mathbf{m} = \rho \mathbf{v}\f$ is the mixture momentum, and + * \f$E\f$ is the mixture total energy. + * + * The mixture density is \f$\rho = \sum_k \alpha_k \rho_k\f$ and the + * mass fractions are \f$Y_k = \alpha_k \rho_k / \rho\f$. The internal + * energy density is \f$\varepsilon = E - \frac{1}{2}\rho|\mathbf{v}|^2\f$. + * + * The mixture pressure is given by an ideal mixture equation of state + * (see Eq. (2.4) in [ClaytonDzanicTovar-2025]): + * \f[ + * p = (\bar{\gamma}(\mathbf{Y}) - 1) \varepsilon, + * \f] + * where \f$\bar{\gamma} = \bar{c}_p / \bar{c}_v\f$ with + * \f$\bar{c}_p = \sum_k Y_k c_{p,k}\f$ and + * \f$\bar{c}_v = \sum_k Y_k c_{v,k}\f$. + * + * The mixture specific entropy is (see Eq. (2.8) in @cite + * ClaytonDzanicTovar-2025): + * \f[ + * s(\mathbf{u}) = \bar{c}_v \log\left(\frac{\rho e}{\rho^{\bar{\gamma}}} + * \right) + K(\mathbf{Y}), + * \f] + * where \f$e = \varepsilon / \rho\f$ is the specific internal energy and + * \f$K(\mathbf{Y})\f$ is a mixing term depending only on mass fractions. + * + * The invariant domain preserved by the numerical scheme is (see Eq. + * (2.11) in [ClaytonDzanicTovar-2025]): + * \f[ + * \mathcal{A} = \{ \mathbf{u} : \alpha_k \rho_k \geq 0 \;\forall k,\; + * \varepsilon(\mathbf{u}) > 0,\; s(\mathbf{u}) \geq s_{\min} \}. + * \f] + * + * @note Currently hardcoded for 2 species (\f$n_s = 2\f$). + * + * @ingroup MultiSpeciesEulerEquations + */ + class HyperbolicSystem final : public dealii::ParameterAcceptor + { + public: + /** + * The name of the hyperbolic system as a string. + */ + static inline std::string problem_name = + "Compressible multi species Euler equations (ideal mixture)"; + + /** + * Constructor. + */ + HyperbolicSystem(const std::string &subsection = "/HyperbolicSystem"); + + /** + * Return a view on the Hyperbolic System for a given dimension @p + * dim and choice of number type @p Number (which can be a scalar + * float, or double, as well as a VectorizedArray holding packed + * scalars. + */ + template + auto view() const + { + return HyperbolicSystemView{*this}; + } + + /** + * Part of step 1 of the hyperbolic update step: Compute "precomputed + * values" and fill into the state vector. + * + * @note The method does not update the ghost range of the state + * vector. The precomputed part has to be synchronized by explicitly + * calling the update ghost values function. + */ + template + void fill_precomputed_values( + const OfflineData &offline_data, + typename HyperbolicSystemView::StateVector + &state_vector, + const bool skip_constrained_dofs = true) const; + + private: + /** + * @name Runtime parameters, internal fields, methods, and friends + */ + //@{ + + dealii::Tensor<1, /* num_species = */ 2, double> cp_for_each_species_; + dealii::Tensor<1, /* num_species = */ 2, double> cv_for_each_species_; + dealii::Tensor<1, /* num_species = */ 2, double> r_for_each_species_; + dealii::Tensor<1, /* num_species = */ 2, double> gamma_for_each_species_; + double reference_density_; + double vacuum_state_relaxation_small_; + double vacuum_state_relaxation_large_; + bool compute_strict_bounds_; + + template + friend class HyperbolicSystemView; + //@} + }; /* HyperbolicSystem */ + + + /** + * A view of the HyperbolicSystem that makes methods available for a + * given dimension @p dim and choice of number type @p Number (which + * can be a scalar float, or double, as well as a VectorizedArray + * holding packed scalars. + * + * Intended usage: + * ``` + * HyperbolicSystem hyperbolic_system; + * const auto view = hyperbolic_system.template view(); + * const auto flux_i = view.flux_contribution(...); + * const auto flux_j = view.flux_contribution(...); + * const auto flux_ij = view.flux_divergence(flux_i, flux_j, c_ij); + * // etc. + * ``` + */ + template + class HyperbolicSystemView + { + public: + /** + * Constructor taking a reference to the underlying + * HyperbolicSystem + */ + HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system) + : hyperbolic_system_(hyperbolic_system) + { + } + + /** + * Create a modified view from the current one: + */ + template + auto view() const + { + return HyperbolicSystemView{hyperbolic_system_}; + } + + /** + * The underlying scalar number type. + */ + using ScalarNumber = typename get_value_type::type; + + /** + * @name Access to runtime parameters + */ + //@{ + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 2, double> + cp_for_each_species() const + { + return hyperbolic_system_.cp_for_each_species_; + } + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 2, double> + cv_for_each_species() const + { + return hyperbolic_system_.cv_for_each_species_; + } + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 2, double> + r_for_each_species() const + { + return hyperbolic_system_.r_for_each_species_; + } + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 2, double> + gamma_for_each_species() const + { + return hyperbolic_system_.gamma_for_each_species_; + } + + + DEAL_II_ALWAYS_INLINE inline ScalarNumber reference_density() const + { + return hyperbolic_system_.reference_density_; + } + + DEAL_II_ALWAYS_INLINE inline ScalarNumber + vacuum_state_relaxation_small() const + { + return hyperbolic_system_.vacuum_state_relaxation_small_; + } + + DEAL_II_ALWAYS_INLINE inline ScalarNumber + vacuum_state_relaxation_large() const + { + return hyperbolic_system_.vacuum_state_relaxation_large_; + } + + DEAL_II_ALWAYS_INLINE inline bool compute_strict_bounds() const + { + return hyperbolic_system_.compute_strict_bounds_; + } + + /** + * constexpr boolean used in the MultiSpeciesEulerInitialStates namespace + */ + static constexpr bool have_gamma = false; + + /** + * constexpr boolean used in the MultiSpeciesEulerInitialStates namespace + */ + static constexpr bool have_eos_interpolation_b = false; + + /** + * constexpr boolean for energy equation presence + */ + static constexpr bool have_energy_equation = true; + + + //@} + /** + * @name Internal data + */ + //@{ + + private: + const HyperbolicSystem &hyperbolic_system_; + + public: + //@} + /** + * @name Types and constexpr constants + */ + //@{ + + /** + * The dimension of the state space. + * FIXME: Generalize for user-defined number of species. + * This is assuming 2 species only. + */ + static constexpr unsigned int problem_dimension = 3 + dim; + + /** + * Storage type for a (conserved) state vector \f$\boldsymbol U\f$. + */ + using state_type = dealii::Tensor<1, problem_dimension, Number>; + + /** + * Storage type for the flux \f$\mathbf{f}\f$. + */ + using flux_type = + dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>; + + /** + * Storage type for the "total mixture" state vector \f$\boldsymbol U\f$. + */ + using mixture_state_type = dealii::Tensor<1, 2 + dim, Number>; + + /** + * Storage type for the "total mixture" flux \f$\mathbf{f^{mix}}\f$. + */ + using mixture_flux_type = + dealii::Tensor<1, 2 + dim, dealii::Tensor<1, dim, Number>>; + + /** + * The storage type used for flux contributions. + */ + using flux_contribution_type = flux_type; + + /** + * An array holding all component names of the conserved state as a + * string. + */ + static inline const auto component_names = + []() -> std::array { + if constexpr (dim == 1) + return {"alpha_rho_0", "alpha_rho_1", "m", "E"}; + else if constexpr (dim == 2) + return {"alpha_rho_0", "alpha_rho_1", "m_1", "m_2", "E"}; + else if constexpr (dim == 3) + return {"alpha_rho_0", "alpha_rho_1", "m_1", "m_2", "m_3", "E"}; + __builtin_trap(); + }(); + + /** + * An array holding all component names of the primitive state as a + * string. + */ + static inline const auto primitive_component_names = + []() -> std::array { + if constexpr (dim == 1) + return {"alpha_rho_0", "alpha_rho_1", "v", "p"}; + else if constexpr (dim == 2) + return {"alpha_rho_0", "alpha_rho_1", "v_1", "v_2", "p"}; + else if constexpr (dim == 3) + return {"alpha_rho_0", "alpha_rho_1", "v_1", "v_2", "v_3", "p"}; + __builtin_trap(); + }(); + + /** + * The number of precomputed values. + */ + static constexpr unsigned int n_precomputed_values = 5; + + /** + * Array type used for precomputed values. + */ + using precomputed_type = std::array; + + /** + * An array holding all component names of the precomputed values. + */ + static inline const auto precomputed_names = + std::array{ + {"rho", + "p", + "surrogate_gamma_min", + "specific_entropy", + "surrogate_harten_entropy"}}; + + /** + * The number of precomputed initial values. + */ + static constexpr unsigned int n_initial_precomputed_values = 0; + + /** + * Array type used for precomputed initial values. + */ + using initial_precomputed_type = + std::array; + + /** + * An array holding all component names of the precomputed values. + */ + static inline const auto initial_precomputed_names = + std::array{}; + + /** + * A compound state vector. + */ + using StateVector = Vectors:: + StateVector; + + /** + * MulticomponentVector for storing the hyperbolic state vector: + */ + using HyperbolicVector = + Vectors::MultiComponentVector; + + /** + * MulticomponentVector for storing a vector of precomputed states: + */ + using PrecomputedVector = + Vectors::MultiComponentVector; + + /** + * MulticomponentVector for storing a vector of precomputed initial + * states: + */ + using InitialPrecomputedVector = + Vectors::MultiComponentVector; + + //@} + /** + * @name Computing derived physical quantities + */ + //@{ + + /** + * For a given (3+dim dimensional) state vector U, return + * the partial density for the first species U[0]. + */ + static Number partial_density_0(const state_type &U); + + /** + * For a given (3+dim dimensional) state vector U, return + * the partial density for the second species U[1]. + */ + static Number partial_density_1(const state_type &U); + + /** + * For a given (3+dim dimensional) state vector U, return + * the total mixture density obtained by summing partial densities. + */ + static Number density(const state_type &U); + + /** + * Given a density @p rho this function returns 0 if the magnitude + * of rho is smaller or equal than relaxation_large * rho_cutoff. + * Otherwise rho is returned unmodified. Here, rho_cutoff is the + * reference density multiplied by eps. + */ + Number filter_vacuum_density(const Number &rho) const; + + /** + * For a given (3+dim dimensional) state vector U, return + * the momentum vector [U[2], ..., U[1+dim]]. + */ + static dealii::Tensor<1, dim, Number> momentum(const state_type &U); + + /** + * For a given (3+dim dimensional) state vector U, return + * the total energy U[2+dim] + */ + static Number total_energy(const state_type &U); + + /** + * For a given (3+dim dimensional) state vector U, compute + * and return the internal energy \f$\varepsilon = (\rho e)\f$. + */ + static Number internal_energy(const state_type &U); + + /** + * For a given (3+dim dimensional) state vector U, compute + * and return the mixture gamma for the mixture EOS: + * \f[ + * \overline{\gamma} = (sum of (alpha_k rho_k) * c_{p, k}) / (sum of + * (alpha_k rho_k) * c_{v, k}) + * \f] + */ + Number gamma_mixture(const state_type &U) const; + + /** + * For a given (3+dim dimensional) state vector U, compute + * and return the mixture specific heat capacity at constant pressure: + * \f[ + * \overline{c}_p = (sum of (alpha_k rho_k) / rho * c_{p, k}) + * \f] + */ + Number cp_mixture(const state_type &U) const; + + /** + * For a given (3+dim dimensional) state vector U, compute + * and return the mixture specific heat capacity at constant volume: + * \f[ + * \overline{c}_v = (sum of (alpha_k rho_k) / rho * c_{v, k}) + * \f] + */ + Number cv_mixture(const state_type &U) const; + + /** + * For a given (3+dim dimensional) state vector U, compute + * and return the pressure \f$p\f$. + * + * We assume that the pressure is given by a mixture ideal EOS + * \f[ + * p = (\overline{\gamma} - 1)\;(\rho e) + * \f] + */ + Number pressure(const state_type &U) const; + + /** + * For a given (3+dim dimensional) state vector U, + * compute and return the speed of sound \f$c\f$: + * + * We assume that the pressure is given by a mixture ideal EOS + * \f[ + * c = sqrt(\overline{\gamma} p / rho) + * \f] + */ + Number speed_of_sound(const state_type &U) const; + + //@} + /** + * @name Surrogate functions for computing various interpolatory + * physical quantities that are needed for Riemann solver, + * indicator and limiter. + */ + //@{ + + /** + * For a given (3+dim dimensional) state vector U, compute + * and return the physical specific entropy. Following Eq. (2.8) in + * [ClaytonDzanicTovar-2025]: + * \f[ + * s(\mathbf{u}) = \bar{r} \log(\rho^{-1}) + \bar{c}_v \log(e) + * + K(\mathbf{Y}), + * \f] + * where \f$\bar{r} = \sum_k Y_k r_k\f$, \f$\bar{c}_v = \sum_k Y_k + * c_{v,k}\f$, \f$e = \varepsilon / \rho\f$ is the specific internal + * energy, and \f$K(\mathbf{Y})\f$ is a mixing term (Eq. 2.9). + */ + Number specific_entropy(const state_type &U) const; + + /** + * For a given (3+dim dimensional) state vector U, compute + * and return a surrogate Harten-type entropy. Following Section 4 of + * [ClaytonDzanicTovar-2025], we use: + * \f[ + * \eta(\mathbf{u}; \gamma_{\min}) = (\rho \varepsilon)^{1/(\gamma_{\min}+1)}, + * \f] + * where \f$\gamma_{\min}\f$ is the minimum surrogate gamma over the + * stencil. This entropy is chosen to ensure convexity properties + * required by the entropy-viscosity indicator. + */ + Number surrogate_harten_entropy(const state_type &U, + const Number &gamma_min) const; + + /** + * For a given (3+dim dimensional) state vector U, compute + * and return the derivative \f$\eta'\f$ of the Harten-type entropy. + */ + mixture_state_type + surrogate_harten_entropy_derivative(const state_type &U, + const Number &eta, + const Number &gamma_min) const; + + /** + * For a given (3+dim dimensional) state vector U and + * pressure p, compute a surrogate gamma. Following + * Section 3 of [ClaytonDzanicTovar-2025]: + * \f[ + * \gamma(\mathbf{u}, p) = 1 + \frac{p}{\varepsilon(\mathbf{u})}. + * \f] + * This surrogate gamma allows us to use the single-species Riemann + * solver machinery with a locally linearized equation of state. It + * is computed point-wise and the minimum over the stencil is used + * for the Riemann solver and limiter. + */ + Number surrogate_gamma(const state_type &U, const Number &p) const; + + /** + * For a given (3+dim dimensional) state vector U and + * gamma gamma, compute a surrogate pressure: + * \f[ + * p(\rho, e, \gamma) = (\gamma - 1) (\rho e) + * \f] + * + * This function is the complementary function to surrogate_gamma(). + */ + Number surrogate_pressure(const state_type &U, const Number &gamma) const; + + /** + * For a given (3+dim dimensional) state vector U and + * gamma gamma, compute a surrogate speed of sound. + */ + Number surrogate_speed_of_sound(const state_type &U, + const Number &gamma) const; + + /** + * Returns whether the state @p U is admissible. If @p U is a + * vectorized state then @p U is admissible if all vectorized values + * are admissible. + */ + bool is_admissible(const state_type &U) const; + + //@} + /** + * @name Special functions for boundary states + */ + //@{ + + /** + * Decomposes a given state @p U into Riemann invariants and then + * replaces the first or second Riemann characteristic from the one + * taken from @p U_bar state. Note that the @p U_bar state is just the + * prescribed dirichlet values. + */ + template + state_type prescribe_riemann_characteristic( + const state_type &U, + const Number &p, + const state_type &U_bar, + const Number &p_bar, + const dealii::Tensor<1, dim, Number> &normal) const; + + /** + * Apply boundary conditions. + * + * For the compressible Euler equations we have: + * + * - Dirichlet boundary conditions by prescribing the return value of + * get_dirichlet_data() as is. + * + * - Slip boundary conditions where we remove the normal component of + * the momentum. + * + * - No slip boundary conditions where we set the momentum to 0. + * + * - "Dynamic boundary" conditions that prescribe different Riemann + * invariants from the return value of get_dirichlet_data() + * depending on the flow state (supersonic versus subsonic, outflow + * versus inflow). + */ + template + state_type + apply_boundary_conditions(const dealii::types::boundary_id id, + const state_type &U, + const dealii::Tensor<1, dim, Number> &normal, + const Lambda &get_dirichlet_data) const; + + //@} + /** + * @name Flux computations + */ + //@{ + + /** + * Given a state @p U and a pressure @p p compute the flux + * \f[ + * \begin{pmatrix} + * \textbf v \alpha_1 \rho_1 \\ + * \textbf v \alpha_2 \rho_2 + * \textbf v\otimes \textbf m + p\mathbb{I}_d \\ + * \textbf v(E+p) + * \end{pmatrix}, + * \f] + */ + flux_type f(const state_type &U, const Number &p) const; + + /** + * Given a state @p U and a pressure @p p compute the "summed system" + * mixture flux needed for indicator: + * \f[ + * \begin{pmatrix} + * \textbf m \\ + * \textbf v\otimes \textbf m + p\mathbb{I}_d \\ + * \textbf v(E+p) + * \end{pmatrix}, + * \f] + */ + mixture_flux_type mixture_f(const state_type &U, const Number &p) const; + + /** + * Given a state @p U_i and an index @p i compute flux contributions. + * + * Intended usage: + * ``` + * Indicator indicator; + * for (unsigned int i = n_internal; i < n_owned; ++i) { + * // ... + * const auto flux_i = flux_contribution(precomputed..., i, U_i); + * for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) { + * // ... + * const auto flux_j = flux_contribution(precomputed..., js, U_j); + * const auto flux_ij = flux_divergence(flux_i, flux_j, c_ij); + * } + * } + * ``` + * + * For the Euler equations we simply compute f(U_i). + */ + flux_contribution_type + flux_contribution(const PrecomputedVector &pv, + const InitialPrecomputedVector &piv, + const unsigned int i, + const state_type &U_i) const; + + flux_contribution_type + flux_contribution(const PrecomputedVector &pv, + const InitialPrecomputedVector &piv, + const unsigned int *js, + const state_type &U_j) const; + + /** + * Given flux contributions @p flux_i and @p flux_j compute the flux + * (-f(U_i) - f(U_j) + */ + 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; + + /** + * The low-order and high-order fluxes are the same: + */ + static constexpr bool have_high_order_flux = false; + + state_type high_order_flux_divergence( + const flux_contribution_type &flux_i, + const flux_contribution_type &flux_j, + const dealii::Tensor<1, dim, Number> &c_ij) const = delete; + + /** + * @name Computing stencil source terms + */ + //@{ + + /** We do not have source terms */ + static constexpr bool have_source_terms = false; + + state_type nodal_source(const PrecomputedVector &pv, + const unsigned int i, + const state_type &U_i, + const ScalarNumber tau) const = delete; + + state_type nodal_source(const PrecomputedVector &pv, + const unsigned int *js, + const state_type &U_j, + const ScalarNumber tau) const = delete; + + //@} + /** + * @name State transformations + */ + //@{ + + /** + * Given a state vector associated with a different spatial + * dimensions than the current one, return an "expanded" version of + * the state vector associated with @a dim spatial dimensions where + * the momentum vector of the conserved state @p state is expaned + * with zeros to a total length of @a dim entries. + * + * @note @a dim has to be larger or equal than the dimension of the + * @a ST vector. + */ + template + state_type expand_state(const ST &state) const; + + /** + * Given an initial state [alpha_0 rho_0, alpha_1 rho_1, u_1, ..., u_d, p] + * return a conserved state [alpha_0 rho_0, alpha_1 rho_1, m_1, ..., m_d, + * E]. + * + * @note This function is used to conveniently convert (user + * provided) primitive initial states with pressure values to a + * conserved state in the MultiSpeciesEulerInitialStateLibrary. + */ + template + state_type from_initial_state(const ST &initial_state) const; + + /** + * Given a primitive state [rho, u_1, ..., u_d, e] return a conserved + * state. + */ + state_type from_primitive_state(const state_type &primitive_state) const; + + /** + * Given a conserved state return a primitive state [rho, u_1, ..., u_d, + * e] + */ + state_type to_primitive_state(const state_type &state) const; + + /** + * Transform the current state according to a given operator + * @p lambda acting on a @a dim dimensional momentum (or velocity) + * vector. + */ + template + state_type apply_galilei_transform(const state_type &state, + const Lambda &lambda) const; + //@} + }; /* HyperbolicSystemView */ + + + /* + * ------------------------------------------------------------------------- + * Inline definitions + * ------------------------------------------------------------------------- + */ + + + inline HyperbolicSystem::HyperbolicSystem( + const std::string &subsection /*= "HyperbolicSystem"*/) + : ParameterAcceptor(subsection) + { + + compute_strict_bounds_ = true; + add_parameter( + "compute strict bounds", + compute_strict_bounds_, + "Compute strict, but significantly more expensive bounds at various " + "places: (a) an expensive, but better upper wavespeed estimate in " + "the approximate RiemannSolver; (b) entropy viscosity-commutator " + "with correct gamma_min over the stencil; (c) mathematically correct " + "surrogate specific entropy minimum with gamma_min over the " + "stencil."); + + reference_density_ = 1.; + add_parameter("reference density", + reference_density_, + "Problem specific density reference"); + + cp_for_each_species_[0] = 1.4; + cv_for_each_species_[0] = 1.0; + + cp_for_each_species_[1] = 1.67; + cv_for_each_species_[1] = 1.0; + + add_parameter( + "c_p for each species", + cp_for_each_species_, + "Specific heat capacity at constant pressure for each species"); + + add_parameter( + "c_v for each species", + cv_for_each_species_, + "Specific heat capacity at constant volume for each species"); + + vacuum_state_relaxation_small_ = 1.e2; + add_parameter("vacuum state relaxation small", + vacuum_state_relaxation_small_, + "Problem specific vacuum relaxation parameter"); + + vacuum_state_relaxation_large_ = 1.e4; + add_parameter("vacuum state relaxation large", + vacuum_state_relaxation_large_, + "Problem specific vacuum relaxation parameter"); + + /* + * And finally populate the r and gamma values. + */ + const auto populate_values = [this]() { + r_for_each_species_[0] = + cp_for_each_species_[0] - cv_for_each_species_[0]; + r_for_each_species_[1] = + cp_for_each_species_[1] - cv_for_each_species_[1]; + + gamma_for_each_species_[0] = + cp_for_each_species_[0] / cv_for_each_species_[0]; + gamma_for_each_species_[1] = + cp_for_each_species_[1] / cv_for_each_species_[1]; + }; + + ParameterAcceptor::parse_parameters_call_back.connect(populate_values); + populate_values(); + } + + + template + inline void HyperbolicSystem::fill_precomputed_values( + const OfflineData &offline_data, + typename HyperbolicSystemView::StateVector + &state_vector, + const bool skip_constrained_dofs) const + { + const unsigned int n_internal = offline_data.n_locally_internal(); + const unsigned int n_owned = offline_data.n_locally_owned(); + const auto &sparsity_simd = offline_data.sparsity_pattern_simd(); + + const auto &U = std::get<0>(state_vector); + auto &precomputed = std::get<1>(state_vector); + + /* Compute values over the diagonal: */ + + const auto body = [&](auto sentinel, unsigned int i) { + using T = decltype(sentinel); + using View = HyperbolicSystemView; + using precomputed_type = typename View::precomputed_type; + + const unsigned int row_length = sparsity_simd.row_length(i); + if (skip_constrained_dofs && row_length == 1) + return; + + const auto U_i = U.template get_tensor(i); + const auto view = this->view(); + const auto rho_i = view.density(U_i); + const auto p_i = view.pressure(U_i); + const auto s_i = view.specific_entropy(U_i); + + const auto gamma_i = view.surrogate_gamma(U_i, p_i); + using PT = precomputed_type; + const PT prec_i{rho_i, p_i, gamma_i, s_i, T(0.)}; + precomputed.template write_tensor(prec_i, i); + }; + + cpu_simd_loop("time_step_1", body, 0, n_internal, n_owned); + precomputed.update_ghost_values(); + + /* Compute gamma_min over the stencil: */ + + const auto body_stencil = [&](auto sentinel, unsigned int i) { + using T = decltype(sentinel); + using View = HyperbolicSystemView; + using PT = typename View::precomputed_type; + + const unsigned int row_length = sparsity_simd.row_length(i); + if (skip_constrained_dofs && row_length == 1) + return; + + const auto U_i = U.template get_tensor(i); + auto prec_i = precomputed.template get_tensor(i); + auto &[rho_i, p_i, gamma_min_i, s_i, harten_i] = prec_i; + + const auto view = this->view(); + + constexpr unsigned int stride_size = get_stride_size; + const unsigned int *js = sparsity_simd.columns(i) + stride_size; + for (unsigned int col_idx = 1; col_idx < row_length; + ++col_idx, js += stride_size) { + + const auto U_j = U.template get_tensor(js); + const auto prec_j = precomputed.template get_tensor(js); + const auto p_j = std::get<1>(prec_j); + const auto gamma_j = view.surrogate_gamma(U_j, p_j); + gamma_min_i = std::min(gamma_min_i, gamma_j); + } + + harten_i = view.surrogate_harten_entropy(U_i, gamma_min_i); + precomputed.template write_tensor(prec_i, i); + }; + + cpu_simd_loop( + "time_step_1", body_stencil, 0, n_internal, n_owned); + } + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::partial_density_0(const state_type &U) + { + return U[0]; + } + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::partial_density_1(const state_type &U) + { + return U[1]; + } + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::density(const state_type &U) + { + return partial_density_0(U) + partial_density_1(U); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::filter_vacuum_density( + const Number &rho) const + { + constexpr ScalarNumber eps = std::numeric_limits::epsilon(); + const Number rho_cutoff_large = + reference_density() * vacuum_state_relaxation_large() * eps; + + return dealii::compare_and_apply_mask( + std::abs(rho), rho_cutoff_large, Number(0.), rho); + } + + + template + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number> + HyperbolicSystemView::momentum(const state_type &U) + { + dealii::Tensor<1, dim, Number> result; + for (unsigned int i = 0; i < dim; ++i) + result[i] = U[2 + i]; + return result; + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::total_energy(const state_type &U) + { + return U[2 + dim]; + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::internal_energy(const state_type &U) + { + /* + * rho e = (E - 1/2*m^2/rho) + */ + const Number rho_inverse = ScalarNumber(1.) / density(U); + const auto m = momentum(U); + const Number E = total_energy(U); + return E - ScalarNumber(0.5) * m.norm_square() * rho_inverse; + } + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::cp_mixture(const state_type &U) const + { + Number result = 0.; + const Number rho = density(U); + + for (unsigned int i = 0; i < 2; ++i) { + result += U[i] / rho * ScalarNumber(cp_for_each_species()[i]); + } + return result; + } + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::cv_mixture(const state_type &U) const + { + Number result = 0.; + const Number rho = density(U); + + for (unsigned int i = 0; i < 2; ++i) { + result += U[i] / rho * ScalarNumber(cv_for_each_species()[i]); + } + return result; + } + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::gamma_mixture(const state_type &U) const + { + return cp_mixture(U) / cv_mixture(U); + } + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::pressure(const state_type &U) const + { + /* p = (\overline{gamma} - 1) * (rho e) */ + return (gamma_mixture(U) - Number(1.)) * internal_energy(U); + } + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::speed_of_sound(const state_type &U) const + { + /* c = sqrt(\overline{gamma} * p / rho) */ + const auto gamma = gamma_mixture(U); + const auto p = pressure(U); + const auto rho = density(U); + return std::sqrt(gamma * p / rho); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::specific_entropy( + const state_type &U) const + { + using ScalarNumber = typename get_value_type::type; + + const auto rho = density(U); + const auto rho_inverse = ScalarNumber(1.) / rho; + const auto e = internal_energy(U) * rho_inverse; + + Number cv_bar = 0.; + Number r_bar = 0.; + + for (unsigned int i = 0; i < 2; ++i) { + const auto Y_i = U[i] * rho_inverse; + r_bar += Y_i * Number(r_for_each_species()[i]); + cv_bar += Y_i * Number(cv_for_each_species()[i]); + } + + Number K_factor = 0.; + for (unsigned int i = 0; i < 2; ++i) { + const auto Y_i = U[i] * rho_inverse; + const auto cv_i = Number(cv_for_each_species()[i]); + const auto r_i = Number(r_for_each_species()[i]); + const auto gm1 = Number(gamma_for_each_species()[i] - 1.); + + const auto K_i = + cv_i * std::log(cv_i / cv_bar * ryujin::pow(r_i / r_bar, gm1)); + K_factor += Y_i * K_i; + } + + const auto s_bar = + r_bar * std::log(rho_inverse) + cv_bar * std::log(e) + K_factor; + + return s_bar; + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::surrogate_harten_entropy( + const state_type &U, const Number &gamma_min) const + { + const Number rho = density(U); + const auto m = momentum(U); + const Number E = total_energy(U); + const Number rho_rho_e = rho * E - ScalarNumber(0.5) * m.norm_square(); + + const Number exponent = ScalarNumber(1.) / (gamma_min + Number(1.)); + return ryujin::pow(rho_rho_e, exponent); + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::surrogate_harten_entropy_derivative( + const state_type &U, const Number &eta, const Number &gamma_min) const + -> mixture_state_type + { + const Number rho = density(U); + const auto m = momentum(U); + const Number E = total_energy(U); + + const auto factor = + ryujin::pow(eta, -gamma_min) * Number(1.) / (gamma_min + Number(1.)); + + mixture_state_type result; + + result[0] = factor * E; + + for (unsigned int i = 0; i < dim; ++i) + result[1 + i] = -factor * m[i]; + + result[dim + 1] = factor * rho; + + return result; + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::surrogate_gamma(const state_type &U, + const Number &p) const + { + return Number(1.) + p / internal_energy(U); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::surrogate_pressure( + const state_type &U, const Number &gamma) const + { + return (gamma - Number(1.)) * internal_energy(U); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::surrogate_speed_of_sound( + const state_type &U, const Number &gamma) const + { + const auto rho = density(U); + const auto rho_e = internal_energy(U); + + auto radicand = gamma * (gamma - Number(1.)) * rho_e / rho; + return std::sqrt(positive_part(radicand)); + } + + + template + DEAL_II_ALWAYS_INLINE inline bool + HyperbolicSystemView::is_admissible(const state_type &U) const + { + const auto rho = density(U); + const auto rho_e = internal_energy(U); + + constexpr auto gt = dealii::SIMDComparison::greater_than; + using T = Number; + const auto test = + dealii::compare_and_apply_mask(rho, T(0.), T(0.), T(-1.)) + // + dealii::compare_and_apply_mask(rho_e, T(0.), T(0.), T(-1.)) + // + dealii::compare_and_apply_mask(pressure(U), T(0.), T(0.), T(-1.)); + +#ifdef DEBUG_OUTPUT + if (!(test == Number(0.))) { + std::cout << std::fixed << std::setprecision(16); + std::cout << "Bounds violation: Negative state [rho, e] detected!\n"; + std::cout << "\t\trho: " << rho << "\n"; + std::cout << "\t\tint: " << rho_e << "\n"; + std::cout << "\t\tp " << pressure(U) << "\n"; + std::cout << "\t\tent: " << specific_entropy(U) << std::endl; + } +#endif + + return (test == Number(0.)); + } + + + template + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::prescribe_riemann_characteristic( + const state_type &U, + const Number &p, + const state_type &U_bar, + const Number &p_bar, + const dealii::Tensor<1, dim, Number> &normal) const -> state_type + { + static_assert(component == 1 || component == 2, + "component has to be 1 or 2"); + + const auto b = Number(0.); + const auto pinf = Number(0.); + const auto q = Number(0.); + + /* + * The "four" Riemann characteristics are formed under the assumption + * of a locally isentropic flow. We further assume that the mass fractions + * are constant along the locally isentropic flow. For this, we first + * transform both states into {rho, vn, vperp, gamma, a}, where we use the + * NASG EOS interpolation to derive a surrogate gamma and speed of sound + * a. + * + * See, e.g., https://arxiv.org/pdf/2004.08750, "Compressible flow in + * a NOble-Abel Stiffened-Gas fluid", M. I. Radulescu. + */ + + const auto rho = density(U); + const auto Y_0 = partial_density_0(U) / rho; + const auto Y_1 = partial_density_1(U) / rho; + + const auto m = momentum(U); + const auto vn = m * normal / rho; + + const auto gamma = surrogate_gamma(U, p); + const auto a = surrogate_speed_of_sound(U, gamma); + const auto covolume = 1. - b * rho; + + const auto rho_bar = density(U_bar); + const auto m_bar = momentum(U_bar); + const auto vn_bar = m_bar * normal / rho_bar; + + const auto gamma_bar = surrogate_gamma(U_bar, p_bar); + const auto a_bar = surrogate_speed_of_sound(U_bar, gamma_bar); + const auto covolume_bar = 1. - b * rho_bar; + + /* + * Now compute the Riemann characteristics {R_1, R_2, vperp, s}: + * R_1 = v * n - 2 / (gamma - 1) * a * (1 - b * rho) + * R_2 = v * n + 2 / (gamma - 1) * a * (1 - b * rho) + * vperp + * S = (p + p_infty) / rho^gamma * (1 - b * rho)^gamma + * + * Here, we replace either R_1, or R_2 with values coming from U_bar: + */ + + const auto R_1 = + component == 1 ? vn_bar - 2. * a_bar / (gamma_bar - 1.) * covolume_bar + : vn - 2. * a / (gamma - 1.) * covolume; + + const auto R_2 = + component == 2 ? vn_bar + 2. * a_bar / (gamma_bar - 1.) * covolume_bar + : vn + 2. * a / (gamma - 1.) * covolume; + + /* + * Note that we are really hoping for the best here... We require + * that R_2 >= R_1 so that we can extract a valid sound speed... + */ + + Assert( + R_2 >= R_1, + dealii::ExcMessage("Encountered R_2 < R_1 in dynamic boundary value " + "enforcement. This implies that the interpolation " + "with Riemann characteristics failed.")); + + const auto vperp = m / rho - vn * normal; + + const auto S = (p + pinf) * ryujin::pow(Number(1.) / rho - b, gamma); + + /* + * Now, we have to reconstruct the actual conserved state U from the + * Riemann characteristics R_1, R_2, vperp, and s. We first set up + * {vn_new, vperp_new, a_new, S} and then solve for {rho_new, p_new} + * with the help of the NASG EOS surrogate formulas: + * + * S = (p + p_infty) / rho^gamma * (1 - b * rho)^gamma + * + * a^2 = gamma * (p + p_infty) / (rho * cov) + * + * This implies: + * + * a^2 / (gamma * S) = rho^{gamma - 1} / (1 - b * rho)^{1 + gamma} + */ + + const auto vn_new = Number(0.5) * (R_1 + R_2); + + /* + * Technically, we would need to solve for rho subject to a number of + * nonlinear relationships: + * + * a = (gamma - 1) * (R_2 - R_1) / (4. * (1 - b * rho)) + * + * a^2 / (gamma * S) = rho^{gamma - 1} / (1 - b * rho)^{gamma + 1} + * + * This seems to be a bit expensive for the fact that our dynamic + * boundary conditions are already terribly heuristic... + * + * So instead, we rewrite this system as: + * + * a * (1 - b * rho) = (gamma - 1) * (R_2 - R_1) / 4. + * + * a^2 / (gamma * S) (1 - b * rho)^2 + * = (rho / (1 - b * rho))^{gamma - 1} + * + * And compute the terms on the left simply with the old covolume and + * solving an easier easier nonlinear equation for the density. The + * resulting system reads: + * + * a = (gamma - 1) * (R_2 - R_1) / (4. * (1 - b * rho_old)) + * A = {a^2 / (gamma * S) (1 - b * rho_old)^{2 gamma}}^{1/(gamma - 1)} + * + * rho = A / (1 + b * A) + */ + + const auto a_new_square = + ryujin::fixed_power<2>((gamma - 1.) * (R_2 - R_1) / (4. * covolume)); + + auto term = ryujin::pow(a_new_square / (gamma * S), 1. / (gamma - 1.)); + if (b != ScalarNumber(0.)) { + term *= std::pow(covolume, 2. / (gamma - 1.)); + } + + const auto rho_new = term / (1. + b * term); + + const auto covolume_new = (1. - b * rho_new); + const auto p_new = a_new_square / gamma * rho_new * covolume_new - pinf; + + /* + * And translate back into conserved quantities: + */ + + const auto rho_e_new = + rho_new * q + (p_new + gamma * pinf) * covolume_new / (gamma - 1.); + + state_type U_new; + + U_new[0] = Y_0 * rho_new; + U_new[1] = Y_1 * rho_new; + for (unsigned int d = 0; d < dim; ++d) { + U_new[2 + d] = rho_new * (vn_new * normal + vperp)[d]; + } + + U_new[2 + dim] = + rho_e_new + 0.5 * rho_new * (vn_new * vn_new + vperp.norm_square()); + + return U_new; + } + + + template + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::apply_boundary_conditions( + dealii::types::boundary_id id, + const state_type &U, + const dealii::Tensor<1, dim, Number> &normal, + const Lambda &get_dirichlet_data) const -> state_type + { + state_type result = U; + + if (id == Boundary::dirichlet) { + result = get_dirichlet_data(); + + } else if (id == Boundary::dirichlet_momentum) { + /* Only enforce Dirichlet conditions on the momentum: */ + auto m_dirichlet = momentum(get_dirichlet_data()); + for (unsigned int k = 0; k < dim; ++k) + result[k + 2] = m_dirichlet[k]; + + } else if (id == Boundary::slip) { + auto m = momentum(U); + m -= 1. * (m * normal) * normal; + for (unsigned int k = 0; k < dim; ++k) + result[k + 2] = m[k]; + + } else if (id == Boundary::no_slip) { + for (unsigned int k = 0; k < dim; ++k) + result[k + 2] = Number(0.); + + } else if (id == Boundary::dynamic) { + /* + * On dynamic boundary conditions, we distinguish four cases: + * + * - supersonic inflow: prescribe full state + * - subsonic inflow: + * decompose into Riemann invariants and leave R_2 + * characteristic untouched. + * - supersonic outflow: do nothing + * - subsonic outflow: + * decompose into Riemann invariants and prescribe incoming + * R_1 characteristic. + */ + const auto m = momentum(U); + const auto rho = density(U); + + /* + * We do not have precomputed values available. Thus, simply query + * the pressure oracle and compute a surrogate speed of sound from + * there: + */ + const auto p = pressure(U); + const auto gamma = surrogate_gamma(U, p); + const auto a = surrogate_speed_of_sound(U, gamma); + const auto vn = m * normal / rho; + + /* Supersonic inflow: */ + if (vn < -a) { + result = get_dirichlet_data(); + } + + /* Subsonic inflow: */ + if (vn >= -a && vn <= 0.) { + const auto U_dirichlet = get_dirichlet_data(); + const auto p_dirichlet = pressure(U_dirichlet); + + result = prescribe_riemann_characteristic<2>( + U_dirichlet, p_dirichlet, U, p, normal); + } + + /* Subsonic outflow: */ + if (vn > 0. && vn <= a) { + const auto U_dirichlet = get_dirichlet_data(); + const auto p_dirichlet = pressure(U_dirichlet); + + result = prescribe_riemann_characteristic<1>( + U, p, U_dirichlet, p_dirichlet, normal); + } + /* Supersonic outflow: do nothing, i.e., keep U as is */ + + } else { + AssertThrow(false, dealii::ExcNotImplemented()); + } + + return result; + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::f(const state_type &U, + const Number &p) const -> flux_type + { + const auto rho_inverse = ScalarNumber(1.) / density(U); + const auto m = momentum(U); + const auto E = total_energy(U); + + flux_type result; + + for (unsigned int i = 0; i < dim; ++i) { + result[0][i] = U[0] * (m[i] * rho_inverse); + result[1][i] = U[1] * (m[i] * rho_inverse); + result[2 + i] = m * (m[i] * rho_inverse); + result[2 + i][i] += p; + } + result[dim + 2] = m * (rho_inverse * (E + p)); + + return result; + } + + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::mixture_f(const state_type &U, + const Number &p) const + -> mixture_flux_type + { + const auto rho_inverse = ScalarNumber(1.) / density(U); + const auto m = momentum(U); + const auto E = total_energy(U); + + mixture_flux_type result; + + result[0] = m; + for (unsigned int i = 0; i < dim; ++i) { + result[1 + i] = m * (m[i] * rho_inverse); + result[1 + i][i] += p; + } + result[dim + 1] = m * (rho_inverse * (E + p)); + + return result; + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::flux_contribution( + const PrecomputedVector &pv, + const InitialPrecomputedVector & /*piv*/, + const unsigned int i, + const state_type &U_i) const -> flux_contribution_type + { + const auto &[rho_i, p_i, surrogate_gamma_i, s_i, surrogate_harten_i] = + pv.template get_tensor(i); + return f(U_i, p_i); + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::flux_contribution( + const PrecomputedVector &pv, + const InitialPrecomputedVector & /*piv*/, + const unsigned int *js, + const state_type &U_j) const -> flux_contribution_type + { + const auto &[rho_j, p_j, surrogate_gamma_j, s_j, surrogate_harten_j] = + pv.template get_tensor(js); + return f(U_j, p_j); + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::flux_divergence( + const flux_contribution_type &flux_i, + const flux_contribution_type &flux_j, + const dealii::Tensor<1, dim, Number> &c_ij) const -> state_type + { + return -contract(add(flux_i, flux_j), c_ij); + } + + + template + template + auto HyperbolicSystemView::expand_state(const ST &state) const + -> state_type + { + using T = typename ST::value_type; + static_assert(std::is_same_v, "template mismatch"); + + constexpr auto dim2 = ST::dimension - 3; + static_assert(dim >= dim2, + "the space dimension of the argument state must not be " + "larger than the one of the target state"); + + state_type result; + result[0] = state[0]; + result[1] = state[1]; + result[dim + 2] = state[dim2 + 2]; + for (unsigned int i = 2; i < dim2 + 2; ++i) + result[i] = state[i]; + + return result; + } + + + template + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::from_initial_state( + const ST &initial_state) const -> state_type + { + const auto primitive_state = expand_state(initial_state); + return from_primitive_state(primitive_state); + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::from_primitive_state( + const state_type &primitive_state) const -> state_type + { + const auto rho = density(primitive_state); + /* extract velocity: */ + const auto u = /*SIC!*/ momentum(primitive_state); + /* extract specific internal energy: */ + const auto &p = primitive_state[dim + 2]; + + auto state = primitive_state; + /* Fix up momentum: */ + for (unsigned int i = 2; i < dim + 2; ++i) + state[i] *= rho; + + /* Compute total energy: */ + const Number gamma_bar = gamma_mixture(primitive_state); + state[dim + 2] = + p / (gamma_bar - Number(1.)) + ScalarNumber(0.5) * rho * u * u; + + return state; + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + HyperbolicSystemView::to_primitive_state( + const state_type &state) const -> state_type + { + const auto rho = density(state); + const auto rho_inverse = Number(1.) / rho; + const auto p = pressure(state); + + auto primitive_state = state; + /* Fix up velocity: */ + for (unsigned int i = 2; i < dim + 2; ++i) + primitive_state[i] *= rho_inverse; + /* Set pressure: */ + primitive_state[dim + 2] = p; + + return primitive_state; + } + + + template + template + auto HyperbolicSystemView::apply_galilei_transform( + const state_type &state, const Lambda &lambda) const -> state_type + { + auto result = state; + const auto M = lambda(momentum(state)); + for (unsigned int d = 0; d < dim; ++d) + result[2 + d] = M[d]; + return result; + } + } // namespace MultiSpeciesEuler +} // namespace ryujin diff --git a/source/multi_euler_ideal/indicator.h b/source/multi_euler_ideal/indicator.h new file mode 100644 index 00000000..a84cc50b --- /dev/null +++ b/source/multi_euler_ideal/indicator.h @@ -0,0 +1,276 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2020 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include "hyperbolic_system.h" + +#include +#include +#include + +#include +#include + + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + template + class IndicatorParameters : public dealii::ParameterAcceptor + { + public: + IndicatorParameters(const std::string &subsection = "/Indicator") + : ParameterAcceptor(subsection) + { + evc_factor_ = ScalarNumber(1.); + add_parameter("evc factor", + evc_factor_, + "Factor for scaling the entropy viscocity commuator"); + } + + ACCESSOR_READ_ONLY(evc_factor); + + private: + ScalarNumber evc_factor_; + }; + + + /** + * This class implements an indicator strategy used to form the + * preliminary high-order update. + * + * The indicator is an entropy-viscosity commutator as described in + * @cite GuermondEtAl2011 and @cite GuermondEtAl2018. For the + * multi-species Euler equations, we use a Harten-type entropy + * \f$\eta(\mathbf{u}) = (\rho \varepsilon)^{1/(\gamma_{\min}+1)}\f$ + * with a surrogate gamma \f$\gamma_{\min}\f$ computed as the minimum + * over the stencil (see Section 4 of [ClaytonDzanicTovar-2025]). + * + * We let \f$\eta'\f$ denote the derivative with respect to the + * "summed system" mixture state variables \f$[\rho, \mathbf{m}, E]\f$. + * We then compute a normalized entropy viscosity ratio + * \f$\alpha_i^n\f$ for the state \f$\mathbf{U}_i^n\f$ as follows: + * \f{align} + * \alpha_i^n\;=\;\frac{N_i^n}{D_i^n}, + * \quad + * N_i^n\;:=\;\left|a_i^n- \eta'(\mathbf{U}^n_i)\cdot\mathbf{b}_i^n + * +\frac{\eta(\mathbf{U}^n_i)}{\rho_i^n}\big(\mathbf{b}_i^n\big)_1 + * \right|, + * \quad + * D_i^n\;:=\;\left|a_i^n\right| + + * \sum_{k=1}^{d+1}\left|\big(\eta'(\mathbf{U}^n_i)\big)_k- + * \delta_{1k}\frac{\eta(\mathbf{U}^n_i)}{\rho_i^n}\right| + * \,\left|\big(\mathbf{b}_i^n\big)_k\right|, + * \f} + * where \f$\big(\,.\,\big)_k\f$ denotes the \f$k\f$-th component of a + * vector, \f$\delta_{ij}\f$ is Kronecker's delta, and where we have set + * \f{align} + * a_i^n \;:=\; + * \sum_{j\in\mathcal{I}_i}\left(\frac{\eta(\mathbf{U}_j^n)}{\rho_j^n} + * -\frac{\eta(\mathbf{U}_i^n)}{\rho_i^n}\right)\, + * \mathbf{m}_j^n\cdot\mathbf{c}_{ij}, + * \qquad + * \mathbf{b}_i^n \;:=\; + * \sum_{j\in\mathcal{I}_i}\left(\mathbf{f}^{\text{mix}}(\mathbf{U}_j^n)- + * \mathbf{f}^{\text{mix}}(\mathbf{U}_i^n)\right)\cdot\mathbf{c}_{ij}, + * \f} + * where \f$\mathbf{f}^{\text{mix}}\f$ is the "summed system" mixture + * flux operating on \f$[\rho, \mathbf{m}, E]\f$. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class Indicator + { + public: + /** + * @name Typedefs and constexpr constants + */ + //@{ + + using View = HyperbolicSystemView; + + using ScalarNumber = typename View::ScalarNumber; + + static constexpr auto problem_dimension = View::problem_dimension; + + using state_type = typename View::state_type; + + using flux_type = typename View::flux_type; + + using mixture_state_type = typename View::mixture_state_type; + + using mixture_flux_type = typename View::mixture_flux_type; + + using precomputed_type = typename View::precomputed_type; + + using PrecomputedVector = typename View::PrecomputedVector; + + using Parameters = IndicatorParameters; + + //@} + /** + * @name Stencil-based computation of indicators + * + * Intended usage: + * ``` + * Indicator indicator; + * for (unsigned int i = n_internal; i < n_owned; ++i) { + * // ... + * indicator.reset(i, U_i); + * for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) { + * // ... + * indicator.accumulate(js, U_j, c_ij); + * } + * indicator.alpha(hd_i); + * } + * ``` + */ + //@{ + + /** + * Constructor taking a HyperbolicSystem instance as argument + */ + Indicator(const HyperbolicSystem &hyperbolic_system, + const Parameters ¶meters, + const PrecomputedVector &precomputed_values) + : hyperbolic_system(hyperbolic_system) + , parameters(parameters) + , precomputed_values(precomputed_values) + { + } + + /** + * Reset temporary storage and initialize for a new row corresponding + * to state vector U_i. + */ + void reset(const unsigned int i, const state_type &U_i); + + /** + * When looping over the sparsity row, add the contribution associated + * with the neighboring state U_j. + */ + void accumulate(const unsigned int *js, + const state_type &U_j, + const dealii::Tensor<1, dim, Number> &c_ij); + + /** + * Return the computed alpha_i value. + */ + Number alpha(const Number h_i) const; + + //@} + + private: + /** + * @name + */ + //@{ + + const HyperbolicSystem &hyperbolic_system; + const Parameters ¶meters; + const PrecomputedVector &precomputed_values; + + Number rho_i_inverse = 0.; + Number eta_i = 0.; + mixture_flux_type f_i; + mixture_state_type d_eta_i; + Number gamma_min; + + Number left = 0.; + mixture_state_type right; + + //@} + }; + + + /* + * ------------------------------------------------------------------------- + * Inline definitions + * ------------------------------------------------------------------------- + */ + + + template + DEAL_II_ALWAYS_INLINE inline void + Indicator::reset(const unsigned int i, const state_type &U_i) + { + /* Entropy viscosity commutator: */ + + const auto view = hyperbolic_system.view(); + + const auto &[rho_i, p_i, gamma_min_i, s_i, new_eta_i] = + precomputed_values.template get_tensor(i); + + gamma_min = gamma_min_i; + + rho_i_inverse = Number(1.) / rho_i; + eta_i = new_eta_i; + + d_eta_i = view.surrogate_harten_entropy_derivative(U_i, eta_i, gamma_min); + d_eta_i[0] -= eta_i * rho_i_inverse; + + const auto surrogate_p_i = view.surrogate_pressure(U_i, gamma_min); + f_i = view.mixture_f(U_i, surrogate_p_i); + + left = 0.; + right = 0.; + } + + + template + DEAL_II_ALWAYS_INLINE inline void Indicator::accumulate( + const unsigned int *js, + const state_type &U_j, + const dealii::Tensor<1, dim, Number> &c_ij) + { + /* Entropy viscosity commutator: */ + const auto &[rho_j, p_j, gamma_min_j, s_j, dont_use_j] = + precomputed_values.template get_tensor(js); + + const auto view = hyperbolic_system.view(); + + const auto eta_j = view.surrogate_harten_entropy(U_j, gamma_min); + const auto rho_j_inverse = Number(1.) / rho_j; + + const auto m_j = view.momentum(U_j); + + const auto surrogate_p_j = view.surrogate_pressure(U_j, gamma_min); + const auto f_j = view.mixture_f(U_j, surrogate_p_j); + + const auto entropy_flux = + (eta_j * rho_j_inverse - eta_i * rho_i_inverse) * (m_j * c_ij); + + left += entropy_flux; + for (unsigned int k = 0; k < problem_dimension - 1; ++k) { + const auto component = (f_j[k] - f_i[k]) * c_ij; + right[k] += component; + } + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + Indicator::alpha(const Number hd_i) const + { + /* Entropy viscosity commutator: */ + + Number numerator = left; + Number denominator = std::abs(left); + for (unsigned int k = 0; k < problem_dimension - 1; ++k) { + numerator -= d_eta_i[k] * right[k]; + denominator += std::abs(d_eta_i[k] * right[k]); + } + + const auto quotient = safe_division(std::abs(numerator), + denominator + hd_i * std::abs(eta_i)); + + return std::min(Number(1.), parameters.evc_factor() * quotient); + } + } // namespace MultiSpeciesEuler +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_contrast.h b/source/multi_euler_ideal/initial_state_contrast.h new file mode 100644 index 00000000..71d2bdf5 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_contrast.h @@ -0,0 +1,96 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2022 - 2024 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + /** + * An initial state formed by a contrast of a given "left" and "right" + * primitive (initial) state. + * + * @note The @p t argument is ignored. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class Contrast : public InitialState + { + public: + using HyperbolicSystem = typename Description::HyperbolicSystem; + using View = + typename Description::template HyperbolicSystemView; + using state_type = typename View::state_type; + + Contrast(const HyperbolicSystem &hyperbolic_system, + const std::string subsection) + : InitialState("contrast", subsection) + , hyperbolic_system_(hyperbolic_system) + { + temp_left_[0] = 0.5; + temp_left_[1] = 1.4; + temp_left_[2] = 0.0; + temp_left_[3] = 1.0; + this->add_parameter( + "primitive state left", + temp_left_, + "Initial 1d primitive state (Y_0, rho, u, p) on the left"); + + temp_right_[0] = 0.5; + temp_right_[1] = 1.4; + temp_right_[2] = 0.0; + temp_right_[3] = 1.0; + this->add_parameter( + "primitive state right", + temp_right_, + "Initial 1d primitive state (Y_0, rho, u, p) on the right"); + + const auto convert_states = [&]() { + const auto view = hyperbolic_system_.template view(); + + const auto primitive_left_ = extend(temp_left_); + const auto primitive_right_ = extend(temp_right_); + + state_left_ = view.from_initial_state(primitive_left_); + state_right_ = view.from_initial_state(primitive_right_); + }; + this->parse_parameters_call_back.connect(convert_states); + convert_states(); + } + + state_type compute(const dealii::Point &point, Number /*t*/) final + { + return (point[0] > 0. ? state_right_ : state_left_); + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + dealii::Tensor<1, 4, Number> temp_left_; + dealii::Tensor<1, 4, Number> temp_right_; + + state_type state_left_; + state_type state_right_; + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> + extend(dealii::Tensor<1, 4, Number> &temp_in) const + { + dealii::Tensor<1, 4, Number> result; + result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; + result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; + + for (unsigned int i = 2; i < 4; ++i) + result[i] = temp_in[i]; + + return result; + } + }; + } // namespace MultiSpeciesEulerInitialStates +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_exact_riemann_solution.h b/source/multi_euler_ideal/initial_state_exact_riemann_solution.h new file mode 100644 index 00000000..298123d1 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_exact_riemann_solution.h @@ -0,0 +1,536 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2024 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include +#include + +#include + +// #define DEBUG_SOLUTION + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + /** + * The exact Riemann solution. + * + * This initial class computes the analytic solution for the + * compressible multi-species Euler equations with an ideal gas equation of + * state for each species. + * + * @note This class returns the analytic solution as a function of time + * @p t and position @p x. + * + * @ingroup MultiSpeciesEulerEquations + */ + + template + class ExactRiemannSolution : public InitialState + { + public: + //@{ + + using HyperbolicSystem = typename Description::HyperbolicSystem; + using View = + typename Description::template HyperbolicSystemView; + using state_type = typename View::state_type; + + using ScalarNumber = typename View::ScalarNumber; + + + ExactRiemannSolution(const HyperbolicSystem &hyperbolic_system, + const std::string subsection) + : InitialState("exact riemann solution", + subsection) + , hyperbolic_system_(hyperbolic_system) + { + primitive_left_[0] = 0.5; + primitive_left_[1] = 1.4; + primitive_left_[2] = 0.0; + primitive_left_[3] = 1.0; + this->add_parameter( + "primitive state left", + primitive_left_, + "Initial 1d primitive state (Y_0, rho, u, p) on the left"); + + primitive_right_[0] = 0.5; + primitive_right_[1] = 1.4; + primitive_right_[2] = 0.0; + primitive_right_[3] = 1.0; + this->add_parameter( + "primitive state right", + primitive_right_, + "Initial 1d primitive state (Y_0, rho, u, p) on the right"); + + // Prepare the Riemann data + const auto prepare_riemann_data = [&]() { + // Convert to alpha_0 rho_0, alpha_1 rho_1, u, p + primitive_left_ = transform(primitive_left_); + primitive_right_ = transform(primitive_right_); + + const Number p_L = primitive_left_[3]; + const Number p_R = primitive_right_[3]; + + p_star_ = compute_pstar(p_L, p_R, primitive_left_, primitive_right_); + + const Number u_L = primitive_left_[2]; + u_star_ = u_L - fZofP(p_star_, primitive_left_); + +#ifdef DEBUG_SOLUTION + const Number u_R = primitive_right_[2]; + std::cout << "left data = " << primitive_left_ + << "\nright data = " << primitive_right_ + << "\np_star = " << p_star_ + << "\nu_star = " << u_star_ + << "\nVerifying u_star = " + << u_R + fZofP(p_star_, primitive_right_) << std::endl; +#endif + + lambda_left_minus_ = lambda(p_star_, primitive_left_, -1.); + lambda_left_plus_ = + lambda_intermediate(p_star_, primitive_left_, -1.); + lambda_right_minus_ = + lambda_intermediate(p_star_, primitive_right_, 1.); + lambda_right_plus_ = lambda(p_star_, primitive_right_, 1.); + + +#ifdef DEBUG_SOLUTION + std::cout << "lambda_left_minus = " << lambda_left_minus_ + << "\nlambda_left_plus = " << lambda_left_plus_ + << "\nlambda_right_minus = " << lambda_right_minus_ + << "\nlambda_right_plus = " << lambda_right_plus_ + << std::endl; +#endif + }; + + this->parse_parameters_call_back.connect(prepare_riemann_data); + prepare_riemann_data(); + } + + + state_type compute(const dealii::Point &point, Number t) final + { + const auto view = hyperbolic_system_.template view(); + + const double &x = point[0]; + + const Number xi = x / t; + + dealii::Tensor<1, 4, Number> primitive_state; + + if (t < 1.e-14 && x < 0.) { + primitive_state = primitive_left_; +#ifdef DEBUG_SOLUTION + std::cout << "Left primitive state: " << primitive_state << std::endl; +#endif + + } else if (t < 1.e-14 && x > 0.) { + primitive_state = primitive_right_; +#ifdef DEBUG_SOLUTION + std::cout << "Right primitive state: " << primitive_state + << std::endl; +#endif + + } else if (xi < lambda_left_minus_) { + /* Left state: */ + primitive_state = primitive_left_; +#ifdef DEBUG_SOLUTION + std::cout << "Left primitive state: " << primitive_state << std::endl; +#endif + + } else if (xi < lambda_left_plus_) { + const auto c_LL = + expansion_solution(p_star_, xi, primitive_left_, -1.); + primitive_state = c_LL; +#ifdef DEBUG_SOLUTION + std::cout << "Left expansion state: " << primitive_state << std::endl; +#endif + + } else if (xi < u_star_) { + primitive_state = cstar_solution(p_star_, u_star_, primitive_left_); + + const Number p_L = primitive_left_[3]; + if (p_star_ < p_L) + primitive_state = expansion_solution( + p_star_, lambda_left_plus_, primitive_left_, -1.); +#ifdef DEBUG_SOLUTION + std::cout << "Left cstar state: " << primitive_state << std::endl; +#endif + + } else if (xi < lambda_right_minus_) { + primitive_state = cstar_solution(p_star_, u_star_, primitive_right_); + + const Number p_R = primitive_right_[3]; + if (p_star_ < p_R) + primitive_state = expansion_solution( + p_star_, lambda_right_minus_, primitive_right_, 1.); +#ifdef DEBUG_SOLUTION + std::cout << "Right cstar state: " << primitive_state << std::endl; +#endif + + } else if (xi < lambda_right_plus_) { + primitive_state = + expansion_solution(p_star_, xi, primitive_right_, 1.); +#ifdef DEBUG_SOLUTION + std::cout << "Right expansion state: " << primitive_state + << std::endl; +#endif + + } else { + /* Right state: */ + primitive_state = primitive_right_; +#ifdef DEBUG_SOLUTION + std::cout << "Right primitive state: " << primitive_state + << std::endl; +#endif + } + + return view.from_initial_state(primitive_state); + } + + private: + //@} + /** + * Runtime parameters + */ + //@{ + + dealii::Tensor<1, 4, Number> primitive_left_; + dealii::Tensor<1, 4, Number> primitive_right_; + + //@} + /** + * Private fields + */ + //@{ + + const HyperbolicSystem &hyperbolic_system_; + + Number p_star_; + Number u_star_; + Number lambda_left_minus_; + Number lambda_left_plus_; + Number lambda_right_minus_; + Number lambda_right_plus_; + + //@} + /** + * Internal helper functions for solving the exact Riemann problem + */ + //@{ + + // Used to transform data to (alpha_rho_0, alpha_rho_1, u, p) + dealii::Tensor<1, 4, Number> + transform(dealii::Tensor<1, 4, Number> &temp_in) const + { + dealii::Tensor<1, 4, Number> result; + result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; + result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; + + for (unsigned int i = 2; i < 4; ++i) + result[i] = temp_in[i]; + + return result; + } + + Number fZofP(const Number &p_in, + const dealii::Tensor<1, 4, Number> &data_in) const + { + // Need hyperbolic system for particular quantities + const auto view = hyperbolic_system_.template view(); + const auto conserved_in = view.from_initial_state(data_in); + + // Get left/right data + const Number rho_Z = view.density(conserved_in); + const Number p_Z = data_in[3]; + const Number gamma_ = view.gamma_mixture(conserved_in); + + const Number c_Z = std::sqrt(gamma_ * p_Z / rho_Z); + + const Number A_Z = 2. / (gamma_ + 1.) / rho_Z; + const Number B_Z = (gamma_ - 1.) / (gamma_ + 1.) * p_Z; + + const Number exp = 0.5 * (gamma_ - 1.) / gamma_; + Number left_brach = 2. * c_Z / (gamma_ - 1.); + left_brach *= (std::pow(p_in / p_Z, exp) - 1.); + + Number f_of_p = (p_in - p_Z) * std::sqrt(A_Z / (p_in + B_Z)); + + if (p_in <= p_Z) + f_of_p = left_brach; + + return f_of_p; + } + + + Number dfZofP(const Number &p_in, + const dealii::Tensor<1, 4, Number> &data_in) const + { + // Need hyperbolic system for particular quantities + const auto view = hyperbolic_system_.template view(); + const auto conserved_in = view.from_initial_state(data_in); + + // Get left/right data + const Number rho_Z = view.density(conserved_in); + const Number p_Z = data_in[3]; + const Number gamma_Z = view.gamma_mixture(conserved_in); + + const Number c_Z = std::sqrt(gamma_Z * p_Z / rho_Z); + + const Number A_Z = 2. / (gamma_Z + 1.) / rho_Z; + const Number B_Z = (gamma_Z - 1.) / (gamma_Z + 1.) * p_Z; + + Number exp = 0.5 * (gamma_Z - 1.) / gamma_Z; + Number left_brach = 2. * c_Z / (gamma_Z - 1.) * exp; + exp -= 1.; + + left_brach *= std::pow(p_in / p_Z, exp - 1.) / p_Z; + + Number right_branch = std::pow(A_Z / (p_in + B_Z), 1.5); + right_branch *= (2. * B_Z + p_in + p_Z) / (2. * A_Z); + + Number df_of_p = right_branch; + + if (p_in <= p_Z) + df_of_p = left_brach; + + return df_of_p; + } + + + Number dphi(const Number &p_in, + const dealii::Tensor<1, 4, Number> &data_left, + const dealii::Tensor<1, 4, Number> &data_right) const + { + return dfZofP(p_in, data_left) + dfZofP(p_in, data_right); + } + + + Number phi(const Number &p_in, + const dealii::Tensor<1, 4, Number> &data_left, + const dealii::Tensor<1, 4, Number> &data_right) const + { + const Number u_L = data_left[2]; + const Number u_R = data_right[2]; + + return fZofP(p_in, data_right) + fZofP(p_in, data_left) + u_R - u_L; + } + + + Number lambda(const Number &p_in, + const dealii::Tensor<1, 4, Number> &data_in, + const Number &sign) const + { + // Need hyperbolic system for particular quantities + const auto view = hyperbolic_system_.template view(); + const auto conserved_in = view.from_initial_state(data_in); + + // Get left/right data + const Number rho_Z = view.density(conserved_in); + const Number u_Z = data_in[2]; + const Number p_Z = data_in[3]; + const Number gamma_Z = view.gamma_mixture(conserved_in); + + const Number c_Z = std::sqrt(gamma_Z * p_Z / rho_Z); + + const Number radicand = + 1. + 0.5 * (gamma_Z + 1.) / gamma_Z * std::max(p_in / p_Z - 1., 0.); + + return u_Z + sign * c_Z * std::sqrt(radicand); + } + + + Number lambda_intermediate(const Number &p_in, + const dealii::Tensor<1, 4, Number> &data_in, + const Number &sign) const + { + // Need hyperbolic system for particular quantities + const auto view = hyperbolic_system_.template view(); + const auto conserved_in = view.from_initial_state(data_in); + + // Get left/right data + const Number rho_Z = view.density(conserved_in); + const Number u_Z = data_in[2]; + const Number p_Z = data_in[3]; + const Number gamma_Z = view.gamma_mixture(conserved_in); + + const Number c_Z = std::sqrt(gamma_Z * p_Z / rho_Z); + + const auto lambda_value = lambda(p_in, data_in, sign); + + const Number f_of_p = fZofP(p_in, data_in); + + const Number exp = 0.5 * (gamma_Z - 1.) / gamma_Z; + const Number expansion_speed = + u_Z + sign * (f_of_p + c_Z * std::pow(p_in / p_Z, exp)); + + Number result = lambda_value; + if (p_in < p_Z) + result = expansion_speed; + + return result; + } + + + dealii::Tensor<1, 4, Number> + cstar_solution(const Number &p_star, + const Number &u_star, + const dealii::Tensor<1, 4, Number> &data_in) const + { + // Need hyperbolic system for particular quantities + const auto view = hyperbolic_system_.template view(); + const auto conserved_in = view.from_initial_state(data_in); + + // Get left/right data + const Number rho_Z = view.density(conserved_in); + const Number p_Z = data_in[3]; + const Number gamma_Z = view.gamma_mixture(conserved_in); + + // Define rho_star + const Number p_ratio = p_star / p_Z; + const Number gamma_ratio = (gamma_Z - 1.) / (gamma_Z + 1.); + + const Number numerator = rho_Z * (p_ratio + gamma_ratio); + const Number denominator = gamma_ratio * p_ratio + 1.; + + Number rho_star = numerator / denominator; + + auto result = data_in; + result[0] = data_in[0] / rho_Z; + result[1] = rho_star; + result[2] = u_star; + result[3] = p_star; + + return transform(result); + } + + + dealii::Tensor<1, 4, Number> + expansion_solution(const Number & /*p_star*/, + const Number &xi, + const dealii::Tensor<1, 4, Number> &data_in, + const Number &sign) const + { + // Need hyperbolic system for particular quantities + const auto view = hyperbolic_system_.template view(); + const auto conserved_in = view.from_initial_state(data_in); + + // Get left/right data + const Number rho_Z = view.density(conserved_in); + const Number u_Z = data_in[2]; + const Number p_Z = data_in[3]; + const Number gamma_Z = view.gamma_mixture(conserved_in); + + const Number c_Z = std::sqrt(gamma_Z * p_Z / rho_Z); + + // Define rho_expansion + const Number gamma_ratio = (gamma_Z - 1.) / (gamma_Z + 1.); + + const Number first = 2. / (gamma_Z + 1.); + const Number second = gamma_ratio / c_Z * (u_Z - xi); + const Number exp = 2. / (gamma_Z - 1.); + + Number rho_expansion = rho_Z * std::pow(first - sign * second, exp); + + // Define p_expansion + const Number factor = p_Z / std::pow(rho_Z, gamma_Z); + const Number p_expansion = factor * std::pow(rho_expansion, gamma_Z); + + // Define u_expansion + const Number u_expansion = u_Z + sign * fZofP(p_expansion, data_in); + + auto result = data_in; + result[0] = data_in[0] / rho_Z; + result[1] = rho_expansion; + result[2] = u_expansion; + result[3] = p_expansion; + + return transform(result); + } + + + /** + * Compute pstar using the quadratic_newton_step() + */ + double compute_pstar(double p_1, + double p_2, + dealii::Tensor<1, 4, Number> data_1, + dealii::Tensor<1, 4, Number> data_2) + { + constexpr Number eps = std::numeric_limits::epsilon(); + + // Ensure that p_1 <= p_2 + + if (p_1 > p_2) { + std::swap(p_1, p_2); + std::swap(data_1, data_2); + } + +#ifdef DEBUG + const double phi_1 = phi(p_1, data_1, data_2); + const double phi_2 = phi(p_2, data_1, data_2); + Assert(phi_1 * phi_2 <= 0., + dealii::ExcMessage( + "Euler::ExactRiemannSolver: failed to compute p_star.")); +#endif + + // + // We simply compute the root of phi with a bisection method down + // to machine precision. This is not terribly efficient but luckily + // happens only once during initialization. + // + +#ifdef DEBUG_SOLUTION + std::cout << "Computing p_star with a bisection method." << std::endl; +#endif + + unsigned int iter = 0; + for (; iter < 200; ++iter) { + + // Check for convergence: + if (std::abs(p_2 - p_1) < 10. * eps * std::max(p_1, p_2)) { + break; + } + + const double phi_2 = phi(p_2, data_1, data_2); + +#ifdef DEBUG_SOLUTION + const double phi_1 = phi(p_1, data_1, data_2); + + std::cout << "\niter: " << iter << "\n"; + std::cout << "p_1: " << p_1 << "\n"; + std::cout << "p_2: " << p_2 << "\n"; + std::cout << "phi_1: " << phi_1 << "\n"; + std::cout << "phi_2: " << phi_2 << "\n"; +#endif + + const auto p_m = 0.5 * (p_2 + p_1); + const double phi_m = phi(p_m, data_1, data_2); + + if (phi_m * phi_2 >= 0.) { + p_2 = p_m; + } else { + p_1 = p_m; + } + } + +#ifdef DEBUG_SOLUTION + const double phi_2 = phi(p_2, data_1, data_2); + std::cout << "After " << iter << " iterations:" + << "\np_star = " << p_2 << "\nphi(p_star) = " << phi_2 + << "\n|p_2 - p_1| = " << std::abs(p_2 - p_1) << std::endl; +#endif + + return p_2; + } + + //@} + }; + } // namespace MultiSpeciesEulerInitialStates +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_function.h b/source/multi_euler_ideal/initial_state_function.h new file mode 100644 index 00000000..8d8bd802 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_function.h @@ -0,0 +1,149 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2024 by the ryujin authors +// + +#pragma once + +#include + +#include + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + /** + * Returns an initial state defined by a set of user specified function. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class Function : public InitialState + { + public: + using HyperbolicSystem = typename Description::HyperbolicSystem; + using View = + typename Description::template HyperbolicSystemView; + using state_type = typename View::state_type; + + Function(const HyperbolicSystem &hyperbolic_system, + const std::string subsection) + : InitialState("function", subsection) + , hyperbolic_system_(hyperbolic_system) + { + Y0_expression_ = "1.0"; + this->add_parameter("Y_0 expression", + Y0_expression_, + "A function expression describing the mass " + "fraction for species 1"); + + density_expression_ = "1.4"; + this->add_parameter("density expression", + density_expression_, + "A function expression describing the density"); + + velocity_x_expression_ = "3.0"; + this->add_parameter( + "velocity x expression", + velocity_x_expression_, + "A function expression describing the x-component of the velocity"); + + if constexpr (dim > 1) { + velocity_y_expression_ = "0.0"; + this->add_parameter("velocity y expression", + velocity_y_expression_, + "A function expression describing the " + "y-component of the velocity"); + } + + if constexpr (dim > 2) { + velocity_z_expression_ = "0.0"; + this->add_parameter("velocity z expression", + velocity_z_expression_, + "A function expression describing the " + "z-component of the velocity"); + } + + pressure_expression_ = "1.0"; + this->add_parameter("pressure expression", + pressure_expression_, + "A function expression describing the pressure"); + + /* + * Set up the muparser object with the final flux description from + * the parameter file: + */ + const auto set_up_muparser = [this] { + using FP = dealii::FunctionParser; + /* + * This variant of the constructor initializes the function + * parser with support for a time-dependent description involving + * a variable »t«: + */ + Y0_function_ = std::make_unique(Y0_expression_); + density_function_ = std::make_unique(density_expression_); + velocity_x_function_ = std::make_unique(velocity_x_expression_); + if constexpr (dim > 1) + velocity_y_function_ = std::make_unique(velocity_y_expression_); + if constexpr (dim > 2) + velocity_z_function_ = std::make_unique(velocity_z_expression_); + pressure_function_ = std::make_unique(pressure_expression_); + }; + + set_up_muparser(); + this->parse_parameters_call_back.connect(set_up_muparser); + } + + state_type compute(const dealii::Point &point, Number t) final + { + const auto view = hyperbolic_system_.template view(); + + state_type full_primitive_state; + + Y0_function_->set_time(t); + density_function_->set_time(t); + + full_primitive_state[0] = Y0_function_->value(point); + full_primitive_state[0] *= density_function_->value(point); + + full_primitive_state[1] = 1. - Y0_function_->value(point); + full_primitive_state[1] *= density_function_->value(point); + + velocity_x_function_->set_time(t); + full_primitive_state[2] = velocity_x_function_->value(point); + + if constexpr (dim > 1) { + velocity_y_function_->set_time(t); + full_primitive_state[3] = velocity_y_function_->value(point); + } + if constexpr (dim > 2) { + velocity_z_function_->set_time(t); + full_primitive_state[4] = velocity_z_function_->value(point); + } + + pressure_function_->set_time(t); + full_primitive_state[2 + dim] = pressure_function_->value(point); + + return view.from_primitive_state(full_primitive_state); + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + std::string Y0_expression_; + std::string density_expression_; + std::string velocity_x_expression_; + std::string velocity_y_expression_; + std::string velocity_z_expression_; + std::string pressure_expression_; + + std::unique_ptr> Y0_function_; + std::unique_ptr> density_function_; + std::unique_ptr> velocity_x_function_; + std::unique_ptr> velocity_y_function_; + std::unique_ptr> velocity_z_function_; + std::unique_ptr> pressure_function_; + }; + } // namespace MultiSpeciesEulerInitialStates +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_icf_like.h b/source/multi_euler_ideal/initial_state_icf_like.h new file mode 100644 index 00000000..3bf1338b --- /dev/null +++ b/source/multi_euler_ideal/initial_state_icf_like.h @@ -0,0 +1,207 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// [LANL Copyright Statement] +// Copyright (C) 2024 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include +#include + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + /** + * An initial state that simulates an "intertial confinement fusion" like + * problem. The set up consists of three regions: (i) a low density state + * inside a perturbed interface; (ii) a high density state outside the + * interface; (iii) an incoming shock wave characterized by its Mach number + * and the state outside the interface as well as starting location (given + * by a radius). The perturbed interface is characterized by the number of + * modes and an amplitude. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class ICFLike : public InitialState + { + public: + using HyperbolicSystem = typename Description::HyperbolicSystem; + using View = + typename Description::template HyperbolicSystemView; + using state_type = typename View::state_type; + + ICFLike(const HyperbolicSystem &hyperbolic_system, + const std::string subsection) + : InitialState("icf like", subsection) + , hyperbolic_system_(hyperbolic_system) + { + + temp_inside_[0] = 0.5; + temp_inside_[1] = 1.0; + temp_inside_[2] = 0.0; + temp_inside_[3] = 1.0; + this->add_parameter("primitive state inside", + temp_inside_, + "Initial primitive state (Y_0, rho, u, p) inside " + "perturbed interface"); + + temp_outside_[0] = 0.5; + temp_outside_[1] = 1.0; + temp_outside_[2] = 0.0; + temp_outside_[3] = 1.0; + this->add_parameter("primitive state outside", + temp_outside_, + "Initial primitive state (Y_0, rho, u, p) outside " + "perturbed interface"); + + interface_radius_ = 1.0; + this->add_parameter( + "interface radius", interface_radius_, "Radius of interface"); + + num_modes_ = 8.0; + this->add_parameter("number of modes", + num_modes_, + "Number of modes for pertburation of interface"); + + amplitude_ = 0.02; + this->add_parameter( + "amplitude", amplitude_, "Amplitude for interface pertburation"); + + mach_number_ = 3.0; + this->add_parameter( + "mach number", mach_number_, "Mach number of incoming shock front"); + + shock_radius_ = 1.2; + this->add_parameter("shock radius", + shock_radius_, + "Radial location of incoming shock front"); + + const auto convert_states = [&]() { + const auto prim_inside = extend(temp_inside_); + const auto prim_outside = extend(temp_outside_); + + const auto view = hyperbolic_system_.template view(); + state_inside_ = view.from_initial_state(prim_inside); + state_outside_ = view.from_initial_state(prim_outside); + }; + this->parse_parameters_call_back.connect(convert_states); + convert_states(); + }; + + state_type compute(const dealii::Point &point, Number) final + { + const auto view = hyperbolic_system_.template view(); + + /* Compute incoming shock state */ + state_type conserved_shock_state; + const auto r_hat = point / point.norm(); + { + const auto gamma_ = view.gamma_mixture(state_outside_); + + const auto &rho_R = temp_outside_[1]; + const auto &u_R = temp_outside_[2]; + const auto &p_R = temp_outside_[3]; + + const Number a_R = std::sqrt(gamma_ * p_R / rho_R); + const Number mach_R = u_R / a_R; + + auto S3_ = mach_number_ * a_R; + const Number delta_mach = mach_R - mach_number_; + + const Number rho_L = + rho_R * (gamma_ + Number(1.)) * delta_mach * delta_mach / + ((gamma_ - Number(1.)) * delta_mach * delta_mach + Number(2.)); + const Number u_L = + (Number(1.) - rho_R / rho_L) * S3_ + rho_R / rho_L * u_R; + const Number p_L = p_R * + (Number(2.) * gamma_ * delta_mach * delta_mach - + (gamma_ - Number(1.))) / + (gamma_ + Number(1.)); + + state_type primitive_shock_state; + const auto Y0_outside = temp_outside_[0]; + + primitive_shock_state[0] = Y0_outside * rho_L; + primitive_shock_state[1] = (1. - Y0_outside) * rho_L; + + for (unsigned int i = 0; i < dim; ++i) { + primitive_shock_state[i + 2] = 0.; + } + + if (point.norm() > 0.) { + for (unsigned int i = 0; i < dim; ++i) { + primitive_shock_state[i + 2] = -u_L * r_hat[i]; + } + } + primitive_shock_state[2 + dim] = p_L; + + conserved_shock_state = + view.from_initial_state(primitive_shock_state); + } + + /* Compute polar (and potentially azimuthal) angle */ + const auto x = point[0]; + const auto y = dim > 1 ? point[1] : 0.; + + const double theta = std::atan2(y, x); + + double phi = 0.; + if constexpr (dim == 3) + phi = std::atan2(point[2], std::sqrt(x * x + y * y)); + + /* Compute perturbation for interface */ + const auto omega = num_modes_; + const double perturbation = + amplitude_ * std::cos(omega * theta) * std::cos(omega * phi); + + /* Compute state depending on location */ + auto full_state = + (point.norm() > interface_radius_ + perturbation ? state_outside_ + : state_inside_); + + if (point.norm() > shock_radius_) { + full_state = conserved_shock_state; + } + + /* Set final state */ + return full_state; + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + Number gamma_; + + dealii::Tensor<1, 4, Number> temp_inside_; + dealii::Tensor<1, 4, Number> temp_outside_; + + state_type state_inside_; + state_type state_outside_; + + double interface_radius_; + double num_modes_; + double amplitude_; + double shock_radius_; + double mach_number_; + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> + extend(dealii::Tensor<1, 4, Number> &temp_in) const + { + dealii::Tensor<1, 4, Number> result; + result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; + result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; + + for (unsigned int i = 2; i < 4; ++i) + result[i] = temp_in[i]; + + return result; + } + }; + + + } // namespace MultiSpeciesEulerInitialStates +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_library.cc b/source/multi_euler_ideal/initial_state_library.cc new file mode 100644 index 00000000..88651650 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_library.cc @@ -0,0 +1,14 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#include "initial_state_library.template.h" + +namespace ryujin +{ + template class InitialStateLibrary; + template class InitialStateLibrary; + template class InitialStateLibrary; +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_library.template.h b/source/multi_euler_ideal/initial_state_library.template.h new file mode 100644 index 00000000..3e301341 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_library.template.h @@ -0,0 +1,41 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2024 by the ryujin authors +// + +#pragma once + +#include "description.h" + +#include "initial_state_library_multi_euler.h" +#include + +namespace ryujin +{ + using Description = MultiSpeciesEuler::Description; + + template + class InitialStateLibrary + { + public: + using HyperbolicSystem = typename Description::HyperbolicSystem; + using ParabolicSystem = typename Description::ParabolicSystem; + + using View = + typename Description::template HyperbolicSystemView; + + using initial_state_list_type = + std::set>>; + + static void + populate_initial_state_list(initial_state_list_type &initial_state_list, + const HyperbolicSystem &h, + const ParabolicSystem & /*p*/, + const std::string &s) + { + MultiSpeciesEulerInitialStates:: + populate_initial_state_list( + initial_state_list, h, s); + } + }; +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_library_multi_euler.h b/source/multi_euler_ideal/initial_state_library_multi_euler.h new file mode 100644 index 00000000..059e56b3 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_library_multi_euler.h @@ -0,0 +1,46 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 by the ryujin authors +// + +#pragma once + +#include + +#include "initial_state_contrast.h" +#include "initial_state_exact_riemann_solution.h" +#include "initial_state_function.h" +#include "initial_state_icf_like.h" +#include "initial_state_radial_contrast.h" +#include "initial_state_shock_bubble.h" +#include "initial_state_smooth_wave.h" +#include "initial_state_three_state_contrast.h" + + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + template + void populate_initial_state_list( + typename ryujin::InitialStateLibrary:: + initial_state_list_type &initial_state_list, + const typename Description::HyperbolicSystem &h, + const std::string &s) + { + auto add = [&](auto &&object) { + initial_state_list.emplace(std::move(object)); + }; + + add(std::make_unique>(h, s)); + add(std::make_unique>(h, + s)); + add(std::make_unique>(h, s)); + add(std::make_unique>(h, s)); + add(std::make_unique>(h, s)); + add(std::make_unique>(h, s)); + add(std::make_unique>(h, s)); + add(std::make_unique>(h, s)); + } + } // namespace MultiSpeciesEulerInitialStates +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_radial_contrast.h b/source/multi_euler_ideal/initial_state_radial_contrast.h new file mode 100644 index 00000000..c7bd2676 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_radial_contrast.h @@ -0,0 +1,127 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2022 - 2024 by the ryujin authors +// + +#pragma once + +#include + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + /** + * A modification of the "contrast" initial state. Now, we have an + * initial state formed by a contrast of a given "left" and "right" + * primitive state where the "left" state is "inside" the radius, R, + * and the "right" state is outside. + * + * @note The @p t argument is ignored. This class always returns the + * initial configuration. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class RadialContrast : public InitialState + { + public: + using HyperbolicSystem = typename Description::HyperbolicSystem; + using View = + typename Description::template HyperbolicSystemView; + using state_type = typename View::state_type; + + RadialContrast(const HyperbolicSystem &hyperbolic_system, + const std::string &subsection) + : InitialState("radial contrast", + subsection) + , hyperbolic_system_(hyperbolic_system) + { + use_radial_velocity_ = false; + this->add_parameter( + "use radial velocity", + use_radial_velocity_, + "If set to true, we transform a non-zero velocity into a radial " + "velocity with scaling 1 / r^(dim - 1)"); + + temp_left_[0] = 0.5; + temp_left_[1] = 1.4; + temp_left_[2] = 0.0; + temp_left_[3] = 1.0; + this->add_parameter( + "primitive state left", + temp_left_, + "Initial 1d primitive state (Y_0, rho, u, p) inside radial area"); + + temp_right_[0] = 0.5; + temp_right_[1] = 1.4; + temp_right_[2] = 0.0; + temp_right_[3] = 1.0; + this->add_parameter( + "primitive state right", + temp_right_, + "Initial 1d primitive state (Y_0, rho, u, p) outside radial area"); + + radius_ = 1.0; + this->add_parameter("radius", radius_, "Radius of radial area"); + + const auto convert_states = [&]() { + const auto view = hyperbolic_system_.template view(); + + const auto primitive_left_ = extend(temp_left_); + const auto primitive_right_ = extend(temp_right_); + + state_left_ = view.from_initial_state(primitive_left_); + state_right_ = view.from_initial_state(primitive_right_); + }; + + this->parse_parameters_call_back.connect(convert_states); + convert_states(); + }; + + state_type compute(const dealii::Point &point, Number /*t*/) final + { + if (point.norm() > 0. && use_radial_velocity_) { + state_left_[2] = temp_left_[2] * point[0] / point.norm(); + state_left_[3] = temp_left_[2] * point[1] / point.norm(); + + state_left_[2] = temp_right_[2] * point[0] / point.norm(); + state_left_[3] = temp_right_[2] * point[1] / point.norm(); + } + + auto final_state = + (point.norm() > radius_ ? state_right_ : state_left_); + + return final_state; + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + dealii::Tensor<1, 4, Number> temp_left_; + dealii::Tensor<1, 4, Number> temp_right_; + + state_type state_left_; + state_type state_right_; + + double radius_; + bool use_radial_velocity_; + + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> + extend(dealii::Tensor<1, 4, Number> &temp_in) const + { + dealii::Tensor<1, 4, Number> result; + result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; + result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; + + for (unsigned int i = 2; i < 4; ++i) + result[i] = temp_in[i]; + + return result; + } + }; + + + } // namespace MultiSpeciesEulerInitialStates +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_shock_bubble.h b/source/multi_euler_ideal/initial_state_shock_bubble.h new file mode 100644 index 00000000..37979ed7 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_shock_bubble.h @@ -0,0 +1,183 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2022 - 2024 by the ryujin authors +// // Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + /** + * A shock-bubble initial state. The shock HAS to be initiated for when + * either the mass fraction is 0 or 1 ONLY. The implementation assumes we + * are working in the direction (1, ..., 0). + * + * @note The @p t argument is ignored. This class always returns the + * initial configuration. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class ShockBubble : public InitialState + { + public: + using HyperbolicSystem = typename Description::HyperbolicSystem; + using View = + typename Description::template HyperbolicSystemView; + using state_type = typename View::state_type; + + ShockBubble(const HyperbolicSystem &hyperbolic_system, + const std::string subsection) + : InitialState("shock bubble", subsection) + , hyperbolic_system_(hyperbolic_system) + { + temp_inside_[0] = 0.5; + temp_inside_[1] = 1.0; + temp_inside_[2] = 0.0; + temp_inside_[3] = 1.0; + this->add_parameter( + "primitive state bubble", + temp_inside_, + "Initial 1d primitive state (Y_0, rho, u, p) inside bubble"); + + temp_ambient_[0] = 0.5; + temp_ambient_[1] = 1.0; + temp_ambient_[2] = 0.0; + temp_ambient_[3] = 1.0; + this->add_parameter( + "primitive state pre-shock", + temp_ambient_, + "Initial 1d primitive state (Y_0, rho, u, p) ambient pre-shock"); + + radius_ = 1.0; + this->add_parameter("radius", radius_, "Bubble radius"); + + bubble_center_[0] = 0.; + if constexpr (dim > 1) + bubble_center_[1] = 0.; + + this->add_parameter( + "bubble center", bubble_center_, "The dim sized bubble center"); + + shock_location_ = 0.0; + this->add_parameter("shock distance from bubble", + shock_location_, + "The distance from the bubble center. Negative " + "value to the left, positive value to the right. "); + + mach_number_ = 2.0; + this->add_parameter( + "mach number", mach_number_, "Mach number of shock front "); + + const auto convert_states = [&]() { + const auto view = hyperbolic_system_.template view(); + + primitive_inside_ = extend(temp_inside_); + primitive_ambient_ = extend(temp_ambient_); + + state_bubble_ = view.from_initial_state(primitive_inside_); + state_ambient_ = view.from_initial_state(primitive_ambient_); + }; + this->parse_parameters_call_back.connect(convert_states); + convert_states(); + } + + + state_type compute(const dealii::Point &point, Number t) final + { + const auto view = hyperbolic_system_.template view(); + + const auto &x = point[0]; + const auto xbar = point - bubble_center_; + const auto &xc = bubble_center_[0]; + auto sign = 1.; + if (xc < shock_location_) + sign = -1.; + + /* Compute shocked state outside */ + { + const auto gamma_ = view.gamma_mixture(state_ambient_); + + const auto &rho_R = temp_ambient_[1]; + const auto &u_R = temp_ambient_[2]; + const auto &p_R = temp_ambient_[3]; + /* a_R^2 = gamma * p / rho / (1 - b * rho) */ + const Number a_R = std::sqrt(gamma_ * p_R / rho_R); + const Number mach_R = u_R / a_R; + + auto S3_ = mach_number_ * a_R; + const Number delta_mach = mach_R - mach_number_; + + const Number rho_L = + rho_R * (gamma_ + Number(1.)) * delta_mach * delta_mach / + ((gamma_ - Number(1.)) * delta_mach * delta_mach + Number(2.)); + const Number u_L = + (Number(1.) - rho_R / rho_L) * S3_ + rho_R / rho_L * u_R; + const Number p_L = p_R * + (Number(2.) * gamma_ * delta_mach * delta_mach - + (gamma_ - Number(1.))) / + (gamma_ + Number(1.)); + + const auto Y0_outside = temp_ambient_[0]; + + primitive_shock_[0] = Y0_outside; + primitive_shock_[1] = rho_L; + primitive_shock_[2] = sign * u_L; + primitive_shock_[3] = p_L; + + primitive_shock_ = extend(primitive_shock_); + state_shock_ = view.from_initial_state(primitive_shock_); + } + + auto final_state = + (xbar.norm() <= radius_ ? state_bubble_ : state_ambient_); + + if (x <= xc + shock_location_) + final_state = state_shock_; + + /* Set final state */ + return final_state; + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + dealii::Tensor<1, 4, Number> temp_inside_; + dealii::Tensor<1, 4, Number> temp_ambient_; + + dealii::Tensor<1, 4, Number> primitive_inside_; + dealii::Tensor<1, 4, Number> primitive_ambient_; + dealii::Tensor<1, 4, Number> primitive_shock_; + + state_type state_bubble_; + state_type state_ambient_; + state_type state_shock_; + + dealii::Point bubble_center_; + + double radius_; + double mach_number_; + double shock_location_; + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> + extend(dealii::Tensor<1, 4, Number> &temp_in) const + { + dealii::Tensor<1, 4, Number> result; + result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; + result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; + + for (unsigned int i = 2; i < 4; ++i) + result[i] = temp_in[i]; + + return result; + } + }; + + + } // namespace MultiSpeciesEulerInitialStates +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_smooth_wave.h b/source/multi_euler_ideal/initial_state_smooth_wave.h new file mode 100644 index 00000000..626cd41a --- /dev/null +++ b/source/multi_euler_ideal/initial_state_smooth_wave.h @@ -0,0 +1,107 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2024 by the ryujin authors +// + +#pragma once + +#include +#include + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + /** + * This is a generalization of the "Smooth traveling wave" problem first + * proposed in Section 5.2 of @cite GuermondEtAl2018 for ideal gas EOS. + * The details for extending to multi-species Euler Equations with each + * species following an ideal gas is soon to be documented. + * + * @note This class returns the analytic solution as a function of time + * @p t and position @p x. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class SmoothWave : public InitialState + { + public: + using HyperbolicSystem = typename Description::HyperbolicSystem; + using View = + typename Description::template HyperbolicSystemView; + using state_type = typename View::state_type; + + using ScalarNumber = typename View::ScalarNumber; + + SmoothWave(const HyperbolicSystem &hyperbolic_system, + const std::string subsection) + : InitialState("smooth wave", subsection) + , hyperbolic_system_(hyperbolic_system) + { + mass_fraction_ = 0.5; + this->add_parameter("mass fraction", + mass_fraction_, + "The mass fraction of first species"); + + density_ref_ = 1.; + this->add_parameter("reference density", + density_ref_, + "The material reference density"); + + pressure_ref_ = 1.; + this->add_parameter("reference pressure", + pressure_ref_, + "The material reference pressure"); + + mach_number_ = 1.0; + this->add_parameter("mach number", + mach_number_, + "Mach number of traveling smooth wave"); + + /* These are the x_0 and x_1 parameters from references above. */ + left_ = 0.1; + right_ = 0.3; + } + + state_type compute(const dealii::Point &point, Number t) final + { + const auto view = hyperbolic_system_.template view(); + + auto point_bar = point; + point_bar[0] = point_bar[0] - mach_number_ * t; + const auto x = Number(point_bar[0]); + + const Number polynomial = Number(64) * + ryujin::fixed_power<3>(x - left_) * + ryujin::fixed_power<3>(right_ - x) / + ryujin::fixed_power<6>(right_ - left_); + + /* Define density profile */ + Number rho = density_ref_; + if (left_ <= point_bar[0] && point_bar[0] <= right_) + rho = density_ref_ + polynomial; + + state_type initial_state; + { + initial_state[0] = mass_fraction_ * rho; // alpha_1 rho_1 + initial_state[1] = (1. - mass_fraction_) * rho; // alpha_2 rho_2 + initial_state[2] = mach_number_; + initial_state[3] = pressure_ref_; + } + + return view.from_initial_state(initial_state); + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + Number mass_fraction_; + Number density_ref_; + Number pressure_ref_; + Number mach_number_; + Number left_; + Number right_; + }; + } // namespace MultiSpeciesEulerInitialStates +} // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_three_state_contrast.h b/source/multi_euler_ideal/initial_state_three_state_contrast.h new file mode 100644 index 00000000..088cc934 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_three_state_contrast.h @@ -0,0 +1,131 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2024 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + /** + * An initial state formed by two contrasts of given "left", "middle" + * and "right" primitive states. The user defines the lengths of the left + * and middle regions. The rest of the domain is populated with the right + * region. For single species, this initial state (default values) can be + * used to replicate the classical Woodward-Colella colliding blast wave + * problem described in @cite Woodward1984 + * + * @note The @p t argument is ignored. This class always returns the + * initial configuration. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class ThreeStateContrast : public InitialState + { + public: + using HyperbolicSystem = typename Description::HyperbolicSystem; + using View = + typename Description::template HyperbolicSystemView; + using state_type = typename View::state_type; + + ThreeStateContrast(const HyperbolicSystem &hyperbolic_system, + const std::string &subsection) + : InitialState("three state contrast", + subsection) + , hyperbolic_system_(hyperbolic_system) + { + + temp_left_[0] = 0.5; + temp_left_[1] = 1.; + temp_left_[2] = 0.; + temp_left_[3] = 1.e3; + this->add_parameter( + "primitive state left", + temp_left_, + "Initial 1d primitive state (Y_0, rho, u, p) on the left"); + + left_length_ = 0.1; + this->add_parameter("left region length", + left_length_, + "The length of the left region"); + + temp_middle_[0] = 0.5; + temp_middle_[1] = 1.; + temp_middle_[2] = 0.; + temp_middle_[3] = 1.e-2; + this->add_parameter( + "primitive state middle", + temp_middle_, + "Initial 1d primitive state (Y_0, rho, u, p) in the middle"); + + middle_length_ = 0.8; + this->add_parameter("middle region length", + middle_length_, + "The length of the middle region"); + + temp_right_[0] = 0.5; + temp_right_[1] = 1.; + temp_right_[2] = 0.; + temp_right_[3] = 1.e2; + this->add_parameter( + "primitive state right", + temp_right_, + "Initial 1d primitive state (Y_0, rho, u, p) on the right"); + + const auto convert_states = [&]() { + const auto view = hyperbolic_system_.template view(); + + const auto primitive_left_ = extend(temp_left_); + const auto primitive_middle_ = extend(temp_middle_); + const auto primitive_right_ = extend(temp_right_); + + state_left_ = view.from_initial_state(primitive_left_); + state_middle_ = view.from_initial_state(primitive_middle_); + state_right_ = view.from_initial_state(primitive_right_); + }; + this->parse_parameters_call_back.connect(convert_states); + convert_states(); + } + + state_type compute(const dealii::Point &point, Number /*t*/) final + { + return point[0] >= left_length_ + middle_length_ ? state_right_ + : point[0] >= left_length_ ? state_middle_ + : state_left_; + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + Number left_length_; + Number middle_length_; + + dealii::Tensor<1, 4, Number> temp_left_; + dealii::Tensor<1, 4, Number> temp_middle_; + dealii::Tensor<1, 4, Number> temp_right_; + + state_type state_left_; + state_type state_middle_; + state_type state_right_; + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> + extend(dealii::Tensor<1, 4, Number> &temp_in) const + { + dealii::Tensor<1, 4, Number> result; + result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; + result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; + + for (unsigned int i = 2; i < 4; ++i) + result[i] = temp_in[i]; + + return result; + } + }; + } // namespace MultiSpeciesEulerInitialStates +} // namespace ryujin diff --git a/source/multi_euler_ideal/instantiate.h b/source/multi_euler_ideal/instantiate.h new file mode 100644 index 00000000..b74f65df --- /dev/null +++ b/source/multi_euler_ideal/instantiate.h @@ -0,0 +1,18 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#ifndef RYUJIN_INCLUDE_INSTANTIATION_ONCE +#define RYUJIN_INCLUDE_INSTANTIATION_ONCE +#else +#error Instantiation files can only be included once. +#endif + +#include "description.h" + +namespace ryujin +{ + using MultiSpeciesEuler::Description; +} // namespace ryujin diff --git a/source/multi_euler_ideal/limiter.cc b/source/multi_euler_ideal/limiter.cc new file mode 100644 index 00000000..56968e89 --- /dev/null +++ b/source/multi_euler_ideal/limiter.cc @@ -0,0 +1,25 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2020 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#include "limiter.template.h" + +using namespace dealii; + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + /* instantiations */ + + template class Limiter<1, NUMBER>; + template class Limiter<2, NUMBER>; + template class Limiter<3, NUMBER>; + + template class Limiter<1, dealii::VectorizedArray>; + template class Limiter<2, dealii::VectorizedArray>; + template class Limiter<3, dealii::VectorizedArray>; + } // namespace MultiSpeciesEuler +} // namespace ryujin diff --git a/source/multi_euler_ideal/limiter.h b/source/multi_euler_ideal/limiter.h new file mode 100644 index 00000000..6bd20f9e --- /dev/null +++ b/source/multi_euler_ideal/limiter.h @@ -0,0 +1,549 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2020 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include "hyperbolic_system.h" + +#include +#include +#include +#include + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + template + class LimiterParameters : public dealii::ParameterAcceptor + { + public: + LimiterParameters(const std::string &subsection = "/Limiter") + : ParameterAcceptor(subsection) + { + iterations_ = 2; + add_parameter( + "iterations", iterations_, "Number of limiter iterations"); + + if constexpr (std::is_same::value) + newton_tolerance_ = 1.e-10; + else + newton_tolerance_ = 1.e-4; + add_parameter("newton tolerance", + newton_tolerance_, + "Tolerance for the quadratic newton stopping criterion"); + + relaxation_factor_ = ScalarNumber(1.); + add_parameter("relaxation factor", + relaxation_factor_, + "Factor for scaling the relaxation window with r_i = " + "factor * (m_i/|Omega|)^(1.5/d)."); + } + + ACCESSOR_READ_ONLY(iterations); + ACCESSOR_READ_ONLY(newton_tolerance); + ACCESSOR_READ_ONLY(relaxation_factor); + + private: + unsigned int iterations_; + ScalarNumber newton_tolerance_; + ScalarNumber relaxation_factor_; + }; + + + /** + * The convex limiter for the multi-species Euler equations. + * + * The class implements a convex limiting technique as described in + * @cite GuermondEtAl2018, @cite ryujin-2021-1, @cite ryujin-2023-4, + * and extended to multi-species flows in [ClaytonDzanicTovar-2025]. + * + * For the multi-species case, we enforce bounds on (see Section 5 in + * [ClaytonDzanicTovar-2025]): + * - Partial densities: \f$(\alpha_k \rho_k)_{\min} \le \alpha_k + * \rho_k \le (\alpha_k \rho_k)_{\max}\f$ for each species \f$k\f$ + * - Internal energy: \f$\varepsilon \ge \varepsilon_{\min}\f$ + * - Specific entropy: \f$s(\mathbf{u}) \ge s_{\min}\f$ + * + * Given a computed set of bounds and an update direction + * \f$\mathbf{P}_{ij}\f$ one determines a candidate \f$\tilde l_{ij}\f$ + * by computing: + * \f{align} + * \tilde l_{ij} = \max_{l\,\in\,[0,1]} + * \,\Big\{(\alpha_k \rho_k)_{\min}\,\le\,\alpha_k \rho_k\,(\mathbf{U}_i + * + \tilde l_{ij}\mathbf{P}_{ij}) \,\le\,(\alpha_k \rho_k)_{\max},\; + * s_{\min}\,\le\,s\,(\mathbf{U}_{i}+\tilde l_{ij}\mathbf{P}_{ij})\Big\}. + * \f} + * + * Algorithmically this is accomplished as follows: Given an initial + * interval \f$[t_L,t_R]\f$, where \f$t_L\f$ is a good state, we first + * make the interval smaller ensuring the bounds on the partial + * densities are fulfilled. We then perform a quadratic Newton + * iteration (updating \f$[t_L,t_R]\f$) solving for the root of the + * 3-convex function: + * \f{align} + * \Psi(\mathbf{U})\;=\;\rho^{\bar{\gamma}+1}(\mathbf{U})\, + * \big(s(\mathbf{U})-s_{\min}\big), + * \f} + * where we use a surrogate gamma \f$\bar{\gamma} = \gamma_{\min}\f$ + * to ensure convexity of \f$\Psi\f$. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class Limiter + { + public: + /** + * @name Typedefs and constexpr constants + */ + //@{ + + using View = HyperbolicSystemView; + + using ScalarNumber = typename View::ScalarNumber; + + static constexpr auto problem_dimension = View::problem_dimension; + + using state_type = typename View::state_type; + + using flux_contribution_type = typename View::flux_contribution_type; + + using precomputed_type = typename View::precomputed_type; + + using PrecomputedVector = typename View::PrecomputedVector; + + using Parameters = LimiterParameters; + + //@} + /** + * @name Computation and manipulation of bounds + */ + // + //@{ + /** + * The number of stored entries in the bounds array. + */ + static constexpr unsigned int n_bounds = 6; + + /** + * Array type used to store accumulated bounds. + */ + using Bounds = std::array; + + /** + * Constructor taking a HyperbolicSystem instance as argument + */ + Limiter(const HyperbolicSystem &hyperbolic_system, + const Parameters ¶meters, + const PrecomputedVector &precomputed_values) + : hyperbolic_system(hyperbolic_system) + , parameters(parameters) + , precomputed_values(precomputed_values) + { + } + + /** + * Given a state @p U_i and an index @p i return "strict" bounds, + * i.e., a minimal convex set containing the state. + */ + Bounds projection_bounds_from_state(const unsigned int i, + const state_type &U_i) const; + + /** + * Given two bounds bounds_left, bounds_right, this function computes + * a larger, combined set of bounds that this is a (convex) superset + * of the two. + */ + Bounds combine_bounds(const Bounds &bounds_left, + const Bounds &bounds_right) const; + + /** + * This function applies a relaxation to a given a (strict) bound @p + * bounds using a non dimensionalized measure @p hd (that should + * scale as $h^d$, where $h$ is the local mesh size). This is done + * for the case of the Euler equations by multiplying maximum bounds + * with $(1+r)$ and minimum bounds with $(1-r)$, while ensuring that + * the bounds still describe an admissible state. + */ + Bounds fully_relax_bounds(const Bounds &bounds, const Number &hd) const; + + //@} + /** + * @name Stencil-based computation of bounds + * + * Intended usage: + * ``` + * Limiter limiter; + * for (unsigned int i = n_internal; i < n_owned; ++i) { + * // ... + * limiter.reset(i, U_i, flux_i); + * for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx) { + * // ... + * limiter.accumulate(js, U_j, flux_j, scaled_c_ij, affine_shift); + * } + * limiter.bounds(hd_i); + * } + * ``` + */ + //@{ + + /** + * Reset temporary storage + */ + void reset(const unsigned int i, + const state_type &U_i, + const flux_contribution_type &flux_i); + + /** + * When looping over the sparsity row, add the contribution associated + * with the neighboring state U_j. + */ + void accumulate(const unsigned int *js, + const state_type &U_j, + const flux_contribution_type &flux_j, + const dealii::Tensor<1, dim, Number> &scaled_c_ij, + const state_type &affine_shift); + + /** + * Return the computed bounds (with relaxation applied). + */ + Bounds bounds(const Number hd_i) const; + + //*} + /** @name Convex limiter */ + //@{ + + /** + * Given a state \f$\mathbf U\f$ and an update \f$\mathbf P\f$ this + * function computes and returns the maximal coefficient \f$t\f$, + * obeying \f$t_{\text{min}} < t < t_{\text{max}}\f$, such that the + * selected local minimum principles are obeyed. + * + * The returned boolean is set to true if the original low-order + * update was within bounds. + * + * @note If the debug option `DEBUG_EXPENSIVE_BOUNDS_CHECK` is set to + * true, then the boolean is set to true if the low-order and the + * resulting high-order update are within bounds. The latter might be + * violated due to round-off errors when computing the limiter + * bounds. + */ + std::tuple limit(const Bounds &bounds, + const state_type &U, + const state_type &P, + const Number t_min = Number(0.), + const Number t_max = Number(1.)) const; + + private: + //@} + /** @name Arguments and internal fields */ + //@{ + + const HyperbolicSystem &hyperbolic_system; + const Parameters ¶meters; + const PrecomputedVector &precomputed_values; + + state_type U_i; + flux_contribution_type flux_i; + + Bounds bounds_; + + Number alpha_rho_0_relaxation_numerator; + Number alpha_rho_1_relaxation_numerator; + Number rho_relaxation_denominator; + Number rho_e_interp_max; + Number s_interp_max; + + Number cv_i; + + //@} + }; + + + /* + * ------------------------------------------------------------------------- + * Inline definitions + * ------------------------------------------------------------------------- + */ + + + template + DEAL_II_ALWAYS_INLINE inline auto + Limiter::projection_bounds_from_state( + const unsigned int i, const state_type &U_i) const -> Bounds + { + const auto view = hyperbolic_system.view(); + + const auto &[rho_i, p_i, gamma_min_i, s_i, harten_i] = + precomputed_values.template get_tensor(i); + + const auto alpha_rho_0 = view.partial_density_0(U_i); + const auto alpha_rho_1 = view.partial_density_1(U_i); + + const auto rho_e = view.internal_energy(U_i); + + return {/*alpha_rho_0_min*/ + alpha_rho_0, + /*alpha_rho_0_max*/ alpha_rho_0, + /*alpha_rho_1_min*/ alpha_rho_1, + /*alpha_rho_1_min*/ alpha_rho_1, + /*rho_e*/ rho_e, + /*s_min*/ s_i}; + } + + + template + DEAL_II_ALWAYS_INLINE inline auto Limiter::combine_bounds( + const Bounds &bounds_left, const Bounds &bounds_right) const -> Bounds + { + const auto &[alpha_rho_0_min_l, + alpha_rho_0_max_l, + alpha_rho_1_min_l, + alpha_rho_1_max_l, + rho_e_min_l, + s_min_l] = bounds_left; + const auto &[alpha_rho_0_min_r, + alpha_rho_0_max_r, + alpha_rho_1_min_r, + alpha_rho_1_max_r, + rho_e_min_r, + s_min_r] = bounds_right; + + return {std::min(alpha_rho_0_min_l, alpha_rho_0_min_r), + std::max(alpha_rho_0_max_l, alpha_rho_0_max_r), + std::min(alpha_rho_1_min_l, alpha_rho_1_min_r), + std::max(alpha_rho_1_max_l, alpha_rho_1_max_r), + std::min(rho_e_min_l, rho_e_min_r), + std::min(s_min_l, s_min_r)}; + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + Limiter::fully_relax_bounds(const Bounds &bounds, + const Number &hd) const -> Bounds + { + auto relaxed_bounds = bounds; + auto &[alpha_rho_0_min_relaxed, + alpha_rho_0_max_relaxed, + alpha_rho_1_min_relaxed, + alpha_rho_1_max_relaxed, + rho_e_min_relaxed, + s_min_relaxed] = relaxed_bounds; + + /* Use r = factor * (m_i / |Omega|) ^ (1.5 / d): */ + + Number r = std::sqrt(hd); // in 3D: ^ 3/6 + if constexpr (dim == 2) // + r = dealii::Utilities::fixed_power<3>(std::sqrt(r)); // in 2D: ^ 3/4 + else if constexpr (dim == 1) // + r = dealii::Utilities::fixed_power<3>(r); // in 1D: ^ 3/2 + r *= parameters.relaxation_factor(); + + constexpr ScalarNumber eps = std::numeric_limits::epsilon(); + + alpha_rho_0_min_relaxed *= std::max(Number(1.) - r, Number(eps)); + alpha_rho_0_max_relaxed *= (Number(1.) + r); + + alpha_rho_1_min_relaxed *= std::max(Number(1.) - r, Number(eps)); + alpha_rho_1_max_relaxed *= (Number(1.) + r); + + rho_e_min_relaxed *= std::max(Number(1.) - r, Number(eps)); + + const Number scaled_ent = std::exp(s_min_relaxed / cv_i); + s_min_relaxed = cv_i * std::log((Number(1.) - r) * scaled_ent); + + return relaxed_bounds; + } + + + template + DEAL_II_ALWAYS_INLINE inline void + Limiter::reset(const unsigned int i, + const state_type &new_U_i, + const flux_contribution_type &new_flux_i) + { + const auto view = hyperbolic_system.view(); + + U_i = new_U_i; + flux_i = new_flux_i; + + cv_i = view.cv_mixture(new_U_i); + + /* Bounds: */ + + auto &[alpha_rho_0_min, + alpha_rho_0_max, + alpha_rho_1_min, + alpha_rho_1_max, + rho_e_min, + s_min] = bounds_; + + alpha_rho_0_min = Number(std::numeric_limits::max()); + alpha_rho_0_max = Number(0.); + + alpha_rho_1_min = Number(std::numeric_limits::max()); + alpha_rho_1_max = Number(0.); + + rho_e_min = Number(std::numeric_limits::max()); + s_min = Number(std::numeric_limits::max()); + + + /* Relaxation: */ + + alpha_rho_0_relaxation_numerator = Number(0.); + alpha_rho_1_relaxation_numerator = Number(0.); + rho_relaxation_denominator = Number(0.); + + rho_e_interp_max = Number(0.); + s_interp_max = Number(0.); + } + + + template + DEAL_II_ALWAYS_INLINE inline void Limiter::accumulate( + const unsigned int *js, + const state_type &U_j, + const flux_contribution_type &flux_j, + const dealii::Tensor<1, dim, Number> &scaled_c_ij, + const state_type &affine_shift) + { + // TODO: Currently we only apply the affine_shift to U_ij_bar (which + // then enters all bounds), but we do not modify s_interp and + // rho_relaxation. When actually adding a source term to the Euler + // equations verify that this does the right thing. + Assert(std::max(affine_shift.norm(), Number(0.)) == Number(0.), + dealii::ExcNotImplemented()); + + const auto view = hyperbolic_system.view(); + + /* Bounds: */ + auto &[alpha_rho_0_min, + alpha_rho_0_max, + alpha_rho_1_min, + alpha_rho_1_max, + rho_e_min, + s_min] = bounds_; + + const auto alpha_rho_0_i = view.partial_density_0(U_i); + const auto alpha_rho_0_j = view.partial_density_0(U_j); + + const auto alpha_rho_1_i = view.partial_density_1(U_i); + const auto alpha_rho_1_j = view.partial_density_1(U_j); + + const auto rho_e_j = view.internal_energy(U_j); + + const auto [rho_j, p_j, gamma_min_j, s_j, harten_j] = + precomputed_values.template get_tensor(js); + + /* bar state shifted by an affine shift: */ + const auto U_ij_bar = + ScalarNumber(0.5) * (U_i + U_j) - + ScalarNumber(0.5) * contract(add(flux_j, -flux_i), scaled_c_ij) + + affine_shift; + + const auto alpha_rho_0_ij_bar = view.partial_density_0(U_ij_bar); + const auto alpha_rho_1_ij_bar = view.partial_density_1(U_ij_bar); + + const auto rho_e_ij_bar = view.internal_energy(U_ij_bar); + + /* Density bounds: */ + + alpha_rho_0_min = std::min(alpha_rho_0_min, alpha_rho_0_ij_bar); + alpha_rho_0_max = std::max(alpha_rho_0_max, alpha_rho_0_ij_bar); + + alpha_rho_1_min = std::min(alpha_rho_1_min, alpha_rho_1_ij_bar); + alpha_rho_1_max = std::max(alpha_rho_1_max, alpha_rho_1_ij_bar); + + /* Density relaxation: */ + + /* Use a uniform weight. */ + const auto beta_ij = Number(1.); + alpha_rho_0_relaxation_numerator += + beta_ij * (alpha_rho_0_i + alpha_rho_0_j); + alpha_rho_1_relaxation_numerator += + beta_ij * (alpha_rho_1_i + alpha_rho_1_j); + rho_relaxation_denominator += std::abs(beta_ij); + + /* Internal energy bounds and relaxation */ + rho_e_min = std::min(rho_e_min, rho_e_j); + rho_e_min = std::min(rho_e_min, rho_e_ij_bar); + + const auto rho_e_interp = + view.internal_energy((U_i + U_j) * ScalarNumber(.5)); + rho_e_interp_max = std::max(rho_e_interp_max, rho_e_interp); + + /* Entropy bounds and relaxation: */ + + const auto s_ij_bar = view.specific_entropy(U_ij_bar); + + const Number s_interp = + view.specific_entropy((U_i + U_j) * ScalarNumber(.5)); + + s_min = std::min(s_min, s_j); + s_min = std::min(s_min, s_ij_bar); + s_interp_max = std::max(s_interp_max, s_interp); + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + Limiter::bounds(const Number hd_i) const -> Bounds + { + const auto &[alpha_rho_0_min, + alpha_rho_0_max, + alpha_rho_1_min, + alpha_rho_1_max, + rho_e_min, + s_min] = bounds_; + + auto relaxed_bounds = fully_relax_bounds(bounds_, hd_i); + auto &[alpha_rho_0_min_relaxed, + alpha_rho_0_max_relaxed, + alpha_rho_1_min_relaxed, + alpha_rho_1_max_relaxed, + rho_e_min_relaxed, + s_min_relaxed] = relaxed_bounds; + + /* Apply a stricter window: */ + + constexpr ScalarNumber eps = std::numeric_limits::epsilon(); + + const auto alpha_rho_0_relaxation = + ScalarNumber(2. * parameters.relaxation_factor()) * + std::abs(alpha_rho_0_relaxation_numerator) / + (std::abs(rho_relaxation_denominator) + Number(eps)); + + const auto alpha_rho_1_relaxation = + ScalarNumber(2. * parameters.relaxation_factor()) * + std::abs(alpha_rho_1_relaxation_numerator) / + (std::abs(rho_relaxation_denominator) + Number(eps)); + + const auto int_relaxation = + parameters.relaxation_factor() * (rho_e_interp_max - rho_e_min); + + const auto entropy_relaxation = + parameters.relaxation_factor() * (s_interp_max - s_min); + + alpha_rho_0_min_relaxed = std::max( + alpha_rho_0_min_relaxed, alpha_rho_0_min - alpha_rho_0_relaxation); + alpha_rho_0_max_relaxed = std::min( + alpha_rho_0_max_relaxed, alpha_rho_0_max + alpha_rho_0_relaxation); + + alpha_rho_1_min_relaxed = std::max( + alpha_rho_1_min_relaxed, alpha_rho_1_min - alpha_rho_1_relaxation); + alpha_rho_1_max_relaxed = std::min( + alpha_rho_1_max_relaxed, alpha_rho_1_max + alpha_rho_1_relaxation); + + rho_e_min_relaxed = + std::max(rho_e_min_relaxed, rho_e_min - int_relaxation); + s_min_relaxed = std::max(s_min_relaxed, s_min - entropy_relaxation); + + return relaxed_bounds; + } + } // namespace MultiSpeciesEuler +} // namespace ryujin diff --git a/source/multi_euler_ideal/limiter.template.h b/source/multi_euler_ideal/limiter.template.h new file mode 100644 index 00000000..f502faad --- /dev/null +++ b/source/multi_euler_ideal/limiter.template.h @@ -0,0 +1,502 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2020 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include "limiter.h" +// #define DEBUG_OUTPUT_LIMITER + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + template + std::tuple + Limiter::limit(const Bounds &bounds, + const state_type &U, + const state_type &P, + const Number t_min /* = Number(0.) */, + const Number t_max /* = Number(1.) */) const + { + const auto view = hyperbolic_system.view(); + + bool success = true; + Number t_r = t_max; + + constexpr ScalarNumber min = std::numeric_limits::min(); + constexpr ScalarNumber eps = std::numeric_limits::epsilon(); + const auto small = view.vacuum_state_relaxation_small(); + const auto large = view.vacuum_state_relaxation_large(); + const ScalarNumber relax_small = ScalarNumber(1. + small * eps); + const ScalarNumber relax = ScalarNumber(1. + large * eps); + + /* + * First limit the fractional density alpha_rho_0. + * + * See [Guermond, Nazarov, Popov, Thomas] (4.8): + */ + + { + /* Be careful with the variable names. Not really rho here but + * alpha_0 * rho_0. + */ + const auto &rho_U = view.partial_density_0(U); + const auto &rho_P = view.partial_density_0(P); + + const auto &rho_min = std::get<0>(bounds); + const auto &rho_max = std::get<1>(bounds); + + /* + * Verify that rho_U is within bounds. This property might be + * violated for relative CFL numbers larger than 1. + */ + const auto test_min = view.filter_vacuum_density( + std::max(Number(0.), rho_U - relax * rho_max)); + const auto test_max = view.filter_vacuum_density( + std::max(Number(0.), rho_min - relax * rho_U)); + if (!(test_min == Number(0.) && test_max == Number(0.))) { +#ifdef DEBUG_OUTPUT + std::cout << std::fixed << std::setprecision(16); + std::cout + << "Bounds violation: low-order [species 0] density (critical)!" + << "\n\t\trho min: " << rho_min + << "\n\t\trho min (delta): " << negative_part(rho_U - rho_min) + << "\n\t\trho: " << rho_U + << "\n\t\trho max (delta): " << positive_part(rho_U - rho_max) + << "\n\t\trho max: " << rho_max << "\n" + << std::endl; +#endif + success = false; + } + + const Number denominator = + ScalarNumber(1.) / (std::abs(rho_P) + eps * rho_max + min); + + t_r = dealii::compare_and_apply_mask( + rho_max, + rho_U + t_r * rho_P, + /* + * rho_P is positive. + * + * Note: Do not take an absolute value here. If we are out of + * bounds we have to ensure that t_r is set to t_min. + */ + (rho_max - rho_U) * denominator, + t_r); + + t_r = dealii::compare_and_apply_mask( + rho_U + t_r * rho_P, + rho_min, + /* + * rho_P is negative. + * + * Note: Do not take an absolute value here. If we are out of + * bounds we have to ensure that t_r is set to t_min. + */ + (rho_U - rho_min) * denominator, + t_r); + + /* + * Ensure that t_min <= t <= t_max. This might not be the case if + * rho_U is outside the interval [rho_min, rho_max]. Furthermore, + * the quotient we take above is prone to numerical cancellation in + * particular in the second pass of the limiter when rho_P might be + * small. + */ + t_r = std::min(t_r, t_max); + t_r = std::max(t_r, t_min); + +#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK + /* + * Verify that the new state is within bounds: + */ + const auto rho_new = view.partial_density_0(U + t_r * P); + const auto test_new_min = view.filter_vacuum_density( + std::max(Number(0.), rho_new - relax * rho_max)); + const auto test_new_max = view.filter_vacuum_density( + std::max(Number(0.), rho_min - relax * rho_new)); + if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) { +#ifdef DEBUG_OUTPUT + std::cout << std::fixed << std::setprecision(16); + std::cout << "Bounds violation: high-order [species 0] density!" + << "\n\t\trho min: " << rho_min + << "\n\t\trho min (delta): " + << negative_part(rho_new - rho_min) + << "\n\t\trho: " << rho_new + << "\n\t\trho max (delta): " + << positive_part(rho_new - rho_max) + << "\n\t\trho max: " << rho_max << "\n" + << std::endl; +#endif + success = false; + } +#endif + } + + /* + * Then limit the fractional density alpha_rho_1. + * + * See [Guermond, Nazarov, Popov, Thomas] (4.8): + */ + { + /* Be careful with the variable names. Not really rho here but + * alpha_1 * rho_1. + */ + auto t_q = t_r; + + const auto &rho_U = view.partial_density_1(U); + const auto &rho_P = view.partial_density_1(P); + + const auto &rho_min = std::get<2>(bounds); + const auto &rho_max = std::get<3>(bounds); + + /* + * Verify that rho_U is within bounds. This property might be + * violated for relative CFL numbers larger than 1. + */ + const auto test_min = view.filter_vacuum_density( + std::max(Number(0.), rho_U - relax * rho_max)); + const auto test_max = view.filter_vacuum_density( + std::max(Number(0.), rho_min - relax * rho_U)); + if (!(test_min == Number(0.) && test_max == Number(0.))) { +#ifdef DEBUG_OUTPUT + std::cout << std::fixed << std::setprecision(16); + std::cout + << "Bounds violation: low-order [species 1] density (critical)!" + << "\n\t\trho min: " << rho_min + << "\n\t\trho min (delta): " << negative_part(rho_U - rho_min) + << "\n\t\trho: " << rho_U + << "\n\t\trho max (delta): " << positive_part(rho_U - rho_max) + << "\n\t\trho max: " << rho_max << "\n" + << std::endl; +#endif + success = false; + } + + const Number denominator = + ScalarNumber(1.) / (std::abs(rho_P) + eps * rho_max + min); + + t_q = dealii::compare_and_apply_mask( + rho_max, + rho_U + t_q * rho_P, + /* + * rho_P is positive. + * + * Note: Do not take an absolute value here. If we are out of + * bounds we have to ensure that t_r is set to t_min. + */ + (rho_max - rho_U) * denominator, + t_q); + + t_q = dealii::compare_and_apply_mask( + rho_U + t_q * rho_P, + rho_min, + /* + * rho_P is negative. + * + * Note: Do not take an absolute value here. If we are out of + * bounds we have to ensure that t_r is set to t_min. + */ + (rho_U - rho_min) * denominator, + t_q); + + /* + * Ensure that t_min <= t <= t_max. This might not be the case if + * rho_U is outside the interval [rho_min, rho_max]. Furthermore, + * the quotient we take above is prone to numerical cancellation in + * particular in the second pass of the limiter when rho_P might be + * small. + */ + t_r = std::min(t_q, t_r); + t_r = std::min(t_r, t_max); + t_r = std::max(t_r, t_min); + +#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK + /* + * Verify that the new state is within bounds: + */ + const auto rho_new = view.partial_density_1(U + t_r * P); + const auto test_new_min = view.filter_vacuum_density( + std::max(Number(0.), rho_new - relax * rho_max)); + const auto test_new_max = view.filter_vacuum_density( + std::max(Number(0.), rho_min - relax * rho_new)); + if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) { +#ifdef DEBUG_OUTPUT + std::cout << std::fixed << std::setprecision(16); + std::cout << "Bounds violation: high-order [species 1] density!" + << "\n\t\trho min: " << rho_min + << "\n\t\trho min (delta): " + << negative_part(rho_new - rho_min) + << "\n\t\trho: " << rho_new + << "\n\t\trho max (delta): " + << positive_part(rho_new - rho_max) + << "\n\t\trho max: " << rho_max << "\n" + << std::endl; +#endif + success = false; + } +#endif + } + + /* + * Then limit the internal energy. Unfortunately, limiting on the + * "mixture" specific entropy does not guarantee positivity of the mixture + * internal energy, so have to do this. + * + * Note that t_r is given by limiting on partial densities. This + * implementation might not be the most efficient, but we can optimize + * later. + */ + + { + Number t_l = t_min; // good state + Number t_q = t_r; // potentially bad state + + const auto &rho_e_min = std::get<4>(bounds); + + const auto U_l = U + t_l * P; + const auto rho_e_l = view.internal_energy(U_l); + auto psi_l = relax_small * rho_e_l - rho_e_min; + + /* + * Verify that the left state is within bounds. This property might + * be violated for relative CFL numbers larger than 1. The small number + * of 1e-10 is sort of arbitrary, but dealing with round off is annoying + * anyway. + */ + const auto lower_bound = (ScalarNumber(1.) - relax) * rho_e_min; + if (!(std::min(Number(0.), psi_l - lower_bound + Number(1.e-10)) == + Number(0.))) { +#ifdef DEBUG_OUTPUT + std::cout << std::fixed << std::setprecision(16); + std::cout + << "Bounds violation: low-order internal energy (critical)!\n"; + std::cout << "\t\tPsi left: 0 <= " << psi_l << "\n" << std::endl; +#endif + success = false; + } + + const auto U_r = U + t_q * P; + const auto rho_e_r = view.internal_energy(U_r); + + auto psi_r = relax_small * rho_e_r - rho_e_min; + + // Compute simple "linear" limiter + const auto new_t_q = dealii::compare_and_apply_mask< + dealii::SIMDComparison::greater_than>(std::abs(psi_l - psi_r), + Number(0.), + (-psi_r * t_l + psi_l * t_r) / + (psi_l - psi_r), + t_r); + + t_q = dealii::compare_and_apply_mask< + dealii::SIMDComparison::greater_than>( + psi_r, Number(0.), t_r, new_t_q); + + t_r = t_q; + t_r = std::min(t_r, t_max); + t_r = std::max(t_r, t_min); + +#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK + /* + * Verify that the new state is within bounds: + */ + { + const auto U_new = U + t_r * P; + const auto rho_e_new = view.internal_energy(U_new); + const bool rho_e_valid = + std::min(Number(0.), rho_e_new) == Number(0.); + + if (!rho_e_valid) { +#ifdef DEBUG_OUTPUT + std::cout << std::fixed << std::setprecision(16); + std::cout << "Bounds violation: high-order internal energy!\n"; + std::cout << "\t\trho e: 0 <= " << rho_e_new << "\n"; +#endif + success = false; + } + } +#endif + + + } // end limiting on internal energy + + /* + * Then limit the concave entropy functional. This approach is not + * documented anywhere. It is a simple approach that works only under the + * assumptions that: + * 1. The entropy (rho * sbar) is concave. + * 2. The low-order method satisfies the minimum principle (which we can + * prove) + * We do this because we can't take derivatives of the true psi + * functional. It's a mess. + * The idea is as follows: Given the concave functional: + * psi(u) = rho s - rho s_min + * we can show that: psi(t) = rho(U + t P) s(U + t P) - rho(U + t P) s_min + * is also concave. Then, we we find a root of psi(t) = 0, with a simple + * linear search. That is: + * given g(t) = psi(t_L) + (t - t_L) (psi(t_R) - psi(t_L)) / (t_R - t_L) + * we get g(t) = 0 from t->(-psi_r * t_l + psi_l * t_r) / (psi_l - psi_r) + */ + + Number t_l = t_min; // good state + + { + const auto &s_min = std::get<5>(bounds); + +#ifdef DEBUG_OUTPUT_LIMITER + std::cout << std::endl; + std::cout << std::fixed << std::setprecision(16); + std::cout << "t_l: (start) " << t_l << std::endl; + std::cout << "t_r: (start) " << t_r << std::endl; +#endif + + const auto U_r = U + t_r * P; + const auto rho_r = view.density(U_r); + const auto s_r = view.specific_entropy(U_r); + + auto psi_r = relax_small * s_r * rho_r - s_min * rho_r; + +#ifndef DEBUG_EXPENSIVE_BOUNDS_CHECK + /* + * If psi_r > 0 the right state is fine, force returning t_r by + * setting t_l = t_r: + */ + t_l = dealii::compare_and_apply_mask< + dealii::SIMDComparison::greater_than>(psi_r, Number(0.), t_r, t_l); + + /* + * If we have set t_l = t_r everywhere then all states state U_r + * with t_r obey the specific entropy inequality and we can + * return. + * + * This is a very important optimization: Only for 1 in (25 to + * 50) cases do we actually need to limit on the specific entropy + * because one of the right states failed. So we can skip + * constructing the left state U_l, which is expensive. + * + * This implies unfortunately that we might not accurately report + * whether the low_order update U itself obeyed bounds because + * U_r = U + t_r * P pushed us back into bounds. We thus skip + * this shortcut if `DEBUG_EXPENSIVE_BOUNDS_CHECK` is set. + */ + if (t_l == t_r) { +#ifdef DEBUG_OUTPUT_LIMITER + std::cout << "shortcut: t_l == t_r" << std::endl; + std::cout << "psi_r: " << psi_r << std::endl; + std::cout << "t_l: " << t_l << std::endl; + std::cout << "t_r: " << t_r << std::endl; +#endif + return {t_l, success}; + } +#endif + + const auto U_l = U + t_l * P; + const auto rho_l = view.density(U_l); + const auto s_l = view.specific_entropy(U_l); + + auto psi_l = relax_small * s_l * rho_l - s_min * rho_l; + + /* + * Verify that the left state is within bounds. This property might + * be violated for relative CFL numbers larger than 1. The small number + * of 1e-10 is sort of arbitrary, but dealing with round off is annoying + * anyway. + */ + const auto lower_bound = (ScalarNumber(1.) - relax) * s_min * rho_l; + if (!(std::min(Number(0.), psi_l - lower_bound + Number(1.e-10)) == + Number(0.))) { +#ifdef DEBUG_OUTPUT + std::cout << std::fixed << std::setprecision(16); + std::cout + << "Bounds violation: low-order specific entropy (critical)!\n"; + std::cout << "\t\tPsi left: 0 <= " << psi_l << "\n" << std::endl; +#endif + success = false; + } + +#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK + /* + * If psi_r > 0 the right state is fine, force returning t_r by + * setting t_l = t_r: + */ + t_l = dealii::compare_and_apply_mask< + dealii::SIMDComparison::greater_than>(psi_r, Number(0.), t_r, t_l); +#endif + + /* + * Return if the window between t_l and t_r is within the prescribed + * tolerance: + */ + const Number tolerance(parameters.newton_tolerance()); + if (std::max(Number(0.), t_r - t_l - tolerance) == Number(0.)) { +#ifdef DEBUG_OUTPUT_LIMITER + std::cout << "break: t_l and t_r within tolerance" << std::endl; + std::cout << "psi_l: " << psi_l << std::endl; + std::cout << "psi_r: " << psi_r << std::endl; + std::cout << "t_l: " << t_l << std::endl; + std::cout << "t_r: " << t_r << std::endl; +#endif + return {t_l, success}; + } + + /* We got unlucky and have to set t_l with linear formula: */ + + t_l = dealii::compare_and_apply_mask< + dealii::SIMDComparison::greater_than>(std::abs(psi_l - psi_r), + Number(0.), + (-psi_r * t_l + psi_l * t_r) / + (psi_l - psi_r), + t_r); + +#ifdef DEBUG_OUTPUT_LIMITER + std::cout << "s_min: " << s_min << std::endl; + std::cout << "psi_l: " << psi_l << std::endl; + std::cout << "psi_r: " << psi_r << std::endl; + std::cout << "t_l (new): " << t_l << std::endl; + std::cout << "t_r: " << t_r << std::endl; +#endif + + +#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK + /* + * Verify that the new state is within bounds: + */ + { + const auto U_new = U + t_l * P; + const auto rho_new = view.density(U_new); + const auto s_new = view.specific_entropy(U_new); + const auto rho_e_new = view.internal_energy(U_new); + + auto psi_new = relax_small * s_new * rho_new - s_min * rho_new; + + const auto lower_bound = (ScalarNumber(1.) - relax) * s_min * rho_new; + + const bool e_valid = std::min(Number(0.), rho_e_new) == Number(0.); + const bool psi_valid = + std::min(Number(0.), psi_new - lower_bound) == Number(0.); + + if (!e_valid || !psi_valid) { +#ifdef DEBUG_OUTPUT + std::cout << std::fixed << std::setprecision(16); + std::cout << "Bounds violation: high-order specific entropy!\n"; + std::cout << "\t\trho e: 0 <= " << rho_e_new << "\n"; + std::cout << "\t\tPsi: 0 <= " << psi_new << "\n" << std::endl; + std::cout << "\t\tLow: " << lower_bound << "\n" << std::endl; + std::cout << "\t\tDiff: " << psi_new - lower_bound << "\n" + << std::endl; +#endif + success = false; + } + } +#endif + } + + return {t_l, success}; + } + + } // namespace MultiSpeciesEuler +} // namespace ryujin diff --git a/source/multi_euler_ideal/riemann_solver.cc b/source/multi_euler_ideal/riemann_solver.cc new file mode 100644 index 00000000..a4a4ac33 --- /dev/null +++ b/source/multi_euler_ideal/riemann_solver.cc @@ -0,0 +1,25 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2020 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#include "riemann_solver.template.h" + +#include + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + /* instantiations */ + + template class RiemannSolver<1, NUMBER>; + template class RiemannSolver<2, NUMBER>; + template class RiemannSolver<3, NUMBER>; + + template class RiemannSolver<1, dealii::VectorizedArray>; + template class RiemannSolver<2, dealii::VectorizedArray>; + template class RiemannSolver<3, dealii::VectorizedArray>; + } // namespace MultiSpeciesEuler +} // namespace ryujin diff --git a/source/multi_euler_ideal/riemann_solver.h b/source/multi_euler_ideal/riemann_solver.h new file mode 100644 index 00000000..6cc87b86 --- /dev/null +++ b/source/multi_euler_ideal/riemann_solver.h @@ -0,0 +1,273 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2020 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include "hyperbolic_system.h" + +#include + +#include +#include + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + template + class RiemannSolverParameters : public dealii::ParameterAcceptor + { + public: + RiemannSolverParameters(const std::string &subsection = "/RiemannSolver") + : ParameterAcceptor(subsection) + { + } + }; + + + /** + * A fast approximate solver for the 1D Riemann problem. The solver + * ensures that the estimate \f$\lambda_{\text{max}}\f$ that is returned + * for the maximal wavespeed is a strict upper bound. + * + * The solver is based on the guaranteed upper bound wavespeed estimates + * derived in @cite ClaytonGuermondPopov-2022. The key idea is to + * compute an upper bound on the intermediate pressure \f$p^*\f$ in the + * Riemann problem and then use this to derive guaranteed upper bounds + * on the wavespeeds \f$\lambda_1^-\f$ and \f$\lambda_3^+\f$. + * + * For the multi-species case, we use a surrogate gamma + * \f$\gamma_{\min}\f$ computed as the minimum over the stencil to + * ensure the wavespeed bounds remain valid. The maximum wavespeed is + * then (see Eq. (3.3) in [ClaytonDzanicTovar-2025]): + * \f[ + * \lambda_{\max} = \max(|\lambda_1^-|, |\lambda_3^+|), + * \f] + * where \f$\lambda_1^-\f$ and \f$\lambda_3^+\f$ are computed using the + * surrogate EOS with \f$\gamma_{\min}\f$. + * + * @ingroup MultiSpeciesEulerEquations + */ + template + class RiemannSolver + { + public: + /** + * @name Typedefs and constexpr constants + */ + //@{ + + using View = HyperbolicSystemView; + + using ScalarNumber = typename View::ScalarNumber; + + static constexpr auto problem_dimension = View::problem_dimension; + + using state_type = typename View::state_type; + + /** + * Number of components in a "riemann data" state, we store \f$[\rho, v, + * p, a, gamma]\f$, thus, 5. + */ + static constexpr unsigned int riemann_data_size = 5; + + /** + * The array type to store the expanded "riemann data" state for the + * Riemann solver \f$[\rho, v, p, a]\f$ + */ + using primitive_type = typename std::array; + + using precomputed_type = typename View::precomputed_type; + + using PrecomputedVector = typename View::PrecomputedVector; + + using Parameters = RiemannSolverParameters; + + //@} + /** + * @name Compute wavespeed estimates + */ + //@{ + + /** + * Constructor taking a HyperbolicSystem instance as argument + */ + RiemannSolver(const HyperbolicSystem &hyperbolic_system, + const Parameters ¶meters, + const PrecomputedVector &precomputed_values) + : hyperbolic_system(hyperbolic_system) + , parameters(parameters) + , precomputed_values(precomputed_values) + { + } + + /** + * For two given 1D primitive states riemann_data_i and + * riemann_data_j, compute an estimate for an upper bound of the + * maximum wavespeed lambda. + */ + Number compute(const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const; + + /** + * For two given states U_i a U_j and a (normalized) "direction" n_ij + * compute an estimate for an upper bound of the maximum wavespeed + * lambda. + */ + Number compute(const state_type &U_i, + const state_type &U_j, + const unsigned int i, + const unsigned int *js, + const dealii::Tensor<1, dim, Number> &n_ij) const; + + //@} + + protected: + /** @name Internal functions used in the Riemann solver */ + //@{ + + /** + * Continuous and monotonic function c(gamma) as defined in (A.3) + * on page A469 of @cite ClaytonGuermondPopov-2022. + * + * Cost: 0x pow, 1x division, 1x sqrt + */ + Number c(const Number &gamma_Z) const; + + /** + * Compute alpha = 2a / (gamma - 1). + * + * Cost: 0x pow, 1x division, 0x sqrt + */ + Number alpha(const Number &gamma, const Number &a) const; + + /** + * Compute the best available, but expensive, upper bound on the + * expansion-shock case as described in §5.4, Eqn. (5.7) and (5.8) in + * @cite ClaytonGuermondPopov-2022 + * + * Cost: 5x pow, 11x division, 1x sqrt + */ + Number p_star_RS_full(const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const; + + /** + * Compute the best available, but expensive, upper bound on the + * shock-shock case as described in §5.5, Eqn. (5.10) and (5.12) in + * @cite ClaytonGuermondPopov-2022 + * + * Cost: 2x pow, 9x division, 3x sqrt + */ + Number p_star_SS_full(const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const; + + /* + * Compute only the failsafe the failsafe bound for \f$\tilde + * p_2^\ast\f$ (5.11) in @cite ClaytonGuermondPopov-2022 + * + * Cost: 0x pow, 3x division, 3x sqrt + */ + Number p_star_failsafe(const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const; + + /* + * Compute a simultaneous upper bound on (5.7) second formula for + * \tilde p_2^\ast (5.8) first formula for \tilde p_1^\ast (5.11) + * formula for \tilde p_2^\ast in @cite ClaytonGuermondPopov-2022 + * + * Cost: 3x pow, 9x division, 2x sqrt + */ + Number p_star_interpolated(const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const; + + +#ifndef DOXYGEN + /* + * Helper function f for Riemann solver. + */ + Number f(const primitive_type &riemann_data, const Number p_star) const; + + + /* + * Helper function phi for Riemann solver. + */ + Number phi(const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j, + const Number p_in) const; +#endif + + + /** + * See @cite ClaytonGuermondPopov-2022 + * + * The approximate Riemann solver is based on a function phi(p) that is + * monotone increasing in p, concave down and whose (weak) third + * derivative is non-negative and locally bounded. Because we + * actually do not perform any iteration for computing our wavespeed + * estimate we can get away by only implementing a specialized + * variant of the phi function that computes phi(p_max). It inlines + * the implementation of the "f" function and eliminates all + * unnecessary branches in "f". + * + * Cost: 0x pow, 2x division, 2x sqrt + */ + Number phi_of_p_max(const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const; + + + /** + * See @cite GuermondPopov2016 page 912, (3.7) + * + * Cost: 0x pow, 1x division, 1x sqrt + */ + Number lambda1_minus(const primitive_type &riemann_data, + const Number p_star) const; + + + /** + * See @cite GuermondPopov2016 page 912, (3.8) + * + * Cost: 0x pow, 1x division, 1x sqrt + */ + Number lambda3_plus(const primitive_type &primitive_state, + const Number p_star) const; + + + /** + * See @cite GuermondPopov2016 page 912, (3.9) + * + * For two given primitive states riemann_data_i and + * riemann_data_j, and a guess p_2, compute an upper bound + * for lambda. + * + * Cost: 0x pow, 2x division, 2x sqrt (inclusive) + */ + Number compute_lambda(const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j, + const Number p_star) const; + + + /** + * For a given (2+dim dimensional) state vector U, and a + * (normalized) "direction" n_ij, first compute the corresponding + * projected state in the corresponding 1D Riemann problem, and then + * compute and return the Riemann data [rho, u, p, a] (used in the + * approximative Riemann solver). + */ + primitive_type + riemann_data_from_state(const state_type &U, + const Number &p, + const Number &rho, + const dealii::Tensor<1, dim, Number> &n_ij) const; + + private: + const HyperbolicSystem &hyperbolic_system; + const Parameters ¶meters; + const PrecomputedVector &precomputed_values; + //@} + }; + } // namespace MultiSpeciesEuler +} /* namespace ryujin */ diff --git a/source/multi_euler_ideal/riemann_solver.template.h b/source/multi_euler_ideal/riemann_solver.template.h new file mode 100644 index 00000000..cc1a38a6 --- /dev/null +++ b/source/multi_euler_ideal/riemann_solver.template.h @@ -0,0 +1,661 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2020 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +#include + +#include "riemann_solver.h" + +#include +#include + +// #define DEBUG_RIEMANN_SOLVER + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + /* + * The RiemannSolver is a guaranteed maximal wavespeed (GMS) estimate + * for the extended Riemann problem outlined in + * @cite ClaytonGuermondPopov-2022. + * + * In contrast to the algorithm outlined in above reference the + * algorithm takes a couple of shortcuts to significantly decrease the + * computational footprint. These simplifications still guarantee that + * we have an upper bound on the maximal wavespeed - but the number + * bound might be larger. In particular: + * + * - We do not check and treat the case phi(p_min) > 0. This + * corresponds to two expansion waves, see §5.2 in the reference. In + * this case we have + * + * 0 < p_star < p_min <= p_max. + * + * And due to the fact that p_star < p_min the wavespeeds reduce to + * a left wavespeed v_L - a_L and right wavespeed v_R + a_R. This + * implies that it is sufficient to set p_2 to ANY value provided + * that p_2 <= p_min hold true in order to compute the correct + * wavespeed. + * + * If p_2 > p_min then a more pessimistic bound is computed. + * + * - FIXME: Simplification in p_star_RS + */ + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::c(const Number &gamma) const + { + /* + * We implement the continuous and monotonic function c(gamma) as + * defined in (A.3) on page A469 of @cite ClaytonGuermondPopov-2022. + * But with a simplified quick cut-off for the case gamma > 3: + * + * c(gamma)^2 = 1 for gamma <= 5 / 3 + * c(gamma)^2 = (3 * gamma + 11) / (6 * gamma + 6) in between + * c(gamma)^2 = max(1/2, 5 / 6 - slope (gamma - 3)) for gamma > 3 + * + * Due to the fact that the function is monotonic we can simply clip + * the values without checking the conditions: + */ + + constexpr ScalarNumber slope = + ScalarNumber(-0.34976871477801828189920753948709); + + const Number first_radicand = (ScalarNumber(3.) * gamma + Number(11.)) / + (ScalarNumber(6.) * gamma + Number(6.)); + + const Number second_radicand = + Number(5. / 6.) + slope * (gamma - Number(3.)); + + Number radicand = std::min(first_radicand, second_radicand); + radicand = std::min(Number(1.), radicand); + radicand = std::max(Number(1. / 2.), radicand); + + return std::sqrt(radicand); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::alpha(const Number &gamma, + const Number &a) const + { + const Number numerator = ScalarNumber(2.) * a; + const Number denominator = gamma - Number(1.); + + return safe_division(numerator, denominator); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::p_star_RS_full( + const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const + { + const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i; + const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j; + const auto alpha_i = alpha(gamma_i, a_i); + const auto alpha_j = alpha(gamma_j, a_j); + + /* + * First get p_min, p_max. + * + * Then, we get gamma_min/max, and alpha_min/max. Note that the + * *_min/max values are associated with p_min/max and are not + * necessarily the minimum/maximum of *_i vs *_j. + */ + + const Number p_min = std::min(p_i, p_j); + const Number p_max = std::max(p_i, p_j); + + const Number gamma_min = + dealii::compare_and_apply_mask( + p_i, p_j, gamma_i, gamma_j); + + const Number alpha_min = + dealii::compare_and_apply_mask( + p_i, p_j, alpha_i, alpha_j); + + const Number alpha_hat_min = c(gamma_min) * alpha_min; + + const Number alpha_max = dealii::compare_and_apply_mask< + dealii::SIMDComparison::greater_than_or_equal>( + p_i, p_j, alpha_i, alpha_j); + + const Number gamma_m = std::min(gamma_i, gamma_j); + const Number gamma_M = std::max(gamma_i, gamma_j); + + const Number numerator = + dealii::compare_and_apply_mask( + p_max, + Number(0.), + Number(0.), + positive_part(alpha_hat_min + alpha_max - (u_j - u_i))); + + /* + * The admissible set is p_min >= 0. But numerically let's avoid + * division by zero and ensure positivity: + */ + const Number p_ratio = safe_division(p_min, p_max); + + /* + * Here, we use a trick: The r-factor only shows up in the formula + * for the case \gamma_min = \gamma_m, otherwise the r-factor + * vanishes. We can accomplish this by using the following modified + * exponent (where we substitute gamma_m by gamma_min): + */ + const Number r_exponent = + (gamma_M - gamma_min) / (ScalarNumber(2.) * gamma_min * gamma_M); + + /* + * Compute (5.7) first formula for \tilde p_1^\ast and (5.8) + * second formula for \tilde p_2^\ast at the same time: + */ + + const Number first_exponent = + (gamma_M - Number(1.)) / (ScalarNumber(2.) * gamma_M); + + const Number first_exponent_inverse = + safe_division(Number(1.), first_exponent); + + const Number first_denom = + alpha_hat_min * ryujin::pow(p_ratio, r_exponent - first_exponent) + + alpha_max; + + const Number p_1_tilde = + p_max * ryujin::pow(safe_division(numerator, first_denom), + first_exponent_inverse); + +#ifdef DEBUG_RIEMANN_SOLVER + std::cout << "RS p_1_tilde = " << p_1_tilde << "\n"; +#endif + + /* + * Compute (5.7) second formula for \tilde p_2^\ast and (5.8) first + * formula for \tilde p_1^\ast at the same time: + */ + + const Number second_exponent = + (gamma_m - Number(1.)) / (ScalarNumber(2.) * gamma_m); + + const Number second_exponent_inverse = + safe_division(Number(1.), second_exponent); + + Number second_denom = + alpha_hat_min * ryujin::pow(p_ratio, -second_exponent) + + alpha_max * ryujin::pow(p_ratio, r_exponent); + + const Number p_2_tilde = + p_max * ryujin::pow(safe_division(numerator, second_denom), + second_exponent_inverse); + +#ifdef DEBUG_RIEMANN_SOLVER + std::cout << "RS p_2_tilde = " << p_2_tilde << "\n"; +#endif + + return std::min(p_1_tilde, p_2_tilde); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::p_star_SS_full( + const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const + { + const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i; + const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j; + + const Number gamma_m = std::min(gamma_i, gamma_j); + + const Number alpha_hat_i = c(gamma_i) * alpha(gamma_i, a_i); + const Number alpha_hat_j = c(gamma_j) * alpha(gamma_j, a_j); + + /* + * Compute (5.10) formula for \tilde p_1^\ast: + * + * Cost: 2x pow, 4x division, 0x sqrt + */ + + const Number exponent = + (gamma_m - Number(1.)) / (ScalarNumber(2.) * gamma_m); + const Number exponent_inverse = Number(1.) / exponent; + + const Number numerator = + dealii::compare_and_apply_mask( + p_j, + Number(0.), + Number(0.), + positive_part(alpha_hat_i + alpha_hat_j - (u_j - u_i))); + + const Number denominator = + alpha_hat_i * ryujin::pow(safe_division(p_i, p_j), -exponent) + + alpha_hat_j; + + const Number p_1_tilde = + p_j * + ryujin::pow(safe_division(numerator, denominator), exponent_inverse); + +#ifdef DEBUG_RIEMANN_SOLVER + std::cout << "SS p_1_tilde = " << p_1_tilde << "\n"; +#endif + + const auto p_2_tilde = p_star_failsafe(riemann_data_i, riemann_data_j); + + return std::min(p_1_tilde, p_2_tilde); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::p_star_failsafe( + const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const + { + const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i; + const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j; + + /* + * Compute (5.11) formula for \tilde p_2^\ast: + * + * Cost: 0x pow, 3x division, 3x sqrt + */ + + const Number p_max = std::max(p_i, p_j); + + const Number radicand_i = + safe_division(ScalarNumber(2.) * p_max, + rho_i * ((gamma_i + Number(1.)) * p_max + + (gamma_i - Number(1.)) * p_i)); + + const Number x_i = std::sqrt(radicand_i); + + const Number radicand_j = + safe_division(ScalarNumber(2.) * p_max, + rho_j * ((gamma_j + Number(1.)) * p_max + + (gamma_j - Number(1.)) * p_j)); + + const Number x_j = std::sqrt(radicand_j); + + const Number a = x_i + x_j; + const Number b = + dealii::compare_and_apply_mask( + a, Number(0.), Number(0.), u_j - u_i); + + const Number c = -p_i * x_i - p_j * x_j; + + const Number base = safe_division( + std::abs(-b + + std::sqrt(positive_part(b * b - ScalarNumber(4.) * a * c))), + std::abs(ScalarNumber(2.) * a)); + + const Number p_2_tilde = base * base; + +#ifdef DEBUG_RIEMANN_SOLVER + std::cout << "SS p_2_tilde = " << p_2_tilde << "\n"; +#endif + return p_2_tilde; + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::p_star_interpolated( + const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const + { + const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i; + const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j; + const auto alpha_i = alpha(gamma_i, a_i); + const auto alpha_j = alpha(gamma_j, a_j); + + /* + * First get p_min, p_max. + * + * Then, we get gamma_min/max, and alpha_min/max. Note that the + * *_min/max values are associated with p_min/max and are not + * necessarily the minimum/maximum of *_i vs *_j. + */ + + const Number p_min = std::min(p_i, p_j); + const Number p_max = std::max(p_i, p_j); + + const Number gamma_min = + dealii::compare_and_apply_mask( + p_i, p_j, gamma_i, gamma_j); + + const Number alpha_min = + dealii::compare_and_apply_mask( + p_i, p_j, alpha_i, alpha_j); + + const Number alpha_hat_min = c(gamma_min) * alpha_min; + + const Number gamma_max = dealii::compare_and_apply_mask< + dealii::SIMDComparison::greater_than_or_equal>( + p_i, p_j, gamma_i, gamma_j); + + const Number alpha_max = dealii::compare_and_apply_mask< + dealii::SIMDComparison::greater_than_or_equal>( + p_i, p_j, alpha_i, alpha_j); + + const Number alpha_hat_max = c(gamma_max) * alpha_max; + + const Number gamma_m = std::min(gamma_i, gamma_j); + const Number gamma_M = std::max(gamma_i, gamma_j); + + const Number p_ratio = safe_division(p_min, p_max); + + /* + * Here, we use a trick: The r-factor only shows up in the formula + * for the case \gamma_min = \gamma_m, otherwise the r-factor + * vanishes. We can accomplish this by using the following modified + * exponent (where we substitute gamma_m by gamma_min): + */ + const Number r_exponent = + (gamma_M - gamma_min) / (ScalarNumber(2.) * gamma_min * gamma_M); + + /* + * Compute a simultaneous upper bound on + * (5.7) second formula for \tilde p_2^\ast + * (5.8) first formula for \tilde p_1^\ast + * (5.11) formula for \tilde p_2^\ast + */ + + const Number exponent = + (gamma_m - Number(1.)) / (ScalarNumber(2.) * gamma_m); + const Number exponent_inverse = Number(1.) / exponent; + + const Number numerator = + positive_part(alpha_hat_min + /*SIC!*/ alpha_max - (u_j - u_i)); + + Number denominator = alpha_hat_min * ryujin::pow(p_ratio, -exponent) + + alpha_hat_max * ryujin::pow(p_ratio, r_exponent); + + const auto temp = safe_division(numerator, denominator); + + const Number p_tilde = p_max * ryujin::pow(temp, exponent_inverse); + +#ifdef DEBUG_RIEMANN_SOLVER + std::cout << "IN p_*_tilde = " << p_tilde << "\n"; +#endif + + return p_tilde; + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::f(const primitive_type &riemann_data, + const Number p_star) const + { + constexpr ScalarNumber min = std::numeric_limits::min(); + + + const auto &[rho, u, p, gamma, a] = riemann_data; + + const Number gamma_minus_one = gamma - Number(1.); + + const Number Az = ScalarNumber(2.) / (rho * (gamma + Number(1.))); + + const Number Bz = (gamma - Number(1.)) / (gamma + Number(1.)) * p; + + const Number radicand = safe_division(Az, p_star + Bz); + + /* true_value is shock case */ + const Number true_value = (p_star - p) * std::sqrt(radicand); + + const auto exponent = ScalarNumber(0.5) * gamma_minus_one / gamma; + + const Number ratio = safe_division(p_star, p); + const Number factor = ryujin::pow(ratio, exponent) - Number(1.); + + /* false_value is rarefaction case */ + const auto false_value = ScalarNumber(2.) * a * factor / + std::max(gamma_minus_one, Number(min)); + + return dealii::compare_and_apply_mask< + dealii::SIMDComparison::greater_than_or_equal>( + p_star, p, true_value, false_value); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::phi(const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j, + const Number p_in) const + { + const Number &u_i = riemann_data_i[1]; + const Number &u_j = riemann_data_j[1]; + + return f(riemann_data_i, p_in) + f(riemann_data_j, p_in) + u_j - u_i; + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::phi_of_p_max( + const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const + { + const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i; + const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j; + + const Number p_max = std::max(p_i, p_j); + + const Number radicand_inverse_i = + ScalarNumber(0.5) * rho_i * + ((gamma_i + Number(1.)) * p_max + (gamma_i - Number(1.)) * p_i); + + const Number value_i = + safe_division(p_max - p_i, std::sqrt(radicand_inverse_i)); + + const Number radicand_inverse_j = + ScalarNumber(0.5) * rho_j * + ((gamma_j + Number(1.)) * p_max + (gamma_j - Number(1.)) * p_j); + + const Number value_j = + safe_division(p_max - p_j, std::sqrt(radicand_inverse_j)); + + return value_i + value_j + u_j - u_i; + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::lambda1_minus( + const primitive_type &riemann_data, const Number p_star) const + { + const auto &[rho, u, p, gamma, a] = riemann_data; + + const auto factor = + ScalarNumber(0.5) * (gamma + ScalarNumber(1.)) / gamma; + + const Number tmp = safe_division(positive_part(p_star - p), p); + + return u - a * std::sqrt(Number(1.) + factor * tmp); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::lambda3_plus(const primitive_type &riemann_data, + const Number p_star) const + { + const auto &[rho, u, p, gamma, a] = riemann_data; + + const auto factor = + ScalarNumber(0.5) * (gamma + ScalarNumber(1.)) / gamma; + + const Number tmp = safe_division(positive_part(p_star - p), p); + + return u + a * std::sqrt(Number(1.) + factor * tmp); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + RiemannSolver::compute_lambda( + const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j, + const Number p_star) const + { + const Number nu_11 = lambda1_minus(riemann_data_i, p_star); + const Number nu_32 = lambda3_plus(riemann_data_j, p_star); + + return std::max(positive_part(nu_32), negative_part(nu_11)); + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + RiemannSolver::riemann_data_from_state( + const state_type &U, + const Number &p, + const Number &rho, + const dealii::Tensor<1, dim, Number> &n_ij) const -> primitive_type + { + const auto view = hyperbolic_system.view(); + + const auto rho_inverse = ScalarNumber(1.0) / rho; + + const auto m = view.momentum(U); + const auto proj_m = n_ij * m; + + const auto gamma = view.gamma_mixture(U); + const auto a = std::sqrt(gamma * p / rho); + +#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK + AssertThrowSIMD( + Number(p), + [](auto val) { return val >= ScalarNumber(0.); }, + dealii::ExcMessage("Internal error: p < 0.")); + + AssertThrowSIMD( + gamma, + [](auto val) { return val >= ScalarNumber(1.); }, + dealii::ExcMessage("Internal error: gamma <= 1.")); +#endif + + return {{rho, proj_m * rho_inverse, p, gamma, a}}; + } + + + template + Number RiemannSolver::compute( + const primitive_type &riemann_data_i, + const primitive_type &riemann_data_j) const + { + const auto view = hyperbolic_system.view(); + + const auto &[rho_i, u_i, p_i, gamma_i, a_i] = riemann_data_i; + const auto &[rho_j, u_j, p_j, gamma_j, a_j] = riemann_data_j; + +#ifdef DEBUG_RIEMANN_SOLVER + std::cout << "rho_left: " << rho_i << std::endl; + std::cout << "u_left: " << u_i << std::endl; + std::cout << "p_left: " << p_i << std::endl; + std::cout << "gamma_left: " << gamma_i << std::endl; + std::cout << "a_left: " << a_i << std::endl; + std::cout << "rho_right: " << rho_j << std::endl; + std::cout << "u_right: " << u_j << std::endl; + std::cout << "p_right: " << p_j << std::endl; + std::cout << "gamma_right: " << gamma_j << std::endl; + std::cout << "a_right: " << a_j << std::endl; +#endif + + const Number p_max = std::max(p_i, p_j); + const Number phi_p_max = phi_of_p_max(riemann_data_i, riemann_data_j); + + if (!view.compute_strict_bounds()) { +#ifdef DEBUG_RIEMANN_SOLVER + const Number p_star_RS = p_star_RS_full(riemann_data_i, riemann_data_j); + const Number p_star_SS = p_star_SS_full(riemann_data_i, riemann_data_j); + const Number p_debug = + dealii::compare_and_apply_mask( + phi_p_max, Number(0.), p_star_SS, std::min(p_max, p_star_RS)); + std::cout << " p^*_debug = " << p_debug << "\n"; + std::cout << " phi(p_*_d) = " + << phi(riemann_data_i, riemann_data_j, p_debug) << "\n"; + std::cout << "-> lambda_deb = " + << compute_lambda(riemann_data_i, riemann_data_j, p_debug) + << std::endl; +#endif + + const Number p_star_tilde = + p_star_interpolated(riemann_data_i, riemann_data_j); + const Number p_star_backup = + p_star_failsafe(riemann_data_i, riemann_data_j); + + const Number p_2 = + dealii::compare_and_apply_mask( + phi_p_max, + Number(0.), + std::min(p_star_tilde, p_star_backup), + std::min(p_max, p_star_tilde)); + +#ifdef DEBUG_RIEMANN_SOLVER + std::cout << " p^*_tilde = " << p_2 << "\n"; + std::cout << " phi(p_*_t) = " + << phi(riemann_data_i, riemann_data_j, p_2) << "\n"; + std::cout << "-> lambda_max = " + << compute_lambda(riemann_data_i, riemann_data_j, p_2) << "\n" + << std::endl; +#endif + + return compute_lambda(riemann_data_i, riemann_data_j, p_2); + } + + const Number p_star_RS = p_star_RS_full(riemann_data_i, riemann_data_j); + const Number p_star_SS = p_star_SS_full(riemann_data_i, riemann_data_j); + + const Number p_2 = + dealii::compare_and_apply_mask( + phi_p_max, Number(0.), p_star_SS, std::min(p_max, p_star_RS)); + +#ifdef DEBUG_RIEMANN_SOLVER + std::cout << " p^*_tilde = " << p_2 << "\n"; + std::cout << " phi(p_*_t) = " + << phi(riemann_data_i, riemann_data_j, p_2) << "\n"; + std::cout << "-> lambda_max = " + << compute_lambda(riemann_data_i, riemann_data_j, p_2) + << std::endl; +#endif + + return compute_lambda(riemann_data_i, riemann_data_j, p_2); + } + + + template + DEAL_II_ALWAYS_INLINE inline Number RiemannSolver::compute( + const state_type &U_i, + const state_type &U_j, + const unsigned int i, + const unsigned int *js, + const dealii::Tensor<1, dim, Number> &n_ij) const + { + const auto &[rho_i, p_i, dont_use_i, s_i, harten_i] = + precomputed_values.template get_tensor(i); + + const auto &[rho_j, p_j, dont_use_j, s_j, harten_j] = + precomputed_values.template get_tensor(js); + + const auto riemann_data_i = + riemann_data_from_state(U_i, p_i, rho_i, n_ij); + const auto riemann_data_j = + riemann_data_from_state(U_j, p_j, rho_j, n_ij); + + return compute(riemann_data_i, riemann_data_j); + } + + + } // namespace MultiSpeciesEuler +} // namespace ryujin From b9b31408e8b35de7f3eb7e54cbd4d106f722502e Mon Sep 17 00:00:00 2001 From: Eric Tovar Date: Sun, 1 Mar 2026 12:12:41 -0700 Subject: [PATCH 2/5] MultiSpeciesEuler: clean up and small bug fixes --- source/multi_euler_ideal/hyperbolic_system.h | 308 ++++++++++-------- source/multi_euler_ideal/indicator.h | 6 +- .../initial_state_contrast.h | 58 ++-- .../initial_state_exact_riemann_solution.h | 54 +-- .../initial_state_function.h | 64 ++-- .../initial_state_icf_like.h | 87 +++-- .../initial_state_library_multi_euler.h | 13 +- .../initial_state_radial_contrast.h | 76 ++--- .../initial_state_shock_bubble.h | 117 +++---- .../initial_state_smooth_wave.h | 36 +- .../initial_state_three_state_contrast.h | 70 ++-- source/multi_euler_ideal/limiter.h | 217 +++++------- source/multi_euler_ideal/limiter.template.h | 151 ++------- source/multi_euler_ideal/riemann_solver.h | 6 +- source/multi_euler_ideal/species_parameters.h | 29 ++ 15 files changed, 630 insertions(+), 662 deletions(-) create mode 100644 source/multi_euler_ideal/species_parameters.h diff --git a/source/multi_euler_ideal/hyperbolic_system.h b/source/multi_euler_ideal/hyperbolic_system.h index 9da92bbc..ef20069f 100644 --- a/source/multi_euler_ideal/hyperbolic_system.h +++ b/source/multi_euler_ideal/hyperbolic_system.h @@ -6,6 +6,8 @@ #pragma once +#include "species_parameters.h" + #include #include @@ -44,6 +46,40 @@ namespace ryujin } + /** + * Convert from user-friendly primitive format + * (Y_0, ..., Y_{n-2}, rho, u, p) + * to internal primitive format + * (alpha_0*rho_0, ..., alpha_{n-1}*rho_{n-1}, u, p) + * + * The last species' partial density is computed from the constraint + * that mass fractions sum to one. + */ + template + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, n_species + 2, Number> + extend_primitive( + const dealii::Tensor<1, n_species + 2, Number> &primitive_in) + { + dealii::Tensor<1, n_species + 2, Number> result; + + const auto rho = primitive_in[n_species - 1]; + Number Y_sum = Number(0.); + + /* Convert mass fractions to partial densities */ + for (unsigned int k = 0; k < n_species - 1; ++k) { + result[k] = primitive_in[k] * rho; + Y_sum += primitive_in[k]; + } + result[n_species - 1] = (Number(1.) - Y_sum) * rho; + + /* Copy velocity and pressure */ + result[n_species] = primitive_in[n_species]; + result[n_species + 1] = primitive_in[n_species + 1]; + + return result; + } + + template class HyperbolicSystemView; @@ -92,7 +128,6 @@ namespace ryujin * \varepsilon(\mathbf{u}) > 0,\; s(\mathbf{u}) \geq s_{\min} \}. * \f] * - * @note Currently hardcoded for 2 species (\f$n_s = 2\f$). * * @ingroup MultiSpeciesEulerEquations */ @@ -143,10 +178,10 @@ namespace ryujin */ //@{ - dealii::Tensor<1, /* num_species = */ 2, double> cp_for_each_species_; - dealii::Tensor<1, /* num_species = */ 2, double> cv_for_each_species_; - dealii::Tensor<1, /* num_species = */ 2, double> r_for_each_species_; - dealii::Tensor<1, /* num_species = */ 2, double> gamma_for_each_species_; + dealii::Tensor<1, n_species, double> cp_for_each_species_; + dealii::Tensor<1, n_species, double> cv_for_each_species_; + dealii::Tensor<1, n_species, double> r_for_each_species_; + dealii::Tensor<1, n_species, double> gamma_for_each_species_; double reference_density_; double vacuum_state_relaxation_small_; double vacuum_state_relaxation_large_; @@ -206,25 +241,25 @@ namespace ryujin */ //@{ - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 2, double> + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, n_species, double> cp_for_each_species() const { return hyperbolic_system_.cp_for_each_species_; } - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 2, double> + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, n_species, double> cv_for_each_species() const { return hyperbolic_system_.cv_for_each_species_; } - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 2, double> + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, n_species, double> r_for_each_species() const { return hyperbolic_system_.r_for_each_species_; } - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 2, double> + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, n_species, double> gamma_for_each_species() const { return hyperbolic_system_.gamma_for_each_species_; @@ -286,11 +321,10 @@ namespace ryujin //@{ /** - * The dimension of the state space. - * FIXME: Generalize for user-defined number of species. - * This is assuming 2 species only. + * The dimension of the state space: n_species partial densities + + * dim momentum components + 1 total energy. */ - static constexpr unsigned int problem_dimension = 3 + dim; + static constexpr unsigned int problem_dimension = n_species + 1 + dim; /** * Storage type for a (conserved) state vector \f$\boldsymbol U\f$. @@ -325,13 +359,17 @@ namespace ryujin */ static inline const auto component_names = []() -> std::array { - if constexpr (dim == 1) - return {"alpha_rho_0", "alpha_rho_1", "m", "E"}; - else if constexpr (dim == 2) - return {"alpha_rho_0", "alpha_rho_1", "m_1", "m_2", "E"}; - else if constexpr (dim == 3) - return {"alpha_rho_0", "alpha_rho_1", "m_1", "m_2", "m_3", "E"}; - __builtin_trap(); + std::array names; + for (unsigned int k = 0; k < n_species; ++k) + names[k] = "alpha_rho_" + std::to_string(k); + if constexpr (dim == 1) { + names[n_species] = "m"; + } else { + for (unsigned int d = 0; d < dim; ++d) + names[n_species + d] = "m_" + std::to_string(d + 1); + } + names[n_species + dim] = "E"; + return names; }(); /** @@ -340,13 +378,17 @@ namespace ryujin */ static inline const auto primitive_component_names = []() -> std::array { - if constexpr (dim == 1) - return {"alpha_rho_0", "alpha_rho_1", "v", "p"}; - else if constexpr (dim == 2) - return {"alpha_rho_0", "alpha_rho_1", "v_1", "v_2", "p"}; - else if constexpr (dim == 3) - return {"alpha_rho_0", "alpha_rho_1", "v_1", "v_2", "v_3", "p"}; - __builtin_trap(); + std::array names; + for (unsigned int k = 0; k < n_species; ++k) + names[k] = "alpha_rho_" + std::to_string(k); + if constexpr (dim == 1) { + names[n_species] = "v"; + } else { + for (unsigned int d = 0; d < dim; ++d) + names[n_species + d] = "v_" + std::to_string(d + 1); + } + names[n_species + dim] = "p"; + return names; }(); /** @@ -420,20 +462,14 @@ namespace ryujin //@{ /** - * For a given (3+dim dimensional) state vector U, return - * the partial density for the first species U[0]. - */ - static Number partial_density_0(const state_type &U); - - /** - * For a given (3+dim dimensional) state vector U, return - * the partial density for the second species U[1]. + * For a given state vector U, return the partial density + * for species k, i.e., U[k]. */ - static Number partial_density_1(const state_type &U); + static Number partial_density(const state_type &U, unsigned int k); /** - * For a given (3+dim dimensional) state vector U, return - * the total mixture density obtained by summing partial densities. + * For a given state vector U, return the total mixture + * density obtained by summing all partial densities. */ static Number density(const state_type &U); @@ -446,25 +482,26 @@ namespace ryujin Number filter_vacuum_density(const Number &rho) const; /** - * For a given (3+dim dimensional) state vector U, return - * the momentum vector [U[2], ..., U[1+dim]]. + * For a given (n_species+1+dim dimensional) state vector U, + * return the momentum vector + * [U[n_species], ..., U[n_species+dim-1]]. */ static dealii::Tensor<1, dim, Number> momentum(const state_type &U); /** - * For a given (3+dim dimensional) state vector U, return - * the total energy U[2+dim] + * For a given (n_species+1+dim dimensional) state vector U, return + * the total energy U[n_species+dim] */ static Number total_energy(const state_type &U); /** - * For a given (3+dim dimensional) state vector U, compute + * For a given (n_species+1+dim dimensional) state vector U, compute * and return the internal energy \f$\varepsilon = (\rho e)\f$. */ static Number internal_energy(const state_type &U); /** - * For a given (3+dim dimensional) state vector U, compute + * For a given (n_species+1+dim dimensional) state vector U, compute * and return the mixture gamma for the mixture EOS: * \f[ * \overline{\gamma} = (sum of (alpha_k rho_k) * c_{p, k}) / (sum of @@ -474,7 +511,7 @@ namespace ryujin Number gamma_mixture(const state_type &U) const; /** - * For a given (3+dim dimensional) state vector U, compute + * For a given (n_species+1+dim dimensional) state vector U, compute * and return the mixture specific heat capacity at constant pressure: * \f[ * \overline{c}_p = (sum of (alpha_k rho_k) / rho * c_{p, k}) @@ -483,7 +520,7 @@ namespace ryujin Number cp_mixture(const state_type &U) const; /** - * For a given (3+dim dimensional) state vector U, compute + * For a given (n_species+1+dim dimensional) state vector U, compute * and return the mixture specific heat capacity at constant volume: * \f[ * \overline{c}_v = (sum of (alpha_k rho_k) / rho * c_{v, k}) @@ -492,7 +529,7 @@ namespace ryujin Number cv_mixture(const state_type &U) const; /** - * For a given (3+dim dimensional) state vector U, compute + * For a given (n_species+1+dim dimensional) state vector U, compute * and return the pressure \f$p\f$. * * We assume that the pressure is given by a mixture ideal EOS @@ -503,7 +540,7 @@ namespace ryujin Number pressure(const state_type &U) const; /** - * For a given (3+dim dimensional) state vector U, + * For a given (n_species+1+dim dimensional) state vector U, * compute and return the speed of sound \f$c\f$: * * We assume that the pressure is given by a mixture ideal EOS @@ -522,7 +559,7 @@ namespace ryujin //@{ /** - * For a given (3+dim dimensional) state vector U, compute + * For a given (n_species+1+dim dimensional) state vector U, compute * and return the physical specific entropy. Following Eq. (2.8) in * [ClaytonDzanicTovar-2025]: * \f[ @@ -536,21 +573,21 @@ namespace ryujin Number specific_entropy(const state_type &U) const; /** - * For a given (3+dim dimensional) state vector U, compute + * For a given (n_species+1+dim dimensional) state vector U, compute * and return a surrogate Harten-type entropy. Following Section 4 of * [ClaytonDzanicTovar-2025], we use: * \f[ - * \eta(\mathbf{u}; \gamma_{\min}) = (\rho \varepsilon)^{1/(\gamma_{\min}+1)}, - * \f] - * where \f$\gamma_{\min}\f$ is the minimum surrogate gamma over the - * stencil. This entropy is chosen to ensure convexity properties - * required by the entropy-viscosity indicator. + * \eta(\mathbf{u}; \gamma_{\min}) = (\rho + * \varepsilon)^{1/(\gamma_{\min}+1)}, \f] where \f$\gamma_{\min}\f$ is + * the minimum surrogate gamma over the stencil. This entropy is chosen to + * ensure convexity properties required by the entropy-viscosity + * indicator. */ Number surrogate_harten_entropy(const state_type &U, const Number &gamma_min) const; /** - * For a given (3+dim dimensional) state vector U, compute + * For a given (n_species+1+dim dimensional) state vector U, compute * and return the derivative \f$\eta'\f$ of the Harten-type entropy. */ mixture_state_type @@ -559,7 +596,7 @@ namespace ryujin const Number &gamma_min) const; /** - * For a given (3+dim dimensional) state vector U and + * For a given (n_species+1+dim dimensional) state vector U and * pressure p, compute a surrogate gamma. Following * Section 3 of [ClaytonDzanicTovar-2025]: * \f[ @@ -573,7 +610,7 @@ namespace ryujin Number surrogate_gamma(const state_type &U, const Number &p) const; /** - * For a given (3+dim dimensional) state vector U and + * For a given (n_species+1+dim dimensional) state vector U and * gamma gamma, compute a surrogate pressure: * \f[ * p(\rho, e, \gamma) = (\gamma - 1) (\rho e) @@ -584,7 +621,7 @@ namespace ryujin Number surrogate_pressure(const state_type &U, const Number &gamma) const; /** - * For a given (3+dim dimensional) state vector U and + * For a given (n_species+1+dim dimensional) state vector U and * gamma gamma, compute a surrogate speed of sound. */ Number surrogate_speed_of_sound(const state_type &U, @@ -762,9 +799,10 @@ namespace ryujin state_type expand_state(const ST &state) const; /** - * Given an initial state [alpha_0 rho_0, alpha_1 rho_1, u_1, ..., u_d, p] - * return a conserved state [alpha_0 rho_0, alpha_1 rho_1, m_1, ..., m_d, - * E]. + * Given an initial state + * [alpha_0 rho_0, ..., alpha_{n-1} rho_{n-1}, u_1, ..., u_d, p] + * return a conserved state + * [alpha_0 rho_0, ..., alpha_{n-1} rho_{n-1}, m_1, ..., m_d, E]. * * @note This function is used to conveniently convert (user * provided) primitive initial states with pressure values to a @@ -825,11 +863,11 @@ namespace ryujin reference_density_, "Problem specific density reference"); - cp_for_each_species_[0] = 1.4; - cv_for_each_species_[0] = 1.0; - - cp_for_each_species_[1] = 1.67; - cv_for_each_species_[1] = 1.0; + /* Set default values for species parameters */ + for (unsigned int k = 0; k < n_species; ++k) { + cp_for_each_species_[k] = 1.4 + 0.27 * k; /* Default: 1.4, 1.67, ... */ + cv_for_each_species_[k] = 1.0; + } add_parameter( "c_p for each species", @@ -855,15 +893,12 @@ namespace ryujin * And finally populate the r and gamma values. */ const auto populate_values = [this]() { - r_for_each_species_[0] = - cp_for_each_species_[0] - cv_for_each_species_[0]; - r_for_each_species_[1] = - cp_for_each_species_[1] - cv_for_each_species_[1]; - - gamma_for_each_species_[0] = - cp_for_each_species_[0] / cv_for_each_species_[0]; - gamma_for_each_species_[1] = - cp_for_each_species_[1] / cv_for_each_species_[1]; + for (unsigned int k = 0; k < n_species; ++k) { + r_for_each_species_[k] = + cp_for_each_species_[k] - cv_for_each_species_[k]; + gamma_for_each_species_[k] = + cp_for_each_species_[k] / cv_for_each_species_[k]; + } }; ParameterAcceptor::parse_parameters_call_back.connect(populate_values); @@ -950,23 +985,21 @@ namespace ryujin template DEAL_II_ALWAYS_INLINE inline Number - HyperbolicSystemView::partial_density_0(const state_type &U) + HyperbolicSystemView::partial_density(const state_type &U, + unsigned int k) { - return U[0]; - } - - template - DEAL_II_ALWAYS_INLINE inline Number - HyperbolicSystemView::partial_density_1(const state_type &U) - { - return U[1]; + AssertIndexRange(k, n_species); + return U[k]; } template DEAL_II_ALWAYS_INLINE inline Number HyperbolicSystemView::density(const state_type &U) { - return partial_density_0(U) + partial_density_1(U); + auto result = Number(0.); + for (unsigned int k = 0; k < n_species; ++k) + result += partial_density(U, k); + return result; } @@ -990,7 +1023,7 @@ namespace ryujin { dealii::Tensor<1, dim, Number> result; for (unsigned int i = 0; i < dim; ++i) - result[i] = U[2 + i]; + result[i] = U[n_species + i]; return result; } @@ -999,7 +1032,7 @@ namespace ryujin DEAL_II_ALWAYS_INLINE inline Number HyperbolicSystemView::total_energy(const state_type &U) { - return U[2 + dim]; + return U[n_species + dim]; } @@ -1023,8 +1056,8 @@ namespace ryujin Number result = 0.; const Number rho = density(U); - for (unsigned int i = 0; i < 2; ++i) { - result += U[i] / rho * ScalarNumber(cp_for_each_species()[i]); + for (unsigned int k = 0; k < n_species; ++k) { + result += U[k] / rho * ScalarNumber(cp_for_each_species()[k]); } return result; } @@ -1036,8 +1069,8 @@ namespace ryujin Number result = 0.; const Number rho = density(U); - for (unsigned int i = 0; i < 2; ++i) { - result += U[i] / rho * ScalarNumber(cv_for_each_species()[i]); + for (unsigned int k = 0; k < n_species; ++k) { + result += U[k] / rho * ScalarNumber(cv_for_each_species()[k]); } return result; } @@ -1083,22 +1116,22 @@ namespace ryujin Number cv_bar = 0.; Number r_bar = 0.; - for (unsigned int i = 0; i < 2; ++i) { - const auto Y_i = U[i] * rho_inverse; - r_bar += Y_i * Number(r_for_each_species()[i]); - cv_bar += Y_i * Number(cv_for_each_species()[i]); + for (unsigned int k = 0; k < n_species; ++k) { + const auto Y_k = U[k] * rho_inverse; + r_bar += Y_k * Number(r_for_each_species()[k]); + cv_bar += Y_k * Number(cv_for_each_species()[k]); } Number K_factor = 0.; - for (unsigned int i = 0; i < 2; ++i) { - const auto Y_i = U[i] * rho_inverse; - const auto cv_i = Number(cv_for_each_species()[i]); - const auto r_i = Number(r_for_each_species()[i]); - const auto gm1 = Number(gamma_for_each_species()[i] - 1.); - - const auto K_i = - cv_i * std::log(cv_i / cv_bar * ryujin::pow(r_i / r_bar, gm1)); - K_factor += Y_i * K_i; + for (unsigned int k = 0; k < n_species; ++k) { + const auto Y_k = U[k] * rho_inverse; + const auto cv_k = Number(cv_for_each_species()[k]); + const auto r_k = Number(r_for_each_species()[k]); + const auto gm1 = Number(gamma_for_each_species()[k] - 1.); + + const auto K_k = + cv_k * std::log(cv_k / cv_bar * ryujin::pow(r_k / r_bar, gm1)); + K_factor += Y_k * K_k; } const auto s_bar = @@ -1239,8 +1272,10 @@ namespace ryujin */ const auto rho = density(U); - const auto Y_0 = partial_density_0(U) / rho; - const auto Y_1 = partial_density_1(U) / rho; + /* Store mass fractions Y_k = (alpha_k rho_k) / rho */ + dealii::Tensor<1, n_species, Number> Y; + for (unsigned int k = 0; k < n_species; ++k) + Y[k] = partial_density(U, k) / rho; const auto m = momentum(U); const auto vn = m * normal / rho; @@ -1357,13 +1392,13 @@ namespace ryujin state_type U_new; - U_new[0] = Y_0 * rho_new; - U_new[1] = Y_1 * rho_new; + for (unsigned int k = 0; k < n_species; ++k) + U_new[k] = Y[k] * rho_new; for (unsigned int d = 0; d < dim; ++d) { - U_new[2 + d] = rho_new * (vn_new * normal + vperp)[d]; + U_new[n_species + d] = rho_new * (vn_new * normal + vperp)[d]; } - U_new[2 + dim] = + U_new[n_species + dim] = rho_e_new + 0.5 * rho_new * (vn_new * vn_new + vperp.norm_square()); return U_new; @@ -1387,18 +1422,18 @@ namespace ryujin } else if (id == Boundary::dirichlet_momentum) { /* Only enforce Dirichlet conditions on the momentum: */ auto m_dirichlet = momentum(get_dirichlet_data()); - for (unsigned int k = 0; k < dim; ++k) - result[k + 2] = m_dirichlet[k]; + for (unsigned int d = 0; d < dim; ++d) + result[n_species + d] = m_dirichlet[d]; } else if (id == Boundary::slip) { auto m = momentum(U); m -= 1. * (m * normal) * normal; - for (unsigned int k = 0; k < dim; ++k) - result[k + 2] = m[k]; + for (unsigned int d = 0; d < dim; ++d) + result[n_species + d] = m[d]; } else if (id == Boundary::no_slip) { - for (unsigned int k = 0; k < dim; ++k) - result[k + 2] = Number(0.); + for (unsigned int d = 0; d < dim; ++d) + result[n_species + d] = Number(0.); } else if (id == Boundary::dynamic) { /* @@ -1470,12 +1505,12 @@ namespace ryujin flux_type result; for (unsigned int i = 0; i < dim; ++i) { - result[0][i] = U[0] * (m[i] * rho_inverse); - result[1][i] = U[1] * (m[i] * rho_inverse); - result[2 + i] = m * (m[i] * rho_inverse); - result[2 + i][i] += p; + for (unsigned int k = 0; k < n_species; ++k) + result[k][i] = U[k] * (m[i] * rho_inverse); + result[n_species + i] = m * (m[i] * rho_inverse); + result[n_species + i][i] += p; } - result[dim + 2] = m * (rho_inverse * (E + p)); + result[n_species + dim] = m * (rho_inverse * (E + p)); return result; } @@ -1550,17 +1585,20 @@ namespace ryujin using T = typename ST::value_type; static_assert(std::is_same_v, "template mismatch"); - constexpr auto dim2 = ST::dimension - 3; + constexpr auto dim2 = ST::dimension - n_species - 1; static_assert(dim >= dim2, "the space dimension of the argument state must not be " "larger than the one of the target state"); state_type result; - result[0] = state[0]; - result[1] = state[1]; - result[dim + 2] = state[dim2 + 2]; - for (unsigned int i = 2; i < dim2 + 2; ++i) - result[i] = state[i]; + /* Copy partial densities */ + for (unsigned int k = 0; k < n_species; ++k) + result[k] = state[k]; + /* Copy total energy */ + result[n_species + dim] = state[n_species + dim2]; + /* Copy momentum components (and zero-fill extra dimensions) */ + for (unsigned int i = 0; i < dim2; ++i) + result[n_species + i] = state[n_species + i]; return result; } @@ -1585,17 +1623,17 @@ namespace ryujin const auto rho = density(primitive_state); /* extract velocity: */ const auto u = /*SIC!*/ momentum(primitive_state); - /* extract specific internal energy: */ - const auto &p = primitive_state[dim + 2]; + /* extract pressure: */ + const auto &p = primitive_state[n_species + dim]; auto state = primitive_state; /* Fix up momentum: */ - for (unsigned int i = 2; i < dim + 2; ++i) + for (unsigned int i = n_species; i < n_species + dim; ++i) state[i] *= rho; /* Compute total energy: */ const Number gamma_bar = gamma_mixture(primitive_state); - state[dim + 2] = + state[n_species + dim] = p / (gamma_bar - Number(1.)) + ScalarNumber(0.5) * rho * u * u; return state; @@ -1613,10 +1651,10 @@ namespace ryujin auto primitive_state = state; /* Fix up velocity: */ - for (unsigned int i = 2; i < dim + 2; ++i) + for (unsigned int i = n_species; i < n_species + dim; ++i) primitive_state[i] *= rho_inverse; /* Set pressure: */ - primitive_state[dim + 2] = p; + primitive_state[n_species + dim] = p; return primitive_state; } @@ -1630,7 +1668,7 @@ namespace ryujin auto result = state; const auto M = lambda(momentum(state)); for (unsigned int d = 0; d < dim; ++d) - result[2 + d] = M[d]; + result[n_species + d] = M[d]; return result; } } // namespace MultiSpeciesEuler diff --git a/source/multi_euler_ideal/indicator.h b/source/multi_euler_ideal/indicator.h index a84cc50b..a561bbb8 100644 --- a/source/multi_euler_ideal/indicator.h +++ b/source/multi_euler_ideal/indicator.h @@ -30,7 +30,7 @@ namespace ryujin evc_factor_ = ScalarNumber(1.); add_parameter("evc factor", evc_factor_, - "Factor for scaling the entropy viscocity commuator"); + "Factor for scaling the entropy viscosity commutator"); } ACCESSOR_READ_ONLY(evc_factor); @@ -247,7 +247,7 @@ namespace ryujin (eta_j * rho_j_inverse - eta_i * rho_i_inverse) * (m_j * c_ij); left += entropy_flux; - for (unsigned int k = 0; k < problem_dimension - 1; ++k) { + for (unsigned int k = 0; k < 2 + dim; ++k) { const auto component = (f_j[k] - f_i[k]) * c_ij; right[k] += component; } @@ -262,7 +262,7 @@ namespace ryujin Number numerator = left; Number denominator = std::abs(left); - for (unsigned int k = 0; k < problem_dimension - 1; ++k) { + for (unsigned int k = 0; k < 2 + dim; ++k) { numerator -= d_eta_i[k] * right[k]; denominator += std::abs(d_eta_i[k] * right[k]); } diff --git a/source/multi_euler_ideal/initial_state_contrast.h b/source/multi_euler_ideal/initial_state_contrast.h index 71d2bdf5..9210e751 100644 --- a/source/multi_euler_ideal/initial_state_contrast.h +++ b/source/multi_euler_ideal/initial_state_contrast.h @@ -6,16 +6,25 @@ #pragma once +#include "hyperbolic_system.h" + #include namespace ryujin { namespace MultiSpeciesEulerInitialStates { + using namespace MultiSpeciesEuler; + /** * An initial state formed by a contrast of a given "left" and "right" * primitive (initial) state. * + * The primitive state format is: + * (Y_0, Y_1, ..., Y_{n-2}, rho, u, p) + * where Y_k are mass fractions for k = 0 to n_species-2, and the last + * mass fraction Y_{n-1} = 1 - sum(Y_k) is computed automatically. + * * @note The @p t argument is ignored. * * @ingroup MultiSpeciesEulerEquations @@ -29,34 +38,43 @@ namespace ryujin typename Description::template HyperbolicSystemView; using state_type = typename View::state_type; + /* 1D primitive state: (Y_0, ..., Y_{n-2}, rho, u, p) */ + static constexpr unsigned int primitive_dim = n_species + 2; + using primitive_state_type = dealii::Tensor<1, primitive_dim, Number>; + Contrast(const HyperbolicSystem &hyperbolic_system, const std::string subsection) : InitialState("contrast", subsection) , hyperbolic_system_(hyperbolic_system) { - temp_left_[0] = 0.5; - temp_left_[1] = 1.4; - temp_left_[2] = 0.0; - temp_left_[3] = 1.0; + /* Default: equal mass fractions */ + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_left_[k] = Number(1.) / Number(n_species); + temp_left_[n_species - 1] = 1.4; /* rho */ + temp_left_[n_species] = 0.0; /* u */ + temp_left_[n_species + 1] = 1.0; /* p */ this->add_parameter( "primitive state left", temp_left_, - "Initial 1d primitive state (Y_0, rho, u, p) on the left"); + "Initial 1d primitive state (Y_0, ..., Y_{n-2}, rho, u, p) on the " + "left"); - temp_right_[0] = 0.5; - temp_right_[1] = 1.4; - temp_right_[2] = 0.0; - temp_right_[3] = 1.0; + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_right_[k] = Number(1.) / Number(n_species); + temp_right_[n_species - 1] = 1.4; /* rho */ + temp_right_[n_species] = 0.0; /* u */ + temp_right_[n_species + 1] = 1.0; /* p */ this->add_parameter( "primitive state right", temp_right_, - "Initial 1d primitive state (Y_0, rho, u, p) on the right"); + "Initial 1d primitive state (Y_0, ..., Y_{n-2}, rho, u, p) on the " + "right"); const auto convert_states = [&]() { const auto view = hyperbolic_system_.template view(); - const auto primitive_left_ = extend(temp_left_); - const auto primitive_right_ = extend(temp_right_); + const auto primitive_left_ = extend_primitive(temp_left_); + const auto primitive_right_ = extend_primitive(temp_right_); state_left_ = view.from_initial_state(primitive_left_); state_right_ = view.from_initial_state(primitive_right_); @@ -73,24 +91,12 @@ namespace ryujin private: const HyperbolicSystem &hyperbolic_system_; - dealii::Tensor<1, 4, Number> temp_left_; - dealii::Tensor<1, 4, Number> temp_right_; + primitive_state_type temp_left_; + primitive_state_type temp_right_; state_type state_left_; state_type state_right_; - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> - extend(dealii::Tensor<1, 4, Number> &temp_in) const - { - dealii::Tensor<1, 4, Number> result; - result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; - result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; - - for (unsigned int i = 2; i < 4; ++i) - result[i] = temp_in[i]; - - return result; - } }; } // namespace MultiSpeciesEulerInitialStates } // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_exact_riemann_solution.h b/source/multi_euler_ideal/initial_state_exact_riemann_solution.h index 298123d1..449ee48b 100644 --- a/source/multi_euler_ideal/initial_state_exact_riemann_solution.h +++ b/source/multi_euler_ideal/initial_state_exact_riemann_solution.h @@ -6,6 +6,8 @@ #pragma once +#include "hyperbolic_system.h" + #include #include @@ -17,6 +19,8 @@ namespace ryujin { namespace MultiSpeciesEulerInitialStates { + using namespace MultiSpeciesEuler; + /** * The exact Riemann solution. * @@ -24,6 +28,9 @@ namespace ryujin * compressible multi-species Euler equations with an ideal gas equation of * state for each species. * + * @note This class is currently only implemented for exactly 2 species. + * A compile-time assertion will fail if n_species != 2. + * * @note This class returns the analytic solution as a function of time * @p t and position @p x. * @@ -43,6 +50,9 @@ namespace ryujin using ScalarNumber = typename View::ScalarNumber; + /* Primitive state for 2 species: (Y_0, rho, u, p) */ + static constexpr unsigned int primitive_dim = 4; + using primitive_state_type = dealii::Tensor<1, primitive_dim, Number>; ExactRiemannSolution(const HyperbolicSystem &hyperbolic_system, const std::string subsection) @@ -50,6 +60,10 @@ namespace ryujin subsection) , hyperbolic_system_(hyperbolic_system) { + Assert(n_species == 2, + dealii::ExcMessage( + "ExactRiemannSolution is only implemented for 2 species")); + primitive_left_[0] = 0.5; primitive_left_[1] = 1.4; primitive_left_[2] = 0.0; @@ -122,7 +136,7 @@ namespace ryujin const Number xi = x / t; - dealii::Tensor<1, 4, Number> primitive_state; + primitive_state_type primitive_state; if (t < 1.e-14 && x < 0.) { primitive_state = primitive_left_; @@ -201,8 +215,8 @@ namespace ryujin */ //@{ - dealii::Tensor<1, 4, Number> primitive_left_; - dealii::Tensor<1, 4, Number> primitive_right_; + primitive_state_type primitive_left_; + primitive_state_type primitive_right_; //@} /** @@ -226,10 +240,10 @@ namespace ryujin //@{ // Used to transform data to (alpha_rho_0, alpha_rho_1, u, p) - dealii::Tensor<1, 4, Number> - transform(dealii::Tensor<1, 4, Number> &temp_in) const + primitive_state_type + transform(primitive_state_type &temp_in) const { - dealii::Tensor<1, 4, Number> result; + primitive_state_type result; result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; @@ -240,7 +254,7 @@ namespace ryujin } Number fZofP(const Number &p_in, - const dealii::Tensor<1, 4, Number> &data_in) const + const primitive_state_type &data_in) const { // Need hyperbolic system for particular quantities const auto view = hyperbolic_system_.template view(); @@ -270,7 +284,7 @@ namespace ryujin Number dfZofP(const Number &p_in, - const dealii::Tensor<1, 4, Number> &data_in) const + const primitive_state_type &data_in) const { // Need hyperbolic system for particular quantities const auto view = hyperbolic_system_.template view(); @@ -305,16 +319,16 @@ namespace ryujin Number dphi(const Number &p_in, - const dealii::Tensor<1, 4, Number> &data_left, - const dealii::Tensor<1, 4, Number> &data_right) const + const primitive_state_type &data_left, + const primitive_state_type &data_right) const { return dfZofP(p_in, data_left) + dfZofP(p_in, data_right); } Number phi(const Number &p_in, - const dealii::Tensor<1, 4, Number> &data_left, - const dealii::Tensor<1, 4, Number> &data_right) const + const primitive_state_type &data_left, + const primitive_state_type &data_right) const { const Number u_L = data_left[2]; const Number u_R = data_right[2]; @@ -324,7 +338,7 @@ namespace ryujin Number lambda(const Number &p_in, - const dealii::Tensor<1, 4, Number> &data_in, + const primitive_state_type &data_in, const Number &sign) const { // Need hyperbolic system for particular quantities @@ -347,7 +361,7 @@ namespace ryujin Number lambda_intermediate(const Number &p_in, - const dealii::Tensor<1, 4, Number> &data_in, + const primitive_state_type &data_in, const Number &sign) const { // Need hyperbolic system for particular quantities @@ -378,10 +392,10 @@ namespace ryujin } - dealii::Tensor<1, 4, Number> + primitive_state_type cstar_solution(const Number &p_star, const Number &u_star, - const dealii::Tensor<1, 4, Number> &data_in) const + const primitive_state_type &data_in) const { // Need hyperbolic system for particular quantities const auto view = hyperbolic_system_.template view(); @@ -411,10 +425,10 @@ namespace ryujin } - dealii::Tensor<1, 4, Number> + primitive_state_type expansion_solution(const Number & /*p_star*/, const Number &xi, - const dealii::Tensor<1, 4, Number> &data_in, + const primitive_state_type &data_in, const Number &sign) const { // Need hyperbolic system for particular quantities @@ -460,8 +474,8 @@ namespace ryujin */ double compute_pstar(double p_1, double p_2, - dealii::Tensor<1, 4, Number> data_1, - dealii::Tensor<1, 4, Number> data_2) + primitive_state_type data_1, + primitive_state_type data_2) { constexpr Number eps = std::numeric_limits::epsilon(); diff --git a/source/multi_euler_ideal/initial_state_function.h b/source/multi_euler_ideal/initial_state_function.h index 8d8bd802..ed296ce1 100644 --- a/source/multi_euler_ideal/initial_state_function.h +++ b/source/multi_euler_ideal/initial_state_function.h @@ -1,20 +1,29 @@ // // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // Copyright (C) 2023 - 2024 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC // #pragma once +#include "hyperbolic_system.h" + #include #include +#include + namespace ryujin { namespace MultiSpeciesEulerInitialStates { + using namespace MultiSpeciesEuler; + /** - * Returns an initial state defined by a set of user specified function. + * Returns an initial state defined by a set of user specified functions. + * Supports n_species mass fractions, where the last species mass fraction + * is computed as Y_{n-1} = 1 - sum(Y_k) for k=0,...,n-2. * * @ingroup MultiSpeciesEulerEquations */ @@ -32,11 +41,16 @@ namespace ryujin : InitialState("function", subsection) , hyperbolic_system_(hyperbolic_system) { - Y0_expression_ = "1.0"; - this->add_parameter("Y_0 expression", - Y0_expression_, - "A function expression describing the mass " - "fraction for species 1"); + /* Initialize mass fraction expressions for species 0 to n-2 */ + for (unsigned int k = 0; k < n_species - 1; ++k) { + Y_expressions_[k] = + std::to_string(Number(1.) / Number(n_species)); /* default */ + this->add_parameter( + "Y_" + std::to_string(k) + " expression", + Y_expressions_[k], + "A function expression describing the mass fraction for species " + + std::to_string(k)); + } density_expression_ = "1.4"; this->add_parameter("density expression", @@ -81,7 +95,9 @@ namespace ryujin * parser with support for a time-dependent description involving * a variable »t«: */ - Y0_function_ = std::make_unique(Y0_expression_); + for (unsigned int k = 0; k < n_species - 1; ++k) + Y_functions_[k] = std::make_unique(Y_expressions_[k]); + density_function_ = std::make_unique(density_expression_); velocity_x_function_ = std::make_unique(velocity_x_expression_); if constexpr (dim > 1) @@ -101,29 +117,36 @@ namespace ryujin state_type full_primitive_state; - Y0_function_->set_time(t); density_function_->set_time(t); - - full_primitive_state[0] = Y0_function_->value(point); - full_primitive_state[0] *= density_function_->value(point); - - full_primitive_state[1] = 1. - Y0_function_->value(point); - full_primitive_state[1] *= density_function_->value(point); + const Number rho = density_function_->value(point); + + /* Compute partial densities for species 0 to n-2 */ + Number Y_sum = Number(0.); + for (unsigned int k = 0; k < n_species - 1; ++k) { + Y_functions_[k]->set_time(t); + const Number Y_k = Y_functions_[k]->value(point); + full_primitive_state[k] = Y_k * rho; + Y_sum += Y_k; + } + /* Last species gets the remaining mass fraction */ + full_primitive_state[n_species - 1] = (Number(1.) - Y_sum) * rho; velocity_x_function_->set_time(t); - full_primitive_state[2] = velocity_x_function_->value(point); + full_primitive_state[n_species] = velocity_x_function_->value(point); if constexpr (dim > 1) { velocity_y_function_->set_time(t); - full_primitive_state[3] = velocity_y_function_->value(point); + full_primitive_state[n_species + 1] = + velocity_y_function_->value(point); } if constexpr (dim > 2) { velocity_z_function_->set_time(t); - full_primitive_state[4] = velocity_z_function_->value(point); + full_primitive_state[n_species + 2] = + velocity_z_function_->value(point); } pressure_function_->set_time(t); - full_primitive_state[2 + dim] = pressure_function_->value(point); + full_primitive_state[n_species + dim] = pressure_function_->value(point); return view.from_primitive_state(full_primitive_state); } @@ -131,14 +154,15 @@ namespace ryujin private: const HyperbolicSystem &hyperbolic_system_; - std::string Y0_expression_; + std::array Y_expressions_; std::string density_expression_; std::string velocity_x_expression_; std::string velocity_y_expression_; std::string velocity_z_expression_; std::string pressure_expression_; - std::unique_ptr> Y0_function_; + std::array>, n_species - 1> + Y_functions_; std::unique_ptr> density_function_; std::unique_ptr> velocity_x_function_; std::unique_ptr> velocity_y_function_; diff --git a/source/multi_euler_ideal/initial_state_icf_like.h b/source/multi_euler_ideal/initial_state_icf_like.h index 3bf1338b..6c85ffb8 100644 --- a/source/multi_euler_ideal/initial_state_icf_like.h +++ b/source/multi_euler_ideal/initial_state_icf_like.h @@ -7,6 +7,8 @@ #pragma once +#include "hyperbolic_system.h" + #include #include @@ -14,8 +16,10 @@ namespace ryujin { namespace MultiSpeciesEulerInitialStates { + using namespace MultiSpeciesEuler; + /** - * An initial state that simulates an "intertial confinement fusion" like + * An initial state that simulates an "inertial confinement fusion" like * problem. The set up consists of three regions: (i) a low density state * inside a perturbed interface; (ii) a high density state outside the * interface; (iii) an incoming shock wave characterized by its Mach number @@ -34,29 +38,33 @@ namespace ryujin typename Description::template HyperbolicSystemView; using state_type = typename View::state_type; + static constexpr unsigned int primitive_dim = n_species + 2; + using primitive_state_type = dealii::Tensor<1, primitive_dim, Number>; + ICFLike(const HyperbolicSystem &hyperbolic_system, const std::string subsection) : InitialState("icf like", subsection) , hyperbolic_system_(hyperbolic_system) { - - temp_inside_[0] = 0.5; - temp_inside_[1] = 1.0; - temp_inside_[2] = 0.0; - temp_inside_[3] = 1.0; + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_inside_[k] = Number(1.) / Number(n_species); + temp_inside_[n_species - 1] = 1.0; + temp_inside_[n_species] = 0.0; + temp_inside_[n_species + 1] = 1.0; this->add_parameter("primitive state inside", temp_inside_, - "Initial primitive state (Y_0, rho, u, p) inside " - "perturbed interface"); - - temp_outside_[0] = 0.5; - temp_outside_[1] = 1.0; - temp_outside_[2] = 0.0; - temp_outside_[3] = 1.0; + "Initial primitive state (Y_0, ..., Y_{n-2}, rho, " + "u, p) inside perturbed interface"); + + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_outside_[k] = Number(1.) / Number(n_species); + temp_outside_[n_species - 1] = 1.0; + temp_outside_[n_species] = 0.0; + temp_outside_[n_species + 1] = 1.0; this->add_parameter("primitive state outside", temp_outside_, - "Initial primitive state (Y_0, rho, u, p) outside " - "perturbed interface"); + "Initial primitive state (Y_0, ..., Y_{n-2}, rho, " + "u, p) outside perturbed interface"); interface_radius_ = 1.0; this->add_parameter( @@ -65,11 +73,11 @@ namespace ryujin num_modes_ = 8.0; this->add_parameter("number of modes", num_modes_, - "Number of modes for pertburation of interface"); + "Number of modes for perturbation of interface"); amplitude_ = 0.02; this->add_parameter( - "amplitude", amplitude_, "Amplitude for interface pertburation"); + "amplitude", amplitude_, "Amplitude for interface perturbation"); mach_number_ = 3.0; this->add_parameter( @@ -81,8 +89,8 @@ namespace ryujin "Radial location of incoming shock front"); const auto convert_states = [&]() { - const auto prim_inside = extend(temp_inside_); - const auto prim_outside = extend(temp_outside_); + const auto prim_inside = extend_primitive(temp_inside_); + const auto prim_outside = extend_primitive(temp_outside_); const auto view = hyperbolic_system_.template view(); state_inside_ = view.from_initial_state(prim_inside); @@ -102,9 +110,9 @@ namespace ryujin { const auto gamma_ = view.gamma_mixture(state_outside_); - const auto &rho_R = temp_outside_[1]; - const auto &u_R = temp_outside_[2]; - const auto &p_R = temp_outside_[3]; + const auto &rho_R = temp_outside_[n_species - 1]; + const auto &u_R = temp_outside_[n_species]; + const auto &p_R = temp_outside_[n_species + 1]; const Number a_R = std::sqrt(gamma_ * p_R / rho_R); const Number mach_R = u_R / a_R; @@ -123,21 +131,26 @@ namespace ryujin (gamma_ + Number(1.)); state_type primitive_shock_state; - const auto Y0_outside = temp_outside_[0]; - primitive_shock_state[0] = Y0_outside * rho_L; - primitive_shock_state[1] = (1. - Y0_outside) * rho_L; + /* Compute partial densities from mass fractions */ + Number Y_sum = Number(0.); + for (unsigned int k = 0; k < n_species - 1; ++k) { + primitive_shock_state[k] = temp_outside_[k] * rho_L; + Y_sum += temp_outside_[k]; + } + primitive_shock_state[n_species - 1] = (Number(1.) - Y_sum) * rho_L; + /* Set radial velocity */ for (unsigned int i = 0; i < dim; ++i) { - primitive_shock_state[i + 2] = 0.; + primitive_shock_state[n_species + i] = Number(0.); } if (point.norm() > 0.) { for (unsigned int i = 0; i < dim; ++i) { - primitive_shock_state[i + 2] = -u_L * r_hat[i]; + primitive_shock_state[n_species + i] = -u_L * r_hat[i]; } } - primitive_shock_state[2 + dim] = p_L; + primitive_shock_state[n_species + dim] = p_L; conserved_shock_state = view.from_initial_state(primitive_shock_state); @@ -174,10 +187,8 @@ namespace ryujin private: const HyperbolicSystem &hyperbolic_system_; - Number gamma_; - - dealii::Tensor<1, 4, Number> temp_inside_; - dealii::Tensor<1, 4, Number> temp_outside_; + primitive_state_type temp_inside_; + primitive_state_type temp_outside_; state_type state_inside_; state_type state_outside_; @@ -188,18 +199,6 @@ namespace ryujin double shock_radius_; double mach_number_; - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> - extend(dealii::Tensor<1, 4, Number> &temp_in) const - { - dealii::Tensor<1, 4, Number> result; - result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; - result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; - - for (unsigned int i = 2; i < 4; ++i) - result[i] = temp_in[i]; - - return result; - } }; diff --git a/source/multi_euler_ideal/initial_state_library_multi_euler.h b/source/multi_euler_ideal/initial_state_library_multi_euler.h index 059e56b3..39b9b14f 100644 --- a/source/multi_euler_ideal/initial_state_library_multi_euler.h +++ b/source/multi_euler_ideal/initial_state_library_multi_euler.h @@ -1,10 +1,13 @@ // // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // Copyright (C) 2023 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC // #pragma once +#include "hyperbolic_system.h" + #include #include "initial_state_contrast.h" @@ -21,6 +24,8 @@ namespace ryujin { namespace MultiSpeciesEulerInitialStates { + using namespace MultiSpeciesEuler; + template void populate_initial_state_list( typename ryujin::InitialStateLibrary:: @@ -33,8 +38,12 @@ namespace ryujin }; add(std::make_unique>(h, s)); - add(std::make_unique>(h, - s)); + + /* ExactRiemannSolution only supports 2 species */ + if constexpr (n_species == 2) + add(std::make_unique>(h, + s)); + add(std::make_unique>(h, s)); add(std::make_unique>(h, s)); add(std::make_unique>(h, s)); diff --git a/source/multi_euler_ideal/initial_state_radial_contrast.h b/source/multi_euler_ideal/initial_state_radial_contrast.h index c7bd2676..bd22b736 100644 --- a/source/multi_euler_ideal/initial_state_radial_contrast.h +++ b/source/multi_euler_ideal/initial_state_radial_contrast.h @@ -5,12 +5,16 @@ #pragma once +#include "hyperbolic_system.h" + #include namespace ryujin { namespace MultiSpeciesEulerInitialStates { + using namespace MultiSpeciesEuler; + /** * A modification of the "contrast" initial state. Now, we have an * initial state formed by a contrast of a given "left" and "right" @@ -31,6 +35,10 @@ namespace ryujin typename Description::template HyperbolicSystemView; using state_type = typename View::state_type; + /* 1D primitive state: (Y_0, ..., Y_{n-2}, rho, u, p) */ + static constexpr unsigned int primitive_dim = n_species + 2; + using primitive_state_type = dealii::Tensor<1, primitive_dim, Number>; + RadialContrast(const HyperbolicSystem &hyperbolic_system, const std::string &subsection) : InitialState("radial contrast", @@ -44,23 +52,28 @@ namespace ryujin "If set to true, we transform a non-zero velocity into a radial " "velocity with scaling 1 / r^(dim - 1)"); - temp_left_[0] = 0.5; - temp_left_[1] = 1.4; - temp_left_[2] = 0.0; - temp_left_[3] = 1.0; + /* Default: equal mass fractions */ + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_left_[k] = Number(1.) / Number(n_species); + temp_left_[n_species - 1] = 1.4; /* rho */ + temp_left_[n_species] = 0.0; /* u */ + temp_left_[n_species + 1] = 1.0; /* p */ this->add_parameter( "primitive state left", temp_left_, - "Initial 1d primitive state (Y_0, rho, u, p) inside radial area"); - - temp_right_[0] = 0.5; - temp_right_[1] = 1.4; - temp_right_[2] = 0.0; - temp_right_[3] = 1.0; + "Initial 1d primitive state (Y_0, ..., Y_{n-2}, rho, u, p) inside " + "radial area"); + + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_right_[k] = Number(1.) / Number(n_species); + temp_right_[n_species - 1] = 1.4; /* rho */ + temp_right_[n_species] = 0.0; /* u */ + temp_right_[n_species + 1] = 1.0; /* p */ this->add_parameter( "primitive state right", temp_right_, - "Initial 1d primitive state (Y_0, rho, u, p) outside radial area"); + "Initial 1d primitive state (Y_0, ..., Y_{n-2}, rho, u, p) outside " + "radial area"); radius_ = 1.0; this->add_parameter("radius", radius_, "Radius of radial area"); @@ -68,8 +81,8 @@ namespace ryujin const auto convert_states = [&]() { const auto view = hyperbolic_system_.template view(); - const auto primitive_left_ = extend(temp_left_); - const auto primitive_right_ = extend(temp_right_); + const auto primitive_left_ = extend_primitive(temp_left_); + const auto primitive_right_ = extend_primitive(temp_right_); state_left_ = view.from_initial_state(primitive_left_); state_right_ = view.from_initial_state(primitive_right_); @@ -81,25 +94,27 @@ namespace ryujin state_type compute(const dealii::Point &point, Number /*t*/) final { - if (point.norm() > 0. && use_radial_velocity_) { - state_left_[2] = temp_left_[2] * point[0] / point.norm(); - state_left_[3] = temp_left_[2] * point[1] / point.norm(); + auto result = + (point.norm() > radius_ ? state_right_ : state_left_); - state_left_[2] = temp_right_[2] * point[0] / point.norm(); - state_left_[3] = temp_right_[2] * point[1] / point.norm(); + if (point.norm() > 0. && use_radial_velocity_) { + const auto view = hyperbolic_system_.template view(); + const auto &temp = + (point.norm() > radius_ ? temp_right_ : temp_left_); + const auto rho = view.density(result); + const auto u = temp[n_species]; /* scalar velocity */ + for (unsigned int d = 0; d < dim; ++d) + result[n_species + d] = rho * u * point[d] / point.norm(); } - auto final_state = - (point.norm() > radius_ ? state_right_ : state_left_); - - return final_state; + return result; } private: const HyperbolicSystem &hyperbolic_system_; - dealii::Tensor<1, 4, Number> temp_left_; - dealii::Tensor<1, 4, Number> temp_right_; + primitive_state_type temp_left_; + primitive_state_type temp_right_; state_type state_left_; state_type state_right_; @@ -107,19 +122,6 @@ namespace ryujin double radius_; bool use_radial_velocity_; - - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> - extend(dealii::Tensor<1, 4, Number> &temp_in) const - { - dealii::Tensor<1, 4, Number> result; - result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; - result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; - - for (unsigned int i = 2; i < 4; ++i) - result[i] = temp_in[i]; - - return result; - } }; diff --git a/source/multi_euler_ideal/initial_state_shock_bubble.h b/source/multi_euler_ideal/initial_state_shock_bubble.h index 37979ed7..db4f95f0 100644 --- a/source/multi_euler_ideal/initial_state_shock_bubble.h +++ b/source/multi_euler_ideal/initial_state_shock_bubble.h @@ -6,12 +6,16 @@ #pragma once +#include "hyperbolic_system.h" + #include namespace ryujin { namespace MultiSpeciesEulerInitialStates { + using namespace MultiSpeciesEuler; + /** * A shock-bubble initial state. The shock HAS to be initiated for when * either the mass fraction is 0 or 1 ONLY. The implementation assumes we @@ -31,28 +35,35 @@ namespace ryujin typename Description::template HyperbolicSystemView; using state_type = typename View::state_type; + static constexpr unsigned int primitive_dim = n_species + 2; + using primitive_state_type = dealii::Tensor<1, primitive_dim, Number>; + ShockBubble(const HyperbolicSystem &hyperbolic_system, const std::string subsection) : InitialState("shock bubble", subsection) , hyperbolic_system_(hyperbolic_system) { - temp_inside_[0] = 0.5; - temp_inside_[1] = 1.0; - temp_inside_[2] = 0.0; - temp_inside_[3] = 1.0; + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_inside_[k] = Number(1.) / Number(n_species); + temp_inside_[n_species - 1] = 1.0; + temp_inside_[n_species] = 0.0; + temp_inside_[n_species + 1] = 1.0; this->add_parameter( "primitive state bubble", temp_inside_, - "Initial 1d primitive state (Y_0, rho, u, p) inside bubble"); - - temp_ambient_[0] = 0.5; - temp_ambient_[1] = 1.0; - temp_ambient_[2] = 0.0; - temp_ambient_[3] = 1.0; + "Initial 1d primitive state (Y_0, ..., Y_{n-2}, rho, u, p) inside " + "bubble"); + + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_ambient_[k] = Number(1.) / Number(n_species); + temp_ambient_[n_species - 1] = 1.0; + temp_ambient_[n_species] = 0.0; + temp_ambient_[n_species + 1] = 1.0; this->add_parameter( "primitive state pre-shock", temp_ambient_, - "Initial 1d primitive state (Y_0, rho, u, p) ambient pre-shock"); + "Initial 1d primitive state (Y_0, ..., Y_{n-2}, rho, u, p) ambient " + "pre-shock"); radius_ = 1.0; this->add_parameter("radius", radius_, "Bubble radius"); @@ -77,36 +88,20 @@ namespace ryujin const auto convert_states = [&]() { const auto view = hyperbolic_system_.template view(); - primitive_inside_ = extend(temp_inside_); - primitive_ambient_ = extend(temp_ambient_); + const auto primitive_inside = extend_primitive(temp_inside_); + const auto primitive_ambient = extend_primitive(temp_ambient_); - state_bubble_ = view.from_initial_state(primitive_inside_); - state_ambient_ = view.from_initial_state(primitive_ambient_); - }; - this->parse_parameters_call_back.connect(convert_states); - convert_states(); - } - - - state_type compute(const dealii::Point &point, Number t) final - { - const auto view = hyperbolic_system_.template view(); + state_bubble_ = view.from_initial_state(primitive_inside); + state_ambient_ = view.from_initial_state(primitive_ambient); - const auto &x = point[0]; - const auto xbar = point - bubble_center_; - const auto &xc = bubble_center_[0]; - auto sign = 1.; - if (xc < shock_location_) - sign = -1.; + /* Precompute shocked state via Rankine-Hugoniot: */ - /* Compute shocked state outside */ - { const auto gamma_ = view.gamma_mixture(state_ambient_); - const auto &rho_R = temp_ambient_[1]; - const auto &u_R = temp_ambient_[2]; - const auto &p_R = temp_ambient_[3]; - /* a_R^2 = gamma * p / rho / (1 - b * rho) */ + const auto &rho_R = temp_ambient_[n_species - 1]; + const auto &u_R = temp_ambient_[n_species]; + const auto &p_R = temp_ambient_[n_species + 1]; + const Number a_R = std::sqrt(gamma_ * p_R / rho_R); const Number mach_R = u_R / a_R; @@ -123,36 +118,42 @@ namespace ryujin (gamma_ - Number(1.))) / (gamma_ + Number(1.)); - const auto Y0_outside = temp_ambient_[0]; + auto sign = 1.; + if (bubble_center_[0] < shock_location_) + sign = -1.; - primitive_shock_[0] = Y0_outside; - primitive_shock_[1] = rho_L; - primitive_shock_[2] = sign * u_L; - primitive_shock_[3] = p_L; + primitive_state_type prim_shock; + for (unsigned int k = 0; k < n_species - 1; ++k) + prim_shock[k] = temp_ambient_[k]; + prim_shock[n_species - 1] = rho_L; + prim_shock[n_species] = sign * u_L; + prim_shock[n_species + 1] = p_L; - primitive_shock_ = extend(primitive_shock_); - state_shock_ = view.from_initial_state(primitive_shock_); - } + state_shock_ = view.from_initial_state(extend_primitive(prim_shock)); + }; + this->parse_parameters_call_back.connect(convert_states); + convert_states(); + } + + + state_type compute(const dealii::Point &point, Number /*t*/) final + { + const auto xbar = point - bubble_center_; auto final_state = (xbar.norm() <= radius_ ? state_bubble_ : state_ambient_); - if (x <= xc + shock_location_) + if (point[0] <= bubble_center_[0] + shock_location_) final_state = state_shock_; - /* Set final state */ return final_state; } private: const HyperbolicSystem &hyperbolic_system_; - dealii::Tensor<1, 4, Number> temp_inside_; - dealii::Tensor<1, 4, Number> temp_ambient_; - - dealii::Tensor<1, 4, Number> primitive_inside_; - dealii::Tensor<1, 4, Number> primitive_ambient_; - dealii::Tensor<1, 4, Number> primitive_shock_; + primitive_state_type temp_inside_; + primitive_state_type temp_ambient_; state_type state_bubble_; state_type state_ambient_; @@ -164,18 +165,6 @@ namespace ryujin double mach_number_; double shock_location_; - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> - extend(dealii::Tensor<1, 4, Number> &temp_in) const - { - dealii::Tensor<1, 4, Number> result; - result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; - result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; - - for (unsigned int i = 2; i < 4; ++i) - result[i] = temp_in[i]; - - return result; - } }; diff --git a/source/multi_euler_ideal/initial_state_smooth_wave.h b/source/multi_euler_ideal/initial_state_smooth_wave.h index 626cd41a..2eaef6cb 100644 --- a/source/multi_euler_ideal/initial_state_smooth_wave.h +++ b/source/multi_euler_ideal/initial_state_smooth_wave.h @@ -5,6 +5,8 @@ #pragma once +#include "hyperbolic_system.h" + #include #include @@ -12,12 +14,17 @@ namespace ryujin { namespace MultiSpeciesEulerInitialStates { + using namespace MultiSpeciesEuler; + /** * This is a generalization of the "Smooth traveling wave" problem first * proposed in Section 5.2 of @cite GuermondEtAl2018 for ideal gas EOS. * The details for extending to multi-species Euler Equations with each * species following an ideal gas is soon to be documented. * + * For n_species > 2, the first species has mass fraction Y_0, and the + * remaining mass (1 - Y_0) is distributed equally among the other species. + * * @note This class returns the analytic solution as a function of time * @p t and position @p x. * @@ -34,6 +41,10 @@ namespace ryujin using ScalarNumber = typename View::ScalarNumber; + /* 1D primitive state for from_initial_state */ + static constexpr unsigned int primitive_dim = n_species + 2; + using primitive_state_type = dealii::Tensor<1, primitive_dim, Number>; + SmoothWave(const HyperbolicSystem &hyperbolic_system, const std::string subsection) : InitialState("smooth wave", subsection) @@ -42,7 +53,9 @@ namespace ryujin mass_fraction_ = 0.5; this->add_parameter("mass fraction", mass_fraction_, - "The mass fraction of first species"); + "The mass fraction of first species (Y_0). " + "Remaining mass is split equally among other " + "species."); density_ref_ = 1.; this->add_parameter("reference density", @@ -82,14 +95,23 @@ namespace ryujin if (left_ <= point_bar[0] && point_bar[0] <= right_) rho = density_ref_ + polynomial; - state_type initial_state; - { - initial_state[0] = mass_fraction_ * rho; // alpha_1 rho_1 - initial_state[1] = (1. - mass_fraction_) * rho; // alpha_2 rho_2 - initial_state[2] = mach_number_; - initial_state[3] = pressure_ref_; + /* Build the primitive state for from_initial_state */ + primitive_state_type initial_state; + + /* Species 0 gets mass_fraction_, the rest share (1 - mass_fraction_) */ + initial_state[0] = mass_fraction_ * rho; + if constexpr (n_species == 2) { + initial_state[1] = (Number(1.) - mass_fraction_) * rho; + } else { + const Number remaining_fraction = + (Number(1.) - mass_fraction_) / Number(n_species - 1); + for (unsigned int k = 1; k < n_species; ++k) + initial_state[k] = remaining_fraction * rho; } + initial_state[n_species] = mach_number_; /* u */ + initial_state[n_species + 1] = pressure_ref_; /* p */ + return view.from_initial_state(initial_state); } diff --git a/source/multi_euler_ideal/initial_state_three_state_contrast.h b/source/multi_euler_ideal/initial_state_three_state_contrast.h index 088cc934..0739e9a2 100644 --- a/source/multi_euler_ideal/initial_state_three_state_contrast.h +++ b/source/multi_euler_ideal/initial_state_three_state_contrast.h @@ -6,12 +6,16 @@ #pragma once +#include "hyperbolic_system.h" + #include namespace ryujin { namespace MultiSpeciesEulerInitialStates { + using namespace MultiSpeciesEuler; + /** * An initial state formed by two contrasts of given "left", "middle" * and "right" primitive states. The user defines the lengths of the left @@ -34,56 +38,66 @@ namespace ryujin typename Description::template HyperbolicSystemView; using state_type = typename View::state_type; + /* 1D primitive state: (Y_0, ..., Y_{n-2}, rho, u, p) */ + static constexpr unsigned int primitive_dim = n_species + 2; + using primitive_state_type = dealii::Tensor<1, primitive_dim, Number>; + ThreeStateContrast(const HyperbolicSystem &hyperbolic_system, const std::string &subsection) : InitialState("three state contrast", subsection) , hyperbolic_system_(hyperbolic_system) { - - temp_left_[0] = 0.5; - temp_left_[1] = 1.; - temp_left_[2] = 0.; - temp_left_[3] = 1.e3; + /* Default: equal mass fractions */ + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_left_[k] = Number(1.) / Number(n_species); + temp_left_[n_species - 1] = 1.; /* rho */ + temp_left_[n_species] = 0.; /* u */ + temp_left_[n_species + 1] = 1.e3; /* p */ this->add_parameter( "primitive state left", temp_left_, - "Initial 1d primitive state (Y_0, rho, u, p) on the left"); + "Initial 1d primitive state (Y_0, ..., Y_{n-2}, rho, u, p) on the " + "left"); left_length_ = 0.1; this->add_parameter("left region length", left_length_, "The length of the left region"); - temp_middle_[0] = 0.5; - temp_middle_[1] = 1.; - temp_middle_[2] = 0.; - temp_middle_[3] = 1.e-2; + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_middle_[k] = Number(1.) / Number(n_species); + temp_middle_[n_species - 1] = 1.; /* rho */ + temp_middle_[n_species] = 0.; /* u */ + temp_middle_[n_species + 1] = 1.e-2; /* p */ this->add_parameter( "primitive state middle", temp_middle_, - "Initial 1d primitive state (Y_0, rho, u, p) in the middle"); + "Initial 1d primitive state (Y_0, ..., Y_{n-2}, rho, u, p) in the " + "middle"); middle_length_ = 0.8; this->add_parameter("middle region length", middle_length_, "The length of the middle region"); - temp_right_[0] = 0.5; - temp_right_[1] = 1.; - temp_right_[2] = 0.; - temp_right_[3] = 1.e2; + for (unsigned int k = 0; k < n_species - 1; ++k) + temp_right_[k] = Number(1.) / Number(n_species); + temp_right_[n_species - 1] = 1.; /* rho */ + temp_right_[n_species] = 0.; /* u */ + temp_right_[n_species + 1] = 1.e2; /* p */ this->add_parameter( "primitive state right", temp_right_, - "Initial 1d primitive state (Y_0, rho, u, p) on the right"); + "Initial 1d primitive state (Y_0, ..., Y_{n-2}, rho, u, p) on the " + "right"); const auto convert_states = [&]() { const auto view = hyperbolic_system_.template view(); - const auto primitive_left_ = extend(temp_left_); - const auto primitive_middle_ = extend(temp_middle_); - const auto primitive_right_ = extend(temp_right_); + const auto primitive_left_ = extend_primitive(temp_left_); + const auto primitive_middle_ = extend_primitive(temp_middle_); + const auto primitive_right_ = extend_primitive(temp_right_); state_left_ = view.from_initial_state(primitive_left_); state_middle_ = view.from_initial_state(primitive_middle_); @@ -106,26 +120,14 @@ namespace ryujin Number left_length_; Number middle_length_; - dealii::Tensor<1, 4, Number> temp_left_; - dealii::Tensor<1, 4, Number> temp_middle_; - dealii::Tensor<1, 4, Number> temp_right_; + primitive_state_type temp_left_; + primitive_state_type temp_middle_; + primitive_state_type temp_right_; state_type state_left_; state_type state_middle_; state_type state_right_; - DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, 4, Number> - extend(dealii::Tensor<1, 4, Number> &temp_in) const - { - dealii::Tensor<1, 4, Number> result; - result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; - result[1] = (1. - temp_in[0]) * temp_in[1]; // = alpha_1 rho_1; - - for (unsigned int i = 2; i < 4; ++i) - result[i] = temp_in[i]; - - return result; - } }; } // namespace MultiSpeciesEulerInitialStates } // namespace ryujin diff --git a/source/multi_euler_ideal/limiter.h b/source/multi_euler_ideal/limiter.h index 6bd20f9e..a6899f86 100644 --- a/source/multi_euler_ideal/limiter.h +++ b/source/multi_euler_ideal/limiter.h @@ -125,9 +125,10 @@ namespace ryujin // //@{ /** - * The number of stored entries in the bounds array. + * The number of stored entries in the bounds array: + * 2 * n_species (min/max per species) + 1 (rho_e_min) + 1 (s_min) */ - static constexpr unsigned int n_bounds = 6; + static constexpr unsigned int n_bounds = 2 * n_species + 2; /** * Array type used to store accumulated bounds. @@ -252,8 +253,7 @@ namespace ryujin Bounds bounds_; - Number alpha_rho_0_relaxation_numerator; - Number alpha_rho_1_relaxation_numerator; + dealii::Tensor<1, n_species, Number> rho_relaxation_numerator; Number rho_relaxation_denominator; Number rho_e_interp_max; Number s_interp_max; @@ -281,18 +281,18 @@ namespace ryujin const auto &[rho_i, p_i, gamma_min_i, s_i, harten_i] = precomputed_values.template get_tensor(i); - const auto alpha_rho_0 = view.partial_density_0(U_i); - const auto alpha_rho_1 = view.partial_density_1(U_i); - - const auto rho_e = view.internal_energy(U_i); + Bounds result; + /* Partial density bounds: [rho_k_min, rho_k_max] for each species */ + for (unsigned int k = 0; k < n_species; ++k) { + const auto alpha_rho_k = view.partial_density(U_i, k); + result[2 * k] = alpha_rho_k; /* min */ + result[2 * k + 1] = alpha_rho_k; /* max */ + } + /* Internal energy and entropy bounds */ + result[2 * n_species] = view.internal_energy(U_i); /* rho_e_min */ + result[2 * n_species + 1] = s_i; /* s_min */ - return {/*alpha_rho_0_min*/ - alpha_rho_0, - /*alpha_rho_0_max*/ alpha_rho_0, - /*alpha_rho_1_min*/ alpha_rho_1, - /*alpha_rho_1_min*/ alpha_rho_1, - /*rho_e*/ rho_e, - /*s_min*/ s_i}; + return result; } @@ -300,25 +300,20 @@ namespace ryujin DEAL_II_ALWAYS_INLINE inline auto Limiter::combine_bounds( const Bounds &bounds_left, const Bounds &bounds_right) const -> Bounds { - const auto &[alpha_rho_0_min_l, - alpha_rho_0_max_l, - alpha_rho_1_min_l, - alpha_rho_1_max_l, - rho_e_min_l, - s_min_l] = bounds_left; - const auto &[alpha_rho_0_min_r, - alpha_rho_0_max_r, - alpha_rho_1_min_r, - alpha_rho_1_max_r, - rho_e_min_r, - s_min_r] = bounds_right; - - return {std::min(alpha_rho_0_min_l, alpha_rho_0_min_r), - std::max(alpha_rho_0_max_l, alpha_rho_0_max_r), - std::min(alpha_rho_1_min_l, alpha_rho_1_min_r), - std::max(alpha_rho_1_max_l, alpha_rho_1_max_r), - std::min(rho_e_min_l, rho_e_min_r), - std::min(s_min_l, s_min_r)}; + Bounds result; + /* Combine partial density bounds for each species */ + for (unsigned int k = 0; k < n_species; ++k) { + result[2 * k] = std::min(bounds_left[2 * k], bounds_right[2 * k]); + result[2 * k + 1] = + std::max(bounds_left[2 * k + 1], bounds_right[2 * k + 1]); + } + /* Combine internal energy and entropy bounds */ + result[2 * n_species] = + std::min(bounds_left[2 * n_species], bounds_right[2 * n_species]); + result[2 * n_species + 1] = std::min(bounds_left[2 * n_species + 1], + bounds_right[2 * n_species + 1]); + + return result; } @@ -328,12 +323,6 @@ namespace ryujin const Number &hd) const -> Bounds { auto relaxed_bounds = bounds; - auto &[alpha_rho_0_min_relaxed, - alpha_rho_0_max_relaxed, - alpha_rho_1_min_relaxed, - alpha_rho_1_max_relaxed, - rho_e_min_relaxed, - s_min_relaxed] = relaxed_bounds; /* Use r = factor * (m_i / |Omega|) ^ (1.5 / d): */ @@ -346,14 +335,17 @@ namespace ryujin constexpr ScalarNumber eps = std::numeric_limits::epsilon(); - alpha_rho_0_min_relaxed *= std::max(Number(1.) - r, Number(eps)); - alpha_rho_0_max_relaxed *= (Number(1.) + r); - - alpha_rho_1_min_relaxed *= std::max(Number(1.) - r, Number(eps)); - alpha_rho_1_max_relaxed *= (Number(1.) + r); + /* Relax partial density bounds for each species */ + for (unsigned int k = 0; k < n_species; ++k) { + relaxed_bounds[2 * k] *= std::max(Number(1.) - r, Number(eps)); + relaxed_bounds[2 * k + 1] *= (Number(1.) + r); + } - rho_e_min_relaxed *= std::max(Number(1.) - r, Number(eps)); + /* Relax internal energy bound */ + relaxed_bounds[2 * n_species] *= std::max(Number(1.) - r, Number(eps)); + /* Relax entropy bound */ + auto &s_min_relaxed = relaxed_bounds[2 * n_species + 1]; const Number scaled_ent = std::exp(s_min_relaxed / cv_i); s_min_relaxed = cv_i * std::log((Number(1.) - r) * scaled_ent); @@ -363,7 +355,7 @@ namespace ryujin template DEAL_II_ALWAYS_INLINE inline void - Limiter::reset(const unsigned int i, + Limiter::reset(const unsigned int /* i */, const state_type &new_U_i, const flux_contribution_type &new_flux_i) { @@ -374,29 +366,18 @@ namespace ryujin cv_i = view.cv_mixture(new_U_i); - /* Bounds: */ - - auto &[alpha_rho_0_min, - alpha_rho_0_max, - alpha_rho_1_min, - alpha_rho_1_max, - rho_e_min, - s_min] = bounds_; - - alpha_rho_0_min = Number(std::numeric_limits::max()); - alpha_rho_0_max = Number(0.); - - alpha_rho_1_min = Number(std::numeric_limits::max()); - alpha_rho_1_max = Number(0.); - - rho_e_min = Number(std::numeric_limits::max()); - s_min = Number(std::numeric_limits::max()); - + /* Initialize bounds for each species */ + for (unsigned int k = 0; k < n_species; ++k) { + bounds_[2 * k] = Number(std::numeric_limits::max()); + bounds_[2 * k + 1] = Number(0.); + } + bounds_[2 * n_species] = Number(std::numeric_limits::max()); + bounds_[2 * n_species + 1] = + Number(std::numeric_limits::max()); /* Relaxation: */ - - alpha_rho_0_relaxation_numerator = Number(0.); - alpha_rho_1_relaxation_numerator = Number(0.); + for (unsigned int k = 0; k < n_species; ++k) + rho_relaxation_numerator[k] = Number(0.); rho_relaxation_denominator = Number(0.); rho_e_interp_max = Number(0.); @@ -421,20 +402,6 @@ namespace ryujin const auto view = hyperbolic_system.view(); - /* Bounds: */ - auto &[alpha_rho_0_min, - alpha_rho_0_max, - alpha_rho_1_min, - alpha_rho_1_max, - rho_e_min, - s_min] = bounds_; - - const auto alpha_rho_0_i = view.partial_density_0(U_i); - const auto alpha_rho_0_j = view.partial_density_0(U_j); - - const auto alpha_rho_1_i = view.partial_density_1(U_i); - const auto alpha_rho_1_j = view.partial_density_1(U_j); - const auto rho_e_j = view.internal_energy(U_j); const auto [rho_j, p_j, gamma_min_j, s_j, harten_j] = @@ -446,30 +413,25 @@ namespace ryujin ScalarNumber(0.5) * contract(add(flux_j, -flux_i), scaled_c_ij) + affine_shift; - const auto alpha_rho_0_ij_bar = view.partial_density_0(U_ij_bar); - const auto alpha_rho_1_ij_bar = view.partial_density_1(U_ij_bar); - const auto rho_e_ij_bar = view.internal_energy(U_ij_bar); - /* Density bounds: */ - - alpha_rho_0_min = std::min(alpha_rho_0_min, alpha_rho_0_ij_bar); - alpha_rho_0_max = std::max(alpha_rho_0_max, alpha_rho_0_ij_bar); - - alpha_rho_1_min = std::min(alpha_rho_1_min, alpha_rho_1_ij_bar); - alpha_rho_1_max = std::max(alpha_rho_1_max, alpha_rho_1_ij_bar); + /* Density bounds and relaxation for each species: */ + const auto beta_ij = Number(1.); + for (unsigned int k = 0; k < n_species; ++k) { + const auto alpha_rho_k_i = view.partial_density(U_i, k); + const auto alpha_rho_k_j = view.partial_density(U_j, k); + const auto alpha_rho_k_ij_bar = view.partial_density(U_ij_bar, k); - /* Density relaxation: */ + bounds_[2 * k] = std::min(bounds_[2 * k], alpha_rho_k_ij_bar); + bounds_[2 * k + 1] = std::max(bounds_[2 * k + 1], alpha_rho_k_ij_bar); - /* Use a uniform weight. */ - const auto beta_ij = Number(1.); - alpha_rho_0_relaxation_numerator += - beta_ij * (alpha_rho_0_i + alpha_rho_0_j); - alpha_rho_1_relaxation_numerator += - beta_ij * (alpha_rho_1_i + alpha_rho_1_j); + rho_relaxation_numerator[k] += + beta_ij * (alpha_rho_k_i + alpha_rho_k_j); + } rho_relaxation_denominator += std::abs(beta_ij); /* Internal energy bounds and relaxation */ + auto &rho_e_min = bounds_[2 * n_species]; rho_e_min = std::min(rho_e_min, rho_e_j); rho_e_min = std::min(rho_e_min, rho_e_ij_bar); @@ -478,7 +440,7 @@ namespace ryujin rho_e_interp_max = std::max(rho_e_interp_max, rho_e_interp); /* Entropy bounds and relaxation: */ - + auto &s_min = bounds_[2 * n_species + 1]; const auto s_ij_bar = view.specific_entropy(U_ij_bar); const Number s_interp = @@ -494,54 +456,37 @@ namespace ryujin DEAL_II_ALWAYS_INLINE inline auto Limiter::bounds(const Number hd_i) const -> Bounds { - const auto &[alpha_rho_0_min, - alpha_rho_0_max, - alpha_rho_1_min, - alpha_rho_1_max, - rho_e_min, - s_min] = bounds_; - auto relaxed_bounds = fully_relax_bounds(bounds_, hd_i); - auto &[alpha_rho_0_min_relaxed, - alpha_rho_0_max_relaxed, - alpha_rho_1_min_relaxed, - alpha_rho_1_max_relaxed, - rho_e_min_relaxed, - s_min_relaxed] = relaxed_bounds; /* Apply a stricter window: */ - constexpr ScalarNumber eps = std::numeric_limits::epsilon(); - const auto alpha_rho_0_relaxation = - ScalarNumber(2. * parameters.relaxation_factor()) * - std::abs(alpha_rho_0_relaxation_numerator) / - (std::abs(rho_relaxation_denominator) + Number(eps)); + /* Apply stricter bounds for each species */ + for (unsigned int k = 0; k < n_species; ++k) { + const auto rho_k_relaxation = + ScalarNumber(2. * parameters.relaxation_factor()) * + std::abs(rho_relaxation_numerator[k]) / + (std::abs(rho_relaxation_denominator) + Number(eps)); + + relaxed_bounds[2 * k] = + std::max(relaxed_bounds[2 * k], bounds_[2 * k] - rho_k_relaxation); + relaxed_bounds[2 * k + 1] = std::min( + relaxed_bounds[2 * k + 1], bounds_[2 * k + 1] + rho_k_relaxation); + } - const auto alpha_rho_1_relaxation = - ScalarNumber(2. * parameters.relaxation_factor()) * - std::abs(alpha_rho_1_relaxation_numerator) / - (std::abs(rho_relaxation_denominator) + Number(eps)); + /* Apply stricter bounds for internal energy and entropy */ + const auto &rho_e_min = bounds_[2 * n_species]; + const auto &s_min = bounds_[2 * n_species + 1]; const auto int_relaxation = parameters.relaxation_factor() * (rho_e_interp_max - rho_e_min); - const auto entropy_relaxation = parameters.relaxation_factor() * (s_interp_max - s_min); - alpha_rho_0_min_relaxed = std::max( - alpha_rho_0_min_relaxed, alpha_rho_0_min - alpha_rho_0_relaxation); - alpha_rho_0_max_relaxed = std::min( - alpha_rho_0_max_relaxed, alpha_rho_0_max + alpha_rho_0_relaxation); - - alpha_rho_1_min_relaxed = std::max( - alpha_rho_1_min_relaxed, alpha_rho_1_min - alpha_rho_1_relaxation); - alpha_rho_1_max_relaxed = std::min( - alpha_rho_1_max_relaxed, alpha_rho_1_max + alpha_rho_1_relaxation); - - rho_e_min_relaxed = - std::max(rho_e_min_relaxed, rho_e_min - int_relaxation); - s_min_relaxed = std::max(s_min_relaxed, s_min - entropy_relaxation); + relaxed_bounds[2 * n_species] = + std::max(relaxed_bounds[2 * n_species], rho_e_min - int_relaxation); + relaxed_bounds[2 * n_species + 1] = std::max( + relaxed_bounds[2 * n_species + 1], s_min - entropy_relaxation); return relaxed_bounds; } diff --git a/source/multi_euler_ideal/limiter.template.h b/source/multi_euler_ideal/limiter.template.h index f502faad..3f5224f7 100644 --- a/source/multi_euler_ideal/limiter.template.h +++ b/source/multi_euler_ideal/limiter.template.h @@ -34,20 +34,17 @@ namespace ryujin const ScalarNumber relax = ScalarNumber(1. + large * eps); /* - * First limit the fractional density alpha_rho_0. + * Limit the partial densities for each species. * * See [Guermond, Nazarov, Popov, Thomas] (4.8): */ - { - /* Be careful with the variable names. Not really rho here but - * alpha_0 * rho_0. - */ - const auto &rho_U = view.partial_density_0(U); - const auto &rho_P = view.partial_density_0(P); + for (unsigned int k = 0; k < n_species; ++k) { + const auto rho_U = view.partial_density(U, k); + const auto rho_P = view.partial_density(P, k); - const auto &rho_min = std::get<0>(bounds); - const auto &rho_max = std::get<1>(bounds); + const auto &rho_min = bounds[2 * k]; + const auto &rho_max = bounds[2 * k + 1]; /* * Verify that rho_U is within bounds. This property might be @@ -60,14 +57,14 @@ namespace ryujin if (!(test_min == Number(0.) && test_max == Number(0.))) { #ifdef DEBUG_OUTPUT std::cout << std::fixed << std::setprecision(16); - std::cout - << "Bounds violation: low-order [species 0] density (critical)!" - << "\n\t\trho min: " << rho_min - << "\n\t\trho min (delta): " << negative_part(rho_U - rho_min) - << "\n\t\trho: " << rho_U - << "\n\t\trho max (delta): " << positive_part(rho_U - rho_max) - << "\n\t\trho max: " << rho_max << "\n" - << std::endl; + std::cout << "Bounds violation: low-order [species " << k + << "] density (critical)!" + << "\n\t\trho min: " << rho_min + << "\n\t\trho min (delta): " << negative_part(rho_U - rho_min) + << "\n\t\trho: " << rho_U + << "\n\t\trho max (delta): " << positive_part(rho_U - rho_max) + << "\n\t\trho max: " << rho_max << "\n" + << std::endl; #endif success = false; } @@ -113,7 +110,7 @@ namespace ryujin /* * Verify that the new state is within bounds: */ - const auto rho_new = view.partial_density_0(U + t_r * P); + const auto rho_new = view.partial_density(U + t_r * P, k); const auto test_new_min = view.filter_vacuum_density( std::max(Number(0.), rho_new - relax * rho_max)); const auto test_new_max = view.filter_vacuum_density( @@ -121,7 +118,8 @@ namespace ryujin if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) { #ifdef DEBUG_OUTPUT std::cout << std::fixed << std::setprecision(16); - std::cout << "Bounds violation: high-order [species 0] density!" + std::cout << "Bounds violation: high-order [species " << k + << "] density!" << "\n\t\trho min: " << rho_min << "\n\t\trho min (delta): " << negative_part(rho_new - rho_min) @@ -134,128 +132,19 @@ namespace ryujin success = false; } #endif - } - - /* - * Then limit the fractional density alpha_rho_1. - * - * See [Guermond, Nazarov, Popov, Thomas] (4.8): - */ - { - /* Be careful with the variable names. Not really rho here but - * alpha_1 * rho_1. - */ - auto t_q = t_r; - - const auto &rho_U = view.partial_density_1(U); - const auto &rho_P = view.partial_density_1(P); - - const auto &rho_min = std::get<2>(bounds); - const auto &rho_max = std::get<3>(bounds); - - /* - * Verify that rho_U is within bounds. This property might be - * violated for relative CFL numbers larger than 1. - */ - const auto test_min = view.filter_vacuum_density( - std::max(Number(0.), rho_U - relax * rho_max)); - const auto test_max = view.filter_vacuum_density( - std::max(Number(0.), rho_min - relax * rho_U)); - if (!(test_min == Number(0.) && test_max == Number(0.))) { -#ifdef DEBUG_OUTPUT - std::cout << std::fixed << std::setprecision(16); - std::cout - << "Bounds violation: low-order [species 1] density (critical)!" - << "\n\t\trho min: " << rho_min - << "\n\t\trho min (delta): " << negative_part(rho_U - rho_min) - << "\n\t\trho: " << rho_U - << "\n\t\trho max (delta): " << positive_part(rho_U - rho_max) - << "\n\t\trho max: " << rho_max << "\n" - << std::endl; -#endif - success = false; - } - - const Number denominator = - ScalarNumber(1.) / (std::abs(rho_P) + eps * rho_max + min); - - t_q = dealii::compare_and_apply_mask( - rho_max, - rho_U + t_q * rho_P, - /* - * rho_P is positive. - * - * Note: Do not take an absolute value here. If we are out of - * bounds we have to ensure that t_r is set to t_min. - */ - (rho_max - rho_U) * denominator, - t_q); - - t_q = dealii::compare_and_apply_mask( - rho_U + t_q * rho_P, - rho_min, - /* - * rho_P is negative. - * - * Note: Do not take an absolute value here. If we are out of - * bounds we have to ensure that t_r is set to t_min. - */ - (rho_U - rho_min) * denominator, - t_q); - - /* - * Ensure that t_min <= t <= t_max. This might not be the case if - * rho_U is outside the interval [rho_min, rho_max]. Furthermore, - * the quotient we take above is prone to numerical cancellation in - * particular in the second pass of the limiter when rho_P might be - * small. - */ - t_r = std::min(t_q, t_r); - t_r = std::min(t_r, t_max); - t_r = std::max(t_r, t_min); - -#ifdef DEBUG_EXPENSIVE_BOUNDS_CHECK - /* - * Verify that the new state is within bounds: - */ - const auto rho_new = view.partial_density_1(U + t_r * P); - const auto test_new_min = view.filter_vacuum_density( - std::max(Number(0.), rho_new - relax * rho_max)); - const auto test_new_max = view.filter_vacuum_density( - std::max(Number(0.), rho_min - relax * rho_new)); - if (!(test_new_min == Number(0.) && test_new_max == Number(0.))) { -#ifdef DEBUG_OUTPUT - std::cout << std::fixed << std::setprecision(16); - std::cout << "Bounds violation: high-order [species 1] density!" - << "\n\t\trho min: " << rho_min - << "\n\t\trho min (delta): " - << negative_part(rho_new - rho_min) - << "\n\t\trho: " << rho_new - << "\n\t\trho max (delta): " - << positive_part(rho_new - rho_max) - << "\n\t\trho max: " << rho_max << "\n" - << std::endl; -#endif - success = false; - } -#endif - } + } /* end loop over species */ /* * Then limit the internal energy. Unfortunately, limiting on the * "mixture" specific entropy does not guarantee positivity of the mixture * internal energy, so have to do this. - * - * Note that t_r is given by limiting on partial densities. This - * implementation might not be the most efficient, but we can optimize - * later. */ { Number t_l = t_min; // good state Number t_q = t_r; // potentially bad state - const auto &rho_e_min = std::get<4>(bounds); + const auto &rho_e_min = bounds[2 * n_species]; const auto U_l = U + t_l * P; const auto rho_e_l = view.internal_energy(U_l); @@ -345,7 +234,7 @@ namespace ryujin Number t_l = t_min; // good state { - const auto &s_min = std::get<5>(bounds); + const auto &s_min = bounds[2 * n_species + 1]; #ifdef DEBUG_OUTPUT_LIMITER std::cout << std::endl; diff --git a/source/multi_euler_ideal/riemann_solver.h b/source/multi_euler_ideal/riemann_solver.h index 6cc87b86..3bfdc3ca 100644 --- a/source/multi_euler_ideal/riemann_solver.h +++ b/source/multi_euler_ideal/riemann_solver.h @@ -69,14 +69,14 @@ namespace ryujin using state_type = typename View::state_type; /** - * Number of components in a "riemann data" state, we store \f$[\rho, v, - * p, a, gamma]\f$, thus, 5. + * Number of components in a "riemann data" state, we store + * \f$[\rho, u, p, \gamma, a]\f$, thus, 5. */ static constexpr unsigned int riemann_data_size = 5; /** * The array type to store the expanded "riemann data" state for the - * Riemann solver \f$[\rho, v, p, a]\f$ + * Riemann solver \f$[\rho, u, p, \gamma, a]\f$. */ using primitive_type = typename std::array; diff --git a/source/multi_euler_ideal/species_parameters.h b/source/multi_euler_ideal/species_parameters.h new file mode 100644 index 00000000..aaea573c --- /dev/null +++ b/source/multi_euler_ideal/species_parameters.h @@ -0,0 +1,29 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2025 by the ryujin authors +// Copyright (C) 2025 by Triad National Security, LLC +// + +#pragma once + +namespace ryujin +{ + namespace MultiSpeciesEuler + { + /** + * Compile-time parameter specifying the number of species in the + * multi-species Euler equations. Change this value and recompile to + * use a different number of species. + * + * The state vector dimension will be: n_species + dim + 1 + * (n_species partial densities + dim momentum components + 1 total energy) + * + * @note n_species must be at most 3. For n_species >= 4 the limiter + * bounds array becomes too large. FIXME. + */ + constexpr unsigned int n_species = 2; + static_assert(n_species >= 1 && n_species <= 3, + "n_species must be between 1 and 3"); + + } // namespace MultiSpeciesEuler +} // namespace ryujin From 018a749ed04abd3e38e5c2c7268c1c17db480b81 Mon Sep 17 00:00:00 2001 From: Eric Tovar Date: Sun, 1 Mar 2026 16:52:54 -0700 Subject: [PATCH 3/5] MultiSpeciesEuler: update files with recent development changes --- source/multi_euler_ideal/hyperbolic_system.h | 28 ++++++++++++------ source/multi_euler_ideal/indicator.h | 4 +-- source/multi_euler_ideal/limiter.h | 4 +-- .../riemann_solver.template.h | 4 +-- source/multi_euler_ideal/species_parameters.h | 29 ------------------- 5 files changed, 25 insertions(+), 44 deletions(-) delete mode 100644 source/multi_euler_ideal/species_parameters.h diff --git a/source/multi_euler_ideal/hyperbolic_system.h b/source/multi_euler_ideal/hyperbolic_system.h index ef20069f..b66b1de9 100644 --- a/source/multi_euler_ideal/hyperbolic_system.h +++ b/source/multi_euler_ideal/hyperbolic_system.h @@ -6,12 +6,11 @@ #pragma once -#include "species_parameters.h" - #include #include #include +#include #include #include #include @@ -26,6 +25,17 @@ namespace ryujin { namespace MultiSpeciesEuler { + /** + * Compile-time parameter specifying the number of species in the + * multi-species Euler equations. Change this value and recompile to + * use a different number of species. + * + * The state vector dimension will be: n_species + dim + 1 + * (n_species partial densities + dim momentum components + 1 total energy) + */ + constexpr unsigned int n_species = 2; + static_assert(n_species >= 1 && n_species <= 3, + "n_species must be between 1 and 3"); /* * For various divisions in the multi-species module we have a * mathematical guarantee that the numerator and denominator are @@ -931,7 +941,7 @@ namespace ryujin if (skip_constrained_dofs && row_length == 1) return; - const auto U_i = U.template get_tensor(i); + const auto U_i = U.template read_tensor(i); const auto view = this->view(); const auto rho_i = view.density(U_i); const auto p_i = view.pressure(U_i); @@ -957,8 +967,8 @@ namespace ryujin if (skip_constrained_dofs && row_length == 1) return; - const auto U_i = U.template get_tensor(i); - auto prec_i = precomputed.template get_tensor(i); + const auto U_i = U.template read_tensor(i); + auto prec_i = precomputed.template read_tensor(i); auto &[rho_i, p_i, gamma_min_i, s_i, harten_i] = prec_i; const auto view = this->view(); @@ -968,8 +978,8 @@ namespace ryujin for (unsigned int col_idx = 1; col_idx < row_length; ++col_idx, js += stride_size) { - const auto U_j = U.template get_tensor(js); - const auto prec_j = precomputed.template get_tensor(js); + const auto U_j = U.template read_tensor(js); + const auto prec_j = precomputed.template read_tensor(js); const auto p_j = std::get<1>(prec_j); const auto gamma_j = view.surrogate_gamma(U_j, p_j); gamma_min_i = std::min(gamma_min_i, gamma_j); @@ -1547,7 +1557,7 @@ namespace ryujin const state_type &U_i) const -> flux_contribution_type { const auto &[rho_i, p_i, surrogate_gamma_i, s_i, surrogate_harten_i] = - pv.template get_tensor(i); + pv.template read_tensor(i); return f(U_i, p_i); } @@ -1561,7 +1571,7 @@ namespace ryujin const state_type &U_j) const -> flux_contribution_type { const auto &[rho_j, p_j, surrogate_gamma_j, s_j, surrogate_harten_j] = - pv.template get_tensor(js); + pv.template read_tensor(js); return f(U_j, p_j); } diff --git a/source/multi_euler_ideal/indicator.h b/source/multi_euler_ideal/indicator.h index a561bbb8..812486d7 100644 --- a/source/multi_euler_ideal/indicator.h +++ b/source/multi_euler_ideal/indicator.h @@ -205,7 +205,7 @@ namespace ryujin const auto view = hyperbolic_system.view(); const auto &[rho_i, p_i, gamma_min_i, s_i, new_eta_i] = - precomputed_values.template get_tensor(i); + precomputed_values.template read_tensor(i); gamma_min = gamma_min_i; @@ -231,7 +231,7 @@ namespace ryujin { /* Entropy viscosity commutator: */ const auto &[rho_j, p_j, gamma_min_j, s_j, dont_use_j] = - precomputed_values.template get_tensor(js); + precomputed_values.template read_tensor(js); const auto view = hyperbolic_system.view(); diff --git a/source/multi_euler_ideal/limiter.h b/source/multi_euler_ideal/limiter.h index a6899f86..480c565a 100644 --- a/source/multi_euler_ideal/limiter.h +++ b/source/multi_euler_ideal/limiter.h @@ -279,7 +279,7 @@ namespace ryujin const auto view = hyperbolic_system.view(); const auto &[rho_i, p_i, gamma_min_i, s_i, harten_i] = - precomputed_values.template get_tensor(i); + precomputed_values.template read_tensor(i); Bounds result; /* Partial density bounds: [rho_k_min, rho_k_max] for each species */ @@ -405,7 +405,7 @@ namespace ryujin const auto rho_e_j = view.internal_energy(U_j); const auto [rho_j, p_j, gamma_min_j, s_j, harten_j] = - precomputed_values.template get_tensor(js); + precomputed_values.template read_tensor(js); /* bar state shifted by an affine shift: */ const auto U_ij_bar = diff --git a/source/multi_euler_ideal/riemann_solver.template.h b/source/multi_euler_ideal/riemann_solver.template.h index cc1a38a6..98ed4fb0 100644 --- a/source/multi_euler_ideal/riemann_solver.template.h +++ b/source/multi_euler_ideal/riemann_solver.template.h @@ -643,10 +643,10 @@ namespace ryujin const dealii::Tensor<1, dim, Number> &n_ij) const { const auto &[rho_i, p_i, dont_use_i, s_i, harten_i] = - precomputed_values.template get_tensor(i); + precomputed_values.template read_tensor(i); const auto &[rho_j, p_j, dont_use_j, s_j, harten_j] = - precomputed_values.template get_tensor(js); + precomputed_values.template read_tensor(js); const auto riemann_data_i = riemann_data_from_state(U_i, p_i, rho_i, n_ij); diff --git a/source/multi_euler_ideal/species_parameters.h b/source/multi_euler_ideal/species_parameters.h deleted file mode 100644 index aaea573c..00000000 --- a/source/multi_euler_ideal/species_parameters.h +++ /dev/null @@ -1,29 +0,0 @@ -// -// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception -// Copyright (C) 2023 - 2025 by the ryujin authors -// Copyright (C) 2025 by Triad National Security, LLC -// - -#pragma once - -namespace ryujin -{ - namespace MultiSpeciesEuler - { - /** - * Compile-time parameter specifying the number of species in the - * multi-species Euler equations. Change this value and recompile to - * use a different number of species. - * - * The state vector dimension will be: n_species + dim + 1 - * (n_species partial densities + dim momentum components + 1 total energy) - * - * @note n_species must be at most 3. For n_species >= 4 the limiter - * bounds array becomes too large. FIXME. - */ - constexpr unsigned int n_species = 2; - static_assert(n_species >= 1 && n_species <= 3, - "n_species must be between 1 and 3"); - - } // namespace MultiSpeciesEuler -} // namespace ryujin From 40d5468b0984153d00c0055f402a92a9f8402095 Mon Sep 17 00:00:00 2001 From: Eric Tovar Date: Sun, 1 Mar 2026 16:56:08 -0700 Subject: [PATCH 4/5] MultiSpeciesEuler: indent --- doc/headers/modules.h | 8 ++ doc/references.bib | 7 ++ source/multi_euler_ideal/hyperbolic_system.h | 83 +++++++++---------- .../initial_state_contrast.h | 13 ++- .../initial_state_exact_riemann_solution.h | 3 +- .../initial_state_function.h | 13 +-- .../initial_state_icf_like.h | 1 - .../initial_state_library_multi_euler.h | 4 +- .../initial_state_radial_contrast.h | 16 ++-- .../initial_state_shock_bubble.h | 1 - .../initial_state_three_state_contrast.h | 13 ++- source/multi_euler_ideal/limiter.template.h | 6 +- 12 files changed, 86 insertions(+), 82 deletions(-) diff --git a/doc/headers/modules.h b/doc/headers/modules.h index 4272cb5e..6366f3d4 100644 --- a/doc/headers/modules.h +++ b/doc/headers/modules.h @@ -85,6 +85,14 @@ */ +/** + * @defgroup MultiSpeciesEulerEquations The Multi-Species Euler Equations + * + * This module contains classes and functions related to solving the + * multi-species compressible Euler equations with ideal gas mixtures. + */ + + /** * @defgroup SkeletonEquations Minimal equation interface * diff --git a/doc/references.bib b/doc/references.bib index d62d3a03..a069b648 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -24,6 +24,13 @@ @article{Brooks1982 doi = {10.1016/0045-7825(82)90071-8} } +@article{ClaytonDzanicTovar-2025, + author = {Bennett Clayton and Tarik Dzanic and Eric Tovar}, + title = {Invariant-domain preserving high-order time stepping for multi-species ideal gas mixtures}, + year = {2025}, + journal = {submitted} +} + @article{Chertock2015, author = {Chertock, A. and Cui, S. and Kurganov, A. and Wu, T.}, title = {Well-balanced positivity preserving central-upwind scheme for diff --git a/source/multi_euler_ideal/hyperbolic_system.h b/source/multi_euler_ideal/hyperbolic_system.h index b66b1de9..949c2039 100644 --- a/source/multi_euler_ideal/hyperbolic_system.h +++ b/source/multi_euler_ideal/hyperbolic_system.h @@ -122,8 +122,8 @@ namespace ryujin * \f$\bar{c}_p = \sum_k Y_k c_{p,k}\f$ and * \f$\bar{c}_v = \sum_k Y_k c_{v,k}\f$. * - * The mixture specific entropy is (see Eq. (2.8) in @cite - * ClaytonDzanicTovar-2025): + * The mixture specific entropy is (see Eq. (2.8) in + * @cite ClaytonDzanicTovar-2025): * \f[ * s(\mathbf{u}) = \bar{c}_v \log\left(\frac{\rho e}{\rho^{\bar{\gamma}}} * \right) + K(\mathbf{Y}), @@ -499,21 +499,20 @@ namespace ryujin static dealii::Tensor<1, dim, Number> momentum(const state_type &U); /** - * For a given (n_species+1+dim dimensional) state vector U, return - * the total energy U[n_species+dim] + * For a given (n_species+1+dim dimensional) state vector U, + * return the total energy U[n_species+dim] */ static Number total_energy(const state_type &U); /** - * For a given (n_species+1+dim dimensional) state vector U, compute - * and return the internal energy \f$\varepsilon = (\rho e)\f$. + * For a given (n_species+1+dim dimensional) state vector U, + * compute and return the internal energy \f$\varepsilon = (\rho e)\f$. */ static Number internal_energy(const state_type &U); /** - * For a given (n_species+1+dim dimensional) state vector U, compute - * and return the mixture gamma for the mixture EOS: - * \f[ + * For a given (n_species+1+dim dimensional) state vector U, + * compute and return the mixture gamma for the mixture EOS: \f[ * \overline{\gamma} = (sum of (alpha_k rho_k) * c_{p, k}) / (sum of * (alpha_k rho_k) * c_{v, k}) * \f] @@ -521,26 +520,24 @@ namespace ryujin Number gamma_mixture(const state_type &U) const; /** - * For a given (n_species+1+dim dimensional) state vector U, compute - * and return the mixture specific heat capacity at constant pressure: - * \f[ - * \overline{c}_p = (sum of (alpha_k rho_k) / rho * c_{p, k}) - * \f] + * For a given (n_species+1+dim dimensional) state vector U, + * compute and return the mixture specific heat capacity at constant + * pressure: \f[ \overline{c}_p = (sum of (alpha_k rho_k) / rho * c_{p, + * k}) \f] */ Number cp_mixture(const state_type &U) const; /** - * For a given (n_species+1+dim dimensional) state vector U, compute - * and return the mixture specific heat capacity at constant volume: - * \f[ - * \overline{c}_v = (sum of (alpha_k rho_k) / rho * c_{v, k}) + * For a given (n_species+1+dim dimensional) state vector U, + * compute and return the mixture specific heat capacity at constant + * volume: \f[ \overline{c}_v = (sum of (alpha_k rho_k) / rho * c_{v, k}) * \f] */ Number cv_mixture(const state_type &U) const; /** - * For a given (n_species+1+dim dimensional) state vector U, compute - * and return the pressure \f$p\f$. + * For a given (n_species+1+dim dimensional) state vector U, + * compute and return the pressure \f$p\f$. * * We assume that the pressure is given by a mixture ideal EOS * \f[ @@ -569,11 +566,10 @@ namespace ryujin //@{ /** - * For a given (n_species+1+dim dimensional) state vector U, compute - * and return the physical specific entropy. Following Eq. (2.8) in - * [ClaytonDzanicTovar-2025]: - * \f[ - * s(\mathbf{u}) = \bar{r} \log(\rho^{-1}) + \bar{c}_v \log(e) + * For a given (n_species+1+dim dimensional) state vector U, + * compute and return the physical specific entropy. Following Eq. (2.8) + * in [ClaytonDzanicTovar-2025]: \f[ s(\mathbf{u}) = \bar{r} + * \log(\rho^{-1}) + \bar{c}_v \log(e) * + K(\mathbf{Y}), * \f] * where \f$\bar{r} = \sum_k Y_k r_k\f$, \f$\bar{c}_v = \sum_k Y_k @@ -583,22 +579,21 @@ namespace ryujin Number specific_entropy(const state_type &U) const; /** - * For a given (n_species+1+dim dimensional) state vector U, compute - * and return a surrogate Harten-type entropy. Following Section 4 of - * [ClaytonDzanicTovar-2025], we use: - * \f[ - * \eta(\mathbf{u}; \gamma_{\min}) = (\rho - * \varepsilon)^{1/(\gamma_{\min}+1)}, \f] where \f$\gamma_{\min}\f$ is - * the minimum surrogate gamma over the stencil. This entropy is chosen to - * ensure convexity properties required by the entropy-viscosity - * indicator. + * For a given (n_species+1+dim dimensional) state vector U, + * compute and return a surrogate Harten-type entropy. Following Section 4 + * of [ClaytonDzanicTovar-2025], we use: \f[ \eta(\mathbf{u}; + * \gamma_{\min}) = (\rho \varepsilon)^{1/(\gamma_{\min}+1)}, \f] where + * \f$\gamma_{\min}\f$ is the minimum surrogate gamma over the stencil. + * This entropy is chosen to ensure convexity properties required by the + * entropy-viscosity indicator. */ Number surrogate_harten_entropy(const state_type &U, const Number &gamma_min) const; /** - * For a given (n_species+1+dim dimensional) state vector U, compute - * and return the derivative \f$\eta'\f$ of the Harten-type entropy. + * For a given (n_species+1+dim dimensional) state vector U, + * compute and return the derivative \f$\eta'\f$ of the Harten-type + * entropy. */ mixture_state_type surrogate_harten_entropy_derivative(const state_type &U, @@ -606,8 +601,8 @@ namespace ryujin const Number &gamma_min) const; /** - * For a given (n_species+1+dim dimensional) state vector U and - * pressure p, compute a surrogate gamma. Following + * For a given (n_species+1+dim dimensional) state vector U + * and pressure p, compute a surrogate gamma. Following * Section 3 of [ClaytonDzanicTovar-2025]: * \f[ * \gamma(\mathbf{u}, p) = 1 + \frac{p}{\varepsilon(\mathbf{u})}. @@ -620,19 +615,17 @@ namespace ryujin Number surrogate_gamma(const state_type &U, const Number &p) const; /** - * For a given (n_species+1+dim dimensional) state vector U and - * gamma gamma, compute a surrogate pressure: - * \f[ - * p(\rho, e, \gamma) = (\gamma - 1) (\rho e) - * \f] + * For a given (n_species+1+dim dimensional) state vector U + * and gamma gamma, compute a surrogate pressure: \f[ p(\rho, + * e, \gamma) = (\gamma - 1) (\rho e) \f] * * This function is the complementary function to surrogate_gamma(). */ Number surrogate_pressure(const state_type &U, const Number &gamma) const; /** - * For a given (n_species+1+dim dimensional) state vector U and - * gamma gamma, compute a surrogate speed of sound. + * For a given (n_species+1+dim dimensional) state vector U + * and gamma gamma, compute a surrogate speed of sound. */ Number surrogate_speed_of_sound(const state_type &U, const Number &gamma) const; diff --git a/source/multi_euler_ideal/initial_state_contrast.h b/source/multi_euler_ideal/initial_state_contrast.h index 9210e751..c5fdfeb8 100644 --- a/source/multi_euler_ideal/initial_state_contrast.h +++ b/source/multi_euler_ideal/initial_state_contrast.h @@ -50,9 +50,9 @@ namespace ryujin /* Default: equal mass fractions */ for (unsigned int k = 0; k < n_species - 1; ++k) temp_left_[k] = Number(1.) / Number(n_species); - temp_left_[n_species - 1] = 1.4; /* rho */ - temp_left_[n_species] = 0.0; /* u */ - temp_left_[n_species + 1] = 1.0; /* p */ + temp_left_[n_species - 1] = 1.4; /* rho */ + temp_left_[n_species] = 0.0; /* u */ + temp_left_[n_species + 1] = 1.0; /* p */ this->add_parameter( "primitive state left", temp_left_, @@ -61,9 +61,9 @@ namespace ryujin for (unsigned int k = 0; k < n_species - 1; ++k) temp_right_[k] = Number(1.) / Number(n_species); - temp_right_[n_species - 1] = 1.4; /* rho */ - temp_right_[n_species] = 0.0; /* u */ - temp_right_[n_species + 1] = 1.0; /* p */ + temp_right_[n_species - 1] = 1.4; /* rho */ + temp_right_[n_species] = 0.0; /* u */ + temp_right_[n_species + 1] = 1.0; /* p */ this->add_parameter( "primitive state right", temp_right_, @@ -96,7 +96,6 @@ namespace ryujin state_type state_left_; state_type state_right_; - }; } // namespace MultiSpeciesEulerInitialStates } // namespace ryujin diff --git a/source/multi_euler_ideal/initial_state_exact_riemann_solution.h b/source/multi_euler_ideal/initial_state_exact_riemann_solution.h index 449ee48b..8a50773e 100644 --- a/source/multi_euler_ideal/initial_state_exact_riemann_solution.h +++ b/source/multi_euler_ideal/initial_state_exact_riemann_solution.h @@ -240,8 +240,7 @@ namespace ryujin //@{ // Used to transform data to (alpha_rho_0, alpha_rho_1, u, p) - primitive_state_type - transform(primitive_state_type &temp_in) const + primitive_state_type transform(primitive_state_type &temp_in) const { primitive_state_type result; result[0] = temp_in[0] * temp_in[1]; // = alpha_0 rho_0; diff --git a/source/multi_euler_ideal/initial_state_function.h b/source/multi_euler_ideal/initial_state_function.h index ed296ce1..30e8a8dc 100644 --- a/source/multi_euler_ideal/initial_state_function.h +++ b/source/multi_euler_ideal/initial_state_function.h @@ -45,11 +45,11 @@ namespace ryujin for (unsigned int k = 0; k < n_species - 1; ++k) { Y_expressions_[k] = std::to_string(Number(1.) / Number(n_species)); /* default */ - this->add_parameter( - "Y_" + std::to_string(k) + " expression", - Y_expressions_[k], - "A function expression describing the mass fraction for species " + - std::to_string(k)); + this->add_parameter("Y_" + std::to_string(k) + " expression", + Y_expressions_[k], + "A function expression describing the mass " + "fraction for species " + + std::to_string(k)); } density_expression_ = "1.4"; @@ -146,7 +146,8 @@ namespace ryujin } pressure_function_->set_time(t); - full_primitive_state[n_species + dim] = pressure_function_->value(point); + full_primitive_state[n_species + dim] = + pressure_function_->value(point); return view.from_primitive_state(full_primitive_state); } diff --git a/source/multi_euler_ideal/initial_state_icf_like.h b/source/multi_euler_ideal/initial_state_icf_like.h index 6c85ffb8..d1b53197 100644 --- a/source/multi_euler_ideal/initial_state_icf_like.h +++ b/source/multi_euler_ideal/initial_state_icf_like.h @@ -198,7 +198,6 @@ namespace ryujin double amplitude_; double shock_radius_; double mach_number_; - }; diff --git a/source/multi_euler_ideal/initial_state_library_multi_euler.h b/source/multi_euler_ideal/initial_state_library_multi_euler.h index 39b9b14f..db6cb428 100644 --- a/source/multi_euler_ideal/initial_state_library_multi_euler.h +++ b/source/multi_euler_ideal/initial_state_library_multi_euler.h @@ -41,8 +41,8 @@ namespace ryujin /* ExactRiemannSolution only supports 2 species */ if constexpr (n_species == 2) - add(std::make_unique>(h, - s)); + add(std::make_unique>( + h, s)); add(std::make_unique>(h, s)); add(std::make_unique>(h, s)); diff --git a/source/multi_euler_ideal/initial_state_radial_contrast.h b/source/multi_euler_ideal/initial_state_radial_contrast.h index bd22b736..34dcd686 100644 --- a/source/multi_euler_ideal/initial_state_radial_contrast.h +++ b/source/multi_euler_ideal/initial_state_radial_contrast.h @@ -55,9 +55,9 @@ namespace ryujin /* Default: equal mass fractions */ for (unsigned int k = 0; k < n_species - 1; ++k) temp_left_[k] = Number(1.) / Number(n_species); - temp_left_[n_species - 1] = 1.4; /* rho */ - temp_left_[n_species] = 0.0; /* u */ - temp_left_[n_species + 1] = 1.0; /* p */ + temp_left_[n_species - 1] = 1.4; /* rho */ + temp_left_[n_species] = 0.0; /* u */ + temp_left_[n_species + 1] = 1.0; /* p */ this->add_parameter( "primitive state left", temp_left_, @@ -66,9 +66,9 @@ namespace ryujin for (unsigned int k = 0; k < n_species - 1; ++k) temp_right_[k] = Number(1.) / Number(n_species); - temp_right_[n_species - 1] = 1.4; /* rho */ - temp_right_[n_species] = 0.0; /* u */ - temp_right_[n_species + 1] = 1.0; /* p */ + temp_right_[n_species - 1] = 1.4; /* rho */ + temp_right_[n_species] = 0.0; /* u */ + temp_right_[n_species + 1] = 1.0; /* p */ this->add_parameter( "primitive state right", temp_right_, @@ -94,8 +94,7 @@ namespace ryujin state_type compute(const dealii::Point &point, Number /*t*/) final { - auto result = - (point.norm() > radius_ ? state_right_ : state_left_); + auto result = (point.norm() > radius_ ? state_right_ : state_left_); if (point.norm() > 0. && use_radial_velocity_) { const auto view = hyperbolic_system_.template view(); @@ -121,7 +120,6 @@ namespace ryujin double radius_; bool use_radial_velocity_; - }; diff --git a/source/multi_euler_ideal/initial_state_shock_bubble.h b/source/multi_euler_ideal/initial_state_shock_bubble.h index db4f95f0..148be0bd 100644 --- a/source/multi_euler_ideal/initial_state_shock_bubble.h +++ b/source/multi_euler_ideal/initial_state_shock_bubble.h @@ -164,7 +164,6 @@ namespace ryujin double radius_; double mach_number_; double shock_location_; - }; diff --git a/source/multi_euler_ideal/initial_state_three_state_contrast.h b/source/multi_euler_ideal/initial_state_three_state_contrast.h index 0739e9a2..8eeb70ad 100644 --- a/source/multi_euler_ideal/initial_state_three_state_contrast.h +++ b/source/multi_euler_ideal/initial_state_three_state_contrast.h @@ -51,9 +51,9 @@ namespace ryujin /* Default: equal mass fractions */ for (unsigned int k = 0; k < n_species - 1; ++k) temp_left_[k] = Number(1.) / Number(n_species); - temp_left_[n_species - 1] = 1.; /* rho */ - temp_left_[n_species] = 0.; /* u */ - temp_left_[n_species + 1] = 1.e3; /* p */ + temp_left_[n_species - 1] = 1.; /* rho */ + temp_left_[n_species] = 0.; /* u */ + temp_left_[n_species + 1] = 1.e3; /* p */ this->add_parameter( "primitive state left", temp_left_, @@ -83,9 +83,9 @@ namespace ryujin for (unsigned int k = 0; k < n_species - 1; ++k) temp_right_[k] = Number(1.) / Number(n_species); - temp_right_[n_species - 1] = 1.; /* rho */ - temp_right_[n_species] = 0.; /* u */ - temp_right_[n_species + 1] = 1.e2; /* p */ + temp_right_[n_species - 1] = 1.; /* rho */ + temp_right_[n_species] = 0.; /* u */ + temp_right_[n_species + 1] = 1.e2; /* p */ this->add_parameter( "primitive state right", temp_right_, @@ -127,7 +127,6 @@ namespace ryujin state_type state_left_; state_type state_middle_; state_type state_right_; - }; } // namespace MultiSpeciesEulerInitialStates } // namespace ryujin diff --git a/source/multi_euler_ideal/limiter.template.h b/source/multi_euler_ideal/limiter.template.h index 3f5224f7..7f7bb8b5 100644 --- a/source/multi_euler_ideal/limiter.template.h +++ b/source/multi_euler_ideal/limiter.template.h @@ -60,9 +60,11 @@ namespace ryujin std::cout << "Bounds violation: low-order [species " << k << "] density (critical)!" << "\n\t\trho min: " << rho_min - << "\n\t\trho min (delta): " << negative_part(rho_U - rho_min) + << "\n\t\trho min (delta): " + << negative_part(rho_U - rho_min) << "\n\t\trho: " << rho_U - << "\n\t\trho max (delta): " << positive_part(rho_U - rho_max) + << "\n\t\trho max (delta): " + << positive_part(rho_U - rho_max) << "\n\t\trho max: " << rho_max << "\n" << std::endl; #endif From 0648bb3c16b00cfea463c0015abf130bb931a71a Mon Sep 17 00:00:00 2001 From: Eric Tovar Date: Sun, 1 Mar 2026 17:28:05 -0700 Subject: [PATCH 5/5] MultiSpeciesEuler: add tests --- tests/multi_euler_ideal/CMakeLists.txt | 20 ++++++++ ...th_wave-1d-erk33.mpirun=2.threads=1.output | 12 +++++ .../verification-smooth_wave-1d-erk33.prm | 50 +++++++++++++++++++ ...tion-smooth_wave-1d-erk33.threads=1.output | 12 +++++ 4 files changed, 94 insertions(+) create mode 100644 tests/multi_euler_ideal/CMakeLists.txt create mode 100644 tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.mpirun=2.threads=1.output create mode 100644 tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.prm create mode 100644 tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.threads=1.output diff --git a/tests/multi_euler_ideal/CMakeLists.txt b/tests/multi_euler_ideal/CMakeLists.txt new file mode 100644 index 00000000..0d84eddd --- /dev/null +++ b/tests/multi_euler_ideal/CMakeLists.txt @@ -0,0 +1,20 @@ +## +## SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +## Copyright (C) 2025 by the ryujin authors +## Copyright (C) 2025 by Triad National Security, LLC +## + +set(EQUATION multi_euler_ideal) + +include_directories( + ${CMAKE_BINARY_DIR}/source/ + ${CMAKE_SOURCE_DIR}/source/${EQUATION} + ${CMAKE_SOURCE_DIR}/source/ + ) + +set(TEST_LIBRARIES obj_common obj_${EQUATION} obj_${EQUATION}_dependent) +set(TEST_TARGET ryujin) + +if(TARGET obj_${EQUATION}) + deal_ii_pickup_tests() +endif() diff --git a/tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.mpirun=2.threads=1.output b/tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.mpirun=2.threads=1.output new file mode 100644 index 00000000..9df9e599 --- /dev/null +++ b/tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.mpirun=2.threads=1.output @@ -0,0 +1,12 @@ +[INFO] initiating flux capacitor +[INFO] dispatching to driver »multi species euler« with dim=1 +[INFO] initializing data structures +[INFO] creating mesh and interpolating initial values +[INFO] preparing compute kernels +[INFO] entering main loop +Normalized consolidated Linf, L1, and L2 errors at final time +#dofs = 1601 +t = 0.6 +Linf = 3.751197956580734e-05 +L1 = 1.877759814590165e-06 +L2 = 6.038551272500228e-06 diff --git a/tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.prm b/tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.prm new file mode 100644 index 00000000..53364875 --- /dev/null +++ b/tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.prm @@ -0,0 +1,50 @@ +subsection A - TimeLoop + set basename = test-smooth-wave + + set enable compute error = true + set error normalize = true + set error quantities = alpha_rho_0, alpha_rho_1, m, E + + set final time = 0.6 + set enforce final time = true + + set timer granularity = 0.6 + set terminal update interval = 0 +end + +subsection B - Equation + set dimension = 1 + set equation = multi species euler + set c_p for each species = 1005., 4041.4 + set c_v for each species = 718., 2420. +end + +subsection C - Discretization + set finite element ansatz = cG Q1 + set geometry = rectangular domain + set mesh refinement = 4 + + subsection rectangular domain + set position bottom left = 0 + set position top right = 1 + set subdivisions x = 100 + set boundary condition left = dirichlet + set boundary condition right = dirichlet + end +end + +subsection E - InitialValues + set configuration = smooth wave + set position = 0. + + subsection smooth wave + set mass fraction = 0.75 + end +end + +subsection H - TimeIntegrator + set cfl min = 0.5 + set cfl max = 0.5 + set cfl recovery strategy = none + set time stepping scheme = erk 33 +end diff --git a/tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.threads=1.output b/tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.threads=1.output new file mode 100644 index 00000000..59350460 --- /dev/null +++ b/tests/multi_euler_ideal/verification-smooth_wave-1d-erk33.threads=1.output @@ -0,0 +1,12 @@ +[INFO] initiating flux capacitor +[INFO] dispatching to driver »multi species euler« with dim=1 +[INFO] initializing data structures +[INFO] creating mesh and interpolating initial values +[INFO] preparing compute kernels +[INFO] entering main loop +Normalized consolidated Linf, L1, and L2 errors at final time +#dofs = 1601 +t = 0.6 +Linf = 3.751197982289476e-05 +L1 = 1.87775980951388e-06 +L2 = 6.038551270344887e-06