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.

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

Public Member Functions
	EosTable (double de, double dnb, double dq, size_t n_e, size_t n_b, size_t n_q)
	Sets up a table p/T/muB/mus/muQ versus (e, nb, nq), where e - energy density, nb - net baryon density, nq - net charge density, p - pressure, T. 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, double nq) const
	Obtain interpolated p/T/muB/muS/muQ from the tabulated equation of state given energy density, net baryon density and net charge density. More...

Private Member Functions
size_t	index (size_t ie, size_t inb, size_t inq) const
	proper index in a 1d vector, where the 3d 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...
double	dq_
	Step in net-charge density. More...
size_t	n_e_
	Number of steps in energy density. More...
size_t	n_nb_
	Number of steps in net-baryon density. More...
size_t	n_q_
	Number of steps in net-charge density. More...

Constructor & Destructor Documentation

EosTable()

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

Sets up a table p/T/muB/mus/muQ versus (e, nb, nq), where e - energy density, nb - net baryon density, nq - net charge density, p - pressure, T.

• temperature, muB - net baryon chemical potential, muS - net strangeness potential, muQ - net charge chemical potential. Net strangeness density and isospin projection density are assumed to be 0 (Note that the corresponding chemical potential is 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]	dq	step in net charge density [GeV/Gev^3]
[in]	n_e	number of steps in energy density
[in]	n_b	number of steps in net baryon density
[in]	n_q	number of steps in net charge density

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

Definition at line 28 of file hadgas_eos.cc.

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

Member Function Documentation

compile_table()

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 (std::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_ >> dq_;
    file >> n_e_ >> n_nb_ >> n_q_;
    table_.resize(n_e_ * n_nb_ * n_q_);
    for (size_t ie = 0; ie < n_e_; ie++) {
      for (size_t inb = 0; inb < n_nb_; inb++) {
        for (size_t iq = 0; iq < n_q_; iq++) {
          double p, T, mub, mus, muq;
          file >> p >> T >> mub >> mus >> muq;
          table_[index(ie, inb, iq)] = {p, T, mub, mus, muq};
        }
      }
    }
    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;
    const size_t iq_step = 1 + n_q_ / 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) {
        for (size_t iq = 0; iq < n_q_; iq += iq_step) {
          const table_element x = table_[index(ie, inb, iq)];
          const bool w = eos.account_for_resonance_widths();
          const double e_comp = eos.energy_density(x.T, x.mub, x.mus, x.muq);
          const double nb_comp =
              eos.net_baryon_density(x.T, x.mub, x.mus, x.muq, w);
          const double ns_comp =
              eos.net_strange_density(x.T, x.mub, x.mus, x.muq, w);
          const double p_comp = eos.pressure(x.T, x.mub, x.mus, x.muq, w);
          const double nq_comp =
              eos.net_charge_density(x.T, x.mub, x.mus, x.muq, w);
          // 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 ||
               std::abs(dq_ * iq - nq_comp) > eps) &&
              (x.T > 0.0)) {
            std::cout << "discrepancy: " << de_ * ie << " = " << e_comp << ", "
                      << dnb_ * inb << " = " << nb_comp << ", " << x.p << " = "
                      << p_comp << ", 0 = " << ns_comp << ", " << dq_ * iq
                      << " = " << nq_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++) {
      std::cout << ie << "/" << n_e_ << "\r" << std::flush;
      const double e = de_ * ie;
      for (size_t inb = 0; inb < n_nb_; inb++) {
        const double nb = dnb_ * inb;
        for (size_t iq = 0; iq < n_q_; iq++) {
          const double q = dq_ * iq;
          // It is physically impossible to have energy density > nucleon
          // mass*nb, therefore eqns have no solutions.
          if (nb >= e || q >= e) {
            table_[index(ie, inb, iq)] = {0.0, 0.0, 0.0, 0.0, 0.0};
            continue;
          }
          // Take extrapolated (T, mub, mus, muq) as initial approximation
          std::array<double, 4> init_approx;
          if (inb >= 2) {
            const table_element y = table_[index(ie, inb - 2, iq)];
            const table_element x = table_[index(ie, inb - 1, iq)];
            init_approx = {2.0 * x.T - y.T, 2.0 * x.mub - y.mub,
                           2.0 * x.mus - y.mus, 2.0 * x.muq - y.muq};
          } else if (iq >= 2) {
            const table_element y = table_[index(ie, inb, iq - 2)];
            const table_element x = table_[index(ie, inb, iq - 1)];
            init_approx = {2.0 * x.T - y.T, 2.0 * x.mub - y.mub,
                           2.0 * x.mus - y.mus, 2.0 * x.muq - y.muq};
          } else {
            init_approx = eos.solve_eos_initial_approximation(e, nb, q);
          }
          const std::array<double, 4> res =
              eos.solve_eos(e, nb, ns, q, init_approx);
          const double T = res[0];
          const double mub = res[1];
          const double mus = res[2];
          const double muq = res[3];
          const bool w = eos.account_for_resonance_widths();
          table_[index(ie, inb, iq)] = {eos.pressure(T, mub, mus, muq, w), T,
                                        mub, mus, muq};
        }
      }
    }
    // 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_ << " " << dq_ << std::endl;
    file << n_e_ << " " << n_nb_ << " " << n_q_ << 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++) {
        for (size_t iq = 0; iq < n_q_; iq++) {
          const EosTable::table_element x = table_[index(ie, inb, iq)];
          file << x.p << " " << x.T << " " << x.mub << " " << x.mus << " "
               << x.muq << std::endl;
        }
      }
    }
  }
}

get()

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

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

Parameters

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

Definition at line 162 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_));
  const size_t iq = static_cast<size_t>(std::floor(q / dq_));
 
  if (ie >= n_e_ - 1 || inb >= n_nb_ - 1 || iq >= n_q_ - 1) {
    res = {-1.0, -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 double aq = q / dq_ - iq;
    const EosTable::table_element s1 = table_[index(ie, inb, iq)];
    const EosTable::table_element s2 = table_[index(ie + 1, inb, iq)];
    const EosTable::table_element s3 = table_[index(ie, inb + 1, iq)];
    const EosTable::table_element s4 = table_[index(ie + 1, inb + 1, iq)];
    const EosTable::table_element s5 = table_[index(ie, inb, iq + 1)];
    const EosTable::table_element s6 = table_[index(ie + 1, inb, iq + 1)];
    const EosTable::table_element s7 = table_[index(ie, inb + 1, iq + 1)];
    const EosTable::table_element s8 = table_[index(ie + 1, inb + 1, iq + 1)];
 
    res.p = interpolate_trilinear(ae, an, aq, s1.p, s2.p, s3.p, s4.p, s5.p,
                                  s6.p, s7.p, s8.p);
    res.T = interpolate_trilinear(ae, an, aq, s1.T, s2.T, s3.T, s4.T, s5.T,
                                  s6.T, s7.T, s8.T);
    res.mub = interpolate_trilinear(ae, an, aq, s1.mub, s2.mub, s3.mub, s4.mub,
                                    s5.mub, s6.mub, s7.mub, s8.mub);
    res.mus = interpolate_trilinear(ae, an, aq, s1.mus, s2.mus, s3.mus, s4.mus,
                                    s5.mus, s6.mus, s7.mus, s8.mus);
    res.muq = interpolate_trilinear(ae, an, aq, s1.muq, s2.muq, s3.muq, s4.muq,
                                    s5.muq, s6.muq, s7.muq, s8.muq);
  }
}

index()

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

inlineprivate

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

Definition at line 93 of file hadgas_eos.h.

                                                        {
    return n_q_ * (ie * n_nb_ + inb) + inq;
  }

Member Data Documentation

table_

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

private

Storage for the tabulated equation of state.

Definition at line 97 of file hadgas_eos.h.

de_

double smash::EosTable::de_

private

Step in energy density.

Definition at line 99 of file hadgas_eos.h.

dnb_

double smash::EosTable::dnb_

private

Step in net-baryon density.

Definition at line 101 of file hadgas_eos.h.

dq_

double smash::EosTable::dq_

private

Step in net-charge density.

Definition at line 103 of file hadgas_eos.h.

n_e_

size_t smash::EosTable::n_e_

private

Number of steps in energy density.

Definition at line 105 of file hadgas_eos.h.

n_nb_

size_t smash::EosTable::n_nb_

private

Number of steps in net-baryon density.

Definition at line 107 of file hadgas_eos.h.

n_q_

size_t smash::EosTable::n_q_

private

Number of steps in net-charge density.

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