42 const int max_iter = 50;
44 double e_previous_step = 0.0;
45 const double tolerance = 5.e-4;
47 for (iter = 0; iter < max_iter; iter++) {
50 if (std::abs(
e_ - e_previous_step) < tolerance) {
53 const double gamma_inv = std::sqrt(1.0 -
v_.
sqr());
70 if (iter == max_iter) {
71 std::cout <<
"Warning from solver: max iterations exceeded." 72 <<
" Accuracy: " << std::abs(
e_ - e_previous_step)
73 <<
" is less than tolerance " << tolerance << std::endl;
91 return out <<
"T[mu,0]: " << node.
Tmu0() <<
", nb: " << node.
nb()
92 <<
", ns: " << node.
ns() <<
", v: " << node.
v()
93 <<
", e: " << node.
e() <<
", p: " << node.
p()
94 <<
", T: " << node.
T() <<
", mub: " << node.
mub()
95 <<
", mus: " << node.
mus();
99 const std::array<int, 3> n_cells,
100 const std::array<double, 3> origin,
101 bool periodicity,
double e_critical,
102 double t_start,
double delta_t,
104 : eos_typelist_(list_eos_particles()),
105 N_sorts_(eos_typelist_.size()),
111 lat_ = make_unique<RectangularLattice<ThermLatticeNode>>(
112 lat_sizes, n_cells, origin, periodicity, upd);
113 const std::array<double, 3> abc =
lat_->cell_sizes();
122 bool ignore_cells_under_treshold) {
126 for (
auto &node : *
lat_) {
130 if (!ignore_cells_under_treshold ||
131 node.Tmu0().x0() + std::abs(node.Tmu0().x1()) +
132 std::abs(node.Tmu0().x2()) + std::abs(node.Tmu0().x3()) >=
134 node.compute_rest_frame_quantities(
eos_);
143 +0.5 *
lat_->cell_sizes()[0]),
145 +0.5 *
lat_->cell_sizes()[1]),
147 +0.5 *
lat_->cell_sizes()[2]));
151 ParticleList &plist,
const FourVector required_total_momentum) {
152 const auto &log = logger<LogArea::GrandcanThermalizer>();
156 log.info(
"Required 4-momentum: ", required_total_momentum);
157 log.info(
"Sampled 4-momentum: ", conserved.momentum());
159 (required_total_momentum.
threevec() - conserved.momentum().threevec()) /
161 log.info(
"Adjusting momenta by ", mom_to_add);
162 for (
auto &particle : plist) {
163 particle.set_4momentum(particle.type().mass(),
164 particle.momentum().threevec() + mom_to_add);
169 const ThreeVector beta_CM_generated = conserved.momentum().velocity();
173 double E_expected = required_total_momentum.
abs();
174 for (
auto &particle : plist) {
175 particle.boost_momentum(beta_CM_generated);
176 E += particle.momentum().x0();
180 double a, a_min, a_max, er;
181 const int max_iter = 50;
183 if (E_expected >= E) {
191 a = 0.5 * (a_min + a_max);
193 for (
const auto &particle : plist) {
194 const double p2 = particle.momentum().threevec().sqr();
195 const double E2 = particle.momentum().x0() * particle.momentum().x0();
196 E += std::sqrt(E2 + a * (a + 2.0) * p2);
204 log.debug(
"Iteration ", iter,
": a = ", a,
", Δ = ", er);
206 }
while (std::abs(er) > tolerance && iter < max_iter);
208 log.info(
"Renormalizing momenta by factor 1+a, a = ", a);
209 for (
auto &particle : plist) {
210 particle.set_4momentum(particle.type().mass(),
211 (1 + a) * particle.momentum().threevec());
212 particle.boost_momentum(-beta_CM_required);
223 for (
size_t i_type = 0; (i_type < N_sorts_) && (N_to_sample > 0); i_type++) {
224 if (
get_class(i_type) != particle_class) {
248 const double gamma = 1.0 / std::sqrt(1.0 - cell.
v().
sqr());
249 const double N_this_cell =
257 for (
int i = 0; i <
mult_int_[type_index]; i++) {
260 double partial_sum = 0.0;
261 int index_only_thermalized = -1;
262 while (partial_sum < r) {
263 index_only_thermalized++;
264 partial_sum +=
N_in_cells_[index_only_thermalized];
279 particle.
set_4momentum(m, phitheta.threevec() * momentum_radial);
283 plist.push_back(particle);
288 double time,
int ntest) {
289 const auto &log = logger<LogArea::GrandcanThermalizer>();
294 const double gamma = 1.0 / std::sqrt(1.0 - cell.
v().
sqr());
295 for (
size_t i = 0; i <
N_sorts_; i++) {
304 for (
size_t i = 0; i <
N_sorts_; i++) {
315 const auto Nbar_antibar = bessel_sampler_B.sample();
322 for (
size_t i = 0; i <
N_sorts_; i++) {
326 std::pair<int, int> NS_antiS;
332 NS_antiS = bessel_sampler_S.
sample();
334 NS_antiS = std::make_pair(
337 if (NS_antiS.first - NS_antiS.second !=
347 for (
size_t i = 0; i <
N_sorts_; i++) {
351 std::pair<int, int> NC_antiC;
356 conserved_initial.
charge() - ch_sampled);
357 NC_antiC = bessel_sampler_C.
sample();
359 NC_antiC = std::make_pair(
362 if (NC_antiC.first - NC_antiC.second !=
363 conserved_initial.
charge() - ch_sampled) {
374 for (
size_t itype = 0; itype <
N_sorts_; itype++) {
378 const double e_init = conserved_initial.
momentum().
x0();
381 e_tot += particle.momentum().x0();
383 if (std::abs(e_tot - e_init) > 0.01 * e_init) {
384 log.debug(
"Rejecting: energy ", e_tot,
" too far from e_init = ", e_init);
394 int S_plus = 0, S_minus = 0, B_plus = 0, B_minus = 0, E_plus = 0, E_minus = 0;
396 auto condition1 = [](int, int, int) {
return true; };
398 while (conserved_initial.
momentum().
x0() > energy ||
401 energy +=
p.momentum().x0();
402 if (
p.pdgcode().strangeness() > 0) {
404 S_plus +=
p.pdgcode().strangeness();
409 auto condition2 = [](
int S, int, int) {
return (S < 0); };
411 while (S_plus + S_minus > conserved_initial.
strangeness()) {
413 const int s_part =
p.pdgcode().strangeness();
415 if (S_plus + S_minus + s_part >= conserved_initial.
strangeness()) {
422 auto condition3 = [](
int S, int, int) {
return (S == 0); };
427 while (conserved_remaining.
momentum().
x0() > energy ||
430 energy +=
p.momentum().x0();
431 if (
p.pdgcode().baryon_number() > 0) {
433 B_plus +=
p.pdgcode().baryon_number();
438 auto condition4 = [](
int S,
int B, int) {
return (S == 0) && (B < 0); };
440 while (B_plus + B_minus > conserved_remaining.
baryon_number()) {
442 const int bar =
p.pdgcode().baryon_number();
443 if (B_plus + B_minus + bar >= conserved_remaining.
baryon_number()) {
450 auto condition5 = [](
int S,
int B, int) {
return (S == 0) && (B == 0); };
454 while (conserved_remaining.
momentum().
x0() > energy ||
455 E_plus < conserved_remaining.
charge()) {
457 energy +=
p.momentum().x0();
458 if (
p.pdgcode().charge() > 0) {
460 E_plus +=
p.pdgcode().charge();
465 auto condition6 = [](
int S,
int B,
int C) {
466 return (S == 0) && (B == 0) && (C < 0);
469 while (E_plus + E_minus > conserved_remaining.
charge()) {
471 const int charge =
p.pdgcode().charge();
472 if (E_plus + E_minus + charge >= conserved_remaining.
charge()) {
479 auto condition7 = [](
int S,
int B,
int C) {
480 return (S == 0) && (B == 0) && (C == 0);
485 while (conserved_remaining.
momentum().
x0() > energy) {
488 energy +=
p.momentum().x0();
494 const auto &log = logger<LogArea::GrandcanThermalizer>();
495 log.info(
"Starting forced thermalization, time ", time,
" fm/c");
502 for (
auto &particle : particles) {
503 const bool is_on_lattice =
504 lat_->value_at(particle.position().threevec(), node);
505 if (is_on_lattice && node.
e() >
e_crit_) {
519 log.info(
"Removed ",
to_remove_.size(),
" particles.");
527 const size_t lattice_total_cells =
lat_->size();
528 for (
size_t i = 0; i < lattice_total_cells; i++) {
533 log.info(
"Number of cells in the thermalization region = ",
547 throw std::invalid_argument(
548 "This thermalization algorithm is" 549 " not yet implemented");
565 struct to_average on_lattice = {0.0, 0.0, 0.0, 0.0, 0.0};
566 struct to_average in_therm_reg = {0.0, 0.0, 0.0, 0.0, 0.0};
567 double e_sum_on_lattice = 0.0, e_sum_in_therm_reg = 0.0;
568 int node_counter = 0;
569 for (
const auto &node : *
lat_) {
570 const double e = node.e();
571 on_lattice.T += node.T() * e;
572 on_lattice.mub += node.mub() * e;
573 on_lattice.mus += node.mus() * e;
574 on_lattice.nb += node.nb() * e;
575 on_lattice.ns += node.ns() * e;
576 e_sum_on_lattice += e;
578 in_therm_reg.T += node.T() * e;
579 in_therm_reg.mub += node.mub() * e;
580 in_therm_reg.mus += node.mus() * e;
581 in_therm_reg.nb += node.nb() * e;
582 in_therm_reg.ns += node.ns() * e;
583 e_sum_in_therm_reg += e;
588 on_lattice.T /= e_sum_on_lattice;
589 on_lattice.mub /= e_sum_on_lattice;
590 on_lattice.mus /= e_sum_on_lattice;
591 on_lattice.nb /= e_sum_on_lattice;
592 on_lattice.ns /= e_sum_on_lattice;
595 in_therm_reg.T /= e_sum_in_therm_reg;
596 in_therm_reg.mub /= e_sum_in_therm_reg;
597 in_therm_reg.mus /= e_sum_in_therm_reg;
598 in_therm_reg.nb /= e_sum_in_therm_reg;
599 in_therm_reg.ns /= e_sum_in_therm_reg;
602 std::cout <<
"Current time [fm/c]: " << clock.
current_time() << std::endl;
603 std::cout <<
"Averages on the lattice - T[GeV], mub[GeV], mus[GeV], " 604 <<
"nb[fm^-3], ns[fm^-3]: " << on_lattice.T <<
" " << on_lattice.mub
605 <<
" " << on_lattice.mus <<
" " << on_lattice.nb <<
" " 606 << on_lattice.ns << std::endl;
607 std::cout <<
"Averages in therm. region - T[GeV], mub[GeV], mus[GeV], " 608 <<
"nb[fm^-3], ns[fm^-3]: " << in_therm_reg.T <<
" " 609 << in_therm_reg.mub <<
" " << in_therm_reg.mus <<
" " 610 << in_therm_reg.nb <<
" " << in_therm_reg.ns << std::endl;
611 std::cout <<
"Volume with e > e_crit [fm^3]: " <<
cell_volume_ * node_counter
double p() const
Get pressure in the rest frame.
const double e_crit_
Critical energy density above which cells are thermalized.
void add_particle(const ParticleData &p, double factor)
Add particle contribution to Tmu0, nb and ns May look like unused at first glance, but it is actually used by update_lattice, where the node type of the lattice is templated.
A class to pre-calculate and store parameters relevant for density calculation.
ThermalizationAlgorithm
Defines the algorithm used for the forced thermalization.
GrandCanThermalizer(const std::array< double, 3 > lat_sizes, const std::array< int, 3 > n_cells, const std::array< double, 3 > origin, bool periodicity, double e_critical, double t_start, double delta_t, ThermalizationAlgorithm algo)
Default constructor for the GranCanThermalizer to allocate the lattice.
const size_t N_sorts_
Number of different species to be sampled.
The ThreeVector class represents a physical three-vector with the components .
void sample_multinomial(HadronClass particle_class, int N)
The sample_multinomial function samples integer numbers n_i distributed according to the multinomial ...
constexpr double really_small
Numerical error tolerance.
ThreeVector velocity() const
Get the velocity (3-vector divided by zero component).
ThermLatticeNode()
Default constructor of thermal quantities on the lattice returning thermodynamic quantities in comput...
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 a...
double mub() const
Get the net baryon chemical potential.
static double net_strange_density(double T, double mub, double mus)
Compute net strangeness density.
Class to handle the equation of state (EoS) of the hadron gas, consisting of all hadrons included int...
FourVector Tmu0_
Four-momentum flow of the cell.
void thermalize(const Particles &particles, double time, int ntest)
Main thermalize function, that chooses the algorithm to follow (BF or mode sampling).
void update_lattice(RectangularLattice< T > *lat, const LatticeUpdate update, const DensityType dens_type, const DensityParameters &par, const Particles &particles, const bool compute_gradient=false)
Updates the contents on the lattice.
void thermalize_BF_algo(QuantumNumbers &conserved_initial, double time, int ntest)
Samples particles according to the BF algorithm by making use of the.
LatticeUpdate
Enumerator option for lattice updates.
const FourVector & momentum() const
Get the particle's 4-momentum.
double ns() const
Get net strangeness density of the cell in the computational frame.
std::pair< int, int > sample()
Sample two numbers from given Poissonians with a fixed difference.
const ThermalizationAlgorithm algorithm_
Algorithm to choose for sampling of particles obeying conservation laws.
void set_formation_time(const double &form_time)
Set the absolute formation time.
double nb() const
Get net baryon density of the cell in the computational frame.
void add_values(const ParticleData &p)
Add the quantum numbers of a single particle to the collection.
void print_statistics(const Clock &clock) const
Generates standard output with information about the thermodynamic properties of the lattice...
void set_rest_frame_quantities(double T0, double mub0, double mus0, const ThreeVector v0)
Set all the rest frame quantities to some values, this is useful for testing.
void update_thermalizer_lattice(const Particles &particles, const DensityParameters &par, bool ignore_cells_under_threshold=true)
Compute all the thermodynamical quantities on the lattice from particles.
ThreeVector v_
Velocity of the rest frame.
double p_
Pressure in the rest frame.
static double net_baryon_density(double T, double mub, double mus)
Compute net baryon density.
ThreeVector threevec() const
HadronClass get_class(size_t typelist_index) const
Defines the class of hadrons by quantum numbers.
double ns_
Net strangeness density of the cell in the computational frame.
int baryon_number() const
double mub_
Net baryon chemical potential.
ParticleList sampled_list_
Newly generated particles by thermalizer.
bool is_tabulated() const
Create an EoS table or not?
void set_4momentum(const FourVector &momentum_vector)
Set the particle's 4-momentum directly.
int baryon_number() const
Non-strange mesons (S = 0) with electric cherge Q < 0.
A container for storing conserved values.
void compute_rest_frame_quantities(HadronGasEos &eos)
Temperature, chemical potentials and rest frame velocity are calculated given the hadron gas equation...
double e() const
Get energy density in the rest frame.
double sample_momenta_from_thermal(const double temperature, const double mass)
Samples a momentum from the Maxwell-Boltzmann (thermal) distribution in a faster way, given by Scott Pratt.
void compute_N_in_cells_mode_algo(F &&condition)
Computes average number of particles in each cell for the mode algorithm.
std::vector< size_t > cells_to_sample_
Cells above critical energy density.
Mesons with strangeness S < 0.
The intention of this class is to efficiently sample from the Bessel distribution ...
double N_total_in_cells_
Total number of particles over all cells in thermalization region.
void thermalize_mode_algo(QuantumNumbers &conserved_initial, double time)
Samples particles to the according to the mode algorithm.
Define the data structure for one element of the table.
void renormalize_momenta(ParticleList &plist, const FourVector required_total_momentum)
Changes energy and momenta of the particles in plist to match the required_total_momentum.
Clock tracks the time in the simulation.
Neutral non-strange mesons.
Non-strange mesons (S = 0) with electric cherge Q > 0.
ParticleData sample_in_random_cell_mode_algo(const double time, F &&condition)
Samples one particle and the species, cell, momentum and coordinate are chosen from the corresponding...
std::vector< double > mult_sort_
Real number multiplicity for each particle type.
ParticleList to_remove_
Particles to be removed after this thermalization step.
Mesons with strangeness S > 0.
std::vector< double > N_in_cells_
Number of particles to be sampled in one cell.
int binomial(const int N, const T &p)
Returns a binomially distributed random number.
double nb_
Net baryon density of the cell in the computational frame.
ThreeVector uniform_in_cell() const
const ParticleType & type() const
Get the type of the particle.
static double partial_density(const ParticleType &ptype, double T, double mub, double mus)
Compute partial density of one hadron sort.
FourVector Tmu0() const
Get Four-momentum flow of the cell.
void set_4position(const FourVector &pos)
Set the particle's 4-position directly.
double e_
Energy density in the rest frame.
The ThermLatticeNode class is intended to compute thermodynamical quantities in a cell given a set of...
HadronClass
Specifier to classify the different hadron species according to their quantum numbers.
std::vector< int > mult_int_
Integer multiplicity for each particle type.
double mus_
Net strangeness chemical potential.
int poisson(const T &lam)
Returns a Poisson distributed random number.
void sample_in_random_cell_BF_algo(ParticleList &plist, const double time, size_t type_index)
The total number of particles of species type_index is defined by mult_int_ array that is returned by...
double T() const
Get the temperature.
double mus() const
Get the net strangeness chemical potential.
std::array< double, 7 > mult_classes_
The different hadron species according to the enum defined in.
Angles provides a common interface for generating directions: i.e., two angles that should be interpr...
std::unique_ptr< RectangularLattice< ThermLatticeNode > > lat_
The lattice on which the thermodynamic quantities are calculated.
void distribute_isotropically()
Populate the object with a new direction.
The Particles class abstracts the storage and manipulation of particles.
double current_time() const
DensityType
Allows to choose which kind of density to calculate.
std::ostream & operator<<(std::ostream &out, const ActionPtr &action)
Convenience: dereferences the ActionPtr to Action.
The FourVector class holds relevant values in Minkowski spacetime with (+, −, −, −) metric signature.
double cell_volume_
Volume of a single cell, necessary to convert thermal densities to actual particle numbers...
double abs() const
calculate the lorentz invariant absolute value
HadronGasEos eos_
Hadron gas equation of state.
FourVector momentum() const
void from_table(EosTable::table_element &res, double e, double nb) const
Get the element of eos table.
ParticleData contains the dynamic information of a certain particle.
void boost_momentum(const ThreeVector &v)
Apply a Lorentz-boost to only the momentum.
const ParticleTypePtrList eos_typelist_
List of particle types from which equation of state is computed.
double mult_class(const HadronClass cl) const
static double pressure(double T, double mub, double mus)
Compute pressure .
ThreeVector v() const
Get 3-velocity of the rest frame.
static double energy_density(double T, double mub, double mus)
Compute energy density.