doc/2.0/grandcan__thermalizer_8h_source.html

/*

 *    Copyright (c) 2016-

 *      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 "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({"Cell_Number"}), origin, periodicity,

            conf.take({"Critical_Edens"}), conf.take({"Start_Time"}),

            conf.take({"Timestep"}),

            conf.take({"Algorithm"}, ThermalizationAlgorithm::BiasedBF),

            conf.take({"Microcanonical"}, false)) {}

  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 Particles& particles,

                                  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) {

    // 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);

    return particle;

  }


  void thermalize_mode_algo(QuantumNumbers& conserved_initial, double time);

  void thermalize(const Particles& particles, double time, int ntest);


  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_