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 into SMASH.

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

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

Definition at line 112 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, 3 >	solve_eos (double e, double nb, double ns, std::array< double, 3 > initial_approximation)
	Compute temperature and chemical potentials given energy-, net baryon-, net strangeness density and an inital approximation. More...
std::array< double, 3 >	solve_eos (double e, double nb, double ns)
	Compute temperature and chemical potentials given energy-, net baryon- and net strangeness density without an inital approximation. More...
std::array< double, 3 >	solve_eos_initial_approximation (double e, double nb)
	Compute a reasonable initial approximation for solve_eos. More...
void	from_table (EosTable::table_element &res, double e, double nb) 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)
	Compute energy density. More...
static double	density (double T, double mub, double mus, bool account_for_resonance_widths=false)
	Compute particle number density. More...
static double	pressure (double T, double mub, double mus, bool account_for_resonance_widths=false)
	Compute pressure \( p = n T \). More...
static double	net_baryon_density (double T, double mub, double mus, bool account_for_resonance_widths=false)
	Compute net baryon density. More...
static double	net_strange_density (double T, double mub, double mus, bool account_for_resonance_widths=false)
	Compute net strangeness density. More...
static double	partial_density (const ParticleType &ptype, double T, double mub, double mus, 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)
	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, 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-2, 1.e-2, 900, 900)
	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_ = 3
	Number of equations in the system of equations to be solved. More...

Constructor & Destructor Documentation

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 161 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_) {
    const auto &log = logger<LogArea::Resonances>();
    log.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:

smash::HadronGasEos::~HadronGasEos

(

)

Definition at line 183 of file hadgas_eos.cc.

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

Member Function Documentation

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

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}{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]

Returns

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

Definition at line 244 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() -
                       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:

double smash::HadronGasEos::density ( double  T, double  mub, double  mus, 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}{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]	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 272 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, account_for_width);
  }
  rho *= prefactor_ * T * T * T;
  return rho;
}

Here is the call graph for this function:

static double smash::HadronGasEos::pressure ( double  T, double  mub, double  mus, 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]	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 180 of file hadgas_eos.h.

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

Here is the call graph for this function:

Here is the caller graph for this function:

double smash::HadronGasEos::net_baryon_density ( double  T, double  mub, double  mus, 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}{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]	account_for_resonance_widths	if false, pole masses are used; if true, then integration over spectral function is included

Returns

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

Definition at line 289 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, 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:

double smash::HadronGasEos::net_strange_density ( double  T, double  mub, double  mus, 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}{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]	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 307 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, 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:

double smash::HadronGasEos::partial_density ( const ParticleType &  ptype, double  T, double  mub, double  mus, 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}{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]	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 234 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, account_for_width);
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 325 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.spectral_function(m0) / ptype.spectral_function_simple(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.spectral_function(m) * thermal_factor /
            ptype.spectral_function_simple(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.spectral_function(m) * thermal_factor /
            ptype.spectral_function_simple(m);
      } while (q < random::uniform(0., max_ratio));
      if (q > max_ratio) {
        const auto &log = logger<LogArea::Resonances>();
        log.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:

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

Compute temperature and chemical potentials given energy-, net baryon-, net strangeness 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]	initial_approximation	(T [GeV], mub [GeV], mus [GeV]) to use as starting point

Returns

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

Definition at line 507 of file hadgas_eos.cc.

                                               {
  int residual_status = GSL_SUCCESS;
  size_t iter = 0;

  struct rparams p = {e, nb, ns, 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_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};
    }

    // 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
                      << ", init. approx.: " << initial_approximation[0] << " "
                      << initial_approximation[1] << " "
                      << initial_approximation[2] << std::endl;
    const auto &log = logger<LogArea::HadronGasEos>();
    log.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)};
}

Here is the call graph for this function:

Here is the caller graph for this function:

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

inline

Compute temperature and chemical potentials given energy-, net baryon- and net strangeness 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}\)]

Returns

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

Definition at line 279 of file hadgas_eos.h.

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

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

Compute a reasonable initial approximation for solve_eos.

Parameters

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

Returns

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

Definition at line 444 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 = 0 and 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);
    }
  }
  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. mus = 0 is typically a good initial approximation

  std::array<double, 3> initial_approximation = {T_init, mub_init, 0.0};
  return initial_approximation;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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

static

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

Parameters

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

Returns

strangeness chemical potential [GeV]

Definition at line 397 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);
    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:

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

inline

Get the element of eos table.

Definition at line 303 of file hadgas_eos.h.

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

Here is the caller graph for this function:

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

inlinestatic

Check if a particle belongs to the EoS.

Definition at line 308 of file hadgas_eos.h.

                                                         {
    return ptype.is_hadron() && ptype.pdgcode().charmness() == 0;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

bool smash::HadronGasEos::is_tabulated ( ) const

inline

Create an EoS table or not?

Definition at line 313 of file hadgas_eos.h.

{ return tabulate_; }

Here is the caller graph for this function:

bool smash::HadronGasEos::account_for_resonance_widths ( ) const

inline

If resonance spectral functions are taken into account.

Definition at line 316 of file hadgas_eos.h.

                                            {
    return account_for_resonance_widths_;
  }

Here is the caller graph for this function:

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 188 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:

double smash::HadronGasEos::scaled_partial_density ( const ParticleType &  ptype, double  beta, double  mub, double  mus, 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]	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 202 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);
  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.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:

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 421 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;
  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);

  gsl_vector_set(f, 0, energy_density(T, mub, mus) - e);
  gsl_vector_set(f, 1, net_baryon_density(T, mub, mus, w) - nb);
  gsl_vector_set(f, 2, net_strange_density(T, mub, mus, w) - ns);

  return GSL_SUCCESS;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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

staticprivate

See also

set_eos_solver_equations()

Definition at line 439 of file hadgas_eos.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 555 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) << ", "
    << "f(x) = " << gsl_vector_get(solver_->f, 0) << " "
                 << gsl_vector_get(solver_->f, 1) << " "
                 << gsl_vector_get(solver_->f, 2) << std::endl;
  // clang-format on
  return s.str();
}

Here is the caller graph for this function:

Member Data Documentation

constexpr double smash::HadronGasEos::prefactor_

staticprivate

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 381 of file hadgas_eos.h.

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

staticprivate

Precision of equation solving.

Definition at line 385 of file hadgas_eos.h.

constexpr size_t smash::HadronGasEos::n_equations_ = 3

staticprivate

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

Definition at line 388 of file hadgas_eos.h.

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

private

EOS Table to be used.

Definition at line 391 of file hadgas_eos.h.

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 399 of file hadgas_eos.h.

gsl_multiroot_fsolver* smash::HadronGasEos::solver_

private

See also

x_

Definition at line 402 of file hadgas_eos.h.

const bool smash::HadronGasEos::tabulate_

private

Create an EoS table or not?

Definition at line 405 of file hadgas_eos.h.

const bool smash::HadronGasEos::account_for_resonance_widths_

private

Use pole masses of resonances or integrate over spectral functions.

Definition at line 408 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