smash::Experiment< Modus > Class Template Reference

#include <experiment.h>

template<typename Modus>
class smash::Experiment< Modus >

The main class, where the simulation of an experiment is executed.

The Experiment class owns all data (maybe indirectly) relevant for the execution of the experiment simulation. The experiment can be conducted in different running modi. Since the abstraction of these differences should not incur any overhead, the design is built around the Policy pattern.

The Policy pattern was defined by Andrei Alexandrescu in his book "Modern C++ Design: Generic Programming and Design Patterns Applied". Addison-Wesley:

A policy defines a class interface or a class template interface. The interface consists of one or all of the following: inner type definitions, member functions, and member variables.

The policy pattern can also be understood as a compile-time variant of the strategy pattern.

The Modus template parameter defines the "policy" of the Experiment class. It determines several aspects of the experiment execution at compile time. The original strategy pattern would select these differences at run time, thus incurring an overhead. This overhead becomes severe in cases where calls to strategy/policy functions are done very frequently. Using the policy pattern, the compiler can fully optimize: It creates a new instance of all functions in Experiment for all different Modus types.

Definition at line 177 of file experiment.h.

Inheritance diagram for smash::Experiment< Modus >:

Collaboration diagram for smash::Experiment< Modus >:

Public Member Functions
void	run () override
	Runs the experiment. More...
	Experiment (Configuration config, const bf::path &output_path)
	Create a new Experiment. More...
void	initialize_new_event ()
	This is called in the beginning of each event. More...
void	run_time_evolution ()
	Runs the time evolution of an event with fixed-sized time steps or without timesteps, from action to actions. More...
void	do_final_decays ()
	Performs the final decays of an event. More...
void	final_output ()
	Output at the end of an event. More...
Particles *	first_ensemble ()
	Provides external access to SMASH particles. More...
std::vector< Particles > *	all_ensembles ()
	Getter for all ensembles. More...
Modus *	modus ()
	Provides external access to SMASH calculation modus. More...
Public Member Functions inherited from smash::ExperimentBase
	ExperimentBase ()=default
virtual	~ExperimentBase ()=default
	The virtual destructor avoids undefined behavior when destroying derived objects. More...

Private Member Functions
bool	perform_action (Action &action, int i_ensemble)
	Perform the given action. More...
void	create_output (const std::string &format, const std::string &content, const bf::path &output_path, const OutputParameters &par)
	Create a list of output files. More...
void	propagate_and_shine (double to_time, Particles &particles)
	Propagate all particles until time to_time without any interactions and shine dileptons. More...
void	run_time_evolution_timestepless (Actions &actions, int i_ensemble, double end_time_propagation)
	Performs all the propagations and actions during a certain time interval neglecting the influence of the potentials. More...
void	intermediate_output ()
	Intermediate output during an event. More...
void	update_potentials ()
	Recompute potentials on lattices if necessary. More...
double	compute_min_cell_length (double dt) const
	Calculate the minimal size for the grid cells such that the ScatterActionsFinder will find all collisions within the maximal transverse distance (which is determined by the maximal cross section). More...
double	next_output_time () const
	Shortcut for next output time. More...

Private Attributes
ExperimentParameters	parameters_
	Struct of several member variables. More...
DensityParameters	density_param_
	Structure to precalculate and hold parameters for density computations. More...
Modus	modus_
	Instance of the Modus template parameter. More...
std::vector< Particles >	ensembles_
	Complete particle list, all ensembles in one vector. More...
std::unique_ptr< Potentials >	potentials_
	An instance of potentials class, that stores parameters of potentials, calculates them and their gradients. More...
std::unique_ptr< PauliBlocker >	pauli_blocker_
	An instance of PauliBlocker class that stores parameters needed for Pauli blocking calculations and computes phase-space density. More...
OutputsList	outputs_
	A list of output formaters. More...
OutputPtr	dilepton_output_
	The Dilepton output. More...
OutputPtr	photon_output_
	The Photon output. More...
std::vector< bool >	projectile_target_interact_
	Whether the projectile and the target collided. More...
std::vector< FourVector >	beam_momentum_ = {}
	The initial nucleons in the ColliderModus propagate with beam_momentum_, if Fermi motion is frozen. More...
std::vector< std::unique_ptr< ActionFinderInterface > >	action_finders_
	The Action finder objects. More...
std::unique_ptr< DecayActionsFinderDilepton >	dilepton_finder_
	The Dilepton Action Finder. More...
std::unique_ptr< ActionFinderInterface >	photon_finder_
	The (Scatter) Actions Finder for Direct Photons. More...
int	n_fractional_photons_
	Number of fractional photons produced per single reaction. More...
std::unique_ptr< DensityLattice >	j_QBS_lat_
	4-current for j_QBS lattice output More...
std::unique_ptr< DensityLattice >	jmu_B_lat_
	Baryon density on the lattice. More...
std::unique_ptr< DensityLattice >	jmu_I3_lat_
	Isospin projection density on the lattice. More...
std::unique_ptr< DensityLattice >	jmu_el_lat_
	Electric charge density on the lattice. More...
std::unique_ptr< FieldsLattice >	fields_lat_
	Mean-field A^mu on the lattice. More...
std::unique_ptr< DensityLattice >	jmu_custom_lat_
	Custom density on the lattices. More...
DensityType	dens_type_lattice_printout_ = DensityType::None
	Type of density for lattice printout. More...
std::unique_ptr< RectangularLattice< FourVector > >	UB_lat_ = nullptr
	Lattices for Skyrme or VDF potentials (evaluated in the local rest frame) times the baryon flow 4-velocity. More...
std::unique_ptr< RectangularLattice< FourVector > >	UI3_lat_ = nullptr
	Lattices for symmetry potentials (evaluated in the local rest frame) times the isospin flow 4-velocity. More...
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > >	FB_lat_
	Lattices for the electric and magnetic components of the Skyrme or VDF force. More...
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > >	FI3_lat_
	Lattices for the electric and magnetic component of the symmetry force. More...
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > >	EM_lat_
	Lattices for electric and magnetic field in fm^-2. More...
std::unique_ptr< RectangularLattice< EnergyMomentumTensor > >	Tmn_
	Lattices of energy-momentum tensors for printout. More...
std::unique_ptr< RectangularLattice< FourVector > >	old_jmu_auxiliary_
	Auxiliary lattice for values of jmu at a time step t0. More...
std::unique_ptr< RectangularLattice< FourVector > >	new_jmu_auxiliary_
	Auxiliary lattice for values of jmu at a time step t0 + dt. More...
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > >	four_gradient_auxiliary_
	Auxiliary lattice for calculating the four-gradient of jmu. More...
std::unique_ptr< RectangularLattice< FourVector > >	old_fields_auxiliary_
	Auxiliary lattice for values of Amu at a time step t0. More...
std::unique_ptr< RectangularLattice< FourVector > >	new_fields_auxiliary_
	Auxiliary lattice for values of Amu at a time step t0 + dt. More...
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > >	fields_four_gradient_auxiliary_
	Auxiliary lattice for calculating the four-gradient of Amu. More...
bool	printout_tmn_ = false
	Whether to print the energy-momentum tensor. More...
bool	printout_tmn_landau_ = false
	Whether to print the energy-momentum tensor in Landau frame. More...
bool	printout_v_landau_ = false
	Whether to print the 4-velocity in Landau frame. More...
bool	printout_j_QBS_ = false
	Whether to print the Q, B, S 4-currents. More...
bool	printout_lattice_td_ = false
	Whether to print the thermodynamics quantities evaluated on the lattices. More...
bool	printout_full_lattice_ascii_td_ = false
	Whether to print the thermodynamics quantities evaluated on the lattices, point by point, in ASCII format. More...
bool	printout_full_lattice_binary_td_ = false
	Whether to print the thermodynamics quantities evaluated on the lattices, point by point, in Binary format. More...
bool	printout_full_lattice_any_td_ = false
	Whether to print the thermodynamics quantities evaluated on the lattices, point by point, in any format. More...
std::unique_ptr< GrandCanThermalizer >	thermalizer_
	Instance of class used for forced thermalization. More...
StringProcess *	process_string_ptr_
	Pointer to the string process class object, which is used to set the random seed for PYTHIA objects in each event. More...
const int	nevents_
	Number of events. More...
int	event_ = 0
	Current event. More...
const double	end_time_
	simulation time at which the evolution is stopped. More...
const double	delta_time_startup_
	The clock's timestep size at start up. More...
const bool	force_decays_
	This indicates whether we force all resonances to decay in the last timestep. More...
const bool	use_grid_
	This indicates whether to use the grid. More...
const ExpansionProperties	metric_
	This struct contains information on the metric to be used. More...
const bool	dileptons_switch_
	This indicates whether dileptons are switched on. More...
const bool	photons_switch_
	This indicates whether photons are switched on. More...
const bool	bremsstrahlung_switch_
	This indicates whether bremsstrahlung is switched on. More...
const bool	IC_output_switch_
	This indicates whether the IC output is enabled. More...
const TimeStepMode	time_step_mode_
	This indicates whether to use time steps. More...
double	max_transverse_distance_sqr_ = std::numeric_limits<double>::max()
	Maximal distance at which particles can interact in case of the geometric criterion, squared. More...
QuantumNumbers	conserved_initial_
	The conserved quantities of the system. More...
double	initial_mean_field_energy_
	The initial total mean field energy in the system. More...
SystemTimePoint	time_start_ = SystemClock::now()
	system starting time of the simulation More...
DensityType	dens_type_ = DensityType::None
	Type of density to be written to collision headers. More...
uint64_t	interactions_total_ = 0
	Total number of interactions for current timestep. More...
uint64_t	previous_interactions_total_ = 0
	Total number of interactions for previous timestep. More...
uint64_t	wall_actions_total_ = 0
	Total number of wall-crossings for current timestep. More...
uint64_t	previous_wall_actions_total_ = 0
	Total number of wall-crossings for previous timestep. More...
uint64_t	total_pauli_blocked_ = 0
	Total number of Pauli-blockings for current timestep. More...
uint64_t	total_hypersurface_crossing_actions_ = 0
	Total number of particles removed from the evolution in hypersurface crossing actions. More...
uint64_t	discarded_interactions_total_ = 0
	Total number of discarded interactions, because they were invalidated before they could be performed. More...
double	total_energy_removed_ = 0.0
	Total energy removed from the system in hypersurface crossing actions. More...
int64_t	seed_ = -1
	random seed for the next event. More...

Friends
class	ExperimentBase
std::ostream &	operator<< (std::ostream &out, const Experiment &e)
	Writes the initial state for the Experiment to the output stream. More...

Additional Inherited Members
Static Public Member Functions inherited from smash::ExperimentBase
static std::unique_ptr< ExperimentBase >	create (Configuration config, const bf::path &output_path)
	Factory method that creates and initializes a new Experiment<Modus>. More...

Constructor & Destructor Documentation

Experiment()

template<typename Modus >

smash::Experiment< Modus >::Experiment ( Configuration  config, const bf::path &  output_path  )

explicit

Create a new Experiment.

This constructor is only called from the ExperimentBase::create factory method.

Parameters

[in]	config	The Configuration object contains all initial setup of the experiment. It is forwarded to the constructors of member variables as needed. Note that the object is passed by non-const reference. This is only necessary for bookkeeping: Values are not only read, but actually taken out of the object. Thus, all values that remain were not used.
[in]	output_path	The directory where the output files are written.

Member Function Documentation

run()

template<typename Modus >

void smash::Experiment< Modus >::run

overridevirtual

Runs the experiment.

The constructor does the setup of the experiment. The run function executes the complete experiment.

Implements smash::ExperimentBase.

Definition at line 2823 of file experiment.h.

                            {
  const auto &mainlog = logg[LMain];
  for (event_ = 0; event_ < nevents_; event_++) {
    mainlog.info() << "Event " << event_;
 
    // Sample initial particles, start clock, some printout and book-keeping
    initialize_new_event();
 
    run_time_evolution();
 
    if (force_decays_) {
      do_final_decays();
    }
 
    // Output at event end
    final_output();
  }
}

initialize_new_event()

template<typename Modus >

void smash::Experiment< Modus >::initialize_new_event

This is called in the beginning of each event.

It initializes particles according to selected modus, resets the clock and saves the initial conserved quantities for subsequent sanity checks.

Definition at line 1840 of file experiment.h.

                                             {
  random::set_seed(seed_);
  logg[LExperiment].info() << "random number seed: " << seed_;
  /* Set seed for the next event. It has to be positive, so it can be entered
   * in the config.
   *
   * We have to be careful about the minimal integer, whose absolute value
   * cannot be represented. */
  int64_t r = random::advance();
  while (r == INT64_MIN) {
    r = random::advance();
  }
  seed_ = std::abs(r);
  /* Set the random seed used in PYTHIA hadronization
   * to be same with the SMASH one.
   * In this way we ensure that the results are reproducible
   * for every event if one knows SMASH random seed. */
  if (process_string_ptr_ != NULL) {
    process_string_ptr_->init_pythia_hadron_rndm();
  }
 
  for (Particles &particles : ensembles_) {
    particles.reset();
  }
 
  // Sample particles according to the initial conditions
  double start_time = -1.0;
 
  // Sample impact parameter only once per all ensembles
  // It should be the same for all ensembles
  if (modus_.is_collider()) {
    modus_.sample_impact();
    logg[LExperiment].info("Impact parameter = ", modus_.impact_parameter(),
                           " fm");
  }
  for (Particles &particles : ensembles_) {
    start_time = modus_.initial_conditions(&particles, parameters_);
  }
  /* For box modus make sure that particles are in the box. In principle, after
   * a correct initialization they should be, so this is just playing it safe.
   */
  for (Particles &particles : ensembles_) {
    modus_.impose_boundary_conditions(&particles, outputs_);
  }
  // Reset the simulation clock
  double timestep = delta_time_startup_;
 
  switch (time_step_mode_) {
    case TimeStepMode::Fixed:
      break;
    case TimeStepMode::None:
      timestep = end_time_ - start_time;
      // Take care of the box modus + timestepless propagation
      const double max_dt = modus_.max_timestep(max_transverse_distance_sqr_);
      if (max_dt > 0. && max_dt < timestep) {
        timestep = max_dt;
      }
      break;
  }
  std::unique_ptr<UniformClock> clock_for_this_event;
  if (modus_.is_list() && (timestep < 0.0)) {
    throw std::runtime_error(
        "Timestep for the given event is negative. \n"
        "This might happen if the formation times of the input particles are "
        "larger than the specified end time of the simulation.");
  }
  clock_for_this_event = make_unique<UniformClock>(start_time, timestep);
  parameters_.labclock = std::move(clock_for_this_event);
 
  // Reset the output clock
  parameters_.outputclock->reset(start_time, true);
  // remove time before starting time in case of custom output times.
  parameters_.outputclock->remove_times_in_past(start_time);
 
  logg[LExperiment].debug(
      "Lab clock: t_start = ", parameters_.labclock->current_time(),
      ", dt = ", parameters_.labclock->timestep_duration());
 
  /* Save the initial conserved quantum numbers and total momentum in
   * the system for conservation checks */
  conserved_initial_ = QuantumNumbers(ensembles_);
  wall_actions_total_ = 0;
  previous_wall_actions_total_ = 0;
  interactions_total_ = 0;
  previous_interactions_total_ = 0;
  discarded_interactions_total_ = 0;
  total_pauli_blocked_ = 0;
  projectile_target_interact_.resize(parameters_.n_ensembles, false);
  total_hypersurface_crossing_actions_ = 0;
  total_energy_removed_ = 0.0;
  // Print output headers
  logg[LExperiment].info() << hline;
  logg[LExperiment].info() << "Time[fm]   Ekin[GeV]   E_MF[GeV]  ETotal[GeV]  "
                           << "ETot/N[GeV]  D(ETot/N)[GeV] Scatt&Decays  "
                           << "Particles     Comp.Time";
  logg[LExperiment].info() << hline;
  double E_mean_field = 0.0;
  if (potentials_) {
    // update_potentials();
    // if (parameters.outputclock->current_time() == 0.0 )
    // using the lattice is necessary
    if ((jmu_B_lat_ != nullptr)) {
      update_lattice(jmu_B_lat_.get(), old_jmu_auxiliary_.get(),
                     new_jmu_auxiliary_.get(), four_gradient_auxiliary_.get(),
                     LatticeUpdate::EveryTimestep, DensityType::Baryon,
                     density_param_, ensembles_,
                     parameters_.labclock->timestep_duration(), true);
      // Because there was no lattice at t=-Delta_t, the time derivatives
      // drho_dt and dj^mu/dt at t=0 are huge, while they shouldn't be; we
      // overwrite the time derivative to zero by hand.
      for (auto &node : *jmu_B_lat_) {
        node.overwrite_drho_dt_to_zero();
        node.overwrite_djmu_dt_to_zero();
      }
      E_mean_field = calculate_mean_field_energy(*potentials_, *jmu_B_lat_,
                                                 EM_lat_.get(), parameters_);
    }
  }
  initial_mean_field_energy_ = E_mean_field;
  logg[LExperiment].info() << format_measurements(
      ensembles_, 0u, conserved_initial_, time_start_,
      parameters_.labclock->current_time(), E_mean_field,
      initial_mean_field_energy_);
 
  // Output at event start
  for (const auto &output : outputs_) {
    for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
      auto event_info =
          fill_event_info(ensembles_, E_mean_field, modus_.impact_parameter(),
                          parameters_, projectile_target_interact_[i_ens]);
      output->at_eventstart(ensembles_[i_ens],
                            // Pretend each ensemble is an independent event
                            event_ * parameters_.n_ensembles + i_ens,
                            event_info);
    }
    // For thermodynamic output
    output->at_eventstart(ensembles_, event_);
    // For thermodynamic lattice output
    if (printout_full_lattice_any_td_) {
      switch (dens_type_lattice_printout_) {
        case DensityType::Baryon:
          output->at_eventstart(event_, ThermodynamicQuantity::EckartDensity,
                                DensityType::Baryon, *jmu_B_lat_);
          break;
        case DensityType::BaryonicIsospin:
          output->at_eventstart(event_, ThermodynamicQuantity::EckartDensity,
                                DensityType::BaryonicIsospin, *jmu_I3_lat_);
          break;
        case DensityType::None:
          break;
        default:
          output->at_eventstart(event_, ThermodynamicQuantity::EckartDensity,
                                DensityType::BaryonicIsospin, *jmu_custom_lat_);
      }
      if (printout_tmn_) {
        output->at_eventstart(event_, ThermodynamicQuantity::Tmn,
                              dens_type_lattice_printout_, *Tmn_);
      }
      if (printout_tmn_landau_) {
        output->at_eventstart(event_, ThermodynamicQuantity::TmnLandau,
                              dens_type_lattice_printout_, *Tmn_);
      }
      if (printout_v_landau_) {
        output->at_eventstart(event_, ThermodynamicQuantity::LandauVelocity,
                              dens_type_lattice_printout_, *Tmn_);
      }
      if (printout_j_QBS_) {
        output->at_eventstart(event_, ThermodynamicQuantity::j_QBS,
                              dens_type_lattice_printout_, *j_QBS_lat_);
      }
    }
  }
 
  /* In the ColliderModus, if Fermi motion is frozen, assign the beam momenta
   * to the nucleons in both the projectile and the target. Every ensemble
   * gets the same beam momenta, so no need to create beam_momenta_ vector
   * for every ensemble.
   */
  if (modus_.is_collider() && modus_.fermi_motion() == FermiMotion::Frozen) {
    for (ParticleData &particle : ensembles_[0]) {
      const double m = particle.effective_mass();
      double v_beam = 0.0;
      if (particle.belongs_to() == BelongsTo::Projectile) {
        v_beam = modus_.velocity_projectile();
      } else if (particle.belongs_to() == BelongsTo::Target) {
        v_beam = modus_.velocity_target();
      }
      const double gamma = 1.0 / std::sqrt(1.0 - v_beam * v_beam);
      beam_momentum_.emplace_back(
          FourVector(gamma * m, 0.0, 0.0, gamma * v_beam * m));
    }  // loop over particles
  }
}

Here is the call graph for this function:

run_time_evolution()

template<typename Modus >

void smash::Experiment< Modus >::run_time_evolution

Runs the time evolution of an event with fixed-sized time steps or without timesteps, from action to actions.

Within one timestep (fixed) evolution from action to action is invoked.

Definition at line 2186 of file experiment.h.

                                           {
  while (parameters_.labclock->current_time() < end_time_) {
    const double t = parameters_.labclock->current_time();
    const double dt =
        std::min(parameters_.labclock->timestep_duration(), end_time_ - t);
    logg[LExperiment].debug("Timestepless propagation for next ", dt, " fm/c.");
 
    // Perform forced thermalization if required
    if (thermalizer_ &&
        thermalizer_->is_time_to_thermalize(parameters_.labclock)) {
      const bool ignore_cells_under_treshold = true;
      // Thermodynamics in thermalizer is computed from all ensembles,
      // but thermalization actions act on each ensemble independently
      thermalizer_->update_thermalizer_lattice(ensembles_, density_param_,
                                               ignore_cells_under_treshold);
      const double current_t = parameters_.labclock->current_time();
      for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
        thermalizer_->thermalize(ensembles_[i_ens], current_t,
                                 parameters_.testparticles);
        ThermalizationAction th_act(*thermalizer_, current_t);
        if (th_act.any_particles_thermalized()) {
          perform_action(th_act, i_ens);
        }
      }
    }
 
    std::vector<Actions> actions(parameters_.n_ensembles);
    for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
      actions[i_ens].clear();
      if (ensembles_[i_ens].size() > 0 && action_finders_.size() > 0) {
        /* (1.a) Create grid. */
        const double min_cell_length = compute_min_cell_length(dt);
        logg[LExperiment].debug("Creating grid with minimal cell length ",
                                min_cell_length);
        const auto &grid =
            use_grid_ ? modus_.create_grid(ensembles_[i_ens], min_cell_length,
                                           dt, parameters_.coll_crit)
                      : modus_.create_grid(ensembles_[i_ens], min_cell_length,
                                           dt, parameters_.coll_crit,
                                           CellSizeStrategy::Largest);
 
        const double gcell_vol = grid.cell_volume();
        /* (1.b) Iterate over cells and find actions. */
        grid.iterate_cells(
            [&](const ParticleList &search_list) {
              for (const auto &finder : action_finders_) {
                actions[i_ens].insert(finder->find_actions_in_cell(
                    search_list, dt, gcell_vol, beam_momentum_));
              }
            },
            [&](const ParticleList &search_list,
                const ParticleList &neighbors_list) {
              for (const auto &finder : action_finders_) {
                actions[i_ens].insert(finder->find_actions_with_neighbors(
                    search_list, neighbors_list, dt, beam_momentum_));
              }
            });
      }
    }
 
    /* \todo (optimizations) Adapt timestep size here */
 
    /* (2) Propagate from action to action until next output or timestep end */
    const double end_timestep_time =
        std::min(parameters_.labclock->next_time(), end_time_);
    while (next_output_time() <= end_timestep_time) {
      for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
        run_time_evolution_timestepless(actions[i_ens], i_ens,
                                        next_output_time());
      }
      ++(*parameters_.outputclock);
 
      // Avoid duplication of final output
      if (parameters_.outputclock->current_time() < end_time_) {
        intermediate_output();
      }
    }
    for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
      run_time_evolution_timestepless(actions[i_ens], i_ens, end_timestep_time);
    }
 
    /* (3) Update potentials (if computed on the lattice) and
     *     compute new momenta according to equations of motion */
    if (potentials_) {
      update_potentials();
      update_momenta(ensembles_, parameters_.labclock->timestep_duration(),
                     *potentials_, FB_lat_.get(), FI3_lat_.get(),
                     EM_lat_.get());
    }
 
    /* (4) Expand universe if non-minkowskian metric; updates
     *     positions and momenta according to the selected expansion */
    if (metric_.mode_ != ExpansionMode::NoExpansion) {
      for (Particles &particles : ensembles_) {
        expand_space_time(&particles, parameters_, metric_);
      }
    }
 
    ++(*parameters_.labclock);
 
    /* (5) Check conservation laws.
     *
     * Check conservation of conserved quantities if potentials and string
     * fragmentation are off.  If potentials are on then momentum is conserved
     * only in average.  If string fragmentation is on, then energy and
     * momentum are only very roughly conserved in high-energy collisions. */
    if (!potentials_ && !parameters_.strings_switch &&
        metric_.mode_ == ExpansionMode::NoExpansion && !IC_output_switch_) {
      std::string err_msg = conserved_initial_.report_deviations(ensembles_);
      if (!err_msg.empty()) {
        logg[LExperiment].error() << err_msg;
        throw std::runtime_error("Violation of conserved quantities!");
      }
    }
  }
 
  if (pauli_blocker_) {
    logg[LExperiment].info(
        "Interactions: Pauli-blocked/performed = ", total_pauli_blocked_, "/",
        interactions_total_ - wall_actions_total_);
  }
}

Here is the call graph for this function:

do_final_decays()

template<typename Modus >

void smash::Experiment< Modus >::do_final_decays

Performs the final decays of an event.

Definition at line 2648 of file experiment.h.

                                        {
  /* At end of time evolution: Force all resonances to decay. In order to handle
   * decay chains, we need to loop until no further actions occur. */
  uint64_t interactions_old = 0;
  do {
    interactions_old = interactions_total_;
 
    for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
      Actions actions;
 
      // Dileptons: shining of remaining resonances
      if (dilepton_finder_ != nullptr) {
        for (const auto &output : outputs_) {
          dilepton_finder_->shine_final(ensembles_[i_ens], output.get(), true);
        }
      }
      // Find actions.
      for (const auto &finder : action_finders_) {
        actions.insert(finder->find_final_actions(ensembles_[i_ens]));
      }
      // Perform actions.
      while (!actions.is_empty()) {
        perform_action(*actions.pop(), i_ens);
      }
    }
    // loop until no more decays occur
  } while (interactions_total_ > interactions_old);
 
  // Dileptons: shining of stable particles at the end
  if (dilepton_finder_ != nullptr) {
    for (const auto &output : outputs_) {
      for (Particles &particles : ensembles_) {
        dilepton_finder_->shine_final(particles, output.get(), false);
      }
    }
  }
}

Here is the call graph for this function:

final_output()

template<typename Modus >

void smash::Experiment< Modus >::final_output

Output at the end of an event.

Definition at line 2687 of file experiment.h.

                                     {
  /* make sure the experiment actually ran (note: we should compare this
   * to the start time, but we don't know that. Therefore, we check that
   * the time is positive, which should heuristically be the same). */
  double E_mean_field = 0.0;
  if (likely(parameters_.labclock > 0)) {
    const uint64_t wall_actions_this_interval =
        wall_actions_total_ - previous_wall_actions_total_;
    const uint64_t interactions_this_interval = interactions_total_ -
                                                previous_interactions_total_ -
                                                wall_actions_this_interval;
    if (potentials_) {
      // using the lattice is necessary
      if ((jmu_B_lat_ != nullptr)) {
        E_mean_field = calculate_mean_field_energy(*potentials_, *jmu_B_lat_,
                                                   EM_lat_.get(), parameters_);
      }
    }
    logg[LExperiment].info() << format_measurements(
        ensembles_, interactions_this_interval, conserved_initial_, time_start_,
        end_time_, E_mean_field, initial_mean_field_energy_);
    int total_particles = 0;
    for (const Particles &particles : ensembles_) {
      total_particles += particles.size();
    }
    if (IC_output_switch_ && (total_particles == 0)) {
      // Verify there is no more energy in the system if all particles were
      // removed when crossing the hypersurface
      const double remaining_energy =
          conserved_initial_.momentum().x0() - total_energy_removed_;
      if (remaining_energy > really_small) {
        throw std::runtime_error(
            "There is remaining energy in the system although all particles "
            "were removed.\n"
            "E_remain = " +
            std::to_string(remaining_energy) + " [GeV]");
      } else {
        logg[LExperiment].info() << hline;
        logg[LExperiment].info()
            << "Time real: " << SystemClock::now() - time_start_;
        logg[LExperiment].info()
            << "Interactions before reaching hypersurface: "
            << interactions_total_ - wall_actions_total_ -
                   total_hypersurface_crossing_actions_;
        logg[LExperiment].info()
            << "Total number of particles removed on hypersurface: "
            << total_hypersurface_crossing_actions_;
      }
    } else {
      const double precent_discarded =
          interactions_total_ > 0
              ? static_cast<double>(discarded_interactions_total_) * 100.0 /
                    interactions_total_
              : 0.0;
      std::stringstream msg_discarded;
      msg_discarded
          << "Discarded interaction number: " << discarded_interactions_total_
          << " (" << precent_discarded
          << "% of the total interaction number including wall crossings)";
 
      logg[LExperiment].info() << hline;
      logg[LExperiment].info()
          << "Time real: " << SystemClock::now() - time_start_;
      logg[LExperiment].debug() << msg_discarded.str();
 
      if (parameters_.coll_crit == CollisionCriterion::Stochastic &&
          precent_discarded > 1.0) {
        // The choosen threshold of 1% is a heuristical value
        logg[LExperiment].warn()
            << msg_discarded.str()
            << "\nThe number of discarded interactions is large, which means "
               "the assumption for the stochastic criterion of\n1 interaction "
               "per particle per timestep is probably violated. Consider "
               "reducing the timestep size.";
      }
 
      logg[LExperiment].info() << "Final interaction number: "
                               << interactions_total_ - wall_actions_total_;
    }
 
    // Check if there are unformed particles
    int unformed_particles_count = 0;
    for (const Particles &particles : ensembles_) {
      for (const ParticleData &particle : particles) {
        if (particle.formation_time() > end_time_) {
          unformed_particles_count++;
        }
      }
    }
    if (unformed_particles_count > 0) {
      logg[LExperiment].warn(
          "End time might be too small. ", unformed_particles_count,
          " unformed particles were found at the end of the evolution.");
    }
  }
 
  for (const auto &output : outputs_) {
    for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
      auto event_info =
          fill_event_info(ensembles_, E_mean_field, modus_.impact_parameter(),
                          parameters_, projectile_target_interact_[i_ens]);
      output->at_eventend(ensembles_[i_ens],
                          // Pretend each ensemble is an independent event
                          event_ * parameters_.n_ensembles + i_ens, event_info);
    }
    // For thermodynamic output
    output->at_eventend(ensembles_, event_);
 
    // For thermodynamic lattice output
    if (dens_type_lattice_printout_ != DensityType::None) {
      output->at_eventend(event_, ThermodynamicQuantity::EckartDensity,
                          dens_type_lattice_printout_);
      output->at_eventend(ThermodynamicQuantity::EckartDensity);
    }
    if (printout_tmn_) {
      output->at_eventend(event_, ThermodynamicQuantity::Tmn,
                          dens_type_lattice_printout_);
      output->at_eventend(ThermodynamicQuantity::Tmn);
    }
    if (printout_tmn_landau_) {
      output->at_eventend(event_, ThermodynamicQuantity::TmnLandau,
                          dens_type_lattice_printout_);
      output->at_eventend(ThermodynamicQuantity::TmnLandau);
    }
    if (printout_v_landau_) {
      output->at_eventend(event_, ThermodynamicQuantity::LandauVelocity,
                          dens_type_lattice_printout_);
      output->at_eventend(ThermodynamicQuantity::LandauVelocity);
    }
    if (printout_j_QBS_) {
      output->at_eventend(ThermodynamicQuantity::j_QBS);
    }
  }
}

Here is the call graph for this function:

first_ensemble()

template<typename Modus >

Particles* smash::Experiment< Modus >::first_ensemble ( )

inline

Provides external access to SMASH particles.

This is helpful if SMASH is used as a 3rd-party library.

Definition at line 231 of file experiment.h.

{ return &ensembles_[0]; }

all_ensembles()

template<typename Modus >

std::vector<Particles>* smash::Experiment< Modus >::all_ensembles ( )

inline

Getter for all ensembles.

Definition at line 233 of file experiment.h.

{ return &ensembles_; }

modus()

template<typename Modus >

Modus* smash::Experiment< Modus >::modus ( )

inline

Provides external access to SMASH calculation modus.

This is helpful if SMASH is used as a 3rd-party library.

Definition at line 239 of file experiment.h.

{ return &modus_; }

perform_action()

template<typename Modus >

bool smash::Experiment< Modus >::perform_action ( Action &  action, int  i_ensemble  )

private

Perform the given action.

Parameters

[in]	action	The action to perform
[in]	i_ensemble	index of ensemble in which action is performed

Returns

False if the action is rejected either due to invalidity or Pauli-blocking, or true if it's accepted and performed.

create_output()

template<typename Modus >

void smash::Experiment< Modus >::create_output ( const std::string &  format, const std::string &  content, const bf::path &  output_path, const OutputParameters &  par  )

private

Create a list of output files.

Parameters

[in]	format	Format of the output file (e.g. Root, Oscar, Vtk)
[in]	content	Content of the output (e.g. particles, collisions)
[in]	output_path	Path of the output file
[in]	par	Output options.(e.g. Extended)

Definition at line 641 of file experiment.h.

                                                                       {
  logg[LExperiment].info() << "Adding output " << content << " of format "
                           << format << std::endl;
 
  if (format == "VTK" && content == "Particles") {
    outputs_.emplace_back(
        make_unique<VtkOutput>(output_path, content, out_par));
  } else if (format == "Root") {
#ifdef SMASH_USE_ROOT
    if (content == "Initial_Conditions") {
      outputs_.emplace_back(
          make_unique<RootOutput>(output_path, "SMASH_IC", out_par));
    } else {
      outputs_.emplace_back(
          make_unique<RootOutput>(output_path, content, out_par));
    }
#else
    logg[LExperiment].error(
        "Root output requested, but Root support not compiled in");
#endif
  } else if (format == "Binary") {
    if (content == "Collisions" || content == "Dileptons" ||
        content == "Photons") {
      outputs_.emplace_back(
          make_unique<BinaryOutputCollisions>(output_path, content, out_par));
    } else if (content == "Particles") {
      outputs_.emplace_back(
          make_unique<BinaryOutputParticles>(output_path, content, out_par));
    } else if (content == "Initial_Conditions") {
      outputs_.emplace_back(make_unique<BinaryOutputInitialConditions>(
          output_path, content, out_par));
    }
  } else if (format == "Oscar1999" || format == "Oscar2013") {
    outputs_.emplace_back(
        create_oscar_output(format, content, output_path, out_par));
  } else if (content == "Thermodynamics" && format == "ASCII") {
    outputs_.emplace_back(
        make_unique<ThermodynamicOutput>(output_path, content, out_par));
  } else if (content == "Thermodynamics" &&
             (format == "Lattice_ASCII" || format == "Lattice_Binary")) {
    printout_full_lattice_any_td_ = true;
    if (format == "Lattice_ASCII") {
      printout_full_lattice_ascii_td_ = true;
    }
    if (format == "Lattice_Binary") {
      printout_full_lattice_binary_td_ = true;
    }
    outputs_.emplace_back(make_unique<ThermodynamicLatticeOutput>(
        output_path, content, out_par, printout_full_lattice_ascii_td_,
        printout_full_lattice_binary_td_));
  } else if (content == "Thermodynamics" && format == "VTK") {
    printout_lattice_td_ = true;
    outputs_.emplace_back(
        make_unique<VtkOutput>(output_path, content, out_par));
  } else if (content == "Initial_Conditions" && format == "ASCII") {
    outputs_.emplace_back(
        make_unique<ICOutput>(output_path, "SMASH_IC", out_par));
  } else if (format == "HepMC") {
#ifdef SMASH_USE_HEPMC
    if (content == "Particles") {
      outputs_.emplace_back(make_unique<HepMcOutput>(
          output_path, "SMASH_HepMC_particles", false));
    } else if (content == "Collisions") {
      outputs_.emplace_back(make_unique<HepMcOutput>(
          output_path, "SMASH_HepMC_collisions", true));
    } else {
      logg[LExperiment].error(
          "HepMC only available for Particles and "
          "Collisions content. Requested for " +
          content + ".");
    }
#else
    logg[LExperiment].error(
        "HepMC output requested, but HepMC support not compiled in");
#endif
  } else if (content == "Coulomb" && format == "VTK") {
    outputs_.emplace_back(
        make_unique<VtkOutput>(output_path, "Fields", out_par));
  } else if (content == "Rivet") {
#ifdef SMASH_USE_RIVET
    // flag to ensure that the Rivet format has not been already assigned
    static bool rivet_format_already_selected = false;
    // if the next check is true, then we are trying to assign the format twice
    if (rivet_format_already_selected) {
      logg[LExperiment].warn(
          "Rivet output format can only be one, either YODA or YODA-full. "
          "Only your first valid choice will be used.");
      return;
    }
    if (format == "YODA") {
      outputs_.emplace_back(
          make_unique<RivetOutput>(output_path, "SMASH_Rivet", false, out_par));
      rivet_format_already_selected = true;
    } else if (format == "YODA-full") {
      outputs_.emplace_back(make_unique<RivetOutput>(
          output_path, "SMASH_Rivet_full", true, out_par));
      rivet_format_already_selected = true;
    } else {
      logg[LExperiment].error("Rivet format " + format +
                              "not one of YODA or YODA-full");
    }
#else
    logg[LExperiment].error(
        "Rivet output requested, but Rivet support not compiled in");
#endif
  } else {
    logg[LExperiment].error()
        << "Unknown combination of format (" << format << ") and content ("
        << content << "). Fix the config.";
  }
}

Here is the call graph for this function:

propagate_and_shine()

template<typename Modus >

void smash::Experiment< Modus >::propagate_and_shine ( double  to_time, Particles &  particles  )

private

Propagate all particles until time to_time without any interactions and shine dileptons.

Parameters

[in]	to_time	Time at the end of propagation [fm/c]
[in,out]	particles	Particles to be propagated

Definition at line 2310 of file experiment.h.

                                                                  {
  const double dt =
      propagate_straight_line(&particles, to_time, beam_momentum_);
  if (dilepton_finder_ != nullptr) {
    for (const auto &output : outputs_) {
      dilepton_finder_->shine(particles, output.get(), dt);
    }
  }
}

Here is the call graph for this function:

run_time_evolution_timestepless()

template<typename Modus >

void smash::Experiment< Modus >::run_time_evolution_timestepless ( Actions &  actions, int  i_ensemble, double  end_time_propagation  )

private

Performs all the propagations and actions during a certain time interval neglecting the influence of the potentials.

This function is called in either the time stepless cases or the cases with time steps.

Parameters

[in,out]	actions	Actions occur during a certain time interval. They provide the ending times of the propagations and are updated during the time interval.
[in]	i_ensemble	index of ensemble to be evolved
[in]	end_time_propagation	time until propagation should be performed

Definition at line 2336 of file experiment.h.

                                                                   {
  Particles &particles = ensembles_[i_ensemble];
  logg[LExperiment].debug(
      "Timestepless propagation: ", "Actions size = ", actions.size(),
      ", end time = ", end_time_propagation);
 
  // iterate over all actions
  while (!actions.is_empty()) {
    if (actions.earliest_time() > end_time_propagation) {
      break;
    }
    // get next action
    ActionPtr act = actions.pop();
    if (!act->is_valid(particles)) {
      discarded_interactions_total_++;
      logg[LExperiment].debug(~einhard::DRed(), "✘ ", act,
                              " (discarded: invalid)");
      continue;
    }
    logg[LExperiment].debug(~einhard::Green(), "✔ ", act,
                            ", action time = ", act->time_of_execution());
 
    /* (1) Propagate to the next action. */
    propagate_and_shine(act->time_of_execution(), particles);
 
    /* (2) Perform action.
     *
     * Update the positions of the incoming particles, because the information
     * in the action object will be outdated as the particles have been
     * propagated since the construction of the action. */
    act->update_incoming(particles);
    const bool performed = perform_action(*act, i_ensemble);
 
    /* No need to update actions for outgoing particles
     * if the action is not performed. */
    if (!performed) {
      continue;
    }
 
    /* (3) Update actions for newly-produced particles. */
 
    const double end_time_timestep =
        std::min(parameters_.labclock->next_time(), end_time_);
    assert(!(end_time_propagation > end_time_timestep));
    // New actions are always search until the end of the current timestep
    const double time_left = end_time_timestep - act->time_of_execution();
    const ParticleList &outgoing_particles = act->outgoing_particles();
    // Grid cell volume set to zero, since there is no grid
    const double gcell_vol = 0.0;
    for (const auto &finder : action_finders_) {
      // Outgoing particles can still decay, cross walls...
      actions.insert(finder->find_actions_in_cell(outgoing_particles, time_left,
                                                  gcell_vol, beam_momentum_));
      // ... and collide with other particles.
      actions.insert(finder->find_actions_with_surrounding_particles(
          outgoing_particles, particles, time_left, beam_momentum_));
    }
 
    check_interactions_total(interactions_total_);
  }
 
  propagate_and_shine(end_time_propagation, particles);
}

Here is the call graph for this function:

intermediate_output()

template<typename Modus >

void smash::Experiment< Modus >::intermediate_output

private

Intermediate output during an event.

Auxiliary variable to communicate the time in the computational frame at the functions printing the thermodynamics lattice output

Definition at line 2402 of file experiment.h.

                                            {
  const uint64_t wall_actions_this_interval =
      wall_actions_total_ - previous_wall_actions_total_;
  previous_wall_actions_total_ = wall_actions_total_;
  const uint64_t interactions_this_interval = interactions_total_ -
                                              previous_interactions_total_ -
                                              wall_actions_this_interval;
  previous_interactions_total_ = interactions_total_;
  double E_mean_field = 0.0;
  double computational_frame_time = 0.0;
  if (potentials_) {
    // using the lattice is necessary
    if ((jmu_B_lat_ != nullptr)) {
      E_mean_field = calculate_mean_field_energy(*potentials_, *jmu_B_lat_,
                                                 EM_lat_.get(), parameters_);
      /*
       * Mean field calculated in a box should remain approximately constant if
       * the system is in equilibrium, and so deviations from its original value
       * may signal a phase transition or other dynamical process. This
       * comparison only makes sense in the Box Modus, hence the condition.
       */
      if (modus_.is_box()) {
        double tmp = (E_mean_field - initial_mean_field_energy_) /
                     (E_mean_field + initial_mean_field_energy_);
        /*
         * This is displayed when the system evolves away from its initial
         * configuration (which is when the total mean field energy in the box
         *  deviates from its initial value).
         */
        if (std::abs(tmp) > 0.01) {
          logg[LExperiment].info()
              << "\n\n\n\t The mean field at t = "
              << parameters_.outputclock->current_time()
              << " [fm/c] differs from the mean field at t = 0:"
              << "\n\t\t                 initial_mean_field_energy_ = "
              << initial_mean_field_energy_ << " [GeV]"
              << "\n\t\t abs[(E_MF - E_MF(t=0))/(E_MF + E_MF(t=0))] = "
              << std::abs(tmp)
              << "\n\t\t                             E_MF/E_MF(t=0) = "
              << E_mean_field / initial_mean_field_energy_ << "\n\n";
        }
      }
    }
  }
 
  logg[LExperiment].info() << format_measurements(
      ensembles_, interactions_this_interval, conserved_initial_, time_start_,
      parameters_.outputclock->current_time(), E_mean_field,
      initial_mean_field_energy_);
  const LatticeUpdate lat_upd = LatticeUpdate::AtOutput;
 
  // save evolution data
  if (!(modus_.is_box() && parameters_.outputclock->current_time() <
                               modus_.equilibration_time())) {
    for (const auto &output : outputs_) {
      if (output->is_dilepton_output() || output->is_photon_output() ||
          output->is_IC_output()) {
        continue;
      }
      for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
        auto event_info =
            fill_event_info(ensembles_, E_mean_field, modus_.impact_parameter(),
                            parameters_, projectile_target_interact_[i_ens]);
 
        output->at_intermediate_time(ensembles_[i_ens], parameters_.outputclock,
                                     density_param_, event_info);
        computational_frame_time = event_info.current_time;
      }
      // For thermodynamic output
      output->at_intermediate_time(ensembles_, parameters_.outputclock,
                                   density_param_);
 
      // Thermodynamic output on the lattice versus time
      switch (dens_type_lattice_printout_) {
        case DensityType::Baryon:
          update_lattice(jmu_B_lat_.get(), lat_upd, DensityType::Baryon,
                         density_param_, ensembles_, false);
          output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
                                        DensityType::Baryon, *jmu_B_lat_);
          output->thermodynamics_lattice_output(*jmu_B_lat_,
                                                computational_frame_time);
          break;
        case DensityType::BaryonicIsospin:
          update_lattice(jmu_I3_lat_.get(), lat_upd,
                         DensityType::BaryonicIsospin, density_param_,
                         ensembles_, false);
          output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
                                        DensityType::BaryonicIsospin,
                                        *jmu_I3_lat_);
          output->thermodynamics_lattice_output(*jmu_I3_lat_,
                                                computational_frame_time);
          break;
        case DensityType::None:
          break;
        default:
          update_lattice(jmu_custom_lat_.get(), lat_upd,
                         dens_type_lattice_printout_, density_param_,
                         ensembles_, false);
          output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
                                        dens_type_lattice_printout_,
                                        *jmu_custom_lat_);
          output->thermodynamics_lattice_output(*jmu_custom_lat_,
                                                computational_frame_time);
      }
      if (printout_tmn_ || printout_tmn_landau_ || printout_v_landau_) {
        update_lattice(Tmn_.get(), lat_upd, dens_type_lattice_printout_,
                       density_param_, ensembles_, false);
        if (printout_tmn_) {
          output->thermodynamics_output(ThermodynamicQuantity::Tmn,
                                        dens_type_lattice_printout_, *Tmn_);
          output->thermodynamics_lattice_output(
              ThermodynamicQuantity::Tmn, *Tmn_, computational_frame_time);
        }
        if (printout_tmn_landau_) {
          output->thermodynamics_output(ThermodynamicQuantity::TmnLandau,
                                        dens_type_lattice_printout_, *Tmn_);
          output->thermodynamics_lattice_output(
              ThermodynamicQuantity::TmnLandau, *Tmn_,
              computational_frame_time);
        }
        if (printout_v_landau_) {
          output->thermodynamics_output(ThermodynamicQuantity::LandauVelocity,
                                        dens_type_lattice_printout_, *Tmn_);
          output->thermodynamics_lattice_output(
              ThermodynamicQuantity::LandauVelocity, *Tmn_,
              computational_frame_time);
        }
      }
      if (EM_lat_) {
        output->fields_output("Efield", "Bfield", *EM_lat_);
      }
      if (printout_j_QBS_) {
        output->thermodynamics_lattice_output(
            *j_QBS_lat_, computational_frame_time, ensembles_, density_param_);
      }
 
      if (thermalizer_) {
        output->thermodynamics_output(*thermalizer_);
      }
    }
  }
}

Here is the call graph for this function:

update_potentials()

template<typename Modus >

void smash::Experiment< Modus >::update_potentials

private

Recompute potentials on lattices if necessary.

Definition at line 2548 of file experiment.h.

                                          {
  if (potentials_) {
    if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
      update_lattice(jmu_I3_lat_.get(), old_jmu_auxiliary_.get(),
                     new_jmu_auxiliary_.get(), four_gradient_auxiliary_.get(),
                     LatticeUpdate::EveryTimestep, DensityType::BaryonicIsospin,
                     density_param_, ensembles_,
                     parameters_.labclock->timestep_duration(), true);
    }
    if ((potentials_->use_skyrme() || potentials_->use_symmetry()) &&
        jmu_B_lat_ != nullptr) {
      update_lattice(jmu_B_lat_.get(), old_jmu_auxiliary_.get(),
                     new_jmu_auxiliary_.get(), four_gradient_auxiliary_.get(),
                     LatticeUpdate::EveryTimestep, DensityType::Baryon,
                     density_param_, ensembles_,
                     parameters_.labclock->timestep_duration(), true);
      const size_t UBlattice_size = UB_lat_->size();
      for (size_t i = 0; i < UBlattice_size; i++) {
        auto jB = (*jmu_B_lat_)[i];
        const FourVector flow_four_velocity_B =
            std::abs(jB.rho()) > very_small_double ? jB.jmu_net() / jB.rho()
                                                   : FourVector();
        double baryon_density = jB.rho();
        ThreeVector baryon_grad_j0 = jB.grad_j0();
        ThreeVector baryon_dvecj_dt = jB.dvecj_dt();
        ThreeVector baryon_curl_vecj = jB.curl_vecj();
        if (potentials_->use_skyrme()) {
          (*UB_lat_)[i] =
              flow_four_velocity_B * potentials_->skyrme_pot(baryon_density);
          (*FB_lat_)[i] =
              potentials_->skyrme_force(baryon_density, baryon_grad_j0,
                                        baryon_dvecj_dt, baryon_curl_vecj);
        }
        if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
          auto jI3 = (*jmu_I3_lat_)[i];
          const FourVector flow_four_velocity_I3 =
              std::abs(jI3.rho()) > very_small_double
                  ? jI3.jmu_net() / jI3.rho()
                  : FourVector();
          (*UI3_lat_)[i] = flow_four_velocity_I3 *
                           potentials_->symmetry_pot(jI3.rho(), baryon_density);
          (*FI3_lat_)[i] = potentials_->symmetry_force(
              jI3.rho(), jI3.grad_j0(), jI3.dvecj_dt(), jI3.curl_vecj(),
              baryon_density, baryon_grad_j0, baryon_dvecj_dt,
              baryon_curl_vecj);
        }
      }
    }
    if (potentials_->use_coulomb()) {
      update_lattice(jmu_el_lat_.get(), LatticeUpdate::EveryTimestep,
                     DensityType::Charge, density_param_, ensembles_, true);
      for (size_t i = 0; i < EM_lat_->size(); i++) {
        ThreeVector electric_field = {0., 0., 0.};
        ThreeVector position = jmu_el_lat_->cell_center(i);
        jmu_el_lat_->integrate_volume(electric_field,
                                      Potentials::E_field_integrand,
                                      potentials_->coulomb_r_cut(), position);
        ThreeVector magnetic_field = {0., 0., 0.};
        jmu_el_lat_->integrate_volume(magnetic_field,
                                      Potentials::B_field_integrand,
                                      potentials_->coulomb_r_cut(), position);
        (*EM_lat_)[i] = std::make_pair(electric_field, magnetic_field);
      }
    }  // if ((potentials_->use_skyrme() || ...
    if (potentials_->use_vdf() && jmu_B_lat_ != nullptr) {
      update_lattice(jmu_B_lat_.get(), old_jmu_auxiliary_.get(),
                     new_jmu_auxiliary_.get(), four_gradient_auxiliary_.get(),
                     LatticeUpdate::EveryTimestep, DensityType::Baryon,
                     density_param_, ensembles_,
                     parameters_.labclock->timestep_duration(), true);
      if (parameters_.field_derivatives_mode == FieldDerivativesMode::Direct) {
        update_fields_lattice(
            fields_lat_.get(), old_fields_auxiliary_.get(),
            new_fields_auxiliary_.get(), fields_four_gradient_auxiliary_.get(),
            jmu_B_lat_.get(), LatticeUpdate::EveryTimestep, *potentials_,
            parameters_.labclock->timestep_duration());
      }
      const size_t UBlattice_size = UB_lat_->size();
      for (size_t i = 0; i < UBlattice_size; i++) {
        auto jB = (*jmu_B_lat_)[i];
        (*UB_lat_)[i] = potentials_->vdf_pot(jB.rho(), jB.jmu_net());
        switch (parameters_.field_derivatives_mode) {
          case FieldDerivativesMode::ChainRule:
            (*FB_lat_)[i] = potentials_->vdf_force(
                jB.rho(), jB.drho_dxnu().x0(), jB.drho_dxnu().threevec(),
                jB.grad_rho_cross_vecj(), jB.jmu_net().x0(), jB.grad_j0(),
                jB.jmu_net().threevec(), jB.dvecj_dt(), jB.curl_vecj());
            break;
          case FieldDerivativesMode::Direct:
            auto Amu = (*fields_lat_)[i];
            (*FB_lat_)[i] = potentials_->vdf_force(
                Amu.grad_A0(), Amu.dvecA_dt(), Amu.curl_vecA());
            break;
        }
      }  // for (size_t i = 0; i < UBlattice_size; i++)
    }    // if potentials_->use_vdf()
  }
}

Here is the call graph for this function:

compute_min_cell_length()

template<typename Modus >

double smash::Experiment< Modus >::compute_min_cell_length ( double  dt ) const

inlineprivate

Calculate the minimal size for the grid cells such that the ScatterActionsFinder will find all collisions within the maximal transverse distance (which is determined by the maximal cross section).

Parameters

[in]

dt

The current time step size [fm/c]

Returns

The minimal required size of cells

Definition at line 301 of file experiment.h.

                                                  {
    if (parameters_.coll_crit == CollisionCriterion::Stochastic) {
      return parameters_.fixed_min_cell_length;
    }
    return std::sqrt(4 * dt * dt + max_transverse_distance_sqr_);
  }

next_output_time()

template<typename Modus >

double smash::Experiment< Modus >::next_output_time ( ) const

inlineprivate

Shortcut for next output time.

Definition at line 309 of file experiment.h.

                                  {
    return parameters_.outputclock->next_time();
  }

Friends And Related Function Documentation

ExperimentBase

template<typename Modus >

friend class ExperimentBase

friend

Definition at line 178 of file experiment.h.

Member Data Documentation

parameters_

template<typename Modus >

ExperimentParameters smash::Experiment< Modus >::parameters_

private

Struct of several member variables.

These variables are combined into a struct for efficient input to functions outside of this class.

Definition at line 318 of file experiment.h.

density_param_

template<typename Modus >

DensityParameters smash::Experiment< Modus >::density_param_

private

Structure to precalculate and hold parameters for density computations.

Definition at line 321 of file experiment.h.

modus_

template<typename Modus >

Modus smash::Experiment< Modus >::modus_

private

Instance of the Modus template parameter.

May store modus-specific data and contains modus-specific function implementations.

Definition at line 327 of file experiment.h.

ensembles_

template<typename Modus >

std::vector<Particles> smash::Experiment< Modus >::ensembles_

private

Complete particle list, all ensembles in one vector.

Definition at line 330 of file experiment.h.

potentials_

template<typename Modus >

std::unique_ptr<Potentials> smash::Experiment< Modus >::potentials_

private

An instance of potentials class, that stores parameters of potentials, calculates them and their gradients.

Definition at line 336 of file experiment.h.

pauli_blocker_

template<typename Modus >

std::unique_ptr<PauliBlocker> smash::Experiment< Modus >::pauli_blocker_

private

An instance of PauliBlocker class that stores parameters needed for Pauli blocking calculations and computes phase-space density.

Definition at line 342 of file experiment.h.

outputs_

template<typename Modus >

OutputsList smash::Experiment< Modus >::outputs_

private

A list of output formaters.

They will be called to write the state of the particles to file.

Definition at line 348 of file experiment.h.

dilepton_output_

template<typename Modus >

OutputPtr smash::Experiment< Modus >::dilepton_output_

private

The Dilepton output.

Definition at line 351 of file experiment.h.

photon_output_

template<typename Modus >

OutputPtr smash::Experiment< Modus >::photon_output_

private

The Photon output.

Definition at line 354 of file experiment.h.

projectile_target_interact_

template<typename Modus >

std::vector<bool> smash::Experiment< Modus >::projectile_target_interact_

private

Whether the projectile and the target collided.

One value for each ensemble.

Definition at line 360 of file experiment.h.

beam_momentum_

template<typename Modus >

std::vector<FourVector> smash::Experiment< Modus >::beam_momentum_ = {}

private

The initial nucleons in the ColliderModus propagate with beam_momentum_, if Fermi motion is frozen.

It's only valid in the ColliderModus, so is set as an empty vector by default.

Definition at line 367 of file experiment.h.

action_finders_

template<typename Modus >

std::vector<std::unique_ptr<ActionFinderInterface> > smash::Experiment< Modus >::action_finders_

private

The Action finder objects.

Definition at line 370 of file experiment.h.

dilepton_finder_

template<typename Modus >

std::unique_ptr<DecayActionsFinderDilepton> smash::Experiment< Modus >::dilepton_finder_

private

The Dilepton Action Finder.

Definition at line 373 of file experiment.h.

photon_finder_

template<typename Modus >

std::unique_ptr<ActionFinderInterface> smash::Experiment< Modus >::photon_finder_

private

The (Scatter) Actions Finder for Direct Photons.

Definition at line 376 of file experiment.h.

n_fractional_photons_

template<typename Modus >

int smash::Experiment< Modus >::n_fractional_photons_

private

Number of fractional photons produced per single reaction.

Definition at line 379 of file experiment.h.

j_QBS_lat_

template<typename Modus >

std::unique_ptr<DensityLattice> smash::Experiment< Modus >::j_QBS_lat_

private

4-current for j_QBS lattice output

Definition at line 382 of file experiment.h.

jmu_B_lat_

template<typename Modus >

std::unique_ptr<DensityLattice> smash::Experiment< Modus >::jmu_B_lat_

private

Baryon density on the lattice.

Definition at line 385 of file experiment.h.

jmu_I3_lat_

template<typename Modus >

std::unique_ptr<DensityLattice> smash::Experiment< Modus >::jmu_I3_lat_

private

Isospin projection density on the lattice.

Definition at line 388 of file experiment.h.

jmu_el_lat_

template<typename Modus >

std::unique_ptr<DensityLattice> smash::Experiment< Modus >::jmu_el_lat_

private

Electric charge density on the lattice.

Definition at line 391 of file experiment.h.

fields_lat_

template<typename Modus >

std::unique_ptr<FieldsLattice> smash::Experiment< Modus >::fields_lat_

private

Mean-field A^mu on the lattice.

Definition at line 394 of file experiment.h.

jmu_custom_lat_

template<typename Modus >

std::unique_ptr<DensityLattice> smash::Experiment< Modus >::jmu_custom_lat_

private

Custom density on the lattices.

In the config user asks for some kind of density for printout. Baryon and isospin projection density are anyway needed for potentials. If user asks for some other density type for printout, it will be handled using jmu_custom variable.

Definition at line 403 of file experiment.h.

dens_type_lattice_printout_

template<typename Modus >

DensityType smash::Experiment< Modus >::dens_type_lattice_printout_ = DensityType::None

private

Type of density for lattice printout.

Definition at line 406 of file experiment.h.

UB_lat_

template<typename Modus >

std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::UB_lat_ = nullptr

private

Lattices for Skyrme or VDF potentials (evaluated in the local rest frame) times the baryon flow 4-velocity.

Definition at line 412 of file experiment.h.

UI3_lat_

template<typename Modus >

std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::UI3_lat_ = nullptr

private

Lattices for symmetry potentials (evaluated in the local rest frame) times the isospin flow 4-velocity.

Definition at line 418 of file experiment.h.

FB_lat_

template<typename Modus >

std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector> > > smash::Experiment< Modus >::FB_lat_

private

Lattices for the electric and magnetic components of the Skyrme or VDF force.

Definition at line 425 of file experiment.h.

FI3_lat_

template<typename Modus >

std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector> > > smash::Experiment< Modus >::FI3_lat_

private

Lattices for the electric and magnetic component of the symmetry force.

Definition at line 429 of file experiment.h.

EM_lat_

template<typename Modus >

std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector> > > smash::Experiment< Modus >::EM_lat_

private

Lattices for electric and magnetic field in fm^-2.

Definition at line 433 of file experiment.h.

Tmn_

template<typename Modus >

std::unique_ptr<RectangularLattice<EnergyMomentumTensor> > smash::Experiment< Modus >::Tmn_

private

Lattices of energy-momentum tensors for printout.

Definition at line 436 of file experiment.h.

old_jmu_auxiliary_

template<typename Modus >

std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::old_jmu_auxiliary_

private

Auxiliary lattice for values of jmu at a time step t0.

Definition at line 439 of file experiment.h.

new_jmu_auxiliary_

template<typename Modus >

std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::new_jmu_auxiliary_

private

Auxiliary lattice for values of jmu at a time step t0 + dt.

Definition at line 441 of file experiment.h.

four_gradient_auxiliary_

template<typename Modus >

std::unique_ptr<RectangularLattice<std::array<FourVector, 4> > > smash::Experiment< Modus >::four_gradient_auxiliary_

private

Auxiliary lattice for calculating the four-gradient of jmu.

Definition at line 444 of file experiment.h.

old_fields_auxiliary_

template<typename Modus >

std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::old_fields_auxiliary_

private

Auxiliary lattice for values of Amu at a time step t0.

Definition at line 447 of file experiment.h.

new_fields_auxiliary_

template<typename Modus >

std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::new_fields_auxiliary_

private

Auxiliary lattice for values of Amu at a time step t0 + dt.

Definition at line 449 of file experiment.h.

fields_four_gradient_auxiliary_

template<typename Modus >

std::unique_ptr<RectangularLattice<std::array<FourVector, 4> > > smash::Experiment< Modus >::fields_four_gradient_auxiliary_

private

Auxiliary lattice for calculating the four-gradient of Amu.

Definition at line 452 of file experiment.h.

printout_tmn_

template<typename Modus >

bool smash::Experiment< Modus >::printout_tmn_ = false

private

Whether to print the energy-momentum tensor.

Definition at line 455 of file experiment.h.

printout_tmn_landau_

template<typename Modus >

bool smash::Experiment< Modus >::printout_tmn_landau_ = false

private

Whether to print the energy-momentum tensor in Landau frame.

Definition at line 458 of file experiment.h.

printout_v_landau_

template<typename Modus >

bool smash::Experiment< Modus >::printout_v_landau_ = false

private

Whether to print the 4-velocity in Landau frame.

Definition at line 461 of file experiment.h.

printout_j_QBS_

template<typename Modus >

bool smash::Experiment< Modus >::printout_j_QBS_ = false

private

Whether to print the Q, B, S 4-currents.

Definition at line 464 of file experiment.h.

printout_lattice_td_

template<typename Modus >

bool smash::Experiment< Modus >::printout_lattice_td_ = false

private

Whether to print the thermodynamics quantities evaluated on the lattices.

Definition at line 467 of file experiment.h.

printout_full_lattice_ascii_td_

template<typename Modus >

bool smash::Experiment< Modus >::printout_full_lattice_ascii_td_ = false

private

Whether to print the thermodynamics quantities evaluated on the lattices, point by point, in ASCII format.

Definition at line 471 of file experiment.h.

printout_full_lattice_binary_td_

template<typename Modus >

bool smash::Experiment< Modus >::printout_full_lattice_binary_td_ = false

private

Whether to print the thermodynamics quantities evaluated on the lattices, point by point, in Binary format.

Definition at line 475 of file experiment.h.

printout_full_lattice_any_td_

template<typename Modus >

bool smash::Experiment< Modus >::printout_full_lattice_any_td_ = false

private

Whether to print the thermodynamics quantities evaluated on the lattices, point by point, in any format.

Definition at line 479 of file experiment.h.

thermalizer_

template<typename Modus >

std::unique_ptr<GrandCanThermalizer> smash::Experiment< Modus >::thermalizer_

private

Instance of class used for forced thermalization.

Definition at line 482 of file experiment.h.

process_string_ptr_

template<typename Modus >

StringProcess* smash::Experiment< Modus >::process_string_ptr_

private

Pointer to the string process class object, which is used to set the random seed for PYTHIA objects in each event.

Definition at line 488 of file experiment.h.

nevents_

template<typename Modus >

const int smash::Experiment< Modus >::nevents_

private

Number of events.

Event is a single simulation of a physical phenomenon: elementary particle or nucleus-nucleus collision. Result of a single SMASH event is random (by construction) as well as result of one collision in nature. To compare simulation with experiment one has to take ensemble averages, i.e. perform simulation and real experiment many times and compare average results.

nevents_ is number of times single phenomenon (particle or nucleus-nucleus collision) will be simulated.

Definition at line 504 of file experiment.h.

event_

template<typename Modus >

int smash::Experiment< Modus >::event_ = 0

private

Current event.

Definition at line 507 of file experiment.h.

end_time_

template<typename Modus >

const double smash::Experiment< Modus >::end_time_

private

simulation time at which the evolution is stopped.

Definition at line 510 of file experiment.h.

delta_time_startup_

template<typename Modus >

const double smash::Experiment< Modus >::delta_time_startup_

private

The clock's timestep size at start up.

Stored here so that the next event will remember this.

Definition at line 517 of file experiment.h.

force_decays_

template<typename Modus >

const bool smash::Experiment< Modus >::force_decays_

private

This indicates whether we force all resonances to decay in the last timestep.

Definition at line 523 of file experiment.h.

use_grid_

template<typename Modus >

const bool smash::Experiment< Modus >::use_grid_

private

This indicates whether to use the grid.

Definition at line 526 of file experiment.h.

metric_

template<typename Modus >

const ExpansionProperties smash::Experiment< Modus >::metric_

private

This struct contains information on the metric to be used.

Definition at line 529 of file experiment.h.

dileptons_switch_

template<typename Modus >

const bool smash::Experiment< Modus >::dileptons_switch_

private

This indicates whether dileptons are switched on.

Definition at line 532 of file experiment.h.

photons_switch_

template<typename Modus >

const bool smash::Experiment< Modus >::photons_switch_

private

This indicates whether photons are switched on.

Definition at line 535 of file experiment.h.

bremsstrahlung_switch_

template<typename Modus >

const bool smash::Experiment< Modus >::bremsstrahlung_switch_

private

This indicates whether bremsstrahlung is switched on.

Definition at line 538 of file experiment.h.

IC_output_switch_

template<typename Modus >

const bool smash::Experiment< Modus >::IC_output_switch_

private

This indicates whether the IC output is enabled.

Definition at line 541 of file experiment.h.

time_step_mode_

template<typename Modus >

const TimeStepMode smash::Experiment< Modus >::time_step_mode_

private

This indicates whether to use time steps.

Definition at line 544 of file experiment.h.

max_transverse_distance_sqr_

template<typename Modus >

double smash::Experiment< Modus >::max_transverse_distance_sqr_ = std::numeric_limits<double>::max()

private

Maximal distance at which particles can interact in case of the geometric criterion, squared.

Definition at line 550 of file experiment.h.

conserved_initial_

template<typename Modus >

QuantumNumbers smash::Experiment< Modus >::conserved_initial_

private

The conserved quantities of the system.

This struct carries the sums of the single particle's various quantities as measured at the beginning of the evolution and can be used to regularly check if they are still good.

Definition at line 559 of file experiment.h.

initial_mean_field_energy_

template<typename Modus >

double smash::Experiment< Modus >::initial_mean_field_energy_

private

The initial total mean field energy in the system.

Note: will only be calculated if lattice is on.

Definition at line 565 of file experiment.h.

time_start_

template<typename Modus >

SystemTimePoint smash::Experiment< Modus >::time_start_ = SystemClock::now()

private

system starting time of the simulation

Definition at line 568 of file experiment.h.

dens_type_

template<typename Modus >

DensityType smash::Experiment< Modus >::dens_type_ = DensityType::None

private

Type of density to be written to collision headers.

Definition at line 571 of file experiment.h.

interactions_total_

template<typename Modus >

uint64_t smash::Experiment< Modus >::interactions_total_ = 0

private

Total number of interactions for current timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 577 of file experiment.h.

previous_interactions_total_

template<typename Modus >

uint64_t smash::Experiment< Modus >::previous_interactions_total_ = 0

private

Total number of interactions for previous timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 583 of file experiment.h.

wall_actions_total_

template<typename Modus >

uint64_t smash::Experiment< Modus >::wall_actions_total_ = 0

private

Total number of wall-crossings for current timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 589 of file experiment.h.

previous_wall_actions_total_

template<typename Modus >

uint64_t smash::Experiment< Modus >::previous_wall_actions_total_ = 0

private

Total number of wall-crossings for previous timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 595 of file experiment.h.

total_pauli_blocked_

template<typename Modus >

uint64_t smash::Experiment< Modus >::total_pauli_blocked_ = 0

private

Total number of Pauli-blockings for current timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 601 of file experiment.h.

total_hypersurface_crossing_actions_

template<typename Modus >

uint64_t smash::Experiment< Modus >::total_hypersurface_crossing_actions_ = 0

private

Total number of particles removed from the evolution in hypersurface crossing actions.

Definition at line 607 of file experiment.h.

discarded_interactions_total_

template<typename Modus >

uint64_t smash::Experiment< Modus >::discarded_interactions_total_ = 0

private

Total number of discarded interactions, because they were invalidated before they could be performed.

Definition at line 613 of file experiment.h.

total_energy_removed_

template<typename Modus >

double smash::Experiment< Modus >::total_energy_removed_ = 0.0

private

Total energy removed from the system in hypersurface crossing actions.

Definition at line 619 of file experiment.h.

seed_

template<typename Modus >

int64_t smash::Experiment< Modus >::seed_ = -1

private

random seed for the next event.

Definition at line 622 of file experiment.h.

---

The documentation for this class was generated from the following file:

• src/include/smash/experiment.h