7 #ifndef SRC_INCLUDE_SMASH_EXPERIMENT_H_
8 #define SRC_INCLUDE_SMASH_EXPERIMENT_H_
40 #ifdef SMASH_USE_HEPMC
63 template <
typename T,
typename Ratio>
65 const chrono::duration<T, Ratio> &seconds) {
66 using Seconds = chrono::duration<double>;
67 using Minutes = chrono::duration<double, std::ratio<60>>;
68 using Hours = chrono::duration<double, std::ratio<60 * 60>>;
69 constexpr Minutes threshold_for_minutes{10};
70 constexpr Hours threshold_for_hours{3};
71 if (seconds < threshold_for_minutes) {
72 return out << Seconds(seconds).count() <<
" [s]";
74 if (seconds < threshold_for_hours) {
75 return out << Minutes(seconds).count() <<
" [min]";
77 return out << Hours(seconds).count() <<
" [h]";
82 static constexpr
int LMain = LogArea::Main::id;
122 const bf::path &output_path);
130 virtual void run() = 0;
138 using std::invalid_argument::invalid_argument;
142 template <
typename Modus>
144 template <
typename Modus>
171 template <
typename Modus>
248 template <
typename Container>
250 const Container &particles_before_actions);
402 std::unique_ptr<RectangularLattice<FourVector>>
UB_lat_ =
nullptr;
408 std::unique_ptr<RectangularLattice<FourVector>>
UI3_lat_ =
nullptr;
411 std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector>>>
415 std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector>>>
419 std::unique_ptr<RectangularLattice<EnergyMomentumTensor>>
Tmn_;
575 friend std::ostream &operator<<<>(std::ostream &out,
const Experiment &e);
579 template <
typename Modus>
581 out <<
"End time: " << e.
end_time_ <<
" fm/c\n";
586 template <
typename Modus>
588 const std::string &content,
589 const bf::path &output_path,
591 logg[
LExperiment].info() <<
"Adding output " << content <<
" of format "
594 if (
format ==
"VTK" && content ==
"Particles") {
595 outputs_.emplace_back(
596 make_unique<VtkOutput>(output_path, content, out_par));
597 }
else if (
format ==
"Root") {
598 #ifdef SMASH_USE_ROOT
599 if (content ==
"Initial_Conditions") {
600 outputs_.emplace_back(
601 make_unique<RootOutput>(output_path,
"SMASH_IC", out_par));
603 outputs_.emplace_back(
604 make_unique<RootOutput>(output_path, content, out_par));
608 "Root output requested, but Root support not compiled in");
610 }
else if (
format ==
"Binary") {
611 if (content ==
"Collisions" || content ==
"Dileptons" ||
612 content ==
"Photons") {
613 outputs_.emplace_back(
614 make_unique<BinaryOutputCollisions>(output_path, content, out_par));
615 }
else if (content ==
"Particles") {
616 outputs_.emplace_back(
617 make_unique<BinaryOutputParticles>(output_path, content, out_par));
618 }
else if (content ==
"Initial_Conditions") {
619 outputs_.emplace_back(make_unique<BinaryOutputInitialConditions>(
620 output_path, content, out_par));
622 }
else if (
format ==
"Oscar1999" ||
format ==
"Oscar2013") {
623 outputs_.emplace_back(
625 }
else if (content ==
"Thermodynamics" &&
format ==
"ASCII") {
626 outputs_.emplace_back(
627 make_unique<ThermodynamicOutput>(output_path, content, out_par));
628 }
else if (content ==
"Thermodynamics" &&
format ==
"VTK") {
629 printout_lattice_td_ =
true;
630 outputs_.emplace_back(
631 make_unique<VtkOutput>(output_path, content, out_par));
632 }
else if (content ==
"Initial_Conditions" &&
format ==
"ASCII") {
633 outputs_.emplace_back(
634 make_unique<ICOutput>(output_path,
"SMASH_IC", out_par));
635 }
else if (content ==
"HepMC" &&
format ==
"ASCII") {
636 #ifdef SMASH_USE_HEPMC
637 outputs_.emplace_back(make_unique<HepMcOutput>(
638 output_path,
"SMASH_HepMC", out_par, modus_.total_N_number(),
639 modus_.proj_N_number()));
642 "HepMC output requested, but HepMC support not compiled in");
646 <<
"Unknown combination of format (" <<
format <<
") and content ("
647 << content <<
"). Fix the config.";
821 template <
typename Modus>
825 modus_(config[
"Modi"], parameters_),
827 nevents_(config.take({
"General",
"Nevents"})),
828 end_time_(config.take({
"General",
"End_Time"})),
829 delta_time_startup_(parameters_.labclock->timestep_duration()),
831 config.take({
"Collision_Term",
"Force_Decays_At_End"},
true)),
832 use_grid_(config.take({
"General",
"Use_Grid"},
true)),
835 config.take({
"General",
"Expansion_Rate"}, 0.1)),
837 config.take({
"Collision_Term",
"Dileptons",
"Decays"},
false)),
838 photons_switch_(config.take(
839 {
"Collision_Term",
"Photons",
"2to2_Scatterings"},
false)),
840 bremsstrahlung_switch_(
841 config.take({
"Collision_Term",
"Photons",
"Bremsstrahlung"},
false)),
842 IC_output_switch_(config.has_value({
"Output",
"Initial_Conditions"})),
849 throw std::invalid_argument(
850 "The stochastic criterion can only be employed for fixed time step "
855 if (dileptons_switch_) {
856 dilepton_finder_ = make_unique<DecayActionsFinderDilepton>();
858 if (photons_switch_ || bremsstrahlung_switch_) {
859 n_fractional_photons_ =
860 config.take({
"Collision_Term",
"Photons",
"Fractional_Photons"}, 100);
862 if (parameters_.two_to_one) {
863 if (parameters_.res_lifetime_factor < 0.) {
864 throw std::invalid_argument(
865 "Resonance lifetime modifier cannot be negative!");
869 "Resonance lifetime set to zero. Make sure resonances cannot "
871 "inelastically (e.g. resonance chains), else SMASH is known to "
874 action_finders_.emplace_back(
875 make_unique<DecayActionsFinder>(parameters_.res_lifetime_factor));
877 bool no_coll = config.take({
"Collision_Term",
"No_Collisions"},
false);
878 if ((parameters_.two_to_one || parameters_.included_2to2.any() ||
879 parameters_.included_multi.any() || parameters_.strings_switch) &&
881 auto scat_finder = make_unique<ScatterActionsFinder>(
882 config, parameters_, nucleon_has_interacted_, modus_.total_N_number(),
883 modus_.proj_N_number());
884 max_transverse_distance_sqr_ =
885 scat_finder->max_transverse_distance_sqr(parameters_.testparticles);
886 process_string_ptr_ = scat_finder->get_process_string_ptr();
887 action_finders_.emplace_back(std::move(scat_finder));
889 max_transverse_distance_sqr_ =
890 parameters_.maximum_cross_section / M_PI *
fm2_mb;
891 process_string_ptr_ = NULL;
893 if (modus_.is_box()) {
894 action_finders_.emplace_back(
895 make_unique<WallCrossActionsFinder>(parameters_.box_length));
897 if (IC_output_switch_) {
898 if (!modus_.is_collider()) {
899 throw std::runtime_error(
900 "Initial conditions can only be extracted in collider modus.");
903 if (config.has_value({
"Output",
"Initial_Conditions",
"Proper_Time"})) {
906 config.take({
"Output",
"Initial_Conditions",
"Proper_Time"});
909 double default_proper_time = modus_.nuclei_passing_time();
911 config.take({
"Output",
"Initial_Conditions",
"Lower_Bound"}, 0.5);
912 if (default_proper_time >= lower_bound) {
913 proper_time = default_proper_time;
916 <<
"Nuclei passing time is too short, hypersurface proper time set "
918 << lower_bound <<
" fm.";
919 proper_time = lower_bound;
922 action_finders_.emplace_back(
923 make_unique<HyperSurfaceCrossActionsFinder>(proper_time));
926 if (config.has_value({
"Collision_Term",
"Pauli_Blocking"})) {
928 pauli_blocker_ = make_unique<PauliBlocker>(
929 config[
"Collision_Term"][
"Pauli_Blocking"], parameters_);
933 config.take({
"Collision_Term",
"Power_Particle_Formation"},
934 modus_.sqrt_s_NN() >= 200. ? -1. : 1.);
971 auto output_conf = config[
"Output"];
1198 <<
"Density type printed to headers: " << dens_type_;
1200 const OutputParameters output_parameters(std::move(output_conf));
1202 std::vector<std::string> output_contents = output_conf.list_upmost_nodes();
1203 for (
const auto &content : output_contents) {
1204 auto this_output_conf = output_conf[content.c_str()];
1205 const std::vector<std::string> formats = this_output_conf.take({
"Format"});
1206 if (output_path ==
"") {
1209 for (
const auto &
format : formats) {
1210 create_output(
format, content, output_path, output_parameters);
1221 if (config.has_value({
"Potentials"})) {
1224 throw std::invalid_argument(
"Can't use potentials without time steps!");
1228 <<
"Potentials don't work with frozen Fermi momenta! "
1229 "Use normal Fermi motion instead.";
1230 throw std::invalid_argument(
1231 "Can't use potentials "
1232 "with frozen Fermi momenta!");
1235 << parameters_.labclock->timestep_duration();
1237 potentials_ = make_unique<Potentials>(config[
"Potentials"], parameters_);
1296 if (config.has_value({
"Lattice"})) {
1298 const std::array<double, 3> l = config.take({
"Lattice",
"Sizes"});
1299 const std::array<int, 3>
n = config.take({
"Lattice",
"Cell_Number"});
1300 const std::array<double, 3> origin = config.take({
"Lattice",
"Origin"});
1301 const bool periodic = config.take({
"Lattice",
"Periodic"});
1303 if (printout_lattice_td_) {
1304 dens_type_lattice_printout_ = output_parameters.td_dens_type;
1305 printout_tmn_ = output_parameters.td_tmn;
1306 printout_tmn_landau_ = output_parameters.td_tmn_landau;
1307 printout_v_landau_ = output_parameters.td_v_landau;
1309 if (printout_tmn_ || printout_tmn_landau_ || printout_v_landau_) {
1310 Tmn_ = make_unique<RectangularLattice<EnergyMomentumTensor>>(
1317 if (potentials_->use_skyrme()) {
1318 jmu_B_lat_ = make_unique<DensityLattice>(l,
n, origin, periodic,
1320 UB_lat_ = make_unique<RectangularLattice<FourVector>>(
1323 RectangularLattice<std::pair<ThreeVector, ThreeVector>>>(
1326 if (potentials_->use_symmetry()) {
1327 jmu_I3_lat_ = make_unique<DensityLattice>(l,
n, origin, periodic,
1329 UI3_lat_ = make_unique<RectangularLattice<FourVector>>(
1332 RectangularLattice<std::pair<ThreeVector, ThreeVector>>>(
1337 jmu_B_lat_ = make_unique<DensityLattice>(l,
n, origin, periodic,
1341 jmu_I3_lat_ = make_unique<DensityLattice>(l,
n, origin, periodic,
1348 jmu_custom_lat_ = make_unique<DensityLattice>(l,
n, origin, periodic,
1351 }
else if (printout_lattice_td_) {
1353 "If you want Thermodynamic VTK output, configure a lattice for it.");
1356 if ((potentials_ !=
nullptr) && (jmu_B_lat_ ==
nullptr)) {
1358 <<
"not going to be calculated.";
1362 if (parameters_.potential_affect_threshold) {
1369 if (config.has_value({
"Forced_Thermalization"})) {
1370 Configuration &&th_conf = config[
"Forced_Thermalization"];
1371 thermalizer_ = modus_.create_grandcan_thermalizer(th_conf);
1376 seed_ = config.take({
"General",
"Randomseed"});
1380 const std::string
hline(113,
'-');
1408 uint64_t scatterings_this_interval,
1409 const QuantumNumbers &conserved_initial,
1411 double E_mean_field,
1412 double E_mean_field_initial);
1426 const Potentials &potentials,
1427 RectangularLattice<smash::DensityOnLattice> &jmu_B_lat,
1428 const ExperimentParameters ¶meters);
1445 EventInfo
fill_event_info(
const Particles &particles,
double E_mean_field,
1446 double modus_impact_parameter,
1447 const ExperimentParameters ¶meters,
1448 bool projectile_target_interact);
1450 template <
typename Modus>
1460 while (r == INT64_MIN) {
1463 seed_ = std::abs(r);
1468 if (process_string_ptr_ != NULL) {
1469 process_string_ptr_->init_pythia_hadron_rndm();
1475 double start_time = modus_.initial_conditions(&particles_, parameters_);
1479 modus_.impose_boundary_conditions(&particles_, outputs_);
1481 double timestep = delta_time_startup_;
1483 switch (time_step_mode_) {
1487 timestep = end_time_ - start_time;
1489 const double max_dt = modus_.max_timestep(max_transverse_distance_sqr_);
1490 if (max_dt > 0. && max_dt < timestep) {
1495 std::unique_ptr<UniformClock> clock_for_this_event;
1496 if (modus_.is_list() && (timestep < 0.0)) {
1497 throw std::runtime_error(
1498 "Timestep for the given event is negative. \n"
1499 "This might happen if the formation times of the input particles are "
1500 "larger than the specified end time of the simulation.");
1502 clock_for_this_event = make_unique<UniformClock>(start_time, timestep);
1503 parameters_.labclock = std::move(clock_for_this_event);
1506 parameters_.outputclock->reset(start_time,
true);
1508 parameters_.outputclock->remove_times_in_past(start_time);
1511 "Lab clock: t_start = ", parameters_.labclock->current_time(),
1512 ", dt = ", parameters_.labclock->timestep_duration());
1517 wall_actions_total_ = 0;
1518 previous_wall_actions_total_ = 0;
1519 interactions_total_ = 0;
1520 previous_interactions_total_ = 0;
1521 discarded_interactions_total_ = 0;
1522 total_pauli_blocked_ = 0;
1523 projectile_target_interact_ =
false;
1524 total_hypersurface_crossing_actions_ = 0;
1525 total_energy_removed_ = 0.0;
1528 logg[
LExperiment].info() <<
"Time[fm] Ekin[GeV] E_MF[GeV] ETotal[GeV] "
1529 <<
"ETot/N[GeV] D(ETot/N)[GeV] Scatt&Decays "
1530 <<
"Particles Comp.Time";
1532 double E_mean_field = 0.0;
1534 update_potentials();
1536 if ((jmu_B_lat_ !=
nullptr)) {
1541 initial_mean_field_energy_ = E_mean_field;
1543 particles_, 0u, conserved_initial_, time_start_,
1544 parameters_.labclock->current_time(), E_mean_field,
1545 initial_mean_field_energy_);
1549 parameters_, projectile_target_interact_);
1552 for (
const auto &output : outputs_) {
1553 output->at_eventstart(particles_, event_number, event_info);
1557 template <
typename Modus>
1558 template <
typename Container>
1560 Action &action,
const Container &particles_before_actions) {
1562 if (!action.
is_valid(particles_)) {
1563 discarded_interactions_total_++;
1565 " (discarded: invalid)");
1570 if (pauli_blocker_ && action.
is_pauli_blocked(particles_, *pauli_blocker_)) {
1571 total_pauli_blocked_++;
1574 if (modus_.is_collider()) {
1577 bool incoming_projectile =
false;
1578 bool incoming_target =
false;
1580 assert(incoming.id() >= 0);
1581 if (incoming.id() < modus_.total_N_number()) {
1582 nucleon_has_interacted_[incoming.id()] =
true;
1584 if (incoming.id() < modus_.proj_N_number()) {
1585 incoming_projectile =
true;
1587 if (incoming.id() >= modus_.proj_N_number() &&
1588 incoming.id() < modus_.total_N_number()) {
1589 incoming_target =
true;
1593 if (incoming_projectile & incoming_target) {
1594 projectile_target_interact_ =
true;
1599 const auto id_process = static_cast<uint32_t>(interactions_total_ + 1);
1600 action.
perform(&particles_, id_process);
1601 interactions_total_++;
1603 wall_actions_total_++;
1606 total_hypersurface_crossing_actions_++;
1613 constexpr
bool compute_grad =
false;
1614 const bool smearing =
true;
1616 particles_before_actions, density_param_,
1617 dens_type_, compute_grad, smearing));
1634 for (
const auto &output : outputs_) {
1635 if (!output->is_dilepton_output() && !output->is_photon_output()) {
1636 if (output->is_IC_output() &&
1638 output->at_interaction(action, rho);
1639 }
else if (!output->is_IC_output()) {
1640 output->at_interaction(action, rho);
1650 if (photons_switch_ &&
1656 constexpr
double action_time = 0.;
1658 n_fractional_photons_,
1674 photon_act.add_single_process();
1676 photon_act.perform_photons(outputs_);
1679 if (bremsstrahlung_switch_ &&
1684 constexpr
double action_time = 0.;
1687 n_fractional_photons_,
1704 brems_act.add_single_process();
1706 brems_act.perform_bremsstrahlung(outputs_);
1713 template <
typename Modus>
1717 while (parameters_.labclock->current_time() < end_time_) {
1718 const double t = parameters_.labclock->current_time();
1720 std::min(parameters_.labclock->timestep_duration(), end_time_ - t);
1721 logg[
LExperiment].debug(
"Timestepless propagation for next ", dt,
" fm/c.");
1725 thermalizer_->is_time_to_thermalize(parameters_.labclock)) {
1726 const bool ignore_cells_under_treshold =
true;
1727 thermalizer_->update_thermalizer_lattice(particles_, density_param_,
1728 ignore_cells_under_treshold);
1729 const double current_t = parameters_.labclock->current_time();
1730 thermalizer_->thermalize(particles_, current_t,
1731 parameters_.testparticles);
1734 perform_action(th_act, particles_);
1738 if (particles_.size() > 0 && action_finders_.size() > 0) {
1740 double min_cell_length = compute_min_cell_length(dt);
1744 use_grid_ ? modus_.create_grid(particles_, min_cell_length, dt)
1745 : modus_.create_grid(particles_, min_cell_length, dt,
1748 const double gcell_vol = grid.cell_volume();
1752 [&](
const ParticleList &search_list) {
1753 for (
const auto &finder : action_finders_) {
1754 actions.
insert(finder->find_actions_in_cell(
1755 search_list, dt, gcell_vol, beam_momentum_));
1758 [&](
const ParticleList &search_list,
1759 const ParticleList &neighbors_list) {
1760 for (
const auto &finder : action_finders_) {
1761 actions.
insert(finder->find_actions_with_neighbors(
1762 search_list, neighbors_list, dt, beam_momentum_));
1770 run_time_evolution_timestepless(actions);
1775 update_potentials();
1776 update_momenta(&particles_, parameters_.labclock->timestep_duration(),
1777 *potentials_, FB_lat_.get(), FI3_lat_.get());
1786 ++(*parameters_.labclock);
1794 if (!potentials_ && !parameters_.strings_switch &&
1796 std::string err_msg = conserved_initial_.report_deviations(particles_);
1797 if (!err_msg.empty()) {
1799 throw std::runtime_error(
"Violation of conserved quantities!");
1804 if (pauli_blocker_) {
1806 "Interactions: Pauli-blocked/performed = ", total_pauli_blocked_,
"/",
1807 interactions_total_ - wall_actions_total_);
1811 template <
typename Modus>
1815 if (dilepton_finder_ !=
nullptr) {
1816 for (
const auto &output : outputs_) {
1817 dilepton_finder_->shine(particles_, output.get(), dt);
1830 constexpr uint64_t max_uint32 = std::numeric_limits<uint32_t>::max();
1831 if (interactions_total >= max_uint32) {
1832 throw std::runtime_error(
"Integer overflow in total interaction number!");
1836 template <
typename Modus>
1838 const double start_time = parameters_.labclock->current_time();
1839 const double end_time =
1840 std::min(parameters_.labclock->next_time(), end_time_);
1841 double time_left = end_time - start_time;
1843 "Timestepless propagation: ",
"Actions size = ", actions.
size(),
1844 ", start time = ", start_time,
", end time = ", end_time);
1849 ActionPtr act = actions.
pop();
1850 if (!act->is_valid(particles_)) {
1851 discarded_interactions_total_++;
1853 " (discarded: invalid)");
1856 if (act->time_of_execution() > end_time) {
1858 act,
" scheduled later than end time: t_action[fm/c] = ",
1859 act->time_of_execution(),
", t_end[fm/c] = ", end_time);
1863 while (next_output_time() <= act->time_of_execution()) {
1865 next_output_time());
1866 propagate_and_shine(next_output_time());
1867 ++(*parameters_.outputclock);
1868 intermediate_output();
1873 ", action time = ", act->time_of_execution());
1874 propagate_and_shine(act->time_of_execution());
1881 act->update_incoming(particles_);
1882 const bool performed = perform_action(*act, particles_);
1892 time_left = end_time - act->time_of_execution();
1893 const ParticleList &outgoing_particles = act->outgoing_particles();
1895 const double gcell_vol = 0.0;
1896 for (
const auto &finder : action_finders_) {
1898 actions.
insert(finder->find_actions_in_cell(outgoing_particles, time_left,
1899 gcell_vol, beam_momentum_));
1901 actions.
insert(finder->find_actions_with_surrounding_particles(
1902 outgoing_particles, particles_, time_left, beam_momentum_));
1908 while (next_output_time() <= end_time) {
1910 next_output_time());
1911 propagate_and_shine(next_output_time());
1912 ++(*parameters_.outputclock);
1914 if (parameters_.outputclock->current_time() < end_time_) {
1915 intermediate_output();
1919 propagate_and_shine(end_time);
1922 template <
typename Modus>
1924 const uint64_t wall_actions_this_interval =
1925 wall_actions_total_ - previous_wall_actions_total_;
1926 previous_wall_actions_total_ = wall_actions_total_;
1927 const uint64_t interactions_this_interval = interactions_total_ -
1928 previous_interactions_total_ -
1929 wall_actions_this_interval;
1930 previous_interactions_total_ = interactions_total_;
1931 double E_mean_field = 0.0;
1934 if ((jmu_B_lat_ !=
nullptr)) {
1943 if (modus_.is_box()) {
1944 double tmp = (E_mean_field - initial_mean_field_energy_) /
1945 (E_mean_field + initial_mean_field_energy_);
1951 if (std::abs(tmp) > 0.01) {
1953 <<
"\n\n\n\t The mean field at t = "
1954 << parameters_.outputclock->current_time()
1955 <<
" [fm/c] differs from the mean field at t = 0:"
1956 <<
"\n\t\t initial_mean_field_energy_ = "
1957 << initial_mean_field_energy_ <<
" [GeV]"
1958 <<
"\n\t\t abs[(E_MF - E_MF(t=0))/(E_MF + E_MF(t=0))] = "
1960 <<
"\n\t\t E_MF/E_MF(t=0) = "
1961 << E_mean_field / initial_mean_field_energy_ <<
"\n\n";
1968 particles_, interactions_this_interval, conserved_initial_, time_start_,
1969 parameters_.outputclock->current_time(), E_mean_field,
1970 initial_mean_field_energy_);
1975 parameters_, projectile_target_interact_);
1977 if (!(modus_.is_box() && parameters_.outputclock->current_time() <
1978 modus_.equilibration_time())) {
1979 for (
const auto &output : outputs_) {
1980 if (output->is_dilepton_output() || output->is_photon_output() ||
1981 output->is_IC_output()) {
1985 output->at_intermediate_time(particles_, parameters_.outputclock,
1986 density_param_, event_info);
1989 switch (dens_type_lattice_printout_) {
1992 density_param_, particles_,
false);
2008 dens_type_lattice_printout_, density_param_,
2011 dens_type_lattice_printout_,
2014 if (printout_tmn_ || printout_tmn_landau_ || printout_v_landau_) {
2016 density_param_, particles_);
2017 if (printout_tmn_) {
2019 dens_type_lattice_printout_, *Tmn_);
2021 if (printout_tmn_landau_) {
2023 dens_type_lattice_printout_, *Tmn_);
2025 if (printout_v_landau_) {
2027 dens_type_lattice_printout_, *Tmn_);
2032 output->thermodynamics_output(*thermalizer_);
2038 template <
typename Modus>
2041 if (potentials_->use_symmetry() && jmu_I3_lat_ !=
nullptr) {
2046 if ((potentials_->use_skyrme() || potentials_->use_symmetry()) &&
2047 jmu_B_lat_ !=
nullptr) {
2050 const size_t UBlattice_size = UB_lat_->size();
2051 for (
size_t i = 0; i < UBlattice_size; i++) {
2052 auto jB = (*jmu_B_lat_)[i];
2054 std::abs(jB.density()) >
really_small ? jB.jmu_net() / jB.density()
2056 double baryon_density = jB.density();
2060 if (potentials_->use_skyrme()) {
2062 flow_four_velocity_B * potentials_->skyrme_pot(baryon_density);
2063 (*FB_lat_)[i] = potentials_->skyrme_force(
2064 baryon_density, baryon_grad_rho, baryon_dj_dt, baryon_rot_j);
2066 if (potentials_->use_symmetry() && jmu_I3_lat_ !=
nullptr) {
2067 auto jI3 = (*jmu_I3_lat_)[i];
2070 ? jI3.jmu_net() / jI3.density()
2073 flow_four_velocity_I3 *
2074 potentials_->symmetry_pot(jI3.density(), baryon_density);
2075 (*FI3_lat_)[i] = potentials_->symmetry_force(
2076 jI3.density(), jI3.grad_rho(), jI3.dj_dt(), jI3.rot_j(),
2077 baryon_density, baryon_grad_rho, baryon_dj_dt, baryon_rot_j);
2084 template <
typename Modus>
2088 uint64_t interactions_old;
2089 const auto particles_before_actions = particles_.copy_to_vector();
2093 interactions_old = interactions_total_;
2096 if (dilepton_finder_ !=
nullptr) {
2097 for (
const auto &output : outputs_) {
2098 dilepton_finder_->shine_final(particles_, output.get(),
true);
2102 for (
const auto &finder : action_finders_) {
2103 actions.
insert(finder->find_final_actions(particles_));
2107 perform_action(*actions.
pop(), particles_before_actions);
2110 }
while (interactions_total_ > interactions_old);
2113 if (dilepton_finder_ !=
nullptr) {
2114 for (
const auto &output : outputs_) {
2115 dilepton_finder_->shine_final(particles_, output.get(),
false);
2120 template <
typename Modus>
2125 double E_mean_field = 0.0;
2126 if (
likely(parameters_.labclock > 0)) {
2127 const uint64_t wall_actions_this_interval =
2128 wall_actions_total_ - previous_wall_actions_total_;
2129 const uint64_t interactions_this_interval = interactions_total_ -
2130 previous_interactions_total_ -
2131 wall_actions_this_interval;
2134 if ((jmu_B_lat_ !=
nullptr)) {
2140 particles_, interactions_this_interval, conserved_initial_, time_start_,
2141 end_time_, E_mean_field, initial_mean_field_energy_);
2142 if (IC_output_switch_ && (particles_.size() == 0)) {
2145 const double remaining_energy =
2146 conserved_initial_.momentum().x0() - total_energy_removed_;
2148 throw std::runtime_error(
2149 "There is remaining energy in the system although all particles "
2152 std::to_string(remaining_energy) +
" [GeV]");
2156 <<
"Time real: " << SystemClock::now() - time_start_;
2158 <<
"Interactions before reaching hypersurface: "
2159 << interactions_total_ - wall_actions_total_ -
2160 total_hypersurface_crossing_actions_;
2162 <<
"Total number of particles removed on hypersurface: "
2163 << total_hypersurface_crossing_actions_;
2166 const double precent_discarded =
2167 interactions_total_ > 0
2168 ? static_cast<double>(discarded_interactions_total_) * 100.0 /
2171 std::stringstream msg_discarded;
2173 <<
"Discarded interaction number: " << discarded_interactions_total_
2174 <<
" (" << precent_discarded
2175 <<
"% of the total interaction number including wall crossings)";
2179 <<
"Time real: " << SystemClock::now() - time_start_;
2183 precent_discarded > 1.0) {
2186 << msg_discarded.str()
2187 <<
"\nThe number of discarded interactions is large, which means "
2188 "the assumption for the stochastic criterion of\n1 interaction"
2189 "per particle per timestep is probably violated. Consider "
2190 "reducing the timestep size.";
2194 << interactions_total_ - wall_actions_total_;
2198 int unformed_particles_count = 0;
2199 for (
const auto &particle : particles_) {
2200 if (particle.formation_time() > end_time_) {
2201 unformed_particles_count++;
2204 if (unformed_particles_count > 0) {
2206 "End time might be too small. ", unformed_particles_count,
2207 " unformed particles were found at the end of the evolution.");
2213 parameters_, projectile_target_interact_);
2215 for (
const auto &output : outputs_) {
2216 output->at_eventend(particles_, evt_num, event_info);
2220 template <
typename Modus>
2223 for (
int j = 0; j < nevents_; j++) {
2224 mainlog.info() <<
"Event " << j;
2227 initialize_new_event(j);
2235 if (modus_.is_collider()) {
2236 if (!modus_.cll_in_nucleus()) {
2237 nucleon_has_interacted_.assign(modus_.total_N_number(),
false);
2239 nucleon_has_interacted_.assign(modus_.total_N_number(),
true);
2245 for (
int i = 0; i < modus_.total_N_number(); i++) {
2246 const auto mass_beam = particles_.copy_to_vector()[i].effective_mass();
2247 const auto v_beam = i < modus_.proj_N_number()
2248 ? modus_.velocity_projectile()
2249 : modus_.velocity_target();
2250 const auto gamma = 1.0 / std::sqrt(1.0 - v_beam * v_beam);
2251 beam_momentum_.emplace_back(
FourVector(gamma * mass_beam, 0.0, 0.0,
2252 gamma * v_beam * mass_beam));
2256 run_time_evolution();
2258 if (force_decays_) {
2269 #endif // SRC_INCLUDE_SMASH_EXPERIMENT_H_