smash::EosTable Class Reference

#include <hadgas_eos.h>

A class to hold, compute and access tabulated EoS.

Definition at line 32 of file hadgas_eos.h.

Collaboration diagram for smash::EosTable:

Classes
struct	table_element
	Define the data structure for one element of the table. More...

Public Member Functions
	EosTable (double de, double dnb, size_t n_e, size_t n_b)
	Sets up a table p/T/muB/mus versus (e, nb), where e is energy density, nb is net baryon density, p - pressure, T - temperature, muB - net baryon chemical potential, muS - net strangeness potential. More...
void	compile_table (HadronGasEos &eos, const std::string &eos_savefile_name="hadgas_eos.dat")
	Computes the actual content of the table (for EosTable description see documentation of the constructor). More...
void	get (table_element &res, double e, double nb) const
	Obtain interpolated p/T/muB/muS from the tabulated equation of state given energy density and net baryon density. More...

Private Member Functions
size_t	index (size_t ie, size_t inb) const
	proper index in a 1d vector, where the 2d table is stored More...

Private Attributes
std::vector< table_element >	table_
	Storage for the tabulated equation of state. More...
double	de_
	Step in energy density. More...
double	dnb_
	Step in net-baryon density. More...
size_t	n_e_
	Number of steps in energy density. More...
size_t	n_nb_
	Number of steos in net-baryon density. More...

Constructor & Destructor Documentation

smash::EosTable::EosTable	(	double	de,
		double	dnb,
		size_t	n_e,
		size_t	n_b
	)

Sets up a table p/T/muB/mus versus (e, nb), where e is energy density, nb is net baryon density, p - pressure, T - temperature, muB - net baryon chemical potential, muS - net strangeness potential.

Net strangeness density and isospin projection density are assumed to be 0 (Note that corresponding chemical potentials are still non-zero, because muB != 0).

After calling this constructor the table is allocated, but it is still empty. To compute values call compile_table.

Parameters

[in]	de	step in energy density [GeV/fm^4]
[in]	dnb	step in net baryon density [GeV/fm^3]
[in]	n_e	number of steps in energy density
[in]	n_b	number of steps in net baryon density

Entry at (ie, inb) corresponds to energy density and net baryon density (e, nb) = (ie*de, inb*dnb) [GeV/fm^4, GeV/fm^3].

Definition at line 29 of file hadgas_eos.cc.

    : de_(de), dnb_(dnb), n_e_(n_e), n_nb_(n_nb) {
  table_.resize(n_e_ * n_nb_);
}

Member Function Documentation

void smash::EosTable::compile_table	(	HadronGasEos &	eos,
		const std::string &	eos_savefile_name = "hadgas_eos.dat"
	)

Computes the actual content of the table (for EosTable description see documentation of the constructor).

Parameters

[in]	eos	equation of state
[in]	eos_savefile_name	name of the file to save tabulated equation of state

Definition at line 34 of file hadgas_eos.cc.

                                                                 {
  bool table_read_success = false, table_consistency = true;
  if (boost::filesystem::exists(eos_savefile_name)) {
    std::cout << "Reading table from file " << eos_savefile_name << std::endl;
    std::ifstream file;
    file.open(eos_savefile_name, std::ios::in);
    file >> de_ >> dnb_;
    file >> n_e_ >> n_nb_;
    table_.resize(n_e_ * n_nb_);
    for (size_t ie = 0; ie < n_e_; ie++) {
      for (size_t inb = 0; inb < n_nb_; inb++) {
        double p, T, mub, mus;
        file >> p >> T >> mub >> mus;
        table_[index(ie, inb)] = {p, T, mub, mus};
      }
    }
    table_read_success = true;
    std::cout << "Table consumed successfully." << std::endl;
  }

  if (table_read_success) {
    // Check if the saved table is consistent with the current particle table
    std::cout << "Checking consistency of the table... " << std::endl;
    constexpr size_t number_of_steps = 50;
    const size_t ie_step = 1 + n_e_ / number_of_steps;
    const size_t inb_step = 1 + n_nb_ / number_of_steps;
    for (size_t ie = 0; ie < n_e_; ie += ie_step) {
      for (size_t inb = 0; inb < n_nb_; inb += inb_step) {
        const table_element x = table_[index(ie, inb)];
        const double e_comp = eos.energy_density(x.T, x.mub, x.mus);
        const double nb_comp = eos.net_baryon_density(x.T, x.mub, x.mus);
        const double ns_comp = eos.net_strange_density(x.T, x.mub, x.mus);
        const double p_comp = eos.pressure(x.T, x.mub, x.mus);
        // Precision is just 10^-3, this is precision of saved data in the file
        const double eps = 1.e-3;
        // Only check the physical region, hence T > 0 condition
        if ((std::abs(de_ * ie - e_comp) > eps ||
             std::abs(dnb_ * inb - nb_comp) > eps || std::abs(ns_comp) > eps ||
             std::abs(x.p - p_comp) > eps) &&
            (x.T > 0.0)) {
          std::cout << "discrepancy: " << de_ * ie << " = " << e_comp << ", "
                    << dnb_ * inb << " = " << nb_comp << ", " << x.p << " = "
                    << p_comp << ", 0 = " << ns_comp << std::endl;
          table_consistency = false;
          goto finish_consistency_check;
        }
      }
    }
  }
finish_consistency_check:

  if (!table_read_success || !table_consistency) {
    std::cout << "Compiling an EoS table..." << std::endl;
    const double ns = 0.0;
    for (size_t ie = 0; ie < n_e_; ie++) {
      const double e = de_ * ie;
      for (size_t inb = 0; inb < n_nb_; inb++) {
        const double nb = dnb_ * inb;
        // It is physically impossible to have energy density > nucleon mass*nb,
        // therefore eqns have no solutions.
        if (nb >= e) {
          table_[index(ie, inb)] = {0.0, 0.0, 0.0, 0.0};
          continue;
        }
        // Take extrapolated (T, mub, mus) as initial approximation
        std::array<double, 3> init_approx;
        if (inb >= 2) {
          const table_element y = table_[index(ie, inb - 2)];
          const table_element x = table_[index(ie, inb - 1)];
          init_approx = {2.0 * x.T - y.T, 2.0 * x.mub - y.mub,
                         2.0 * x.mus - y.mus};
        } else {
          init_approx = eos.solve_eos_initial_approximation(e, nb);
        }
        const std::array<double, 3> res = eos.solve_eos(e, nb, ns, init_approx);
        const double T = res[0];
        const double mub = res[1];
        const double mus = res[2];
        table_[index(ie, inb)] = {eos.pressure(T, mub, mus), T, mub, mus};
      }
    }
    // Save table to file
    std::cout << "Saving table to file " << eos_savefile_name << std::endl;
    std::ofstream file;
    file.open(eos_savefile_name, std::ios::out);
    file << de_ << " " << dnb_ << std::endl;
    file << n_e_ << " " << n_nb_ << std::endl;
    file << std::setprecision(7);
    file << std::fixed;
    for (size_t ie = 0; ie < n_e_; ie++) {
      for (size_t inb = 0; inb < n_nb_; inb++) {
        const EosTable::table_element x = table_[index(ie, inb)];
        file << x.p << " " << x.T << " " << x.mub << " " << x.mus << std::endl;
      }
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

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

Obtain interpolated p/T/muB/muS from the tabulated equation of state given energy density and net baryon density.

Parameters

[in]	e	energy density
[in]	nb	net baryon density
[out]	res	structure, that contains p/T/muB/muS

Definition at line 133 of file hadgas_eos.cc.

                                                                        {
  const size_t ie = static_cast<size_t>(std::floor(e / de_));
  const size_t inb = static_cast<size_t>(std::floor(nb / dnb_));

  if (ie >= n_e_ - 1 || inb >= n_nb_ - 1) {
    res = {-1.0, -1.0, -1.0, -1.0};
  } else {
    // 1st order interpolation
    const double ae = e / de_ - ie;
    const double an = nb / dnb_ - inb;
    const EosTable::table_element s1 = table_[index(ie, inb)];
    const EosTable::table_element s2 = table_[index(ie + 1, inb)];
    const EosTable::table_element s3 = table_[index(ie, inb + 1)];
    const EosTable::table_element s4 = table_[index(ie + 1, inb + 1)];
    res.p = ae * (an * s4.p + (1.0 - an) * s2.p) +
            (1.0 - ae) * (an * s3.p + (1.0 - an) * s1.p);
    res.T = ae * (an * s4.T + (1.0 - an) * s2.T) +
            (1.0 - ae) * (an * s3.T + (1.0 - an) * s1.T);
    res.mub = ae * (an * s4.mub + (1.0 - an) * s2.mub) +
              (1.0 - ae) * (an * s3.mub + (1.0 - an) * s1.mub);
    res.mus = ae * (an * s4.mus + (1.0 - an) * s2.mus) +
              (1.0 - ae) * (an * s3.mus + (1.0 - an) * s1.mus);
  }
}

Here is the call graph for this function:

size_t smash::EosTable::index ( size_t  ie, size_t  inb  ) const

inlineprivate

proper index in a 1d vector, where the 2d table is stored

Definition at line 87 of file hadgas_eos.h.

{ return ie * n_nb_ + inb; }

Here is the caller graph for this function:

Member Data Documentation

std::vector<table_element> smash::EosTable::table_

private

Storage for the tabulated equation of state.

Definition at line 89 of file hadgas_eos.h.

double smash::EosTable::de_

private

Step in energy density.

Definition at line 91 of file hadgas_eos.h.

double smash::EosTable::dnb_

private

Step in net-baryon density.

Definition at line 93 of file hadgas_eos.h.

size_t smash::EosTable::n_e_

private

Number of steps in energy density.

Definition at line 95 of file hadgas_eos.h.

size_t smash::EosTable::n_nb_

private

Number of steos in net-baryon density.

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