#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 (const bool tabulate=false)
	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...

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)
	Compute particle number density. More...

static double	pressure (double T, double mub, double mus)
	Compute pressure \( p = n T \). More...

static double	net_baryon_density (double T, double mub, double mus)
	Compute net baryon density. More...

static double	net_strange_density (double T, double mub, double mus)
	Compute net strangeness density. More...

static double	partial_density (const ParticleType &ptype, double T, double mub, double mus)
	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)
	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...

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

◆ HadronGasEos()

smash::HadronGasEos::HadronGasEos ( const bool tabulate = false )

explicit

Constructor of HadronGasEos.

Definition at line 158 of file hadgas_eos.cc.

     : x_(gsl_vector_alloc(n_equations_)), tabulate_(tabulate) {
   const gsl_multiroot_fsolver_type *solver_type;
   solver_type = gsl_multiroot_fsolver_hybrid;
   solver_ = gsl_multiroot_fsolver_alloc(solver_type, n_equations_);
   if (tabulate_) {
     eos_table_.compile_table(*this);
   }
 }

Here is the call graph for this function:

◆ ~HadronGasEos()

smash::HadronGasEos::~HadronGasEos ( )

Definition at line 168 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
	)

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

◆ density()

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

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]

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

Definition at line 256 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);
   }
   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
	)

inlinestatic

Compute pressure \( p = n T \).

Parameters

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

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

Definition at line 159 of file hadgas_eos.h.

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

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
	)

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]

Returns: net baryon density density \(n_B\) [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 (!ptype.is_baryon() || !is_eos_particle(ptype)) {
       continue;
     }
     rho +=
         scaled_partial_density(ptype, beta, mub, mus) * 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
	)

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]

Returns: net strangeness density density \(n_S\) [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.strangeness() == 0 || !is_eos_particle(ptype)) {
       continue;
     }
     rho += scaled_partial_density(ptype, beta, mub, mus) * 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:

◆ partial_density()

double smash::HadronGasEos::partial_density	(	const ParticleType &	ptype,
		double	T,
		double	mub,
		double	mus
	)

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]

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

Definition at line 219 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);
 }

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

◆ solve_eos() [1/2]

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 486 of file hadgas_eos.cc.

                                                {
   int residual_status = GSL_SUCCESS;
   size_t iter = 0;
 
   struct rparams p = {e, nb, ns};
   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 << "Solver run with "
                       << "e = " << e << ", nb = " << nb << ", ns = " << ns
                       << ", init. approx.: " << initial_approximation[0] << " "
                       << initial_approximation[1] << " "
                       << initial_approximation[2] << std::endl;
     throw std::runtime_error(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:

◆ solve_eos() [2/2]

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

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

Here is the call graph for this function:

◆ solve_eos_initial_approximation()

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

◆ mus_net_strangeness0()

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

◆ from_table()

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

inline

Get the element of eos table.

Definition at line 272 of file hadgas_eos.h.

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

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

◆ is_tabulated()

bool smash::HadronGasEos::is_tabulated ( ) const

inline

Create an EoS table or not?

Definition at line 282 of file hadgas_eos.h.

282 { return tabulate_; }

smash::HadronGasEos::tabulate_

const bool tabulate_

Create an EoS table or not?

Definition: hadgas_eos.h:363

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 173 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
	)

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]

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

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

◆ 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 401 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;
 
   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) - nb);
   gsl_vector_set(f, 2, net_strange_density(T, mub, mus) - ns);
 
   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 418 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:

◆ 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 534 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

◆ prefactor_

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

◆ tolerance_

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

staticprivate

Precision of equation solving.

Definition at line 343 of file hadgas_eos.h.

◆ n_equations_

constexpr size_t smash::HadronGasEos::n_equations_ = 3

staticprivate

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

Definition at line 346 of file hadgas_eos.h.

◆ eos_table_

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

private

EOS Table to be used.

Definition at line 349 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 357 of file hadgas_eos.h.

◆ solver_

gsl_multiroot_fsolver* smash::HadronGasEos::solver_

private

See also: x_

Definition at line 360 of file hadgas_eos.h.

◆ tabulate_

const bool smash::HadronGasEos::tabulate_

private

Create an EoS table or not?

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

Classes

Public Member Functions

Static Public Member Functions

Private Member Functions

Static Private Member Functions

Private Attributes

Static Private Attributes

Constructor & Destructor Documentation

◆ HadronGasEos()

◆ ~HadronGasEos()

Member Function Documentation

◆ energy_density()

◆ density()

◆ pressure()

◆ net_baryon_density()

◆ net_strange_density()

◆ partial_density()

◆ sample_mass_thermal()

◆ solve_eos() [1/2]

◆ solve_eos() [2/2]

◆ solve_eos_initial_approximation()

◆ mus_net_strangeness0()

◆ from_table()

◆ is_eos_particle()

◆ is_tabulated()

◆ scaled_partial_density_auxiliary()

◆ scaled_partial_density()

◆ set_eos_solver_equations()

◆ e_equation()

◆ print_solver_state()

Member Data Documentation

◆ prefactor_

◆ tolerance_

◆ n_equations_

◆ eos_table_

◆ x_

◆ solver_

◆ tabulate_