Version: SMASH-3.2
grandcan_thermalizer.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2016-2020,2022,2024
3  * SMASH Team
4  *
5  * GNU General Public License (GPLv3 or later)
6  */
7 #ifndef SRC_INCLUDE_SMASH_GRANDCAN_THERMALIZER_H_
8 #define SRC_INCLUDE_SMASH_GRANDCAN_THERMALIZER_H_
9 
10 #include <memory>
11 #include <vector>
12 
13 #include "angles.h"
14 #include "clock.h"
15 #include "configuration.h"
16 #include "density.h"
17 #include "distributions.h"
18 #include "forwarddeclarations.h"
19 #include "hadgas_eos.h"
20 #include "input_keys.h"
21 #include "lattice.h"
22 #include "particledata.h"
23 #include "quantumnumbers.h"
24 
25 namespace smash {
26 
50  public:
67  void add_particle(const ParticleData& p, double factor);
92  void set_rest_frame_quantities(double T0, double mub0, double mus0,
93  double muq0, const ThreeVector v0);
95  FourVector Tmu0() const { return Tmu0_; }
97  double nb() const { return nb_; }
99  double ns() const { return ns_; }
101  double nq() const { return nq_; }
103  double e() const { return e_; }
105  double p() const { return p_; }
107  ThreeVector v() const { return v_; }
109  double T() const { return T_; }
111  double mub() const { return mub_; }
113  double mus() const { return mus_; }
115  double muq() const { return muq_; }
116 
117  private:
121  double nb_;
123  double ns_;
125  double nq_;
127  double e_;
129  double p_;
133  double T_;
135  double mub_;
137  double mus_;
139  double muq_;
140 };
141 
148 std::ostream& operator<<(std::ostream& s, const ThermLatticeNode& node);
149 
154 enum class HadronClass {
156  Baryon = 0,
158  Antibaryon = 1,
160  PositiveSMeson = 2,
162  NegativeSMeson = 3,
168  ZeroQZeroSMeson = 6,
169 };
170 
190  public:
207  GrandCanThermalizer(const std::array<double, 3> lat_sizes,
208  const std::array<int, 3> n_cells,
209  const std::array<double, 3> origin, bool periodicity,
210  double e_critical, double t_start, double delta_t,
211  ThermalizationAlgorithm algo, bool BF_microcanonical);
214  const std::array<double, 3> lat_sizes,
215  const std::array<double, 3> origin, bool periodicity)
217  lat_sizes, conf.take(InputKeys::forcedThermalization_cellNumber),
218  origin, periodicity,
219  conf.take(InputKeys::forcedThermalization_criticalEDensity),
220  conf.take(InputKeys::forcedThermalization_startTime),
221  conf.take(InputKeys::forcedThermalization_timestep),
222  conf.take(InputKeys::forcedThermalization_algorithm),
223  conf.take(InputKeys::forcedThermalization_microcanonical)) {}
229  bool is_time_to_thermalize(std::unique_ptr<Clock>& clock) const {
230  const double t = clock->current_time();
231  const int n = static_cast<int>(std::floor((t - t_start_) / period_));
232  return (t > t_start_ &&
233  t < t_start_ + n * period_ + clock->timestep_duration());
234  }
244  void update_thermalizer_lattice(const std::vector<Particles>& ensembles,
245  const DensityParameters& par,
246  bool ignore_cells_under_threshold = true);
256  void renormalize_momenta(ParticleList& plist,
257  const FourVector required_total_momentum);
258 
259  // Functions for BF-sampling algorithm
260 
269  void sample_multinomial(HadronClass particle_class, int N);
281  void sample_in_random_cell_BF_algo(ParticleList& plist, const double time,
282  size_t type_index);
295  void thermalize_BF_algo(QuantumNumbers& conserved_initial, double time,
296  int ntest);
297 
298  // Functions for mode-sampling algorithm
299 
304  template <typename F>
305  void compute_N_in_cells_mode_algo(F&& condition) {
306  N_in_cells_.clear();
307  N_total_in_cells_ = 0.0;
308  for (auto cell_index : cells_to_sample_) {
309  const ThermLatticeNode cell = (*lat_)[cell_index];
310  const double gamma = 1.0 / std::sqrt(1.0 - cell.v().sqr());
311  double N_tot = 0.0;
312  for (ParticleTypePtr i : eos_typelist_) {
313  if (condition(i->strangeness(), i->baryon_number(), i->charge())) {
314  // N_i = n u^mu dsigma_mu = (isochronous hypersurface) n * V * gamma
315  N_tot += lat_cell_volume_ * gamma *
316  HadronGasEos::partial_density(*i, cell.T(), cell.mub(),
317  cell.mus(), 0.0);
318  }
319  }
320  N_in_cells_.push_back(N_tot);
321  N_total_in_cells_ += N_tot;
322  }
323  }
324 
335  template <typename F>
337  F&& condition) {
338  // Choose random cell, probability = N_in_cell/N_total
339  double r = random::uniform(0.0, N_total_in_cells_);
340  double partial_sum = 0.0;
341  int index_only_thermalized = -1;
342  while (partial_sum < r) {
343  index_only_thermalized++;
344  partial_sum += N_in_cells_[index_only_thermalized];
345  }
346  const int cell_index = cells_to_sample_[index_only_thermalized];
347  const ThermLatticeNode cell = (*lat_)[cell_index];
348  const ThreeVector cell_center = lat_->cell_center(cell_index);
349  const double gamma = 1.0 / std::sqrt(1.0 - cell.v().sqr());
350  const double N_in_cell = N_in_cells_[index_only_thermalized];
351  // Which sort to sample - probability N_i/N_tot
352  r = random::uniform(0.0, N_in_cell);
353  double N_sum = 0.0;
354  ParticleTypePtr type_to_sample;
355  for (ParticleTypePtr i : eos_typelist_) {
356  if (!condition(i->strangeness(), i->baryon_number(), i->charge())) {
357  continue;
358  }
359  N_sum += lat_cell_volume_ * gamma *
360  HadronGasEos::partial_density(*i, cell.T(), cell.mub(),
361  cell.mus(), 0.0);
362  if (N_sum >= r) {
363  type_to_sample = i;
364  break;
365  }
366  }
367  ParticleData particle(*type_to_sample);
368  // Note: it's pole mass for resonances!
369  const double m = type_to_sample->mass();
370  // Position
371  particle.set_4position(FourVector(time, cell_center + uniform_in_cell()));
372  // Momentum
373  double momentum_radial = sample_momenta_from_thermal(cell.T(), m);
374  Angles phitheta;
375  phitheta.distribute_isotropically();
376  particle.set_4momentum(m, phitheta.threevec() * momentum_radial);
377  particle.boost_momentum(-cell.v());
378  particle.set_formation_time(time);
379  return particle;
380  }
381 
390  void thermalize_mode_algo(QuantumNumbers& conserved_initial, double time);
398  void thermalize(const Particles& particles, double time, int ntest);
399 
406  void print_statistics(const Clock& clock) const;
410  double e_crit() const { return e_crit_; }
412  ParticleList particles_to_remove() const { return to_remove_; }
414  ParticleList particles_to_insert() const { return sampled_list_; }
415 
416  private:
421  ParticleTypePtrList list_eos_particles() const {
422  ParticleTypePtrList res;
423  for (const ParticleType& ptype : ParticleType::list_all()) {
424  if (HadronGasEos::is_eos_particle(ptype)) {
425  res.push_back(&ptype);
426  }
427  }
428  return res;
429  }
434  HadronClass get_class(size_t typelist_index) const {
435  const int B = eos_typelist_[typelist_index]->baryon_number();
436  const int S = eos_typelist_[typelist_index]->strangeness();
437  const int ch = eos_typelist_[typelist_index]->charge();
438  // clang-format off
439  return (B > 0) ? HadronClass::Baryon :
440  (B < 0) ? HadronClass::Antibaryon :
446  // clang-format on
447  }
449  double mult_class(const HadronClass cl) const {
450  return mult_classes_[static_cast<size_t>(cl)];
451  }
453  std::vector<double> N_in_cells_;
455  std::vector<size_t> cells_to_sample_;
459  std::unique_ptr<RectangularLattice<ThermLatticeNode>> lat_;
461  ParticleList to_remove_;
463  ParticleList sampled_list_;
472  const ParticleTypePtrList eos_typelist_;
474  const size_t N_sorts_;
476  std::vector<double> mult_sort_;
478  std::vector<int> mult_int_;
483  std::array<double, 7> mult_classes_;
492  const double e_crit_;
494  const double t_start_;
496  const double period_;
501 };
502 
503 } // namespace smash
504 
505 #endif // SRC_INCLUDE_SMASH_GRANDCAN_THERMALIZER_H_
Angles provides a common interface for generating directions: i.e., two angles that should be interpr...
Definition: angles.h:59
ThreeVector threevec() const
Definition: angles.h:288
void distribute_isotropically()
Populate the object with a new direction.
Definition: angles.h:199
Clock tracks the time in the simulation.
Definition: clock.h:87
Interface to the SMASH configuration files.
A class to pre-calculate and store parameters relevant for density calculation.
Definition: density.h:92
The FourVector class holds relevant values in Minkowski spacetime with (+, −, −, −) metric signature.
Definition: fourvector.h:33
The GrandCanThermalizer class implements the following functionality:
void thermalize_BF_algo(QuantumNumbers &conserved_initial, double time, int ntest)
Samples particles according to the BF algorithm by making use of the.
void compute_N_in_cells_mode_algo(F &&condition)
Computes average number of particles in each cell for the mode algorithm.
void print_statistics(const Clock &clock) const
Generates standard output with information about the thermodynamic properties of the lattice,...
std::vector< double > mult_sort_
Real number multiplicity for each particle type.
std::vector< int > mult_int_
Integer multiplicity for each particle type.
GrandCanThermalizer(Configuration &conf, const std::array< double, 3 > lat_sizes, const std::array< double, 3 > origin, bool periodicity)
HadronGasEos eos_
Hadron gas equation of state.
ParticleList to_remove_
Particles to be removed after this thermalization step.
const bool BF_enforce_microcanonical_
Enforce energy conservation as part of BF sampling algorithm or not.
ThreeVector uniform_in_cell() const
std::unique_ptr< RectangularLattice< ThermLatticeNode > > lat_
The lattice on which the thermodynamic quantities are calculated.
ParticleList particles_to_insert() const
List of newly created particles to be inserted in the simulation.
double mult_class(const HadronClass cl) const
ParticleList particles_to_remove() const
List of particles to be removed from the simulation.
HadronClass get_class(size_t typelist_index) const
Defines the class of hadrons by quantum numbers.
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, bool BF_microcanonical)
Default constructor for the GranCanThermalizer to allocate the lattice.
const double t_start_
Starting time of the simulation.
const double period_
Defines periodicity of the lattice in fm.
ParticleTypePtrList list_eos_particles() const
Extracts the particles in the hadron gas equation of state from the complete list of particle types i...
const double e_crit_
Critical energy density above which cells are thermalized.
ParticleList sampled_list_
Newly generated particles by thermalizer.
std::array< double, 7 > mult_classes_
The different hadron species according to the enum defined in.
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...
void update_thermalizer_lattice(const std::vector< Particles > &ensembles, const DensityParameters &par, bool ignore_cells_under_threshold=true)
Compute all the thermodynamical quantities on the lattice from particles.
std::vector< size_t > cells_to_sample_
Cells above critical energy density.
const ThermalizationAlgorithm algorithm_
Algorithm to choose for sampling of particles obeying conservation laws.
bool is_time_to_thermalize(std::unique_ptr< Clock > &clock) const
Check that the clock is close to n * period of thermalization, since the thermalization only happens ...
RectangularLattice< ThermLatticeNode > & lattice() const
Getter function for the lattice.
double e_crit() const
Get the critical energy density.
const ParticleTypePtrList eos_typelist_
List of particle types from which equation of state is computed.
std::vector< double > N_in_cells_
Number of particles to be sampled in one cell.
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...
double lat_cell_volume_
Volume of a single lattice cell, necessary to convert thermal densities to actual particle numbers.
const size_t N_sorts_
Number of different species to be sampled.
double N_total_in_cells_
Total number of particles over all cells in thermalization region.
void thermalize(const Particles &particles, double time, int ntest)
Main thermalize function, that chooses the algorithm to follow (BF or mode sampling).
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.
void sample_multinomial(HadronClass particle_class, int N)
The sample_multinomial function samples integer numbers n_i distributed according to the multinomial ...
void thermalize_mode_algo(QuantumNumbers &conserved_initial, double time)
Samples particles to the according to the mode algorithm.
Class to handle the equation of state (EoS) of the hadron gas, consisting of all hadrons included in ...
Definition: hadgas_eos.h:125
static double partial_density(const ParticleType &ptype, double T, double mub, double mus, double muq, bool account_for_resonance_widths=false)
Compute partial density of one hadron sort.
Definition: hadgas_eos.cc:270
static bool is_eos_particle(const ParticleType &ptype)
Check if a particle belongs to the EoS.
Definition: hadgas_eos.h:355
ParticleData contains the dynamic information of a certain particle.
Definition: particledata.h:58
void set_4momentum(const FourVector &momentum_vector)
Set the particle's 4-momentum directly.
Definition: particledata.h:164
void set_4position(const FourVector &pos)
Set the particle's 4-position directly.
Definition: particledata.h:209
void set_formation_time(const double &form_time)
Set the absolute formation time.
Definition: particledata.h:251
void boost_momentum(const ThreeVector &v)
Apply a Lorentz-boost to only the momentum.
Definition: particledata.h:332
A pointer-like interface to global references to ParticleType objects.
Definition: particletype.h:679
Particle type contains the static properties of a particle species.
Definition: particletype.h:98
static const ParticleTypeList & list_all()
Definition: particletype.cc:51
double mass() const
Definition: particletype.h:145
The Particles class abstracts the storage and manipulation of particles.
Definition: particles.h:33
A container for storing conserved values.
A container class to hold all the arrays on the lattice and access them.
Definition: lattice.h:49
The ThermLatticeNode class is intended to compute thermodynamical quantities in a cell given a set of...
void compute_rest_frame_quantities(HadronGasEos &eos)
Temperature, chemical potentials and rest frame velocity are calculated given the hadron gas equation...
double muq() const
Get the net charge chemical potential.
FourVector Tmu0() const
Get Four-momentum flow of the cell.
double ns_
Net strangeness density of the cell in the computational frame.
void set_rest_frame_quantities(double T0, double mub0, double mus0, double muq0, const ThreeVector v0)
Set all the rest frame quantities to some values, this is useful for testing.
double p() const
Get pressure in the rest frame.
double p_
Pressure in the rest frame.
double mus_
Net strangeness chemical potential.
ThreeVector v() const
Get 3-velocity of the rest frame.
double ns() const
Get net strangeness density of the cell in the computational frame.
double e_
Energy density in the rest frame.
double nb_
Net baryon density of the cell in the computational frame.
void add_particle(const ParticleData &p, double factor)
Add particle contribution to Tmu0, nb, ns and nq May look like unused at first glance,...
double muq_
Net charge chemical potential.
void add_particle_for_derivatives(const ParticleData &, double, ThreeVector)
dummy function for update_lattice
double nq() const
Get net charge density of the cell in the computational frame.
double mub_
Net baryon chemical potential.
FourVector Tmu0_
Four-momentum flow of the cell.
double mus() const
Get the net strangeness chemical potential.
ThreeVector v_
Velocity of the rest frame.
ThermLatticeNode()
Default constructor of thermal quantities on the lattice returning thermodynamic quantities in comput...
double mub() const
Get the net baryon chemical potential.
double nb() const
Get net baryon density of the cell in the computational frame.
double nq_
Net charge density of the cell in the computational frame.
double e() const
Get energy density in the rest frame.
double T() const
Get the temperature.
The ThreeVector class represents a physical three-vector with the components .
Definition: threevector.h:31
double sqr() const
Definition: threevector.h:266
ThermalizationAlgorithm
Defines the algorithm used for the forced thermalization.
std::ostream & operator<<(std::ostream &out, const ActionPtr &action)
Convenience: dereferences the ActionPtr to Action.
Definition: action.h:546
constexpr int n
Neutron.
T uniform(T min, T max)
Definition: random.h:88
Definition: action.h:24
double sample_momenta_from_thermal(const double temperature, const double mass)
Samples a momentum from the Maxwell-Boltzmann (thermal) distribution in a faster way,...
HadronClass
Specifier to classify the different hadron species according to their quantum numbers.
@ Antibaryon
All anti-baryons.
@ ZeroQZeroSMeson
Neutral non-strange mesons.
@ NegativeSMeson
Mesons with strangeness S < 0.
@ NegativeQZeroSMeson
Non-strange mesons (S = 0) with electric cherge Q < 0.
@ PositiveSMeson
Mesons with strangeness S > 0.
@ Baryon
All baryons.
@ PositiveQZeroSMeson
Non-strange mesons (S = 0) with electric cherge Q > 0.
#define S(x, n)
Definition: sha256.cc:54
A container to keep track of all ever existed input keys.
Definition: input_keys.h:1067