smash::HadronGasEos Class Reference

#include <hadgas_eos.h>

Class to handle the equation of state (EoS) of the hadron gas, consisting of all hadrons included in SMASH.

This implementation deals with an ideal Boltzmann gas and allows to compute:

• energy density \(\epsilon\), pressure \(p\), density \(n\), net baryon density \(n_B\), net strangeness \(n_S\) and net charge density \(n_Q\) as a function of temperature \(T\), baryon chemical potential \(\mu_B\), strange chemical potential \(\mu_S\) and charge chemical potential \(\mu_Q\).
• Temperature and chemical potentials given energy-, net baryon- and net strangeness density. This requires solving a system of nonlinear equations.

Definition at line 124 of file hadgas_eos.h.

Collaboration diagram for smash::HadronGasEos:

Classes
struct	eparams
	Another structure for passing energy density to the gnu library. More...
struct	rparams
	A structure for passing equation parameters to the gnu library. More...

Public Member Functions
	HadronGasEos (bool tabulate, bool account_for_widths)
	Constructor of HadronGasEos. More...
	~HadronGasEos ()
std::array< double, 4 >	solve_eos (double e, double nb, double ns, double nq, std::array< double, 4 > initial_approximation)
	Compute temperature and chemical potentials given energy-, net baryon-, net strangeness- and net charge density and an inital approximation. More...
std::array< double, 4 >	solve_eos (double e, double nb, double ns, double nq)
	Compute temperature and chemical potentials given energy-, net baryon-, net strangeness- and net charge density without an inital approximation. More...
std::array< double, 4 >	solve_eos_initial_approximation (double e, double nb, double nq)
	Compute a reasonable initial approximation for solve_eos. More...
void	from_table (EosTable::table_element &res, double e, double nb, double nq) const
	Get the element of eos table. More...
bool	is_tabulated () const
	Create an EoS table or not? More...
bool	account_for_resonance_widths () const
	If resonance spectral functions are taken into account. More...

Static Public Member Functions
static double	energy_density (double T, double mub, double mus, double muq)
	Compute energy density. More...
static double	density (double T, double mub, double mus, double muq, bool account_for_resonance_widths=false)
	Compute particle number density. More...
static double	pressure (double T, double mub, double mus, double muq, bool account_for_resonance_widths=false)
	Compute pressure \( p = n T \). More...
static double	net_baryon_density (double T, double mub, double mus, double muq, bool account_for_resonance_widths=false)
	Compute net baryon density. More...
static double	net_strange_density (double T, double mub, double mus, double muq, bool account_for_resonance_widths=false)
	Compute net strangeness density. More...
static double	net_charge_density (double T, double mub, double mus, double muq, bool account_for_resonance_widths=false)
	Compute net charge density. More...
static double	partial_density (const ParticleType &ptype, double T, double mub, double mus, double muq, bool account_for_resonance_widths=false)
	Compute partial density of one hadron sort. More...
static double	sample_mass_thermal (const ParticleType &ptype, double beta)
	Sample resonance mass in a thermal medium. More...
static double	mus_net_strangeness0 (double T, double mub, double muq)
	Compute strangeness chemical potential, requiring that net strangeness = 0. More...
static bool	is_eos_particle (const ParticleType &ptype)
	Check if a particle belongs to the EoS. More...

Private Member Functions
std::string	print_solver_state (size_t iter) const
	Helpful printout, useful for debugging if gnu equation solving goes crazy. More...

Static Private Member Functions
static double	scaled_partial_density_auxiliary (double m_over_T, double mu_over_T)
	Function used to avoid duplications in density calculations. More...
static double	scaled_partial_density (const ParticleType &ptype, double beta, double mub, double mus, double muq, bool account_for_width=false)
	Compute (unnormalized) density of one hadron sort - helper functions used to reduce code duplication. More...
static int	set_eos_solver_equations (const gsl_vector *x, void *params, gsl_vector *f)
	Interface EoS equations to be solved to gnu library. More...
static double	e_equation (double T, void *params)

Private Attributes
EosTable	eos_table_ = EosTable(1.e-1, 1.e-1, 1.e-1, 90, 90, 90)
	EOS Table to be used. More...
gsl_vector *	x_
	Variables used by gnu equation solver. More...
gsl_multiroot_fsolver *	solver_
const bool	tabulate_
	Create an EoS table or not? More...
const bool	account_for_resonance_widths_
	Use pole masses of resonances or integrate over spectral functions. More...

Static Private Attributes
static constexpr double	prefactor_
	Constant factor, that appears in front of many thermodyn. expressions. More...
static constexpr double	tolerance_ = 1.e-8
	Precision of equation solving. More...
static constexpr size_t	n_equations_ = 4
	Number of equations in the system of equations to be solved. More...

Constructor & Destructor Documentation

HadronGasEos()

smash::HadronGasEos::HadronGasEos	(	bool	tabulate,
		bool	account_for_widths
	)

Constructor of HadronGasEos.

Parameters

[in]	tabulate	Whether the equation of state should be tabulated Tabulation takes time once (typically around 5 minutes), but makes the further usage of the class much faster. Tabulated values are saved in a file and loaded at the next run.
[in]	account_for_widths	Whether equation of state should account for resonance spectral functions. Normally one wants to do it, if HadronGasEos is used for density calculations, for example in the box initialization. However, it is not recommended to account for spectral functions if EoS is tabulated (tabulate = true), because this makes tabulation incredibly slow. Therefore, for HadronGasEos to be used in thermalizer, this option has to be false. Also note that presently width account is not implemented for energy density calculation.

Definition at line 200 of file hadgas_eos.cc.

    : x_(gsl_vector_alloc(n_equations_)),
      tabulate_(tabulate),
      account_for_resonance_widths_(account_for_width) {
  const gsl_multiroot_fsolver_type *solver_type;
  solver_type = gsl_multiroot_fsolver_hybrid;
  solver_ = gsl_multiroot_fsolver_alloc(solver_type, n_equations_);
  if (tabulate_ && account_for_resonance_widths_) {
    logg[LResonances].error(
        "Compilation of hadron gas EoS table requested with"
        " account of resonance spectral functions. This is not "
        "advised, as it will likely take a few days to finish."
        " Besides, the effect of resonance widths is currently not "
        "implemented for energy density computation, so the computed"
        " table will be inconsistent anyways.");
  }
  if (tabulate_) {
    eos_table_.compile_table(*this);
  }
}

Here is the call graph for this function:

~HadronGasEos()

smash::HadronGasEos::~HadronGasEos

(

)

Definition at line 221 of file hadgas_eos.cc.

                            {
  gsl_multiroot_fsolver_free(solver_);
  gsl_vector_free(x_);
}

Member Function Documentation

energy_density()

double smash::HadronGasEos::energy_density ( double  T, double  mub, double  mus, double  muq  )

static

Compute energy density.

Grand-canonical Boltzmann ideal gas, consisting of all hadrons in SMASH:

\[ \epsilon = \sum \frac{g_i m_i^2 T^2}{2\pi^2(\hbar c)^3} exp \left(\frac{\mu_B B_i + \mu_S S_i + \mu_Q Q_i}{T} \right) \times \left[ 3 K_2\left( \frac{m_i}{T}\right) + \frac{m_i}{T} K_1\left( \frac{m_i}{T}\right)\right] \]

Parameters

[in]	T	temperature [GeV]
[in]	mub	baryon chemical potential [GeV]
[in]	mus	strangeness chemical potential [GeV]
[in]	muq	charge chemical potential [GeV]

Returns

energy density e [GeV/fm \(^3\)]

Definition at line 284 of file hadgas_eos.cc.

                                                {
  if (T < really_small) {
    return 0.0;
  }
  const double beta = 1.0 / T;
  double e = 0.0;
  for (const ParticleType &ptype : ParticleType::list_all()) {
    if (!is_eos_particle(ptype)) {
      continue;
    }
    const double z = ptype.mass() * beta;
    double x = beta * (mub * ptype.baryon_number() + mus * ptype.strangeness() +
                       muq * ptype.charge() - ptype.mass());
    if (x < -500.0) {
      return 0.0;
    }
    x = std::exp(x);
    const size_t g = ptype.spin() + 1;
    // Small mass case, z*z*K_2(z) -> 2, z*z*z*K_1(z) -> 0 at z->0
    e += (z < really_small) ? 3.0 * g * x
                            : z * z * g * x *
                                  (3.0 * gsl_sf_bessel_Kn_scaled(2, z) +
                                   z * gsl_sf_bessel_K1_scaled(z));
  }
  e *= prefactor_ * T * T * T * T;
  return e;
}

Here is the call graph for this function:

Here is the caller graph for this function:

density()

double smash::HadronGasEos::density ( double  T, double  mub, double  mus, double  muq, bool  account_for_resonance_widths = false  )

static

Compute particle number density.

Grand-canonical Boltzmann ideal gas, consisting of all hadrons in SMASH:

\[ n = \sum \frac{g_i m_i^2 T}{2\pi^2(\hbar c)^3} exp \left(\frac{\mu_B B_i + \mu_S S_i + \mu_Q Q_i}{T} \right) K_2\left( \frac{m_i}{T}\right) \]

Parameters

[in]	T	temperature [GeV]
[in]	mub	baryon chemical potential [GeV]
[in]	mus	strangeness chemical potential [GeV]
[in]	muq	charge chemical potential [GeV]
[in]	account_for_resonance_widths	if false, pole masses are used; if true, then integration over spectral function is included

Returns

particle number density n [fm \(^{-3}\)]

Definition at line 313 of file hadgas_eos.cc.

                                                     {
  if (T < really_small) {
    return 0.0;
  }
  const double beta = 1.0 / T;
  double rho = 0.0;
  for (const ParticleType &ptype : ParticleType::list_all()) {
    if (!is_eos_particle(ptype)) {
      continue;
    }
    rho +=
        scaled_partial_density(ptype, beta, mub, mus, muq, account_for_width);
  }
  rho *= prefactor_ * T * T * T;
  return rho;
}

Here is the call graph for this function:

Here is the caller graph for this function:

pressure()

static double smash::HadronGasEos::pressure ( double  T, double  mub, double  mus, double  muq, bool  account_for_resonance_widths = false  )

inlinestatic

Compute pressure \( p = n T \).

Parameters

[in]	T	temperature [GeV]
[in]	mub	baryon chemical potential [GeV]
[in]	mus	strangeness chemical potential [GeV]
[in]	muq	charge chemical potential [GeV]
[in]	account_for_resonance_widths	if false, pole masses are used; if true, then integration over spectral function is included

Returns

pressure p [GeV/fm \(^{-3}\)]

Definition at line 194 of file hadgas_eos.h.

                                                                    {
    return T * density(T, mub, mus, muq, account_for_resonance_widths);
  }

Here is the call graph for this function:

Here is the caller graph for this function:

net_baryon_density()

double smash::HadronGasEos::net_baryon_density ( double  T, double  mub, double  mus, double  muq, bool  account_for_resonance_widths = false  )

static

Compute net baryon density.

Grand-canonical Boltzmann ideal gas, consisting of all hadrons in SMASH:

\[ n_B = \sum B_i \frac{g_i m_i^2 T}{2\pi^2(\hbar c)^3} exp \left(\frac{\mu_B B_i + \mu_S S_i + \mu_Q Q_i}{T} \right) K_2\left( \frac{m_i}{T}\right) \]

Parameters

[in]	T	temperature [GeV]
[in]	mub	baryon chemical potential [GeV]
[in]	mus	strangeness chemical potential [GeV]
[in]	muq	charge chemical potential [GeV]
[in]	account_for_resonance_widths	if false, pole masses are used; if true, then integration over spectral function is included

Returns

net baryon density \(n_B\) [fm \(^{-3}\)]

Definition at line 331 of file hadgas_eos.cc.

                                                                            {
  if (T < really_small) {
    return 0.0;
  }
  const double beta = 1.0 / T;
  double rho = 0.0;
  for (const ParticleType &ptype : ParticleType::list_all()) {
    if (!ptype.is_baryon() || !is_eos_particle(ptype)) {
      continue;
    }
    rho +=
        scaled_partial_density(ptype, beta, mub, mus, muq, account_for_width) *
        ptype.baryon_number();
  }
  rho *= prefactor_ * T * T * T;
  return rho;
}

Here is the call graph for this function:

Here is the caller graph for this function:

net_strange_density()

double smash::HadronGasEos::net_strange_density ( double  T, double  mub, double  mus, double  muq, bool  account_for_resonance_widths = false  )

static

Compute net strangeness density.

Grand-canonical Boltzmann ideal gas, consisting of all hadrons in SMASH:

\[ n_S = \sum S_i \frac{g_i m_i^2 T}{2\pi^2(\hbar c)^3} exp \left(\frac{\mu_B B_i + \mu_S S_i + \mu_Q Q_i}{T} \right) K_2\left( \frac{m_i}{T}\right) \]

Parameters

[in]	T	temperature [GeV]
[in]	mub	baryon chemical potential [GeV]
[in]	mus	strangeness chemical potential [GeV]
[in]	muq	charge chemical potential [GeV]
[in]	account_for_resonance_widths	if false, pole masses are used; if true, then integration over spectral function is included

Returns

net strangeness density density \(n_S\) [fm \(^{-3}\)]

Definition at line 350 of file hadgas_eos.cc.

                                                                             {
  if (T < really_small) {
    return 0.0;
  }
  const double beta = 1.0 / T;
  double rho = 0.0;
  for (const ParticleType &ptype : ParticleType::list_all()) {
    if (ptype.strangeness() == 0 || !is_eos_particle(ptype)) {
      continue;
    }
    rho +=
        scaled_partial_density(ptype, beta, mub, mus, muq, account_for_width) *
        ptype.strangeness();
  }
  rho *= prefactor_ * T * T * T;
  return rho;
}

Here is the call graph for this function:

Here is the caller graph for this function:

net_charge_density()

double smash::HadronGasEos::net_charge_density ( double  T, double  mub, double  mus, double  muq, bool  account_for_resonance_widths = false  )

static

Compute net charge density.

Grand-canonical Boltzmann ideal gas, consisting of all hadrons in SMASH:

\[ n_Q = \sum Q_i \frac{g_i m_i^2 T}{2\pi^2(\hbar c)^3} exp \left(\frac{\mu_B B_i + \mu_S S_i + \mu_Q Q_i}{T} \right) K_2\left( \frac{m_i}{T}\right) \]

Parameters

[in]	T	temperature [GeV]
[in]	mub	baryon chemical potential [GeV]
[in]	mus	strangeness chemical potential [GeV]
[in]	muq	charge chemical potential [GeV]
[in]	account_for_resonance_widths	if false, pole masses are used; if true, then integration over spectral function is included

Returns

net charge density density \(n_S\) [fm \(^{-3}\)]

Definition at line 369 of file hadgas_eos.cc.

                                                                            {
  if (T < really_small) {
    return 0.0;
  }
  const double beta = 1.0 / T;
  double rho = 0.0;
  for (const ParticleType &ptype : ParticleType::list_all()) {
    if (ptype.charge() == 0 || !is_eos_particle(ptype)) {
      continue;
    }
    rho +=
        scaled_partial_density(ptype, beta, mub, mus, muq, account_for_width) *
        ptype.charge();
  }
  rho *= prefactor_ * T * T * T;
  return rho;
}

Here is the call graph for this function:

Here is the caller graph for this function:

partial_density()

double smash::HadronGasEos::partial_density ( const ParticleType &  ptype, double  T, double  mub, double  mus, double  muq, bool  account_for_resonance_widths = false  )

static

Compute partial density of one hadron sort.

Grand-canonical Boltzmann ideal gas:

\[ n = \frac{g m^2 T}{2\pi^2(\hbar c)^3} exp \left(\frac{\mu_B B + \mu_S S + \mu_Q Q_i}{T} \right) K_2\left( \frac{m}{T}\right) \]

Parameters

[in]	ptype	the hadron sort, for which partial density is computed
[in]	T	temperature [GeV]
[in]	mub	baryon chemical potential [GeV]
[in]	mus	strangeness chemical potential [GeV]
[in]	muq	charge chemical potential [GeV]
[in]	account_for_resonance_widths	if false, pole masses are used; if true, then integration over spectral function is included

Returns

partial density of the given hadron sort \(n\) [fm \(^{-3}\)]

Definition at line 273 of file hadgas_eos.cc.

                                                             {
  if (T < really_small) {
    return 0.0;
  }
  return prefactor_ * T * T * T *
         scaled_partial_density(ptype, 1.0 / T, mub, mus, muq,
                                account_for_width);
}

Here is the call graph for this function:

Here is the caller graph for this function:

sample_mass_thermal()

double smash::HadronGasEos::sample_mass_thermal ( const ParticleType &  ptype, double  beta  )

static

Sample resonance mass in a thermal medium.

Samples mass from the distribution

\[ dN/dm \sim A(m) m^2 K_2\left( \frac{m}{T}\right) \]

For stable particles always returns pole mass.

Parameters

[in]	ptype	the hadron sort, for which mass is sampled
[in]	beta	inverse temperature 1/T [1/GeV]

Returns

sampled mass

Definition at line 388 of file hadgas_eos.cc.

                                                      {
  if (ptype.is_stable()) {
    return ptype.mass();
  }
  // Sampling mass m from A(m) x^2 BesselK_2(x), where x = beta m.
  // Strategy employs the idea of importance sampling:
  // -- Sample mass from the simple Breit-Wigner first, then
  //    reject by the ratio.
  // -- Note that f(x) = x^2 BesselK_2(x) monotonously decreases
  //    and has maximum at x = xmin. That is why instead of f(x)
  //    the ratio f(x)/f(xmin) is used.
 
  const double max_mass = 5.0;  // GeV
  double m, q;
  {
    // Allow underflows in exponentials
    DisableFloatTraps guard(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
    const double w0 = ptype.width_at_pole();
    const double mth = ptype.min_mass_spectral();
    const double m0 = ptype.mass();
    double max_ratio = m0 * m0 * std::exp(-beta * m0) *
                       gsl_sf_bessel_Kn_scaled(2, m0 * beta) *
                       ptype.full_spectral_function(m0) /
                       ptype.breit_wigner_spectral_function(m0);
    // Heuristic adaptive maximum search to find max_ratio
    constexpr int npoints = 31;
    double m_lower = mth, m_upper = max_mass, m_where_max = m0;
 
    for (size_t n_iterations = 0; n_iterations < 2; n_iterations++) {
      const double dm = (m_upper - m_lower) / npoints;
      for (size_t i = 1; i < npoints; i++) {
        m = m_lower + dm * i;
        const double thermal_factor =
            m * m * std::exp(-beta * m) * gsl_sf_bessel_Kn_scaled(2, m * beta);
        q = ptype.full_spectral_function(m) * thermal_factor /
            ptype.breit_wigner_spectral_function(m);
        if (q > max_ratio) {
          max_ratio = q;
          m_where_max = m;
        }
      }
      m_lower = m_where_max - (m_where_max - m_lower) * 0.1;
      m_upper = m_where_max + (m_upper - m_where_max) * 0.1;
    }
    // Safety factor
    max_ratio *= 1.5;
 
    do {
      // sample mass from A(m)
      do {
        m = random::cauchy(m0, 0.5 * w0, mth, max_mass);
        const double thermal_factor =
            m * m * std::exp(-beta * m) * gsl_sf_bessel_Kn_scaled(2, m * beta);
        q = ptype.full_spectral_function(m) * thermal_factor /
            ptype.breit_wigner_spectral_function(m);
      } while (q < random::uniform(0., max_ratio));
      if (q > max_ratio) {
        logg[LResonances].warn(
            ptype.name(), " - maximum increased in",
            " sample_mass_thermal from ", max_ratio, " to ", q, ", mass = ", m,
            " previously assumed maximum at m = ", m_where_max);
        max_ratio = q;
        m_where_max = m;
      } else {
        break;
      }
    } while (true);
  }
  return m;
}

Here is the call graph for this function:

Here is the caller graph for this function:

solve_eos() [1/2]

std::array< double, 4 > smash::HadronGasEos::solve_eos	(	double	e,
		double	nb,
		double	ns,
		double	nq,
		std::array< double, 4 >	initial_approximation
	)

Compute temperature and chemical potentials given energy-, net baryon-, net strangeness- and net charge density and an inital approximation.

Parameters

[in]	e	energy density [GeV/fm \(^3\)]
[in]	nb	net baryon density [fm \(^{-3}\)]
[in]	ns	net strangeness density [fm \(^{-3}\)]
[in]	nq	net charge density [fm \(^{-3}\)]
[in]	initial_approximation	(T [GeV], mub [GeV], mus [GeV]) to use as starting point

Returns

array of 4 values: temperature, baryon chemical potential, strange chemical potential and charge chemical potential

Definition at line 590 of file hadgas_eos.cc.

                                               {
  int residual_status = GSL_SUCCESS;
  size_t iter = 0;
 
  struct rparams p = {e, nb, ns, nq, account_for_resonance_widths_};
  gsl_multiroot_function f = {&HadronGasEos::set_eos_solver_equations,
                              n_equations_, &p};
 
  gsl_vector_set(x_, 0, initial_approximation[0]);
  gsl_vector_set(x_, 1, initial_approximation[1]);
  gsl_vector_set(x_, 2, initial_approximation[2]);
  gsl_vector_set(x_, 3, initial_approximation[3]);
 
  gsl_multiroot_fsolver_set(solver_, &f, x_);
  do {
    iter++;
    const auto iterate_status = gsl_multiroot_fsolver_iterate(solver_);
    // std::cout << print_solver_state(iter);
 
    // Avoiding too low temperature
    if (gsl_vector_get(solver_->x, 0) < 0.015) {
      return {0.0, 0.0, 0.0, 0.0};
    }
 
    // check if solver is stuck
    if (iterate_status) {
      break;
    }
    residual_status = gsl_multiroot_test_residual(solver_->f, tolerance_);
  } while (residual_status == GSL_CONTINUE && iter < 1000);
 
  if (residual_status != GSL_SUCCESS) {
    std::stringstream solver_parameters;
    solver_parameters << "\nSolver run with "
                      << "e = " << e << ", nb = " << nb << ", ns = " << ns
                      << ", nq = " << nq
                      << ", init. approx.: " << initial_approximation[0] << " "
                      << initial_approximation[1] << " "
                      << initial_approximation[2] << " "
                      << initial_approximation[3] << std::endl;
    logg[LResonances].warn(gsl_strerror(residual_status) +
                           solver_parameters.str() + print_solver_state(iter));
  }
 
  return {gsl_vector_get(solver_->x, 0), gsl_vector_get(solver_->x, 1),
          gsl_vector_get(solver_->x, 2), gsl_vector_get(solver_->x, 3)};
}

Here is the call graph for this function:

Here is the caller graph for this function:

solve_eos() [2/2]

std::array<double, 4> smash::HadronGasEos::solve_eos ( double  e, double  nb, double  ns, double  nq  )

inline

Compute temperature and chemical potentials given energy-, net baryon-, net strangeness- and net charge density without an inital approximation.

Parameters

[in]	e	energy density [GeV/fm \(^3\)]
[in]	nb	net baryon density [fm \(^{-3}\)]
[in]	ns	net strangeness density [fm \(^{-3}\)]
[in]	nq	net charge density [fm \(^{-3}\)]

Returns

array of 4 values: temperature, baryon chemical potential and strange chemical potential and charge

Definition at line 321 of file hadgas_eos.h.

                                                                           {
    return solve_eos(e, nb, ns, nq, solve_eos_initial_approximation(e, nb, nq));
  }

Here is the call graph for this function:

solve_eos_initial_approximation()

std::array< double, 4 > smash::HadronGasEos::solve_eos_initial_approximation	(	double	e,
		double	nb,
		double	nq
	)

Compute a reasonable initial approximation for solve_eos.

Parameters

[in]	e	energy density [GeV/fm \(^3\)]
[in]	nb	net baryon density [fm \(^{-3}\)]
[in]	nq	net charge density [fm \(^{-3}\)]

Returns

array of 3 values: temperature, baryon chemical potential and strange chemical potential

Definition at line 510 of file hadgas_eos.cc.

                                                                               {
  assert(e >= 0.0);
  // 1. Get temperature from energy density assuming zero chemical potentials
  int degeneracies_sum = 0.0;
  for (const ParticleType &ptype : ParticleType::list_all()) {
    if (is_eos_particle(ptype)) {
      degeneracies_sum += ptype.spin() + 1;
    }
  }
  // Temperature in case of massless gas. For massive it should be larger.
  const double T_min = std::pow(e / prefactor_ / 6 / degeneracies_sum, 1. / 4.);
  // Simply assume that the temperature is not higher than 2 GeV.
  const double T_max = 2.0;
 
  struct eparams parameters = {e};
  gsl_function F = {&e_equation, &parameters};
  const gsl_root_fsolver_type *T = gsl_root_fsolver_brent;
  gsl_root_fsolver *e_solver;
  e_solver = gsl_root_fsolver_alloc(T);
  gsl_root_fsolver_set(e_solver, &F, T_min, T_max);
 
  int iter = 0, status, max_iter = 100;
  double T_init = 0.0;
 
  do {
    iter++;
    status = gsl_root_fsolver_iterate(e_solver);
    if (status != GSL_SUCCESS) {
      break;
    }
    T_init = gsl_root_fsolver_root(e_solver);
    double x_lo = gsl_root_fsolver_x_lower(e_solver);
    double x_hi = gsl_root_fsolver_x_upper(e_solver);
    status = gsl_root_test_interval(x_lo, x_hi, 0.0, 0.001);
  } while (status == GSL_CONTINUE && iter < max_iter);
 
  if (status != GSL_SUCCESS) {
    std::stringstream err_msg;
    err_msg << "Solver of equation for temperature with e = " << e
            << " failed to converge. Maybe Tmax = " << T_max << " is too small?"
            << std::endl;
    throw std::runtime_error(gsl_strerror(status) + err_msg.str());
  }
 
  gsl_root_fsolver_free(e_solver);
 
  // 2. Get the baryon chemical potential for muS = muQ = 0 with previously
  //    obtained T
  double n_only_baryons = 0.0;
  for (const ParticleType &ptype : ParticleType::list_all()) {
    if (is_eos_particle(ptype) && ptype.baryon_number() == 1) {
      n_only_baryons +=
          scaled_partial_density(ptype, 1.0 / T_init, 0.0, 0.0, 0.0);
    }
  }
  const double nb_scaled = nb / prefactor_ / (T_init * T_init * T_init);
  double mub_init = T_init * std::asinh(nb_scaled / n_only_baryons / 2.0);
 
  // 3. Get the charge chemical potential assuming muB = muS = 0 with previously
  //    obtained T
  double n_only_charge_1_particles = 0.0;
  for (const ParticleType &ptype : ParticleType::list_all()) {
    if (is_eos_particle(ptype) && ptype.charge() == 1) {
      n_only_charge_1_particles +=
          scaled_partial_density(ptype, 1.0 / T_init, 0.0, 0.0, 0.0);
    }
  }
  const double q_scaled = nq / prefactor_ / (T_init * T_init * T_init);
  double muq_init =
      T_init * std::asinh(q_scaled / n_only_charge_1_particles / 2.0);
 
  // 4. Get the strange chemical potential, where mus = 0 is typically a good
  //    initial approximation
  std::array<double, 4> initial_approximation = {T_init, mub_init, 0.0,
                                                 muq_init};
  return initial_approximation;
}

Here is the call graph for this function:

Here is the caller graph for this function:

mus_net_strangeness0()

double smash::HadronGasEos::mus_net_strangeness0 ( double  T, double  mub, double  muq  )

static

Compute strangeness chemical potential, requiring that net strangeness = 0.

Parameters

[in]	T	temperature [GeV]
[in]	mub	baryon chemical potential [GeV]
[in]	muq	charge chemical potential [GeV]

Returns

strangeness chemical potential [GeV]

Definition at line 460 of file hadgas_eos.cc.

                                                                          {
  // Binary search
  double mus_u = mub + T;
  double mus_l = 0.0;
  double mus, rhos;
  size_t iteration = 0;
  // 50 iterations should give precision 2^-50 ~ 10^-15
  const size_t max_iteration = 50;
  do {
    mus = 0.5 * (mus_u + mus_l);
    rhos = net_strange_density(T, mub, mus, muq);
    if (rhos > 0.0) {
      mus_u = mus;
    } else {
      mus_l = mus;
    }
    iteration++;
  } while (std::abs(rhos) > tolerance_ && iteration < max_iteration);
  if (iteration == max_iteration) {
    throw std::runtime_error("Solving rho_s = 0: too many iterations.");
  }
  return mus;
}

Here is the call graph for this function:

from_table()

void smash::HadronGasEos::from_table ( EosTable::table_element &  res, double  e, double  nb, double  nq  ) const

inline

Get the element of eos table.

Definition at line 348 of file hadgas_eos.h.

                                   {
    eos_table_.get(res, e, nb, nq);
  }

Here is the call graph for this function:

Here is the caller graph for this function:

is_eos_particle()

static bool smash::HadronGasEos::is_eos_particle ( const ParticleType &  ptype )

inlinestatic

Check if a particle belongs to the EoS.

Definition at line 354 of file hadgas_eos.h.

                                                         {
    return ptype.is_hadron() && !ptype.pdgcode().is_heavy_flavor();
  }

Here is the call graph for this function:

Here is the caller graph for this function:

is_tabulated()

bool smash::HadronGasEos::is_tabulated ( ) const

inline

Create an EoS table or not?

Definition at line 359 of file hadgas_eos.h.

{ return tabulate_; }

Here is the caller graph for this function:

account_for_resonance_widths()

bool smash::HadronGasEos::account_for_resonance_widths ( ) const

inline

If resonance spectral functions are taken into account.

Definition at line 362 of file hadgas_eos.h.

                                            {
    return account_for_resonance_widths_;
  }

Here is the caller graph for this function:

scaled_partial_density_auxiliary()

double smash::HadronGasEos::scaled_partial_density_auxiliary ( double  m_over_T, double  mu_over_T  )

staticprivate

Function used to avoid duplications in density calculations.

Parameters

[in]	m_over_T	mass to temperature ratio \( m/T \)
[in]	mu_over_T	chemical potential to temperature ratio \( \mu/T \)

Returns

calculated \( (m/T)^2 exp(\mu/T) K_2(m/T) \)

Definition at line 226 of file hadgas_eos.cc.

                                                                        {
  double x = mu_over_T - m_over_T;
  if (x < -500.0) {
    return 0.0;
  }
  x = std::exp(x);
  // In the case of small masses: K_n(z) -> (n-1)!/2 *(2/z)^n, z -> 0,
  // z*z*K_2(z) -> 2
  return (m_over_T < really_small)
             ? 2.0 * x
             : m_over_T * m_over_T * x * gsl_sf_bessel_Kn_scaled(2, m_over_T);
}

Here is the caller graph for this function:

scaled_partial_density()

double smash::HadronGasEos::scaled_partial_density ( const ParticleType &  ptype, double  beta, double  mub, double  mus, double  muq, bool  account_for_width = false  )

staticprivate

Compute (unnormalized) density of one hadron sort - helper functions used to reduce code duplication.

Parameters

[in]	ptype	the hadron sort, for which partial density is computed
[in]	beta	inverse temperature [1/GeV]
[in]	mub	baryon chemical potential [GeV]
[in]	mus	strangeness chemical potential [GeV]
[in]	muq	charge chemical potential [GeV]
[in]	account_for_width	Take hadron spectral functions into account or not. When taken into account, they result in a considerable slow down.

Returns

partial (unnormalized) density of the given hadron sort \(n\) [fm \(^{-3}\)]

Definition at line 240 of file hadgas_eos.cc.

                                                                    {
  const double m_over_T = ptype.mass() * beta;
  double mu_over_T = beta * (ptype.baryon_number() * mub +
                             ptype.strangeness() * mus + ptype.charge() * muq);
  const double g = ptype.spin() + 1;
  if (ptype.is_stable() || !account_for_width) {
    return g * scaled_partial_density_auxiliary(m_over_T, mu_over_T);
  } else {
    // Integral \int_{threshold}^{\infty} A(m) N_{thermal}(m) dm,
    // where A(m) is the spectral function of the resonance.
    const double m0 = ptype.mass();
    const double w0 = ptype.width_at_pole();
    const double mth = ptype.min_mass_spectral();
    const double u_min = std::atan(2.0 * (mth - m0) / w0);
    const double u_max = 0.5 * M_PI;
    Integrator integrate;
    const double result =
        g * integrate(u_min, u_max, [&](double u) {
          // One of many possible variable substitutions. Not clear if it has
          // any advantages, except transforming (m_th, inf) to finite interval.
          const double tanu = std::tan(u);
          const double m = m0 + 0.5 * w0 * tanu;
          const double jacobian = 0.5 * w0 * (1.0 + tanu * tanu);
          return ptype.full_spectral_function(m) * jacobian *
                 scaled_partial_density_auxiliary(m * beta, mu_over_T);
        });
    return result;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

set_eos_solver_equations()

int smash::HadronGasEos::set_eos_solver_equations ( const gsl_vector *  x, void *  params, gsl_vector *  f  )

staticprivate

Interface EoS equations to be solved to gnu library.

Definition at line 484 of file hadgas_eos.cc.

                                                          {
  double e = reinterpret_cast<struct rparams *>(params)->e;
  double nb = reinterpret_cast<struct rparams *>(params)->nb;
  double ns = reinterpret_cast<struct rparams *>(params)->ns;
  double nq = reinterpret_cast<struct rparams *>(params)->nq;
  bool w = reinterpret_cast<struct rparams *>(params)->account_for_width;
 
  const double T = gsl_vector_get(x, 0);
  const double mub = gsl_vector_get(x, 1);
  const double mus = gsl_vector_get(x, 2);
  const double muq = gsl_vector_get(x, 3);
 
  gsl_vector_set(f, 0, energy_density(T, mub, mus, muq) - e);
  gsl_vector_set(f, 1, net_baryon_density(T, mub, mus, muq, w) - nb);
  gsl_vector_set(f, 2, net_strange_density(T, mub, mus, muq, w) - ns);
  gsl_vector_set(f, 3, net_charge_density(T, mub, mus, muq, w) - nq);
 
  return GSL_SUCCESS;
}

Here is the call graph for this function:

Here is the caller graph for this function:

e_equation()

double smash::HadronGasEos::e_equation ( double  T, void *  params  )

staticprivate

See also

set_eos_solver_equations()

Definition at line 505 of file hadgas_eos.cc.

                                                      {
  const double edens = reinterpret_cast<struct eparams *>(params)->edens;
  return edens - energy_density(T, 0.0, 0.0, 0.0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

print_solver_state()

std::string smash::HadronGasEos::print_solver_state ( size_t  iter ) const

private

Helpful printout, useful for debugging if gnu equation solving goes crazy.

Parameters

[in]

iter

current value of iterator

Returns

debug output string with iter, x and f(x) from solver

Definition at line 640 of file hadgas_eos.cc.

                                                            {
  std::stringstream s;
  // clang-format off
  s << "iter = " << iter << ","
    << " x = " << gsl_vector_get(solver_->x, 0) << " "
               << gsl_vector_get(solver_->x, 1) << " "
               << gsl_vector_get(solver_->x, 2) << ", "
               << gsl_vector_get(solver_->x, 3) << ", "
    << "f(x) = " << gsl_vector_get(solver_->f, 0) << " "
                 << gsl_vector_get(solver_->f, 1) << " "
                 << gsl_vector_get(solver_->f, 2) << " "
                 << gsl_vector_get(solver_->f, 3) << std::endl;
  // clang-format on
  return s.str();
}

Here is the caller graph for this function:

Member Data Documentation

prefactor_

constexpr double smash::HadronGasEos::prefactor_

staticconstexprprivate

Initial value:

=
      0.5 * M_1_PI * M_1_PI / (hbarc * hbarc * hbarc)

Constant factor, that appears in front of many thermodyn. expressions.

Definition at line 430 of file hadgas_eos.h.

tolerance_

constexpr double smash::HadronGasEos::tolerance_ = 1.e-8

staticconstexprprivate

Precision of equation solving.

Definition at line 434 of file hadgas_eos.h.

n_equations_

constexpr size_t smash::HadronGasEos::n_equations_ = 4

staticconstexprprivate

Number of equations in the system of equations to be solved.

Definition at line 437 of file hadgas_eos.h.

eos_table_

EosTable smash::HadronGasEos::eos_table_ = EosTable(1.e-1, 1.e-1, 1.e-1, 90, 90, 90)

private

EOS Table to be used.

Definition at line 440 of file hadgas_eos.h.

x_

gsl_vector* smash::HadronGasEos::x_

private

Variables used by gnu equation solver.

They are stored here to allocate and deallocate memory for them only once. It is expected that this class will be used for solving the EoS many times, so multiple allocations and frees are unwanted.

Definition at line 448 of file hadgas_eos.h.

solver_

gsl_multiroot_fsolver* smash::HadronGasEos::solver_

private

See also

x_

Definition at line 451 of file hadgas_eos.h.

tabulate_

const bool smash::HadronGasEos::tabulate_

private

Create an EoS table or not?

Definition at line 454 of file hadgas_eos.h.

account_for_resonance_widths_

const bool smash::HadronGasEos::account_for_resonance_widths_

private

Use pole masses of resonances or integrate over spectral functions.

Definition at line 457 of file hadgas_eos.h.

---

The documentation for this class was generated from the following files:

• src/include/smash/hadgas_eos.h
• src/hadgas_eos.cc