grandcan_thermalizer.h

Go to the documentation of this file.

/*
 *    Copyright (c) 2016-2020,2022,2024-2025
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */
#ifndef SRC_INCLUDE_SMASH_GRANDCAN_THERMALIZER_H_
#define SRC_INCLUDE_SMASH_GRANDCAN_THERMALIZER_H_
 
#include <memory>
#include <vector>
 
#include "angles.h"
#include "clock.h"
#include "configuration.h"
#include "density.h"
#include "distributions.h"
#include "forwarddeclarations.h"
#include "hadgas_eos.h"
#include "input_keys.h"
#include "lattice.h"
#include "particledata.h"
#include "quantumnumbers.h"
 
namespace smash {
 
class ThermLatticeNode {
 public:
  ThermLatticeNode();
  void add_particle(const ParticleData& p, double factor);
  void add_particle_for_derivatives(const ParticleData&, double, ThreeVector) {}
  void compute_rest_frame_quantities(HadronGasEos& eos);
  void set_rest_frame_quantities(double T0, double mub0, double mus0,
                                 double muq0, const ThreeVector v0);
  FourVector Tmu0() const { return Tmu0_; }
  double nb() const { return nb_; }
  double ns() const { return ns_; }
  double nq() const { return nq_; }
  double e() const { return e_; }
  double p() const { return p_; }
  ThreeVector v() const { return v_; }
  double T() const { return T_; }
  double mub() const { return mub_; }
  double mus() const { return mus_; }
  double muq() const { return muq_; }
 
 private:
  FourVector Tmu0_;
  double nb_;
  double ns_;
  double nq_;
  double e_;
  double p_;
  ThreeVector v_;
  double T_;
  double mub_;
  double mus_;
  double muq_;
};
 
std::ostream& operator<<(std::ostream& s, const ThermLatticeNode& node);
 
enum class HadronClass {
  Baryon = 0,
  Antibaryon = 1,
  PositiveSMeson = 2,
  NegativeSMeson = 3,
  PositiveQZeroSMeson = 4,
  NegativeQZeroSMeson = 5,
  ZeroQZeroSMeson = 6,
};
 
class GrandCanThermalizer {
 public:
  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);
  GrandCanThermalizer(Configuration& conf,
                      const std::array<double, 3> lat_sizes,
                      const std::array<double, 3> origin, bool periodicity)
      : GrandCanThermalizer(
            lat_sizes, conf.take(InputKeys::forcedThermalization_cellNumber),
            origin, periodicity,
            conf.take(InputKeys::forcedThermalization_criticalEDensity),
            conf.take(InputKeys::forcedThermalization_startTime),
            conf.take(InputKeys::forcedThermalization_timestep),
            conf.take(InputKeys::forcedThermalization_algorithm),
            conf.take(InputKeys::forcedThermalization_microcanonical)) {}
  bool is_time_to_thermalize(std::unique_ptr<Clock>& clock) const {
    const double t = clock->current_time();
    const int n = static_cast<int>(std::floor((t - t_start_) / period_));
    return (t > t_start_ &&
            t < t_start_ + n * period_ + clock->timestep_duration());
  }
  void update_thermalizer_lattice(const std::vector<Particles>& ensembles,
                                  const DensityParameters& par,
                                  bool ignore_cells_under_threshold = true);
  ThreeVector uniform_in_cell() const;
  void renormalize_momenta(ParticleList& plist,
                           const FourVector required_total_momentum);
 
  // Functions for BF-sampling algorithm
 
  void sample_multinomial(HadronClass particle_class, int N);
  void sample_in_random_cell_BF_algo(ParticleList& plist, const double time,
                                     size_t type_index);
  void thermalize_BF_algo(QuantumNumbers& conserved_initial, double time,
                          int ntest);
 
  // Functions for mode-sampling algorithm
 
  template <typename F>
  void compute_N_in_cells_mode_algo(F&& condition) {
    N_in_cells_.clear();
    N_total_in_cells_ = 0.0;
    for (auto cell_index : cells_to_sample_) {
      const ThermLatticeNode cell = (*lat_)[cell_index];
      const double gamma = 1.0 / std::sqrt(1.0 - cell.v().sqr());
      double N_tot = 0.0;
      for (ParticleTypePtr i : eos_typelist_) {
        if (condition(i->strangeness(), i->baryon_number(), i->charge())) {
          // N_i = n u^mu dsigma_mu = (isochronous hypersurface) n * V * gamma
          N_tot += lat_cell_volume_ * gamma *
                   HadronGasEos::partial_density(*i, cell.T(), cell.mub(),
                                                 cell.mus(), 0.0);
        }
      }
      N_in_cells_.push_back(N_tot);
      N_total_in_cells_ += N_tot;
    }
  }
 
  template <typename F>
  ParticleData sample_in_random_cell_mode_algo(
      const double time, F&& condition,
      SpinInteractionType spin_interaction_type = SpinInteractionType::Off) {
    // Choose random cell, probability = N_in_cell/N_total
    double r = random::uniform(0.0, N_total_in_cells_);
    double partial_sum = 0.0;
    int index_only_thermalized = -1;
    while (partial_sum < r) {
      index_only_thermalized++;
      partial_sum += N_in_cells_[index_only_thermalized];
    }
    const int cell_index = cells_to_sample_[index_only_thermalized];
    const ThermLatticeNode cell = (*lat_)[cell_index];
    const ThreeVector cell_center = lat_->cell_center(cell_index);
    const double gamma = 1.0 / std::sqrt(1.0 - cell.v().sqr());
    const double N_in_cell = N_in_cells_[index_only_thermalized];
    // Which sort to sample - probability N_i/N_tot
    r = random::uniform(0.0, N_in_cell);
    double N_sum = 0.0;
    ParticleTypePtr type_to_sample;
    for (ParticleTypePtr i : eos_typelist_) {
      if (!condition(i->strangeness(), i->baryon_number(), i->charge())) {
        continue;
      }
      N_sum += lat_cell_volume_ * gamma *
               HadronGasEos::partial_density(*i, cell.T(), cell.mub(),
                                             cell.mus(), 0.0);
      if (N_sum >= r) {
        type_to_sample = i;
        break;
      }
    }
    ParticleData particle(*type_to_sample);
    // Note: it's pole mass for resonances!
    const double m = type_to_sample->mass();
    // Position
    particle.set_4position(FourVector(time, cell_center + uniform_in_cell()));
    // Momentum
    double momentum_radial = sample_momenta_from_thermal(cell.T(), m);
    Angles phitheta;
    phitheta.distribute_isotropically();
    particle.set_4momentum(m, phitheta.threevec() * momentum_radial);
    particle.boost_momentum(-cell.v());
    particle.set_formation_time(time);
    if (spin_interaction_type != SpinInteractionType::Off) {
      particle.set_unpolarized_spin_vector();
    }
    return particle;
  }
 
  void thermalize_mode_algo(
      QuantumNumbers& conserved_initial, double time,
      SpinInteractionType spin_interaction_type = SpinInteractionType::Off);
  void thermalize(
      const Particles& particles, double time, int ntest,
      SpinInteractionType spin_interaction_type = SpinInteractionType::Off);
 
  void print_statistics(const Clock& clock) const;
  RectangularLattice<ThermLatticeNode>& lattice() const { return *lat_; }
  double e_crit() const { return e_crit_; }
  ParticleList particles_to_remove() const { return to_remove_; }
  ParticleList particles_to_insert() const { return sampled_list_; }
 
 private:
  ParticleTypePtrList list_eos_particles() const {
    ParticleTypePtrList res;
    for (const ParticleType& ptype : ParticleType::list_all()) {
      if (HadronGasEos::is_eos_particle(ptype)) {
        res.push_back(&ptype);
      }
    }
    return res;
  }
  HadronClass get_class(size_t typelist_index) const {
    const int B = eos_typelist_[typelist_index]->baryon_number();
    const int S = eos_typelist_[typelist_index]->strangeness();
    const int ch = eos_typelist_[typelist_index]->charge();
    // clang-format off
    return (B > 0) ? HadronClass::Baryon :
           (B < 0) ? HadronClass::Antibaryon :
           (S > 0) ? HadronClass::PositiveSMeson :
           (S < 0) ? HadronClass::NegativeSMeson :
           (ch > 0) ? HadronClass::PositiveQZeroSMeson :
           (ch < 0) ? HadronClass::NegativeQZeroSMeson :
                      HadronClass::ZeroQZeroSMeson;
    // clang-format on
  }
  double mult_class(const HadronClass cl) const {
    return mult_classes_[static_cast<size_t>(cl)];
  }
  std::vector<double> N_in_cells_;
  std::vector<size_t> cells_to_sample_;
  HadronGasEos eos_ = HadronGasEos(true, false);
  std::unique_ptr<RectangularLattice<ThermLatticeNode>> lat_;
  ParticleList to_remove_;
  ParticleList sampled_list_;
  const ParticleTypePtrList eos_typelist_;
  const size_t N_sorts_;
  std::vector<double> mult_sort_;
  std::vector<int> mult_int_;
  std::array<double, 7> mult_classes_;
  double N_total_in_cells_;
  double lat_cell_volume_;
  const double e_crit_;
  const double t_start_;
  const double period_;
  const ThermalizationAlgorithm algorithm_;
  const bool BF_enforce_microcanonical_;
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_GRANDCAN_THERMALIZER_H_