diff --git a/doc/headers/modules.h b/doc/headers/modules.h index 4272cb5e1..6366f3d46 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 d62d3a03c..a069b6482 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/CMakeLists.txt b/source/multi_euler_ideal/CMakeLists.txt new file mode 100644 index 000000000..7b169e0d6 --- /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 000000000..9ee0ef15e --- /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 000000000..a17fc589c --- /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 000000000..949c20393 --- /dev/null +++ b/source/multi_euler_ideal/hyperbolic_system.h @@ -0,0 +1,1678 @@ +// +// 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 + +#include + +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 + * 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)); + } + + + /** + * 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; + + /** + * 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] + * + * + * @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, 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_; + 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, n_species, double> + cp_for_each_species() const + { + return hyperbolic_system_.cp_for_each_species_; + } + + 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, n_species, double> + r_for_each_species() const + { + return hyperbolic_system_.r_for_each_species_; + } + + DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, n_species, 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: n_species partial densities + + * dim momentum components + 1 total energy. + */ + static constexpr unsigned int problem_dimension = n_species + 1 + 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 { + 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; + }(); + + /** + * An array holding all component names of the primitive state as a + * string. + */ + static inline const auto primitive_component_names = + []() -> std::array { + 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; + }(); + + /** + * 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 state vector U, return the partial density + * for species k, i.e., U[k]. + */ + static Number partial_density(const state_type &U, unsigned int k); + + /** + * For a given state vector U, return the total mixture + * density obtained by summing all 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 (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 (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$. + */ + 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[ + * \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 (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}) + * \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$. + * + * 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 (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 + * \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 (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 + * 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 (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. + */ + mixture_state_type + surrogate_harten_entropy_derivative(const state_type &U, + const Number &eta, + const Number &gamma_min) const; + + /** + * 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})}. + * \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 (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. + */ + 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_{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 + * 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"); + + /* 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", + 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]() { + 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); + 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 read_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 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(); + + 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 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); + } + + 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(const state_type &U, + unsigned int k) + { + AssertIndexRange(k, n_species); + return U[k]; + } + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::density(const state_type &U) + { + auto result = Number(0.); + for (unsigned int k = 0; k < n_species; ++k) + result += partial_density(U, k); + return result; + } + + + 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[n_species + i]; + return result; + } + + + template + DEAL_II_ALWAYS_INLINE inline Number + HyperbolicSystemView::total_energy(const state_type &U) + { + return U[n_species + 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 k = 0; k < n_species; ++k) { + result += U[k] / rho * ScalarNumber(cp_for_each_species()[k]); + } + 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 k = 0; k < n_species; ++k) { + result += U[k] / rho * ScalarNumber(cv_for_each_species()[k]); + } + 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 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 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 = + 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); + /* 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; + + 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; + + 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[n_species + d] = rho_new * (vn_new * normal + vperp)[d]; + } + + U_new[n_species + 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 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 d = 0; d < dim; ++d) + result[n_species + d] = m[d]; + + } else if (id == Boundary::no_slip) { + for (unsigned int d = 0; d < dim; ++d) + result[n_species + d] = 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) { + 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[n_species + dim] = 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 read_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 read_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 - 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; + /* 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; + } + + + 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 pressure: */ + const auto &p = primitive_state[n_species + dim]; + + auto state = primitive_state; + /* Fix up momentum: */ + 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[n_species + dim] = + 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 = n_species; i < n_species + dim; ++i) + primitive_state[i] *= rho_inverse; + /* Set pressure: */ + primitive_state[n_species + dim] = 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[n_species + 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 000000000..812486d7a --- /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 viscosity commutator"); + } + + 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 read_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 read_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 < 2 + dim; ++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 < 2 + dim; ++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 000000000..c5fdfeb87 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_contrast.h @@ -0,0 +1,101 @@ +// +// 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 "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 + */ + template + class Contrast : public InitialState + { + public: + using HyperbolicSystem = typename Description::HyperbolicSystem; + using View = + 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) + { + /* 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, ..., Y_{n-2}, rho, u, p) on the " + "left"); + + 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, ..., 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_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_); + }; + 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_; + + primitive_state_type temp_left_; + primitive_state_type temp_right_; + + 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 new file mode 100644 index 000000000..8a50773ed --- /dev/null +++ b/source/multi_euler_ideal/initial_state_exact_riemann_solution.h @@ -0,0 +1,549 @@ +// +// 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 + +// #define DEBUG_SOLUTION + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + using namespace MultiSpeciesEuler; + + /** + * 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 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. + * + * @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; + + /* 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) + : InitialState("exact riemann solution", + 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; + 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; + + primitive_state_type 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 + */ + //@{ + + primitive_state_type primitive_left_; + primitive_state_type 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) + 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; + 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 primitive_state_type &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 primitive_state_type &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 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 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]; + + return fZofP(p_in, data_right) + fZofP(p_in, data_left) + u_R - u_L; + } + + + Number lambda(const Number &p_in, + const primitive_state_type &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 primitive_state_type &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; + } + + + primitive_state_type + cstar_solution(const Number &p_star, + const Number &u_star, + const primitive_state_type &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); + } + + + primitive_state_type + expansion_solution(const Number & /*p_star*/, + const Number &xi, + const primitive_state_type &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, + primitive_state_type data_1, + primitive_state_type 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 000000000..30e8a8dc5 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_function.h @@ -0,0 +1,174 @@ +// +// 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 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 + */ + 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) + { + /* 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", + 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«: + */ + 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) + 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; + + density_function_->set_time(t); + 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[n_species] = velocity_x_function_->value(point); + + if constexpr (dim > 1) { + velocity_y_function_->set_time(t); + full_primitive_state[n_species + 1] = + velocity_y_function_->value(point); + } + if constexpr (dim > 2) { + velocity_z_function_->set_time(t); + full_primitive_state[n_species + 2] = + velocity_z_function_->value(point); + } + + pressure_function_->set_time(t); + full_primitive_state[n_species + dim] = + pressure_function_->value(point); + + return view.from_primitive_state(full_primitive_state); + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + 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::array>, n_species - 1> + Y_functions_; + 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 000000000..d1b531978 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_icf_like.h @@ -0,0 +1,205 @@ +// +// 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 "hyperbolic_system.h" + +#include +#include + +namespace ryujin +{ + namespace MultiSpeciesEulerInitialStates + { + using namespace MultiSpeciesEuler; + + /** + * 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 + * 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; + + 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) + { + 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, ..., 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, ..., Y_{n-2}, 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 perturbation of interface"); + + amplitude_ = 0.02; + this->add_parameter( + "amplitude", amplitude_, "Amplitude for interface perturbation"); + + 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_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); + 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_[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; + + 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; + + /* 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[n_species + i] = Number(0.); + } + + if (point.norm() > 0.) { + for (unsigned int i = 0; i < dim; ++i) { + primitive_shock_state[n_species + i] = -u_L * r_hat[i]; + } + } + primitive_shock_state[n_species + 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_; + + primitive_state_type temp_inside_; + primitive_state_type temp_outside_; + + state_type state_inside_; + state_type state_outside_; + + double interface_radius_; + double num_modes_; + double amplitude_; + double shock_radius_; + double mach_number_; + }; + + + } // 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 000000000..886516508 --- /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 000000000..3e301341e --- /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 000000000..db6cb4285 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_library_multi_euler.h @@ -0,0 +1,55 @@ +// +// 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" +#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 + { + using namespace MultiSpeciesEuler; + + 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)); + + /* 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)); + 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 000000000..34dcd6866 --- /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 "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" + * 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; + + /* 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", + 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)"); + + /* 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, ..., 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, ..., Y_{n-2}, 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_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_); + }; + + this->parse_parameters_call_back.connect(convert_states); + convert_states(); + }; + + state_type compute(const dealii::Point &point, Number /*t*/) final + { + auto result = (point.norm() > radius_ ? state_right_ : state_left_); + + 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(); + } + + return result; + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + primitive_state_type temp_left_; + primitive_state_type temp_right_; + + state_type state_left_; + state_type state_right_; + + double radius_; + bool use_radial_velocity_; + }; + + + } // 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 000000000..148be0bdb --- /dev/null +++ b/source/multi_euler_ideal/initial_state_shock_bubble.h @@ -0,0 +1,171 @@ +// +// 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 "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 + * 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; + + 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) + { + 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, ..., 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, ..., Y_{n-2}, 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(); + + 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); + + /* Precompute shocked state via Rankine-Hugoniot: */ + + const auto gamma_ = view.gamma_mixture(state_ambient_); + + 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; + + 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.)); + + auto sign = 1.; + if (bubble_center_[0] < shock_location_) + sign = -1.; + + 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; + + 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 (point[0] <= bubble_center_[0] + shock_location_) + final_state = state_shock_; + + return final_state; + } + + private: + const HyperbolicSystem &hyperbolic_system_; + + primitive_state_type temp_inside_; + primitive_state_type temp_ambient_; + + state_type state_bubble_; + state_type state_ambient_; + state_type state_shock_; + + dealii::Point bubble_center_; + + double radius_; + double mach_number_; + double shock_location_; + }; + + + } // 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 000000000..2eaef6cb9 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_smooth_wave.h @@ -0,0 +1,129 @@ +// +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// Copyright (C) 2023 - 2024 by the ryujin authors +// + +#pragma once + +#include "hyperbolic_system.h" + +#include +#include + +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. + * + * @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; + + /* 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) + , hyperbolic_system_(hyperbolic_system) + { + mass_fraction_ = 0.5; + this->add_parameter("mass fraction", + mass_fraction_, + "The mass fraction of first species (Y_0). " + "Remaining mass is split equally among other " + "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; + + /* 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); + } + + 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 000000000..8eeb70ad2 --- /dev/null +++ b/source/multi_euler_ideal/initial_state_three_state_contrast.h @@ -0,0 +1,132 @@ +// +// 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 + +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 + * 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; + + /* 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) + { + /* 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, ..., 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"); + + 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, ..., 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"); + + 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, ..., 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_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_); + 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_; + + 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_; + }; + } // 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 000000000..b74f65df0 --- /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 000000000..56968e89c --- /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 000000000..480c565aa --- /dev/null +++ b/source/multi_euler_ideal/limiter.h @@ -0,0 +1,494 @@ +// +// 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: + * 2 * n_species (min/max per species) + 1 (rho_e_min) + 1 (s_min) + */ + static constexpr unsigned int n_bounds = 2 * n_species + 2; + + /** + * 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_; + + dealii::Tensor<1, n_species, Number> rho_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 read_tensor(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 result; + } + + + template + DEAL_II_ALWAYS_INLINE inline auto Limiter::combine_bounds( + const Bounds &bounds_left, const Bounds &bounds_right) const -> Bounds + { + 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; + } + + + template + DEAL_II_ALWAYS_INLINE inline auto + Limiter::fully_relax_bounds(const Bounds &bounds, + const Number &hd) const -> Bounds + { + auto relaxed_bounds = 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(); + + /* 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); + } + + /* 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); + + 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); + + /* 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: */ + 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.); + 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(); + + 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 read_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 rho_e_ij_bar = view.internal_energy(U_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); + + 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); + + 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); + + 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: */ + auto &s_min = bounds_[2 * n_species + 1]; + 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 + { + auto relaxed_bounds = fully_relax_bounds(bounds_, hd_i); + + /* Apply a stricter window: */ + constexpr ScalarNumber eps = std::numeric_limits::epsilon(); + + /* 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); + } + + /* 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); + + 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; + } + } // 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 000000000..7f7bb8b5a --- /dev/null +++ b/source/multi_euler_ideal/limiter.template.h @@ -0,0 +1,393 @@ +// +// 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); + + /* + * Limit the partial densities for each species. + * + * See [Guermond, Nazarov, Popov, Thomas] (4.8): + */ + + 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 = bounds[2 * k]; + const auto &rho_max = bounds[2 * k + 1]; + + /* + * 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 " << 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; + } + + 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(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( + 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 " << k + << "] 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. + */ + + { + Number t_l = t_min; // good state + Number t_q = t_r; // potentially bad state + + 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); + 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 = bounds[2 * n_species + 1]; + +#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 000000000..a4a4ac332 --- /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 000000000..3bfdc3ca8 --- /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, 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, u, p, \gamma, 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 000000000..98ed4fb0a --- /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 read_tensor(i); + + const auto &[rho_j, p_j, dont_use_j, s_j, harten_j] = + precomputed_values.template read_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 diff --git a/tests/multi_euler_ideal/CMakeLists.txt b/tests/multi_euler_ideal/CMakeLists.txt new file mode 100644 index 000000000..0d84eddd9 --- /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 000000000..9df9e5998 --- /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 000000000..533648750 --- /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 000000000..59350460a --- /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