equation_of_state_bumpy_barotropic_pressure.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
// Copyright (C) 2026 by the ryujin authors
//
 
#pragma once
 
#include <ryujin/base/compile_time_options.h>
#include <ryujin/euler_aeos/equation_of_state.h>
 
namespace ryujin
{
  namespace EquationsOfState
  {
    class BumpyBarotropicPressure : public EquationOfState
    {
    public:
      using EquationOfState::pressure;
      using EquationOfState::specific_internal_energy;
      using EquationOfState::speed_of_sound;
      using EquationOfState::temperature;
 
      BumpyBarotropicPressure(const std::string &subsection)
          : EquationOfState("bumpy barotropic pressure", subsection)
      {
        /*
         * Parameters of the hydrodynamical (Noble-Abel stiffened gas)
         * constituent:
         */
 
        gamma_ = 7. / 5.;
        this->add_parameter("gamma", gamma_, "The ratio of specific heats");
 
        /*
         * R is the specific gas constant with units [J / (Kg K)]. More details
         * can be found at:
         * https://en.wikipedia.org/wiki/Gas_constant#Specific_gas_constant
         */
        R_ = 287.052874;
        this->add_parameter(
            "gas constant R", R_, "The specific gas constant R");
 
        cv_ = R_ / (gamma_ - 1.);
 
        b_ = 0.;
        this->add_parameter(
            "covolume b", b_, "The maximum compressibility constant");
 
        q_ = 0.;
        this->add_parameter("reference specific internal energy",
                            q_,
                            "The reference specific internal energy");
 
        pinf_ = 0.;
        this->add_parameter(
            "reference pressure", pinf_, "The reference pressure p infinity");
 
        s0_ = 0.;
        this->add_parameter("reference specific entropy",
                            s0_,
                            "The reference specific entropy");
 
        /*
         * Parameters of the barotropic constituent:
         */
 
        rho_0_ = 1.;
        this->add_parameter("barotropic reference density",
                            rho_0_,
                            "The density rho_0 at which the barotropic speed "
                            "of sound is centered");
 
        c_0_ = 1.;
        this->add_parameter("barotropic sound speed",
                            c_0_,
                            "The strength c_0 of the barotropic speed of "
                            "sound");
 
        eps_ = 1.e-2;
        this->add_parameter("barotropic bump width",
                            eps_,
                            "The width epsilon of the barotropic bump");
 
        e_0_ = 0.;
        this->add_parameter("barotropic reference specific internal energy",
                            e_0_,
                            "The (arbitrary) additive constant e_0 of the "
                            "barotropic specific internal energy");
 
        /*
         * Update the EOS interpolation parameters on parameter read in
         * and specific heat at constant volume:
         *
         * @note The interpolation parameters are the ones of the
         * hydrodynamical constituent. They drive the surrogate NASG
         * interpolation performed in HyperbolicSystemView
         * (surrogate_gamma(), surrogate_pressure(),
         * surrogate_specific_entropy(), ...) which is unaware of the
         * barotropic contribution.
         */
        ParameterAcceptor::parse_parameters_call_back.connect([this] {
          this->covolume_constant_ = b_;
          this->interpolation_pinfty_ = pinf_;
          this->interpolation_q_ = q_;
          cv_ = R_ / (gamma_ - 1.);
        });
      }
 
      double pressure(double rho, double e) const final
      {
        const auto e_h = e - barotropic_specific_internal_energy(rho);
        return hydrodynamic_pressure(rho, e_h) + barotropic_pressure(rho);
      }
 
 
      double specific_internal_energy(double rho, double p) const final
      {
        const auto p_h = p - barotropic_pressure(rho);
        return hydrodynamic_specific_internal_energy(rho, p_h) +
               barotropic_specific_internal_energy(rho);
      }
 
      double temperature(double rho, double e) const final
      {
        const auto e_h = e - barotropic_specific_internal_energy(rho);
        return hydrodynamic_temperature(rho, e_h);
      }
 
      double cold_curve_bound(double rho) const final
      {
        return hydrodynamic_cold_curve_bound(rho) +
               barotropic_specific_internal_energy(rho);
      }
 
      double specific_entropy(double rho, double e) const final
      {
        const auto e_h = e - barotropic_specific_internal_energy(rho);
        return hydrodynamic_specific_entropy(rho, e_h);
      }
 
      double speed_of_sound(double rho, double e) const final
      {
        const auto e_h = e - barotropic_specific_internal_energy(rho);
        return std::sqrt(hydrodynamic_sound_speed_squared(rho, e_h) +
                         barotropic_sound_speed_squared(rho));
      }
 
    private:
      double barotropic_sound_speed_squared(double rho) const
      {
        const auto radicand = (rho - rho_0_) * (rho - rho_0_) + eps_;
        return c_0_ * c_0_ / (radicand * std::sqrt(radicand));
      }
 
      double barotropic_pressure(double rho) const
      {
        const auto s = std::sqrt((rho - rho_0_) * (rho - rho_0_) + eps_);
        const auto s_0 = std::sqrt(rho_0_ * rho_0_ + eps_);
        return c_0_ * c_0_ / eps_ * ((rho - rho_0_) / s + rho_0_ / s_0);
      }
 
      double barotropic_specific_internal_energy(double rho) const
      {
        const auto s = std::sqrt((rho - rho_0_) * (rho - rho_0_) + eps_);
        const auto s_0_squared = rho_0_ * rho_0_ + eps_;
        const auto s_0 = std::sqrt(s_0_squared);
 
        const auto first_term =
            c_0_ * c_0_ * rho_0_ * (s - s_0) / (eps_ * s_0_squared * rho);
 
        const auto argument =
            (2. * rho * rho_0_ - s_0_squared) / ((rho + s) * s_0);
        const auto second_term =
            2. * c_0_ * c_0_ / (s_0_squared * s_0) * std::atanh(argument);
 
        return first_term + second_term + e_0_;
      }
 
      double hydrodynamic_pressure(double rho, double e_h) const
      {
        return (gamma_ - 1.) * rho * (e_h - q_) / (1. - b_ * rho) -
               gamma_ * pinf_;
      }
 
      double hydrodynamic_specific_internal_energy(double rho, double p_h) const
      {
        const auto numerator = (p_h + gamma_ * pinf_) * (1. - b_ * rho);
        const auto denominator = rho * (gamma_ - 1.);
        return q_ + numerator / denominator;
      }
 
      double hydrodynamic_temperature(double rho, double e_h) const
      {
        return (e_h - q_ - pinf_ * (1. / rho - b_)) / cv_;
      }
 
      double hydrodynamic_cold_curve_bound(double rho) const
      {
        return q_ + pinf_ * (1. / rho - b_);
      }
 
      double hydrodynamic_specific_entropy(double rho, double e_h) const
      {
        const auto covolume_term = 1. / rho - b_;
        const auto p_plus_pinf = (gamma_ - 1.) *
                                 ((e_h - q_) - pinf_ * covolume_term) /
                                 covolume_term;
        const auto first_term = cv_ * std::log(p_plus_pinf);
        const auto second_term =
            cv_ * gamma_ * std::log((gamma_ - 1.) * cv_ / covolume_term);
        return first_term - second_term + s0_;
      }
 
      double hydrodynamic_sound_speed_squared(double rho, double e_h) const
      {
        const auto covolume = 1. - b_ * rho;
        auto result =
            (rho * (e_h - q_) - pinf_ * covolume) / (covolume * covolume * rho);
        result *= gamma_ * (gamma_ - 1.);
        return result;
      }
 
      double gamma_;
      double R_;
      double cv_;
      double b_;
      double q_;
      double pinf_;
      double s0_;
 
      double rho_0_;
      double c_0_;
      double eps_;
      double e_0_;
    };
  } // namespace EquationsOfState
} // namespace ryujin