hyperbolic_system.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2020 - 2026 by the ryujin authors
//
 
#pragma once
 
#include <compile_time_options.h>
 
#include "barotropic_equation_of_state_library.h"
 
#include <convenience_macros.h>
#include <discretization.h>
#include <loop.h>
#include <multicomponent_vector.h>
#include <patterns_conversion.h>
#include <simd.h>
#include <state_vector.h>
 
#include <deal.II/base/parameter_acceptor.h>
#include <deal.II/base/tensor.h>
 
#include <array>
 
namespace ryujin
{
  namespace EulerBarotropic
  {
    /*
     * For various divisions in the barotropic equation of state 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 <typename Number>
    DEAL_II_ALWAYS_INLINE inline Number safe_division(const Number &numerator,
                                                      const Number &denominator)
    {
      using ScalarNumber = typename get_value_type<Number>::type;
      constexpr ScalarNumber min = std::numeric_limits<ScalarNumber>::min();
 
      return std::max(numerator, Number(0.)) /
             std::max(denominator, Number(min));
    }
 
 
    template <int dim, typename Number>
    class HyperbolicSystemView;
 
    class HyperbolicSystem final : public dealii::ParameterAcceptor
    {
    public:
      static inline std::string problem_name =
          "Compressible Euler equations (barotropic EOS, optimized barotropic)";
 
      HyperbolicSystem(const std::string &subsection = "/HyperbolicSystem");
 
      template <int dim, typename Number = double>
      using View = HyperbolicSystemView<dim, Number>;
 
      template <int dim, typename Number>
      auto view() const
      {
        return View<dim, Number>{*this};
      }
 
      template <int dim, typename ScalarNumber>
      void fill_precomputed_values(
          const OfflineData<dim, ScalarNumber> &offline_data,
          typename HyperbolicSystemView<dim, ScalarNumber>::StateVector
              &state_vector,
          const bool skip_constrained_dofs = true) const;
 
    private:
 
      std::string barotropic_equation_of_state_;
 
      double reference_density_;
      double vacuum_state_relaxation_small_;
      double vacuum_state_relaxation_large_;
 
 
 
      BarotropicEquationOfStateLibrary::equation_of_state_list_type
          barotropic_equation_of_state_list_;
 
      using BarotropicEquationOfState =
          BarotropicEquationOfStateLibrary::BarotropicEquationOfState;
      std::shared_ptr<BarotropicEquationOfState>
          selected_barotropic_equation_of_state_;
 
      template <int dim, typename Number>
      friend class HyperbolicSystemView;
 
    }; /* HyperbolicSystem */
 
 
    template <int dim, typename Number>
    class HyperbolicSystemView
    {
    public:
 
      using ScalarNumber = typename get_value_type<Number>::type;
 
      static constexpr unsigned int problem_dimension = 1 + dim;
 
      using state_type = dealii::Tensor<1, problem_dimension, Number>;
 
      using flux_type =
          dealii::Tensor<1, problem_dimension, dealii::Tensor<1, dim, Number>>;
 
      using flux_contribution_type = flux_type;
 
      static inline const auto component_names =
          []() -> std::array<std::string, problem_dimension> {
        if constexpr (dim == 1)
          return {"rho", "m"};
        else if constexpr (dim == 2)
          return {"rho", "m_1", "m_2"};
        else if constexpr (dim == 3)
          return {"rho", "m_1", "m_2", "m_3"};
        __builtin_trap();
      }();
 
      static inline const auto primitive_component_names =
          []() -> std::array<std::string, problem_dimension> {
        if constexpr (dim == 1)
          return {"rho", "v"};
        else if constexpr (dim == 2)
          return {"rho", "v_1", "v_2"};
        else if constexpr (dim == 3)
          return {"rho", "v_1", "v_2", "v_3"};
        __builtin_trap();
      }();
 
      static constexpr unsigned int n_precomputed_values = 3;
 
      using precomputed_type = std::array<Number, n_precomputed_values>;
 
      static inline const auto precomputed_names =
          std::array<std::string, n_precomputed_values>{{"e", "p", "a"}};
 
      static constexpr unsigned int n_initial_precomputed_values = 0;
 
      using initial_precomputed_type =
          std::array<Number, n_initial_precomputed_values>;
 
      static inline const auto initial_precomputed_names =
          std::array<std::string, n_initial_precomputed_values>{};
 
      using StateVector = Vectors::
          StateVector<ScalarNumber, problem_dimension, n_precomputed_values>;
 
      using HyperbolicVector =
          Vectors::MultiComponentVector<ScalarNumber, problem_dimension>;
 
      using PrecomputedVector =
          Vectors::MultiComponentVector<ScalarNumber, n_precomputed_values>;
 
      using PrecomputedVectorView = Vectors::MultiComponentVectorView<
          ScalarNumber,
          n_precomputed_values,
          dealii::VectorizedArray<ScalarNumber>::size(),
          dealii::MemorySpace::Host,
          /*writable=*/false>;
 
      using InitialPrecomputedVector =
          Vectors::MultiComponentVector<ScalarNumber,
                                        n_initial_precomputed_values>;
 
      using InitialPrecomputedVectorView = Vectors::MultiComponentVectorView<
          ScalarNumber,
          n_initial_precomputed_values,
          dealii::VectorizedArray<ScalarNumber>::size(),
          dealii::MemorySpace::Host,
          /*writable=*/false>;
 
 
 
      HyperbolicSystemView(const HyperbolicSystem &hyperbolic_system)
          : hyperbolic_system_(hyperbolic_system)
      {
      }
 
 
 
      DEAL_II_ALWAYS_INLINE inline const std::string &
      barotropic_equation_of_state() const
      {
        return hyperbolic_system_.barotropic_equation_of_state_;
      }
 
      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 Number
      beos_specific_internal_energy(const Number &rho) const
      {
        const auto &beos =
            hyperbolic_system_.selected_barotropic_equation_of_state_;
 
        if constexpr (std::is_same_v<ScalarNumber, Number>) {
          return ScalarNumber(beos->specific_internal_energy(rho));
        } else {
          Number e;
          for (unsigned int k = 0; k < Number::size(); ++k) {
            e[k] = ScalarNumber(beos->specific_internal_energy(rho[k]));
          }
          return e;
        }
      }
 
      DEAL_II_ALWAYS_INLINE inline Number beos_pressure(const Number &rho) const
      {
        const auto &beos =
            hyperbolic_system_.selected_barotropic_equation_of_state_;
 
        if constexpr (std::is_same_v<ScalarNumber, Number>) {
          return ScalarNumber(beos->pressure(rho));
        } else {
          Number p;
          for (unsigned int k = 0; k < Number::size(); ++k) {
            p[k] = ScalarNumber(beos->pressure(rho[k]));
          }
          return p;
        }
      }
 
      DEAL_II_ALWAYS_INLINE inline Number
      beos_speed_of_sound(const Number &rho) const
      {
        const auto &beos =
            hyperbolic_system_.selected_barotropic_equation_of_state_;
 
        if constexpr (std::is_same_v<ScalarNumber, Number>) {
          return ScalarNumber(beos->speed_of_sound(rho));
        } else {
          Number a;
          for (unsigned int k = 0; k < Number::size(); ++k) {
            a[k] = ScalarNumber(beos->speed_of_sound(rho[k]));
          }
          return a;
        }
      }
 
 
 
      static constexpr bool have_gamma = false;
      static constexpr bool have_covolume_constant = false;
      static constexpr bool have_energy_equation = false;
 
 
 
      static Number density(const state_type &U);
 
      Number filter_vacuum_density(const Number &rho) const;
 
      static dealii::Tensor<1, dim, Number> momentum(const state_type &U);
 
      Number total_energy(const state_type &U,
                          const Number &specific_internal_energy) const;
 
      state_type total_energy_derivative(const state_type &U,
                                         const Number &specific_internal_energy,
                                         const Number &pressure) const;
 
      bool is_admissible(const state_type &U) const;
 
 
 
      template <int component>
      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;
 
      template <typename Lambda>
      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;
 
 
 
      flux_type f(const state_type &U, const Number &p) const;
 
      flux_contribution_type
      flux_contribution(const PrecomputedVectorView &pv,
                        const InitialPrecomputedVectorView &piv,
                        const unsigned int i,
                        const state_type &U_i) const;
 
      flux_contribution_type
      flux_contribution(const PrecomputedVectorView &pv,
                        const InitialPrecomputedVectorView &piv,
                        const unsigned int *js,
                        const state_type &U_j) const;
 
      state_type
      flux_divergence(const flux_contribution_type &flux_i,
                      const flux_contribution_type &flux_j,
                      const dealii::Tensor<1, dim, Number> &c_ij) const;
 
      static constexpr bool have_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;
 
 
 
      static constexpr bool have_source_terms = false;
 
      state_type nodal_source(const PrecomputedVectorView &pv,
                              const unsigned int i,
                              const state_type &U_i,
                              const ScalarNumber tau) const = delete;
 
      state_type nodal_source(const PrecomputedVectorView &pv,
                              const unsigned int *js,
                              const state_type &U_j,
                              const ScalarNumber tau) const = delete;
 
 
 
      template <typename ST>
      state_type expand_state(const ST &state) const;
 
      template <typename ST>
      state_type from_initial_state(const ST &initial_state) const;
 
      state_type from_primitive_state(const state_type &primitive_state) const;
 
      state_type to_primitive_state(const state_type &state) const;
 
      template <typename Lambda>
      state_type apply_galilei_transform(const state_type &state,
                                         const Lambda &lambda) const;
 
    private:
 
 
      const HyperbolicSystem &hyperbolic_system_;
 
    }; /* HyperbolicSystemView */
 
 
    /*
     * -------------------------------------------------------------------------
     * Inline definitions
     * -------------------------------------------------------------------------
     */
 
 
    inline HyperbolicSystem::HyperbolicSystem(
        const std::string &subsection /*= "HyperbolicSystem"*/)
        : ParameterAcceptor(subsection)
    {
      barotropic_equation_of_state_ = "isothermal";
      add_parameter("barotropic equation of state",
                    barotropic_equation_of_state_,
                    "The barotropic equation of state. Valid names are given "
                    "by any of the subsections defined below");
 
      reference_density_ = 1.;
      add_parameter("reference density",
                    reference_density_,
                    "Problem specific density reference");
 
      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 equation of state list with all equation of
       * state configurations defined in the EquationOfState namespace:
       */
      BarotropicEquationOfStateLibrary::populate_equation_of_state_list(
          barotropic_equation_of_state_list_, subsection);
 
      const auto populate_functions = [this]() {
        bool initialized = false;
        for (auto &it : barotropic_equation_of_state_list_)
 
          /* Populate EOS-specific quantities and functions */
          if (it->name() == barotropic_equation_of_state_) {
            selected_barotropic_equation_of_state_ = it;
            problem_name = "Compressible Euler equations (" + it->name() +
                           " EOS, optimized barotropic)";
            initialized = true;
            break;
          }
 
        AssertThrow(initialized,
                    dealii::ExcMessage("Could not find a barotropic equation "
                                       "of state description with name \"" +
                                       barotropic_equation_of_state_ + "\""));
      };
 
      ParameterAcceptor::parse_parameters_call_back.connect(populate_functions);
      populate_functions();
    }
 
 
    template <int dim, typename ScalarNumber>
    inline void HyperbolicSystem::fill_precomputed_values(
        const OfflineData<dim, ScalarNumber> &offline_data,
        typename HyperbolicSystemView<dim, ScalarNumber>::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_view =
          offline_data.sparsity_pattern_simd().view();
      using VA = dealii::VectorizedArray<ScalarNumber>;
 
      const auto U_view = std::get<0>(state_vector).view();
      const auto precomputed_view = std::get<1>(state_vector).view();
 
      const auto body = [&](auto sentinel, unsigned int i) {
        using T = decltype(sentinel);
        using View = HyperbolicSystemView<dim, T>;
        using precomputed_type = typename View::precomputed_type;
 
        const unsigned int row_length = sparsity_simd_view.row_length(i);
        if (skip_constrained_dofs && row_length == 1)
          return;
 
        const auto U_i = U_view.template read_tensor<T>(i);
        const auto view = this->view<dim, T>();
        const auto rho_i = view.density(U_i);
 
        const auto e_i = view.beos_specific_internal_energy(rho_i);
        const auto p_i = view.beos_pressure(rho_i);
        const auto a_i = view.beos_speed_of_sound(rho_i);
 
        const precomputed_type prec_i{e_i, p_i, a_i};
        precomputed_view.template write_tensor<T>(prec_i, i);
      };
 
      cpu_simd_loop<ScalarNumber>("time_step_1", body, 0, n_internal, n_owned);
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline Number
    HyperbolicSystemView<dim, Number>::density(const state_type &U)
    {
      return U[0];
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline Number
    HyperbolicSystemView<dim, Number>::filter_vacuum_density(
        const Number &rho) const
    {
      constexpr ScalarNumber eps = std::numeric_limits<ScalarNumber>::epsilon();
      const Number rho_cutoff_large =
          reference_density() * vacuum_state_relaxation_large() * eps;
 
      return dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
          std::abs(rho), rho_cutoff_large, Number(0.), rho);
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline dealii::Tensor<1, dim, Number>
    HyperbolicSystemView<dim, Number>::momentum(const state_type &U)
    {
      dealii::Tensor<1, dim, Number> result;
      for (unsigned int i = 0; i < dim; ++i)
        result[i] = U[1 + i];
      return result;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline Number
    HyperbolicSystemView<dim, Number>::total_energy(
        const state_type &U, const Number &specific_internal_energy) const
    {
      const auto rho = density(U);
      const auto rho_inverse = ScalarNumber(1.) / rho;
      const auto m = momentum(U);
 
      return rho * specific_internal_energy +
             ScalarNumber(0.5) * rho_inverse * m.norm_square();
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::total_energy_derivative(
        const state_type &U,
        const Number &specific_internal_energy,
        const Number &pressure) const -> state_type
    {
      const auto rho = density(U);
      const auto rho_inverse = ScalarNumber(1.) / rho;
      const auto m = momentum(U);
 
      state_type result;
 
      result[0] =
          specific_internal_energy + rho_inverse * pressure -
          ScalarNumber(0.5) * rho_inverse * rho_inverse * m.norm_square();
      for (unsigned int i = 0; i < dim; ++i)
        result[1 + i] = rho_inverse * m[i];
 
      return result;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline bool
    HyperbolicSystemView<dim, Number>::is_admissible(const state_type &U) const
    {
      const auto rho = density(U);
      constexpr auto gt = dealii::SIMDComparison::greater_than;
      using T = Number;
      const auto test =
          dealii::compare_and_apply_mask<gt>(rho, 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";
      }
#endif
 
      return (test == Number(0.));
    }
 
 
    template <int dim, typename Number>
    template <int component>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::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
    {
      // FIXME
      AssertThrow(false, dealii::ExcNotImplemented());
      __builtin_trap();
      return state_type{};
    }
 
 
    template <int dim, typename Number>
    template <typename Lambda>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::apply_boundary_conditions(
        dealii::types::boundary_id id,
        const state_type &U,
        const dealii::Tensor<1, dim, Number> &normal,
        const Lambda &get_dirichlet_data) const -> state_type
    {
      state_type result = U;
 
      if (id == Boundary::dirichlet) {
        result = get_dirichlet_data();
 
      } else if (id == Boundary::dirichlet_momentum) {
        /* Only enforce Dirichlet conditions on the momentum: */
        auto m_dirichlet = momentum(get_dirichlet_data());
        for (unsigned int k = 0; k < dim; ++k)
          result[k + 1] = m_dirichlet[k];
 
      } else if (id == Boundary::dirichlet_velocity) {
        /* Only enforce Dirichlet conditions on the velocity: */
        const auto U_dirichlet = get_dirichlet_data();
        const auto rho_dirichlet = density(U_dirichlet);
        const auto v_dirichlet = momentum(U_dirichlet) / rho_dirichlet;
        const auto rho = density(result);
        for (unsigned int k = 0; k < dim; ++k)
          result[k + 1] = rho * v_dirichlet[k];
 
      } else if (id == Boundary::slip) {
        auto m = momentum(U);
        m -= 1. * (m * normal) * normal;
        for (unsigned int k = 0; k < dim; ++k)
          result[k + 1] = m[k];
 
      } else if (id == Boundary::no_slip) {
        for (unsigned int k = 0; k < dim; ++k)
          result[k + 1] = 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 and speed of sound oracle:
         */
        const auto p = beos_pressure(rho);
        const auto a = beos_speed_of_sound(rho);
        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 rho_dirichlet = density(U_dirichlet);
          const auto p_dirichlet = beos_pressure(rho_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 rho_dirichlet = density(U_dirichlet);
          const auto p_dirichlet = beos_pressure(rho_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 <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::f(const state_type &U,
                                         const Number &p) const -> flux_type
    {
      const auto rho_inverse = ScalarNumber(1.) / density(U);
      const auto m = momentum(U);
 
      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;
      }
 
      return result;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::flux_contribution(
        const PrecomputedVectorView &pv,
        const InitialPrecomputedVectorView & /*piv*/,
        const unsigned int i,
        const state_type &U_i) const -> flux_contribution_type
    {
      const auto &[e_i, p_i, a_i] =
          pv.template read_tensor<Number, precomputed_type>(i);
      return f(U_i, p_i);
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::flux_contribution(
        const PrecomputedVectorView &pv,
        const InitialPrecomputedVectorView & /*piv*/,
        const unsigned int *js,
        const state_type &U_j) const -> flux_contribution_type
    {
      const auto &[e_j, p_j, a_j] =
          pv.template read_tensor<Number, precomputed_type>(js);
      return f(U_j, p_j);
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::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 <int dim, typename Number>
    template <typename ST>
    auto HyperbolicSystemView<dim, Number>::expand_state(const ST &state) const
        -> state_type
    {
      using T = typename ST::value_type;
      static_assert(std::is_same_v<Number, T>, "template mismatch");
 
      constexpr auto dim2 = ST::dimension - 1;
      static_assert(dim >= dim2,
                    "the space dimension of the argument state must not be "
                    "larger than the one of the target state");
 
      state_type result;
      result[0] = state[0];
      for (unsigned int i = 1; i < dim2 + 1; ++i)
        result[i] = state[i];
 
      return result;
    }
 
 
    template <int dim, typename Number>
    template <typename ST>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::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 <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::from_primitive_state(
        const state_type &primitive_state) const -> state_type
    {
      const auto rho = density(primitive_state);
 
      auto state = primitive_state;
      /* Fix up momentum: */
      for (unsigned int i = 1; i < dim + 1; ++i)
        state[i] *= rho;
 
      return state;
    }
 
 
    template <int dim, typename Number>
    DEAL_II_ALWAYS_INLINE inline auto
    HyperbolicSystemView<dim, Number>::to_primitive_state(
        const state_type &state) const -> state_type
    {
      const auto rho = density(state);
      const auto rho_inverse = Number(1.) / rho;
 
      auto primitive_state = state;
      /* Fix up velocity: */
      for (unsigned int i = 1; i < dim + 1; ++i)
        primitive_state[i] *= rho_inverse;
 
      return primitive_state;
    }
 
 
    template <int dim, typename Number>
    template <typename Lambda>
    auto HyperbolicSystemView<dim, Number>::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[1 + d] = M[d];
      return result;
    }
  } // namespace EulerBarotropic
} // namespace ryujin