diff --git a/include/openmc/event.h b/include/openmc/event.h index 5b60922a4..7428eedd5 100644 --- a/include/openmc/event.h +++ b/include/openmc/event.h @@ -31,8 +31,9 @@ struct EventQueueItem{ // Constructors EventQueueItem() = default; - EventQueueItem(const Particle& p, int64_t buffer_idx) : - idx(buffer_idx), type(p.type_), material(p.material_), E(p.E_) {} + EventQueueItem(const Particle& p, int64_t buffer_idx) + : idx(buffer_idx), type(p.type()), material(p.material()), E(p.E()) + {} // Compare by particle type, then by material type (4.5% fuel/7.0% fuel/cladding/etc), // then by energy. diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 0aa3f4ac5..4b7c8b991 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -187,25 +187,6 @@ public: Particle(); - //========================================================================== - // Methods and accessors - - // Accessors for position in global coordinates - Position& r() { return coord_[0].r; } - const Position& r() const { return coord_[0].r; } - - // Accessors for position in local coordinates - Position& r_local() { return coord_[n_coord_ - 1].r; } - const Position& r_local() const { return coord_[n_coord_ - 1].r; } - - // Accessors for direction in global coordinates - Direction& u() { return coord_[0].u; } - const Direction& u() const { return coord_[0].u; } - - // Accessors for direction in local coordinates - Direction& u_local() { return coord_[n_coord_ - 1].u; } - const Direction& u_local() const { return coord_[n_coord_ - 1].u; } - //! resets all coordinate levels for the particle void clear(); @@ -279,8 +260,16 @@ public: uint64_t* current_seed() {return seeds_ + stream_;} const uint64_t* current_seed() const {return seeds_ + stream_;} + //! Force recalculation of neutron xs by setting last energy to zero + void invalidate_neutron_xs() + { + for (auto& micro : neutron_xs_) + micro.last_E = 0.0; + } + +private: //========================================================================== - // Data members + // Data members (accessor methods are below) // Cross section caches std::vector neutron_xs_; //!< Microscopic neutron cross sections @@ -387,6 +376,138 @@ public: #endif int64_t n_progeny_ {0}; // Number of progeny produced by this particle + +public: + //========================================================================== + // Methods and accessors + + NuclideMicroXS& neutron_xs(const int& i) { return neutron_xs_[i]; } + const NuclideMicroXS& neutron_xs(const int& i) const + { + return neutron_xs_[i]; + } + ElementMicroXS& photon_xs(const int& i) { return photon_xs_[i]; } + MacroXS& macro_xs() { return macro_xs_; } + const MacroXS& macro_xs() const { return macro_xs_; } + + int64_t& id() { return id_; } + Type& type() { return type_; } + const Type& type() const { return type_; } + + int& n_coord() { return n_coord_; } + const int& n_coord() const { return n_coord_; } + int& cell_instance() { return cell_instance_; } + const int& cell_instance() const { return cell_instance_; } + LocalCoord& coord(const int& i) { return coord_[i]; } + const LocalCoord& coord(const int& i) const { return coord_[i]; } + + int& n_coord_last() { return n_coord_last_; } + const int& n_coord_last() const { return n_coord_last_; } + int& cell_last(const int& i) { return cell_last_[i]; } + const int& cell_last(const int& i) const { return cell_last_[i]; } + + double& E() { return E_; } + const double& E() const { return E_; } + double& E_last() { return E_last_; } + const double& E_last() const { return E_last_; } + int& g() { return g_; } + const int& g() const { return g_; } + int& g_last() { return g_last_; } + const int& g_last() const { return g_last_; } + + double& wgt() { return wgt_; } + double& mu() { return mu_; } + const double& mu() const { return mu_; } + bool& alive() { return alive_; } + + Position& r_last_current() { return r_last_current_; } + const Position& r_last_current() const { return r_last_current_; } + Position& r_last() { return r_last_; } + const Position& r_last() const { return r_last_; } + Position& u_last() { return u_last_; } + const Position& u_last() const { return u_last_; } + double& wgt_last() { return wgt_last_; } + const double& wgt_last() const { return wgt_last_; } + double& wgt_absorb() { return wgt_absorb_; } + const double& wgt_absorb() const { return wgt_absorb_; } + + bool& fission() { return fission_; } + TallyEvent& event() { return event_; } + const TallyEvent& event() const { return event_; } + int& event_nuclide() { return event_nuclide_; } + const int& event_nuclide() const { return event_nuclide_; } + int& event_mt() { return event_mt_; } + int& delayed_group() { return delayed_group_; } + + int& n_bank() { return n_bank_; } + int& n_bank_second() { return n_bank_second_; } + double& wgt_bank() { return wgt_bank_; } + int* n_delayed_bank() { return n_delayed_bank_; } + int& n_delayed_bank(const int& i) { return n_delayed_bank_[i]; } + + int& surface() { return surface_; } + const int& surface() const { return surface_; } + int& cell_born() { return cell_born_; } + const int& cell_born() const { return cell_born_; } + int& material() { return material_; } + const int& material() const { return material_; } + int& material_last() { return material_last_; } + + BoundaryInfo& boundary() { return boundary_; } + + double& sqrtkT() { return sqrtkT_; } + const double& sqrtkT() const { return sqrtkT_; } + double& sqrtkT_last() { return sqrtkT_last_; } + + int& n_collision() { return n_collision_; } + + bool& write_track() { return write_track_; } + uint64_t& seeds(const int& i) { return seeds_[i]; } + uint64_t* seeds() { return seeds_; } + int& stream() { return stream_; } + + Particle::Bank& secondary_bank(const int& i) { return secondary_bank_[i]; } + decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; } + int64_t& current_work() { return current_work_; } + decltype(flux_derivs_)& flux_derivs() { return flux_derivs_; } + const decltype(flux_derivs_)& flux_derivs() const { return flux_derivs_; } + decltype(filter_matches_)& filter_matches() { return filter_matches_; } + FilterMatch& filter_matches(const int& i) { return filter_matches_[i]; } + decltype(tracks_)& tracks() { return tracks_; } + decltype(nu_bank_)& nu_bank() { return nu_bank_; } + NuBank& nu_bank(const int& i) { return nu_bank_[i]; } + + double& keff_tally_absorption() { return keff_tally_absorption_; } + double& keff_tally_collision() { return keff_tally_collision_; } + double& keff_tally_tracklength() { return keff_tally_tracklength_; } + double& keff_tally_leakage() { return keff_tally_leakage_; } + + bool& trace() { return trace_; } + double& collision_distance() { return collision_distance_; } + int& n_event() { return n_event_; } + +#ifdef DAGMC + moab::DagMC::RayHistory& rayhistory() { return history_; } + Direction& last_dir() { return last_dir_; } +#endif + + int64_t& n_progeny() { return n_progeny_; } + + // Accessors for position in global coordinates + Position& r() { return coord_[0].r; } + const Position& r() const { return coord_[0].r; } + + // Accessors for position in local coordinates + Position& r_local() { return coord_[n_coord_ - 1].r; } + const Position& r_local() const { return coord_[n_coord_ - 1].r; } + + // Accessors for direction in global coordinates + Direction& u() { return coord_[0].u; } + const Direction& u() const { return coord_[0].u; } + + // Accessors for direction in local coordinates + Direction& u_local() { return coord_[n_coord_ - 1].u; } + const Direction& u_local() const { return coord_[n_coord_ - 1].u; } }; //============================================================================ diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 7b66c6696..44da79e99 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -171,7 +171,7 @@ T PlotBase::get_map() const { Particle p; p.r() = xyz; p.u() = dir; - p.coord_[0].universe = model::root_universe; + p.coord(0).universe = model::root_universe; int level = level_; int j{}; @@ -180,10 +180,10 @@ T PlotBase::get_map() const { p.r()[out_i] = xyz[out_i] - out_pixel * y; for (int x = 0; x < width; x++) { p.r()[in_i] = xyz[in_i] + in_pixel * x; - p.n_coord_ = 1; + p.n_coord() = 1; // local variables bool found_cell = exhaustive_find_cell(p); - j = p.n_coord_ - 1; + j = p.n_coord() - 1; if (level >= 0) { j = level; } if (found_cell) { data.set_value(y, x, p, j); diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp index fccae191b..884b06f81 100644 --- a/src/boundary_condition.cpp +++ b/src/boundary_condition.cpp @@ -113,7 +113,7 @@ void TranslationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const { // TODO: off-by-one on surface indices throughout this function. - int i_particle_surf = std::abs(p.surface_) - 1; + int i_particle_surf = std::abs(p.surface()) - 1; // Figure out which of the two BC surfaces were struck then find the // particle's new location and surface. @@ -121,10 +121,10 @@ TranslationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const int new_surface; if (i_particle_surf == i_surf_) { new_r = p.r() + translation_; - new_surface = p.surface_ > 0 ? j_surf_ + 1 : -(j_surf_ + 1); + new_surface = p.surface() > 0 ? j_surf_ + 1 : -(j_surf_ + 1); } else if (i_particle_surf == j_surf_) { new_r = p.r() - translation_; - new_surface = p.surface_ > 0 ? i_surf_ + 1 : -(i_surf_ + 1); + new_surface = p.surface() > 0 ? i_surf_ + 1 : -(i_surf_ + 1); } else { throw std::runtime_error("Called BoundaryCondition::handle_particle after " "hitting a surface, but that surface is not recognized by the BC."); @@ -222,7 +222,7 @@ void RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const { // TODO: off-by-one on surface indices throughout this function. - int i_particle_surf = std::abs(p.surface_) - 1; + int i_particle_surf = std::abs(p.surface()) - 1; // Figure out which of the two BC surfaces were struck to figure out if a // forward or backward rotation is required. Specify the other surface as @@ -231,10 +231,10 @@ RotationalPeriodicBC::handle_particle(Particle& p, const Surface& surf) const int new_surface; if (i_particle_surf == i_surf_) { theta = angle_; - new_surface = p.surface_ > 0 ? -(j_surf_ + 1) : j_surf_ + 1; + new_surface = p.surface() > 0 ? -(j_surf_ + 1) : j_surf_ + 1; } else if (i_particle_surf == j_surf_) { theta = -angle_; - new_surface = p.surface_ > 0 ? -(i_surf_ + 1) : i_surf_ + 1; + new_surface = p.surface() > 0 ? -(i_surf_ + 1) : i_surf_ + 1; } else { throw std::runtime_error("Called BoundaryCondition::handle_particle after " "hitting a surface, but that surface is not recognized by the BC."); diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index a7fae50c1..f96b60c1e 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -28,20 +28,22 @@ std::vector ttb; void thick_target_bremsstrahlung(Particle& p, double* E_lost) { - if (p.material_ == MATERIAL_VOID) return; + if (p.material() == MATERIAL_VOID) + return; int photon = static_cast(Particle::Type::photon); - if (p.E_ < settings::energy_cutoff[photon]) return; + if (p.E() < settings::energy_cutoff[photon]) + return; // Get bremsstrahlung data for this material and particle type BremsstrahlungData* mat; - if (p.type_ == Particle::Type::positron) { - mat = &model::materials[p.material_]->ttb_->positron; + if (p.type() == Particle::Type::positron) { + mat = &model::materials[p.material()]->ttb_->positron; } else { - mat = &model::materials[p.material_]->ttb_->electron; + mat = &model::materials[p.material()]->ttb_->electron; } - double e = std::log(p.E_); + double e = std::log(p.E()); auto n_e = data::ttb_e_grid.size(); // Find the lower bounding index of the incident electron energy @@ -108,7 +110,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) if (w > settings::energy_cutoff[photon]) { // Create secondary photon - p.create_secondary(p.wgt_, p.u(), w, Particle::Type::photon); + p.create_secondary(p.wgt(), p.u(), w, Particle::Type::photon); *E_lost += w; } } diff --git a/src/event.cpp b/src/event.cpp index d1b4a36a7..6072a7712 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -50,7 +50,8 @@ void free_event_queues(void) void dispatch_xs_event(int64_t buffer_idx) { Particle& p = simulation::particles[buffer_idx]; - if (p.material_ == MATERIAL_VOID || !model::materials[p.material_]->fissionable_) { + if (p.material() == MATERIAL_VOID || + !model::materials[p.material()]->fissionable_) { simulation::calculate_nonfuel_xs_queue.thread_safe_append({p, buffer_idx}); } else { simulation::calculate_fuel_xs_queue.thread_safe_append({p, buffer_idx}); @@ -109,7 +110,7 @@ void process_advance_particle_events() int64_t buffer_idx = simulation::advance_particle_queue[i].idx; Particle& p = simulation::particles[buffer_idx]; p.event_advance(); - if (p.collision_distance_ > p.boundary_.distance) { + if (p.collision_distance() > p.boundary().distance) { simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx}); } else { simulation::collision_queue.thread_safe_append({p, buffer_idx}); @@ -131,7 +132,7 @@ void process_surface_crossing_events() Particle& p = simulation::particles[buffer_idx]; p.event_cross_surface(); p.event_revive_from_secondary(); - if (p.alive_) + if (p.alive()) dispatch_xs_event(buffer_idx); } @@ -150,7 +151,7 @@ void process_collision_events() Particle& p = simulation::particles[buffer_idx]; p.event_collide(); p.event_revive_from_secondary(); - if (p.alive_) + if (p.alive()) dispatch_xs_event(buffer_idx); } diff --git a/src/geometry.cpp b/src/geometry.cpp index ce6dfd65b..7103b0a5f 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -37,25 +37,25 @@ std::vector overlap_check_count; bool check_cell_overlap(Particle& p, bool error) { - int n_coord = p.n_coord_; + int n_coord = p.n_coord(); // Loop through each coordinate level for (int j = 0; j < n_coord; j++) { - Universe& univ = *model::universes[p.coord_[j].universe]; + Universe& univ = *model::universes[p.coord(j).universe]; // Loop through each cell on this level for (auto index_cell : univ.cells_) { Cell& c = *model::cells[index_cell]; - if (c.contains(p.coord_[j].r, p.coord_[j].u, p.surface_)) { - if (index_cell != p.coord_[j].cell) { + if (c.contains(p.coord(j).r, p.coord(j).u, p.surface())) { + if (index_cell != p.coord(j).cell) { if (error) { - fatal_error(fmt::format( - "Overlapping cells detected: {}, {} on universe {}", - c.id_, model::cells[p.coord_[j].cell]->id_, univ.id_)); + fatal_error( + fmt::format("Overlapping cells detected: {}, {} on universe {}", + c.id_, model::cells[p.coord(j).cell]->id_, univ.id_)); } return true; } - #pragma omp atomic +#pragma omp atomic ++model::overlap_check_count[index_cell]; } } @@ -78,15 +78,15 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) i_cell = *it; // Make sure the search cell is in the same universe. - int i_universe = p.coord_[p.n_coord_-1].universe; + int i_universe = p.coord(p.n_coord() - 1).universe; if (model::cells[i_cell]->universe_ != i_universe) continue; // Check if this cell contains the particle. Position r {p.r_local()}; Direction u {p.u_local()}; - auto surf = p.surface_; + auto surf = p.surface(); if (model::cells[i_cell]->contains(r, u, surf)) { - p.coord_[p.n_coord_-1].cell = i_cell; + p.coord(p.n_coord() - 1).cell = i_cell; found = true; break; } @@ -102,7 +102,7 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) } // Check successively lower coordinate levels until finding material fill - for (;;++p.n_coord_) { + for (;;++p.n_coord()) { // If we did not attempt to use neighbor lists, i_cell is still C_NONE. In // that case, we should now do an exhaustive search to find the right value // of i_cell. @@ -112,7 +112,7 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) // code below this conditional, we set i_cell back to C_NONE to indicate // that. if (i_cell == C_NONE) { - int i_universe = p.coord_[p.n_coord_-1].universe; + int i_universe = p.coord(p.n_coord() - 1).universe; const auto& univ {*model::universes[i_universe]}; const auto& cells { !univ.partitioner_ @@ -124,15 +124,15 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) i_cell = *it; // Make sure the search cell is in the same universe. - int i_universe = p.coord_[p.n_coord_-1].universe; + int i_universe = p.coord(p.n_coord() - 1).universe; if (model::cells[i_cell]->universe_ != i_universe) continue; // Check if this cell contains the particle. Position r {p.r_local()}; Direction u {p.u_local()}; - auto surf = p.surface_; + auto surf = p.surface(); if (model::cells[i_cell]->contains(r, u, surf)) { - p.coord_[p.n_coord_-1].cell = i_cell; + p.coord(p.n_coord() - 1).cell = i_cell; found = true; break; } @@ -143,7 +143,7 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) } // Announce the cell that the particle is entering. - if (found && (settings::verbosity >= 10 || p.trace_)) { + if (found && (settings::verbosity >= 10 || p.trace())) { auto msg = fmt::format(" Entering cell {}", model::cells[i_cell]->id_); write_message(msg, 1); } @@ -155,8 +155,8 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) // Find the distribcell instance number. int offset = 0; if (c.distribcell_index_ >= 0) { - for (int i = 0; i < p.n_coord_; i++) { - const auto& c_i {*model::cells[p.coord_[i].cell]}; + for (int i = 0; i < p.n_coord(); i++) { + const auto& c_i {*model::cells[p.coord(i).cell]}; if (c_i.type_ == Fill::UNIVERSE) { offset += c_i.offset_[c.distribcell_index_]; } else if (c_i.type_ == Fill::LATTICE) { @@ -169,20 +169,20 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) } } } - p.cell_instance_ = offset; + p.cell_instance() = offset; // Set the material and temperature. - p.material_last_ = p.material_; + p.material_last() = p.material(); if (c.material_.size() > 1) { - p.material_ = c.material_[p.cell_instance_]; + p.material() = c.material_[p.cell_instance()]; } else { - p.material_ = c.material_[0]; + p.material() = c.material_[0]; } - p.sqrtkT_last_ = p.sqrtkT_; + p.sqrtkT_last() = p.sqrtkT(); if (c.sqrtkT_.size() > 1) { - p.sqrtkT_ = c.sqrtkT_[p.cell_instance_]; + p.sqrtkT() = c.sqrtkT_[p.cell_instance()]; } else { - p.sqrtkT_ = c.sqrtkT_[0]; + p.sqrtkT() = c.sqrtkT_[0]; } return true; @@ -192,7 +192,7 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) //! Found a lower universe, update this coord level then search the next. // Set the lower coordinate level universe. - auto& coord {p.coord_[p.n_coord_]}; + auto& coord {p.coord(p.n_coord())}; coord.universe = c.fill_; // Set the position and direction. @@ -214,7 +214,7 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) Lattice& lat {*model::lattices[c.fill_]}; // Set the position and direction. - auto& coord {p.coord_[p.n_coord_]}; + auto& coord {p.coord(p.n_coord())}; coord.r = p.r_local(); coord.u = p.u_local(); @@ -247,7 +247,7 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) } else { warning(fmt::format("Particle {} is outside lattice {} but the " "lattice has no defined outer universe.", - p.id_, lat.id_)); + p.id(), lat.id_)); return false; } } @@ -265,13 +265,13 @@ bool neighbor_list_find_cell(Particle& p) { // Reset all the deeper coordinate levels. - for (int i = p.n_coord_; i < p.coord_.size(); i++) { - p.coord_[i].reset(); + for (int i = p.n_coord(); i < model::n_coord_levels; i++) { + p.coord(i).reset(); } // Get the cell this particle was in previously. - auto coord_lvl = p.n_coord_ - 1; - auto i_cell = p.coord_[coord_lvl].cell; + auto coord_lvl = p.n_coord() - 1; + auto i_cell = p.coord(coord_lvl).cell; Cell& c {*model::cells[i_cell]}; // Search for the particle in that cell's neighbor list. Return if we @@ -285,21 +285,21 @@ bool neighbor_list_find_cell(Particle& p) // neighboring cell. found = find_cell_inner(p, nullptr); if (found) - c.neighbors_.push_back(p.coord_[coord_lvl].cell); + c.neighbors_.push_back(p.coord(coord_lvl).cell); return found; } bool exhaustive_find_cell(Particle& p) { - int i_universe = p.coord_[p.n_coord_-1].universe; + int i_universe = p.coord(p.n_coord() - 1).universe; if (i_universe == C_NONE) { - p.coord_[0].universe = model::root_universe; - p.n_coord_ = 1; + p.coord(0).universe = model::root_universe; + p.n_coord() = 1; i_universe = model::root_universe; } // Reset all the deeper coordinate levels. - for (int i = p.n_coord_; i < p.coord_.size(); i++) { - p.coord_[i].reset(); + for (int i = p.n_coord(); i < model::n_coord_levels; i++) { + p.coord(i).reset(); } return find_cell_inner(p, nullptr); } @@ -309,10 +309,10 @@ bool exhaustive_find_cell(Particle& p) void cross_lattice(Particle& p, const BoundaryInfo& boundary) { - auto& coord {p.coord_[p.n_coord_ - 1]}; + auto& coord {p.coord(p.n_coord() - 1)}; auto& lat {*model::lattices[coord.lattice]}; - if (settings::verbosity >= 10 || p.trace_) { + if (settings::verbosity >= 10 || p.trace()) { write_message(fmt::format( " Crossing lattice {}. Current position ({},{},{}). r={}", lat.id_, coord.lattice_i[0], coord.lattice_i[1], coord.lattice_i[2], p.r()), 1); @@ -325,7 +325,7 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary) std::array i_xyz {coord.lattice_x, coord.lattice_y, coord.lattice_z}; // Set the new coordinate position. - const auto& upper_coord {p.coord_[p.n_coord_ - 2]}; + const auto& upper_coord {p.coord(p.n_coord() - 2)}; const auto& cell {model::cells[upper_coord.cell]}; Position r = upper_coord.r; r -= cell->translation_; @@ -336,11 +336,12 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary) if (!lat.are_valid_indices(i_xyz)) { // The particle is outside the lattice. Search for it from the base coords. - p.n_coord_ = 1; + p.n_coord() = 1; bool found = exhaustive_find_cell(p); - if (!found && p.alive_) { + if (!found && p.alive()) { p.mark_as_lost(fmt::format("Could not locate particle {} after " - "crossing a lattice boundary", p.id_)); + "crossing a lattice boundary", + p.id())); } } else { @@ -351,11 +352,12 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary) if (!found) { // A particle crossing the corner of a lattice tile may not be found. In // this case, search for it from the base coords. - p.n_coord_ = 1; + p.n_coord() = 1; bool found = exhaustive_find_cell(p); - if (!found && p.alive_) { + if (!found && p.alive()) { p.mark_as_lost(fmt::format("Could not locate particle {} after " - "crossing a lattice boundary", p.id_)); + "crossing a lattice boundary", + p.id())); } } } @@ -372,14 +374,14 @@ BoundaryInfo distance_to_boundary(Particle& p) std::array level_lat_trans {}; // Loop over each coordinate level. - for (int i = 0; i < p.n_coord_; i++) { - const auto& coord {p.coord_[i]}; - Position r {coord.r}; - Direction u {coord.u}; + for (int i = 0; i < p.n_coord(); i++) { + const auto& coord {p.coord(i)}; + const Position& r {coord.r}; + const Direction& u {coord.u}; Cell& c {*model::cells[coord.cell]}; // Find the oncoming surface in this cell and the distance to it. - auto surface_distance = c.distance(r, u, p.surface_, &p); + auto surface_distance = c.distance(r, u, p.surface(), &p); d_surf = surface_distance.first; level_surf_cross = surface_distance.second; @@ -395,8 +397,8 @@ BoundaryInfo distance_to_boundary(Particle& p) lattice_distance = lat.distance(r, u, i_xyz); break; case LatticeType::hex: - auto& cell_above {model::cells[p.coord_[i-1].cell]}; - Position r_hex {p.coord_[i-1].r}; + auto& cell_above {model::cells[p.coord(i - 1).cell]}; + Position r_hex {p.coord(i - 1).r}; r_hex -= cell_above->translation_; if (coord.rotated) { r_hex = r_hex.rotate(cell_above->rotation_); @@ -410,7 +412,7 @@ BoundaryInfo distance_to_boundary(Particle& p) if (d_lat < 0) { p.mark_as_lost(fmt::format( - "Particle {} had a negative distance to a lattice boundary", p.id_)); + "Particle {} had a negative distance to a lattice boundary", p.id())); } } @@ -473,8 +475,8 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance) return OPENMC_E_GEOMETRY; } - *index = p.coord_[p.n_coord_-1].cell; - *instance = p.cell_instance_; + *index = p.coord(p.n_coord() - 1).cell; + *instance = p.cell_instance(); return 0; } diff --git a/src/material.cpp b/src/material.cpp index 89d5d5add..4cc2f7056 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -748,14 +748,14 @@ void Material::init_nuclide_index() void Material::calculate_xs(Particle& p) const { // Set all material macroscopic cross sections to zero - p.macro_xs_.total = 0.0; - p.macro_xs_.absorption = 0.0; - p.macro_xs_.fission = 0.0; - p.macro_xs_.nu_fission = 0.0; + p.macro_xs().total = 0.0; + p.macro_xs().absorption = 0.0; + p.macro_xs().fission = 0.0; + p.macro_xs().nu_fission = 0.0; - if (p.type_ == Particle::Type::neutron) { + if (p.type() == Particle::Type::neutron) { this->calculate_neutron_xs(p); - } else if (p.type_ == Particle::Type::photon) { + } else if (p.type() == Particle::Type::photon) { this->calculate_photon_xs(p); } } @@ -764,7 +764,8 @@ void Material::calculate_neutron_xs(Particle& p) const { // Find energy index on energy grid int neutron = static_cast(Particle::Type::neutron); - int i_grid = std::log(p.E_/data::energy_min[neutron])/simulation::log_spacing; + int i_grid = + std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing; // Determine if this material has S(a,b) tables bool check_sab = (thermal_tables_.size() > 0); @@ -791,7 +792,8 @@ void Material::calculate_neutron_xs(Particle& p) const // If particle energy is greater than the highest energy for the // S(a,b) table, then don't use the S(a,b) table - if (p.E_ > data::thermal_scatt[i_sab]->energy_max_) i_sab = C_NONE; + if (p.E() > data::thermal_scatt[i_sab]->energy_max_) + i_sab = C_NONE; // Increment position in thermal_tables_ ++j; @@ -808,11 +810,9 @@ void Material::calculate_neutron_xs(Particle& p) const int i_nuclide = nuclide_[i]; // Calculate microscopic cross section for this nuclide - const auto& micro {p.neutron_xs_[i_nuclide]}; - if (p.E_ != micro.last_E - || p.sqrtkT_ != micro.last_sqrtkT - || i_sab != micro.index_sab - || sab_frac != micro.sab_frac) { + const auto& micro {p.neutron_xs(i_nuclide)}; + if (p.E() != micro.last_E || p.sqrtkT() != micro.last_sqrtkT || + i_sab != micro.index_sab || sab_frac != micro.sab_frac) { data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p); } @@ -823,19 +823,19 @@ void Material::calculate_neutron_xs(Particle& p) const double atom_density = atom_density_(i); // Add contributions to cross sections - p.macro_xs_.total += atom_density * micro.total; - p.macro_xs_.absorption += atom_density * micro.absorption; - p.macro_xs_.fission += atom_density * micro.fission; - p.macro_xs_.nu_fission += atom_density * micro.nu_fission; + p.macro_xs().total += atom_density * micro.total; + p.macro_xs().absorption += atom_density * micro.absorption; + p.macro_xs().fission += atom_density * micro.fission; + p.macro_xs().nu_fission += atom_density * micro.nu_fission; } } void Material::calculate_photon_xs(Particle& p) const { - p.macro_xs_.coherent = 0.0; - p.macro_xs_.incoherent = 0.0; - p.macro_xs_.photoelectric = 0.0; - p.macro_xs_.pair_production = 0.0; + p.macro_xs().coherent = 0.0; + p.macro_xs().incoherent = 0.0; + p.macro_xs().photoelectric = 0.0; + p.macro_xs().pair_production = 0.0; // Add contribution from each nuclide in material for (int i = 0; i < nuclide_.size(); ++i) { @@ -846,8 +846,8 @@ void Material::calculate_photon_xs(Particle& p) const int i_element = element_[i]; // Calculate microscopic cross section for this nuclide - const auto& micro {p.photon_xs_[i_element]}; - if (p.E_ != micro.last_E) { + const auto& micro {p.photon_xs(i_element)}; + if (p.E() != micro.last_E) { data::elements[i_element]->calculate_xs(p); } @@ -858,11 +858,11 @@ void Material::calculate_photon_xs(Particle& p) const double atom_density = atom_density_(i); // Add contributions to material macroscopic cross sections - p.macro_xs_.total += atom_density * micro.total; - p.macro_xs_.coherent += atom_density * micro.coherent; - p.macro_xs_.incoherent += atom_density * micro.incoherent; - p.macro_xs_.photoelectric += atom_density * micro.photoelectric; - p.macro_xs_.pair_production += atom_density * micro.pair_production; + p.macro_xs().total += atom_density * micro.total; + p.macro_xs().coherent += atom_density * micro.coherent; + p.macro_xs().incoherent += atom_density * micro.incoherent; + p.macro_xs().photoelectric += atom_density * micro.photoelectric; + p.macro_xs().pair_production += atom_density * micro.pair_production; } } diff --git a/src/mgxs.cpp b/src/mgxs.cpp index dec505f3c..02d27f960 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -627,13 +627,13 @@ Mgxs::calculate_xs(Particle& p) #else int tid = 0; #endif - set_temperature_index(p.sqrtkT_); + set_temperature_index(p.sqrtkT()); set_angle_index(p.u_local()); XsData* xs_t = &xs[cache[tid].t]; - p.macro_xs_.total = xs_t->total(cache[tid].a, p.g_); - p.macro_xs_.absorption = xs_t->absorption(cache[tid].a, p.g_); - p.macro_xs_.nu_fission = - fissionable ? xs_t->nu_fission(cache[tid].a, p.g_) : 0.; + p.macro_xs().total = xs_t->total(cache[tid].a, p.g()); + p.macro_xs().absorption = xs_t->absorption(cache[tid].a, p.g()); + p.macro_xs().nu_fission = + fissionable ? xs_t->nu_fission(cache[tid].a, p.g()) : 0.; } //============================================================================== diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 5ecbdc3b4..832677fff 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -519,7 +519,7 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const void Nuclide::calculate_elastic_xs(Particle& p) const { // Get temperature index, grid index, and interpolation factor - auto& micro {p.neutron_xs_[index_]}; + auto& micro {p.neutron_xs(index_)}; int i_temp = micro.index_temp; int i_grid = micro.index_grid; double f = micro.interp_factor; @@ -555,7 +555,7 @@ double Nuclide::elastic_xs_0K(double E) const void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle& p) { - auto& micro {p.neutron_xs_[index_]}; + auto& micro {p.neutron_xs(index_)}; // Initialize cached cross sections to zero micro.elastic = CACHE_INVALID; @@ -565,21 +565,21 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle // Check to see if there is multipole data present at this energy bool use_mp = false; if (multipole_) { - use_mp = (p.E_ >= multipole_->E_min_ && p.E_ <= multipole_->E_max_); + use_mp = (p.E() >= multipole_->E_min_ && p.E() <= multipole_->E_max_); } // Evaluate multipole or interpolate if (use_mp) { // Call multipole kernel double sig_s, sig_a, sig_f; - std::tie(sig_s, sig_a, sig_f) = multipole_->evaluate(p.E_, p.sqrtkT_); + std::tie(sig_s, sig_a, sig_f) = multipole_->evaluate(p.E(), p.sqrtkT()); micro.total = sig_s + sig_a; micro.elastic = sig_s; micro.absorption = sig_a; micro.fission = sig_f; - micro.nu_fission = fissionable_ ? - sig_f * this->nu(p.E_, EmissionMode::total) : 0.0; + micro.nu_fission = + fissionable_ ? sig_f * this->nu(p.E(), EmissionMode::total) : 0.0; if (simulation::need_depletion_rx) { // Only non-zero reaction is (n,gamma) @@ -607,7 +607,7 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle } else { // Find the appropriate temperature index. - double kT = p.sqrtkT_*p.sqrtkT_; + double kT = p.sqrtkT() * p.sqrtkT(); double f; int i_temp = -1; switch (settings::temperature_method) { @@ -644,9 +644,9 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle const auto& xs {xs_[i_temp]}; int i_grid; - if (p.E_ < grid.energy.front()) { + if (p.E() < grid.energy.front()) { i_grid = 0; - } else if (p.E_ > grid.energy.back()) { + } else if (p.E() > grid.energy.back()) { i_grid = grid.energy.size() - 2; } else { // Determine bounding indices based on which equal log-spaced @@ -655,15 +655,16 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle int i_high = grid.grid_index[i_log_union + 1] + 1; // Perform binary search over reduced range - i_grid = i_low + lower_bound_index(&grid.energy[i_low], &grid.energy[i_high], p.E_); + i_grid = i_low + lower_bound_index( + &grid.energy[i_low], &grid.energy[i_high], p.E()); } // check for rare case where two energy points are the same if (grid.energy[i_grid] == grid.energy[i_grid + 1]) ++i_grid; // calculate interpolation factor - f = (p.E_ - grid.energy[i_grid]) / - (grid.energy[i_grid + 1]- grid.energy[i_grid]); + f = (p.E() - grid.energy[i_grid]) / + (grid.energy[i_grid + 1] - grid.energy[i_grid]); micro.index_temp = i_temp; micro.index_grid = i_grid; @@ -750,19 +751,19 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle // probability tables, we need to determine cross sections from the table if (settings::urr_ptables_on && urr_present_ && !use_mp) { int n = urr_data_[micro.index_temp].n_energy_; - if ((p.E_ > urr_data_[micro.index_temp].energy_(0)) && - (p.E_ < urr_data_[micro.index_temp].energy_(n-1))) { + if ((p.E() > urr_data_[micro.index_temp].energy_(0)) && + (p.E() < urr_data_[micro.index_temp].energy_(n - 1))) { this->calculate_urr_xs(micro.index_temp, p); } } - micro.last_E = p.E_; - micro.last_sqrtkT = p.sqrtkT_; + micro.last_E = p.E(); + micro.last_sqrtkT = p.sqrtkT(); } void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p) { - auto& micro {p.neutron_xs_[index_]}; + auto& micro {p.neutron_xs(index_)}; // Set flag that S(a,b) treatment should be used for scattering micro.index_sab = i_sab; @@ -771,7 +772,8 @@ void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p) int i_temp; double elastic; double inelastic; - data::thermal_scatt[i_sab]->calculate_xs(p.E_, p.sqrtkT_, &i_temp, &elastic, &inelastic, p.current_seed()); + data::thermal_scatt[i_sab]->calculate_xs( + p.E(), p.sqrtkT(), &i_temp, &elastic, &inelastic, p.current_seed()); // Store the S(a,b) cross sections. micro.thermal = sab_frac * (elastic + inelastic); @@ -791,7 +793,7 @@ void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p) void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const { - auto& micro = p.neutron_xs_[index_]; + auto& micro = p.neutron_xs(index_); micro.use_ptable = true; // Create a shorthand for the URR data @@ -799,19 +801,21 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const // Determine the energy table int i_energy = 0; - while (p.E_ >= urr.energy_(i_energy + 1)) {++i_energy;}; + while (p.E() >= urr.energy_(i_energy + 1)) { + ++i_energy; + }; // Sample the probability table using the cumulative distribution - // Random nmbers for the xs calculation are sampled from a separate stream. + // Random numbers for the xs calculation are sampled from a separate stream. // This guarantees the randomness and, at the same time, makes sure we // reuse random numbers for the same nuclide at different temperatures, // therefore preserving correlation of temperature in probability tables. - p.stream_ = STREAM_URR_PTABLE; + p.stream() = STREAM_URR_PTABLE; //TODO: to maintain the same random number stream as the Fortran code this //replaces, the seed is set with index_ + 1 instead of index_ double r = future_prn(static_cast(index_ + 1), *p.current_seed()); - p.stream_ = STREAM_TRACKING; + p.stream() = STREAM_TRACKING; int i_low = 0; while (urr.prob_(i_energy, URRTableParam::CUM_PROB, i_low) <= r) {++i_low;}; @@ -827,8 +831,8 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const double f; if (urr.interp_ == Interpolation::lin_lin) { // Determine the interpolation factor on the table - f = (p.E_ - urr.energy_(i_energy)) / - (urr.energy_(i_energy + 1) - urr.energy_(i_energy)); + f = (p.E() - urr.energy_(i_energy)) / + (urr.energy_(i_energy + 1) - urr.energy_(i_energy)); elastic = (1. - f) * urr.prob_(i_energy, URRTableParam::ELASTIC, i_low) + f * urr.prob_(i_energy + 1, URRTableParam::ELASTIC, i_up); @@ -838,8 +842,8 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const f * urr.prob_(i_energy + 1, URRTableParam::N_GAMMA, i_up); } else if (urr.interp_ == Interpolation::log_log) { // Determine interpolation factor on the table - f = std::log(p.E_ / urr.energy_(i_energy)) / - std::log(urr.energy_(i_energy + 1) / urr.energy_(i_energy)); + f = std::log(p.E() / urr.energy_(i_energy)) / + std::log(urr.energy_(i_energy + 1) / urr.energy_(i_energy)); // Calculate the elastic cross section/factor if ((urr.prob_(i_energy, URRTableParam::ELASTIC, i_low) > 0.) && @@ -916,7 +920,7 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const // Determine nu-fission cross-section if (fissionable_) { - micro.nu_fission = nu(p.E_, EmissionMode::total) * micro.fission; + micro.nu_fission = nu(p.E(), EmissionMode::total) * micro.fission; } } diff --git a/src/output.cpp b/src/output.cpp index 89a9070a2..458026766 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -146,62 +146,62 @@ std::string time_stamp() void print_particle(Particle& p) { // Display particle type and ID. - switch (p.type_) { - case Particle::Type::neutron: - fmt::print("Neutron "); - break; - case Particle::Type::photon: - fmt::print("Photon "); - break; - case Particle::Type::electron: - fmt::print("Electron "); - break; - case Particle::Type::positron: - fmt::print("Positron "); - break; - default: - fmt::print("Unknown Particle "); + switch (p.type()) { + case Particle::Type::neutron: + fmt::print("Neutron "); + break; + case Particle::Type::photon: + fmt::print("Photon "); + break; + case Particle::Type::electron: + fmt::print("Electron "); + break; + case Particle::Type::positron: + fmt::print("Positron "); + break; + default: + fmt::print("Unknown Particle "); } - fmt::print("{}\n", p.id_); + fmt::print("{}\n", p.id()); // Display particle geometry hierarchy. - for (auto i = 0; i < p.n_coord_; i++) { + for (auto i = 0; i < p.n_coord(); i++) { fmt::print(" Level {}\n", i); - if (p.coord_[i].cell != C_NONE) { - const Cell& c {*model::cells[p.coord_[i].cell]}; + if (p.coord(i).cell != C_NONE) { + const Cell& c {*model::cells[p.coord(i).cell]}; fmt::print(" Cell = {}\n", c.id_); } - if (p.coord_[i].universe != C_NONE) { - const Universe& u {*model::universes[p.coord_[i].universe]}; + if (p.coord(i).universe != C_NONE) { + const Universe& u {*model::universes[p.coord(i).universe]}; fmt::print(" Universe = {}\n", u.id_); } - if (p.coord_[i].lattice != C_NONE) { - const Lattice& lat {*model::lattices[p.coord_[i].lattice]}; + if (p.coord(i).lattice != C_NONE) { + const Lattice& lat {*model::lattices[p.coord(i).lattice]}; fmt::print(" Lattice = {}\n", lat.id_); fmt::print(" Lattice position = ({},{},{})\n", p.coord_[i].lattice_x, p.coord_[i].lattice_y, p.coord_[i].lattice_z); } - fmt::print(" r = {}\n", p.coord_[i].r); - fmt::print(" u = {}\n", p.coord_[i].u); + fmt::print(" r = {}\n", p.coord(i).r); + fmt::print(" u = {}\n", p.coord(i).u); } // Display miscellaneous info. - if (p.surface_ != 0) { + if (p.surface() != 0) { // Surfaces identifiers are >= 1, but indices are >= 0 so we need -1 - const Surface& surf {*model::surfaces[std::abs(p.surface_)-1]}; - fmt::print(" Surface = {}\n", (p.surface_ > 0) ? surf.id_ : -surf.id_); + const Surface& surf {*model::surfaces[std::abs(p.surface()) - 1]}; + fmt::print(" Surface = {}\n", (p.surface() > 0) ? surf.id_ : -surf.id_); } - fmt::print(" Weight = {}\n", p.wgt_); + fmt::print(" Weight = {}\n", p.wgt()); if (settings::run_CE) { - fmt::print(" Energy = {}\n", p.E_); + fmt::print(" Energy = {}\n", p.E()); } else { - fmt::print(" Energy Group = {}\n", p.g_); + fmt::print(" Energy Group = {}\n", p.g()); } - fmt::print(" Delayed Group = {}\n\n", p.delayed_group_); + fmt::print(" Delayed Group = {}\n\n", p.delayed_group()); } //============================================================================== diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index be9eda77f..9d9be30ef 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -41,28 +41,28 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) } else if (mode == "fixed source") { previous_run_mode = RunMode::FIXED_SOURCE; } - read_dataset(file_id, "id", p.id_); + read_dataset(file_id, "id", p.id()); int type; read_dataset(file_id, "type", type); - p.type_ = static_cast(type); - read_dataset(file_id, "weight", p.wgt_); - read_dataset(file_id, "energy", p.E_); + p.type() = static_cast(type); + read_dataset(file_id, "weight", p.wgt()); + read_dataset(file_id, "energy", p.E()); read_dataset(file_id, "xyz", p.r()); read_dataset(file_id, "uvw", p.u()); // Set energy group and average energy in multi-group mode if (!settings::run_CE) { - p.g_ = p.E_; - p.E_ = data::mg.energy_bin_avg_[p.g_]; + p.g() = p.E(); + p.E() = data::mg.energy_bin_avg_[p.g()]; } // Set particle last attributes - p.wgt_last_ = p.wgt_; - p.r_last_current_ = p.r(); - p.r_last_ = p.r(); - p.u_last_ = p.u(); - p.E_last_ = p.E_; - p.g_last_ = p.g_; + p.wgt_last() = p.wgt(); + p.r_last_current() = p.r(); + p.r_last() = p.r(); + p.u_last() = p.u(); + p.E_last() = p.E(); + p.g_last() = p.g(); // Close hdf5 file file_close(file_id); @@ -84,7 +84,8 @@ void run_particle_restart() read_particle_restart(p, previous_run_mode); // write track if that was requested on command line - if (settings::write_all_tracks) p.write_track_ = true; + if (settings::write_all_tracks) + p.write_track() = true; // Set all tallies to 0 for now (just tracking errors) model::tallies.clear(); @@ -93,33 +94,37 @@ void run_particle_restart() int64_t particle_seed; switch (previous_run_mode) { case RunMode::EIGENVALUE: - particle_seed = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + p.id_; + particle_seed = (simulation::total_gen + overall_generation() - 1) * + settings::n_particles + + p.id(); break; case RunMode::FIXED_SOURCE: - particle_seed = p.id_; + particle_seed = p.id(); break; default: throw std::runtime_error{"Unexpected run mode: " + std::to_string(static_cast(previous_run_mode))}; } - init_particle_seeds(particle_seed, p.seeds_); + init_particle_seeds(particle_seed, p.seeds()); // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { - for (auto& micro : p.neutron_xs_) micro.last_E = 0.0; + p.invalidate_neutron_xs(); } // Prepare to write out particle track. - if (p.write_track_) add_particle_track(p); + if (p.write_track()) + add_particle_track(p); // Every particle starts with no accumulated flux derivative. if (!model::active_tallies.empty()) { - p.flux_derivs_.resize(model::tally_derivs.size(), 0.0); - std::fill(p.flux_derivs_.begin(), p.flux_derivs_.end(), 0.0); + p.flux_derivs().resize(model::tally_derivs.size(), 0.0); + std::fill(p.flux_derivs().begin(), p.flux_derivs().end(), 0.0); } - // Allocate space for tally filter matches - p.filter_matches_.resize(model::tally_filters.size()); + // Allocate space for tally filter matches (TODO shouldn't this be in the + // particle constructor, instead?) + p.filter_matches().resize(model::tally_filters.size()); // Transport neutron transport_history_based_single_particle(p); diff --git a/src/photon.cpp b/src/photon.cpp index 6333f9b29..4722a2fac 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -454,7 +454,7 @@ void PhotonInteraction::calculate_xs(Particle& p) const // Perform binary search on the element energy grid in order to determine // which points to interpolate between int n_grid = energy_.size(); - double log_E = std::log(p.E_); + double log_E = std::log(p.E()); int i_grid; if (log_E <= energy_[0]) { i_grid = 0; @@ -472,7 +472,7 @@ void PhotonInteraction::calculate_xs(Particle& p) const // calculate interpolation factor double f = (log_E - energy_(i_grid)) / (energy_(i_grid+1) - energy_(i_grid)); - auto& xs {p.photon_xs_[index_]}; + auto& xs {p.photon_xs(index_)}; xs.index_grid = i_grid; xs.interp_factor = f; @@ -506,7 +506,7 @@ void PhotonInteraction::calculate_xs(Particle& p) const // Calculate microscopic total cross section xs.total = xs.coherent + xs.incoherent + xs.photoelectric + xs.pair_production; - xs.last_E = p.E_; + xs.last_E = p.E(); } double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* seed) const @@ -669,7 +669,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); double E = shell.binding_energy; - p.create_secondary(p.wgt_, u, E, Particle::Type::photon); + p.create_secondary(p.wgt(), u, E, Particle::Type::photon); return; } @@ -701,7 +701,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Non-radiative transition -- Auger/Coster-Kronig effect // Create auger electron - p.create_secondary(p.wgt_, u, E, Particle::Type::electron); + p.create_secondary(p.wgt(), u, E, Particle::Type::electron); // Fill hole left by emitted auger electron int i_hole = shell_map_.at(secondary); @@ -711,7 +711,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Radiative transition -- get X-ray energy // Create fluorescent photon - p.create_secondary(p.wgt_, u, E, Particle::Type::photon); + p.create_secondary(p.wgt(), u, E, Particle::Type::photon); } // Fill hole created by electron transitioning to the photoelectron hole diff --git a/src/physics.cpp b/src/physics.cpp index 7a8be4ee0..5b58e0aa0 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -36,10 +36,10 @@ namespace openmc { void collision(Particle& p) { // Add to collision counter for particle - ++(p.n_collision_); + ++(p.n_collision()); // Sample reaction for the material the particle is in - switch (p.type_) { + switch (p.type()) { case Particle::Type::neutron: sample_neutron_reaction(p); break; @@ -55,27 +55,27 @@ void collision(Particle& p) } // Kill particle if energy falls below cutoff - int type = static_cast(p.type_); - if (p.E_ < settings::energy_cutoff[type]) { - p.alive_ = false; - p.wgt_ = 0.0; + int type = static_cast(p.type()); + if (p.E() < settings::energy_cutoff[type]) { + p.alive() = false; + p.wgt() = 0.0; } // Display information about collision - if (settings::verbosity >= 10 || p.trace_) { + if (settings::verbosity >= 10 || p.trace()) { std::string msg; - if (p.event_ == TallyEvent::KILL) { - msg = fmt::format(" Killed. Energy = {} eV.", p.E_); - } else if (p.type_ == Particle::Type::neutron) { + if (p.event() == TallyEvent::KILL) { + msg = fmt::format(" Killed. Energy = {} eV.", p.E()); + } else if (p.type() == Particle::Type::neutron) { msg = fmt::format(" {} with {}. Energy = {} eV.", - reaction_name(p.event_mt_), data::nuclides[p.event_nuclide_]->name_, - p.E_); - } else if (p.type_ == Particle::Type::photon) { + reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_, + p.E()); + } else if (p.type() == Particle::Type::photon) { msg = fmt::format(" {} with {}. Energy = {} eV.", - reaction_name(p.event_mt_), - to_element(data::nuclides[p.event_nuclide_]->name_), p.E_); + reaction_name(p.event_mt()), + to_element(data::nuclides[p.event_nuclide()]->name_), p.E()); } else { - msg = fmt::format(" Disappeared. Energy = {} eV.", p.E_); + msg = fmt::format(" Disappeared. Energy = {} eV.", p.E()); } write_message(msg, 1); } @@ -87,7 +87,7 @@ void sample_neutron_reaction(Particle& p) int i_nuclide = sample_nuclide(p); // Save which nuclide particle had collision with - p.event_nuclide_ = i_nuclide; + p.event_nuclide() = i_nuclide; // Create fission bank sites. Note that while a fission reaction is sampled, // it never actually "happens", i.e. the weight of the particle does not @@ -106,7 +106,7 @@ void sample_neutron_reaction(Particle& p) // Make sure particle population doesn't grow out of control for // subcritical multiplication problems. - if (p.secondary_bank_.size() >= 10000) { + if (p.secondary_bank().size() >= 10000) { fatal_error("The secondary particle bank appears to be growing without " "bound. You are likely running a subcritical multiplication problem " "with k-effective close to or greater than one."); @@ -116,36 +116,38 @@ void sample_neutron_reaction(Particle& p) // Create secondary photons if (settings::photon_transport) { - p.stream_ = STREAM_PHOTON; + p.stream() = STREAM_PHOTON; sample_secondary_photons(p, i_nuclide); - p.stream_ = STREAM_TRACKING; + p.stream() = STREAM_TRACKING; } // If survival biasing is being used, the following subroutine adjusts the // weight of the particle. Otherwise, it checks to see if absorption occurs - if (p.neutron_xs_[i_nuclide].absorption > 0.0) { + if (p.neutron_xs(i_nuclide).absorption > 0.0) { absorption(p, i_nuclide); } else { - p.wgt_absorb_ = 0.0; + p.wgt_absorb() = 0.0; } - if (!p.alive_) return; + if (!p.alive()) + return; // Sample a scattering reaction and determine the secondary energy of the // exiting neutron scatter(p, i_nuclide); // Advance URR seed stream 'N' times after energy changes - if (p.E_ != p.E_last_) { - p.stream_ = STREAM_URR_PTABLE; + if (p.E() != p.E_last()) { + p.stream() = STREAM_URR_PTABLE; advance_prn_seed(data::nuclides.size(), p.current_seed()); - p.stream_ = STREAM_TRACKING; + p.stream() = STREAM_TRACKING; } // Play russian roulette if survival biasing is turned on if (settings::survival_biasing) { russian_roulette(p); - if (!p.alive_) return; + if (!p.alive()) + return; } } @@ -157,8 +159,9 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p.wgt_ / simulation::keff * weight * p.neutron_xs_[ - i_nuclide].nu_fission / p.neutron_xs_[i_nuclide].total; + double nu_t = p.wgt() / simulation::keff * weight * + p.neutron_xs(i_nuclide).nu_fission / + p.neutron_xs(i_nuclide).total; // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -173,9 +176,9 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) double nu_d[MAX_DELAYED_GROUPS] = {0.}; // Clear out particle's nu fission bank - p.nu_bank_.clear(); + p.nu_bank().clear(); - p.fission_ = true; + p.fission() = true; int skipped = 0; // Determine whether to place fission sites into the shared fission bank @@ -188,12 +191,12 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) site.r = p.r(); site.particle = Particle::Type::neutron; site.wgt = 1. / weight; - site.parent_id = p.id_; - site.progeny_id = p.n_progeny_++; + site.parent_id = p.id(); + site.progeny_id = p.n_progeny()++; site.surf_id = 0; // Sample delayed group and angle/energy for fission reaction - sample_fission_neutron(i_nuclide, rx, p.E_, &site, p.current_seed()); + sample_fission_neutron(i_nuclide, rx, p.E(), &site, p.current_seed()); // Store fission site in bank if (use_fission_bank) { @@ -205,20 +208,20 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) break; } } else { - p.secondary_bank_.push_back(site); + p.secondary_bank().push_back(site); } // Set the delayed group on the particle as well - p.delayed_group_ = site.delayed_group; + p.delayed_group() = site.delayed_group; // Increment the number of neutrons born delayed - if (p.delayed_group_ > 0) { - nu_d[p.delayed_group_-1]++; + if (p.delayed_group() > 0) { + nu_d[p.delayed_group() - 1]++; } // Write fission particles to nuBank - p.nu_bank_.emplace_back(); - Particle::NuBank* nu_bank_entry = &p.nu_bank_.back(); + p.nu_bank().emplace_back(); + Particle::NuBank* nu_bank_entry = &p.nu_bank().back(); nu_bank_entry->wgt = site.wgt; nu_bank_entry->E = site.E; nu_bank_entry->delayed_group = site.delayed_group; @@ -227,7 +230,7 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // If shared fission bank was full, and no fissions could be added, // set the particle fission flag to false. if (nu == skipped) { - p.fission_ = false; + p.fission() = false; return; } @@ -236,10 +239,10 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) nu -= skipped; // Store the total weight banked for analog fission tallies - p.n_bank_ = nu; - p.wgt_bank_ = nu / weight; + p.n_bank() = nu; + p.wgt_bank() = nu / weight; for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { - p.n_delayed_bank_[d] = nu_d[d]; + p.n_delayed_bank(d) = nu_d[d]; } } @@ -249,19 +252,19 @@ void sample_photon_reaction(Particle& p) // photons with energy below the cutoff may have been produced by neutrons // reactions or atomic relaxation int photon = static_cast(Particle::Type::photon); - if (p.E_ < settings::energy_cutoff[photon]) { - p.E_ = 0.0; - p.alive_ = false; + if (p.E() < settings::energy_cutoff[photon]) { + p.E() = 0.0; + p.alive() = false; return; } // Sample element within material int i_element = sample_element(p); - const auto& micro {p.photon_xs_[i_element]}; + const auto& micro {p.photon_xs(i_element)}; const auto& element {*data::elements[i_element]}; // Calculate photon energy over electron rest mass equivalent - double alpha = p.E_/MASS_ELECTRON_EV; + double alpha = p.E() / MASS_ELECTRON_EV; // For tallying purposes, this routine might be called directly. In that // case, we need to sample a reaction via the cutoff variable @@ -273,8 +276,8 @@ void sample_photon_reaction(Particle& p) if (prob > cutoff) { double mu = element.rayleigh_scatter(alpha, p.current_seed()); p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed()); - p.event_ = TallyEvent::SCATTER; - p.event_mt_ = COHERENT; + p.event() = TallyEvent::SCATTER; + p.event_mt() = COHERENT; return; } @@ -302,7 +305,7 @@ void sample_photon_reaction(Particle& p) double mu_electron = (alpha - alpha_out*mu) / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu); Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed()); - p.create_secondary(p.wgt_, u, E_electron, Particle::Type::electron); + p.create_secondary(p.wgt(), u, E_electron, Particle::Type::electron); } // TODO: Compton subshell data does not match atomic relaxation data @@ -314,10 +317,10 @@ void sample_photon_reaction(Particle& p) } phi += PI; - p.E_ = alpha_out*MASS_ELECTRON_EV; + p.E() = alpha_out * MASS_ELECTRON_EV; p.u() = rotate_angle(p.u(), mu, &phi, p.current_seed()); - p.event_ = TallyEvent::SCATTER; - p.event_mt_ = INCOHERENT; + p.event() = TallyEvent::SCATTER; + p.event_mt() = INCOHERENT; return; } @@ -340,7 +343,7 @@ void sample_photon_reaction(Particle& p) prob += xs; if (prob > cutoff) { - double E_electron = p.E_ - shell.binding_energy; + double E_electron = p.E() - shell.binding_energy; // Sample mu using non-relativistic Sauter distribution. // See Eqns 3.19 and 3.20 in "Implementing a photon physics @@ -363,15 +366,15 @@ void sample_photon_reaction(Particle& p) u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); // Create secondary electron - p.create_secondary(p.wgt_, u, E_electron, Particle::Type::electron); + p.create_secondary(p.wgt(), u, E_electron, Particle::Type::electron); // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons element.atomic_relaxation(shell, p); - p.event_ = TallyEvent::ABSORB; - p.event_mt_ = 533 + shell.index_subshell; - p.alive_ = false; - p.E_ = 0.0; + p.event() = TallyEvent::ABSORB; + p.event_mt() = 533 + shell.index_subshell; + p.alive() = false; + p.E() = 0.0; return; } } @@ -388,16 +391,16 @@ void sample_photon_reaction(Particle& p) // Create secondary electron Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed()); - p.create_secondary(p.wgt_, u, E_electron, Particle::Type::electron); + p.create_secondary(p.wgt(), u, E_electron, Particle::Type::electron); // Create secondary positron u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed()); - p.create_secondary(p.wgt_, u, E_positron, Particle::Type::positron); + p.create_secondary(p.wgt(), u, E_positron, Particle::Type::positron); - p.event_ = TallyEvent::ABSORB; - p.event_mt_ = PAIR_PROD; - p.alive_ = false; - p.E_ = 0.0; + p.event() = TallyEvent::ABSORB; + p.event_mt() = PAIR_PROD; + p.alive() = false; + p.E() = 0.0; } } @@ -410,9 +413,9 @@ void sample_electron_reaction(Particle& p) thick_target_bremsstrahlung(p, &E_lost); } - p.E_ = 0.0; - p.alive_ = false; - p.event_ = TallyEvent::ABSORB; + p.E() = 0.0; + p.alive() = false; + p.event() = TallyEvent::ABSORB; } void sample_positron_reaction(Particle& p) @@ -433,21 +436,21 @@ void sample_positron_reaction(Particle& p) u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); // Create annihilation photon pair traveling in opposite directions - p.create_secondary(p.wgt_, u, MASS_ELECTRON_EV, Particle::Type::photon); - p.create_secondary(p.wgt_, -u, MASS_ELECTRON_EV, Particle::Type::photon); + p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, Particle::Type::photon); + p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, Particle::Type::photon); - p.E_ = 0.0; - p.alive_ = false; - p.event_ = TallyEvent::ABSORB; + p.E() = 0.0; + p.alive() = false; + p.event() = TallyEvent::ABSORB; } int sample_nuclide(Particle& p) { // Sample cumulative distribution function - double cutoff = prn(p.current_seed()) * p.macro_xs_.total; + double cutoff = prn(p.current_seed()) * p.macro_xs().total; // Get pointers to nuclide/density arrays - const auto& mat {model::materials[p.material_]}; + const auto& mat {model::materials[p.material()]}; int n = mat->nuclide_.size(); double prob = 0.0; @@ -457,7 +460,7 @@ int sample_nuclide(Particle& p) double atom_density = mat->atom_density_[i]; // Increment probability to compare to cutoff - prob += atom_density * p.neutron_xs_[i_nuclide].total; + prob += atom_density * p.neutron_xs(i_nuclide).total; if (prob >= cutoff) return i_nuclide; } @@ -469,10 +472,10 @@ int sample_nuclide(Particle& p) int sample_element(Particle& p) { // Sample cumulative distribution function - double cutoff = prn(p.current_seed()) * p.macro_xs_.total; + double cutoff = prn(p.current_seed()) * p.macro_xs().total; // Get pointers to elements, densities - const auto& mat {model::materials[p.material_]}; + const auto& mat {model::materials[p.material()]}; double prob = 0.0; for (int i = 0; i < mat->element_.size(); ++i) { @@ -481,13 +484,13 @@ int sample_element(Particle& p) double atom_density = mat->atom_density_[i]; // Determine microscopic cross section - double sigma = atom_density * p.photon_xs_[i_element].total; + double sigma = atom_density * p.photon_xs(i_element).total; // Increment probability to compare to cutoff prob += sigma; if (prob > cutoff) { // Save which nuclide particle had collision with for tally purpose - p.event_nuclide_ = mat->nuclide_[i]; + p.event_nuclide() = mat->nuclide_[i]; return i_element; } @@ -506,23 +509,23 @@ Reaction& sample_fission(int i_nuclide, Particle& p) // If we're in the URR, by default use the first fission reaction. We also // default to the first reaction if we know that there are no partial fission // reactions - if (p.neutron_xs_[i_nuclide].use_ptable || !nuc->has_partial_fission_) { + if (p.neutron_xs(i_nuclide).use_ptable || !nuc->has_partial_fission_) { return *nuc->fission_rx_[0]; } // Check to see if we are in a windowed multipole range. WMP only supports // the first fission reaction. if (nuc->multipole_) { - if (p.E_ >= nuc->multipole_->E_min_ && p.E_ <= nuc->multipole_->E_max_) { + if (p.E() >= nuc->multipole_->E_min_ && p.E() <= nuc->multipole_->E_max_) { return *nuc->fission_rx_[0]; } } // Get grid index and interpolatoin factor and sample fission cdf - int i_temp = p.neutron_xs_[i_nuclide].index_temp; - int i_grid = p.neutron_xs_[i_nuclide].index_grid; - double f = p.neutron_xs_[i_nuclide].interp_factor; - double cutoff = prn(p.current_seed()) * p.neutron_xs_[i_nuclide].fission; + int i_temp = p.neutron_xs(i_nuclide).index_temp; + int i_grid = p.neutron_xs(i_nuclide).index_grid; + double f = p.neutron_xs(i_nuclide).interp_factor; + double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).fission; double prob = 0.0; // Loop through each partial fission reaction type @@ -546,10 +549,10 @@ Reaction& sample_fission(int i_nuclide, Particle& p) void sample_photon_product(int i_nuclide, Particle& p, int* i_rx, int* i_product) { // Get grid index and interpolation factor and sample photon production cdf - int i_temp = p.neutron_xs_[i_nuclide].index_temp; - int i_grid = p.neutron_xs_[i_nuclide].index_grid; - double f = p.neutron_xs_[i_nuclide].interp_factor; - double cutoff = prn(p.current_seed()) * p.neutron_xs_[i_nuclide].photon_prod; + int i_temp = p.neutron_xs(i_nuclide).index_temp; + int i_grid = p.neutron_xs(i_nuclide).index_grid; + double f = p.neutron_xs(i_nuclide).interp_factor; + double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).photon_prod; double prob = 0.0; // Loop through each reaction type @@ -573,15 +576,15 @@ void sample_photon_product(int i_nuclide, Particle& p, int* i_rx, int* i_product if (settings::delayed_photon_scaling) { if (is_fission(rx->mt_)) { if (nuc->prompt_photons_ && nuc->delayed_photons_) { - double energy_prompt = (*nuc->prompt_photons_)(p.E_); - double energy_delayed = (*nuc->delayed_photons_)(p.E_); + double energy_prompt = (*nuc->prompt_photons_)(p.E()); + double energy_delayed = (*nuc->delayed_photons_)(p.E()); f = (energy_prompt + energy_delayed)/(energy_prompt); } } } // add to cumulative probability - prob += f * (*rx->products_[j].yield_)(p.E_) * xs; + prob += f * (*rx->products_[j].yield_)(p.E()) * xs; *i_rx = i; *i_product = j; @@ -595,31 +598,33 @@ void absorption(Particle& p, int i_nuclide) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p.wgt_absorb_ = p.wgt_ * p.neutron_xs_[i_nuclide].absorption / - p.neutron_xs_[i_nuclide].total; + p.wgt_absorb() = p.wgt() * p.neutron_xs(i_nuclide).absorption / + p.neutron_xs(i_nuclide).total; // Adjust weight of particle by probability of absorption - p.wgt_ -= p.wgt_absorb_; - p.wgt_last_ = p.wgt_; + p.wgt() -= p.wgt_absorb(); + p.wgt_last() = p.wgt(); // Score implicit absorption estimate of keff if (settings::run_mode == RunMode::EIGENVALUE) { - p.keff_tally_absorption_ += p.wgt_absorb_ * p.neutron_xs_[ - i_nuclide].nu_fission / p.neutron_xs_[i_nuclide].absorption; + p.keff_tally_absorption() += p.wgt_absorb() * + p.neutron_xs(i_nuclide).nu_fission / + p.neutron_xs(i_nuclide).absorption; } } else { // See if disappearance reaction happens - if (p.neutron_xs_[i_nuclide].absorption > - prn(p.current_seed()) * p.neutron_xs_[i_nuclide].total) { + if (p.neutron_xs(i_nuclide).absorption > + prn(p.current_seed()) * p.neutron_xs(i_nuclide).total) { // Score absorption estimate of keff if (settings::run_mode == RunMode::EIGENVALUE) { - p.keff_tally_absorption_ += p.wgt_ * p.neutron_xs_[ - i_nuclide].nu_fission / p.neutron_xs_[i_nuclide].absorption; + p.keff_tally_absorption() += p.wgt() * + p.neutron_xs(i_nuclide).nu_fission / + p.neutron_xs(i_nuclide).absorption; } - p.alive_ = false; - p.event_ = TallyEvent::ABSORB; - p.event_mt_ = N_DISAPPEAR; + p.alive() = false; + p.event() = TallyEvent::ABSORB; + p.event_mt() = N_DISAPPEAR; } } } @@ -631,7 +636,7 @@ void scatter(Particle& p, int i_nuclide) // Get pointer to nuclide and grid index/interpolation factor const auto& nuc {data::nuclides[i_nuclide]}; - const auto& micro {p.neutron_xs_[i_nuclide]}; + const auto& micro {p.neutron_xs(i_nuclide)}; int i_temp = micro.index_temp; int i_grid = micro.index_grid; double f = micro.interp_factor; @@ -652,12 +657,12 @@ void scatter(Particle& p, int i_nuclide) // NON-S(A,B) ELASTIC SCATTERING // Determine temperature - double kT = nuc->multipole_ ? p.sqrtkT_*p.sqrtkT_ : nuc->kTs_[i_temp]; + double kT = nuc->multipole_ ? p.sqrtkT() * p.sqrtkT() : nuc->kTs_[i_temp]; // Perform collision physics for elastic scattering elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p); - p.event_mt_ = ELASTIC; + p.event_mt() = ELASTIC; sampled = true; } @@ -668,7 +673,7 @@ void scatter(Particle& p, int i_nuclide) sab_scatter(i_nuclide, micro.index_sab, p); - p.event_mt_ = ELASTIC; + p.event_mt() = ELASTIC; sampled = true; } @@ -700,14 +705,14 @@ void scatter(Particle& p, int i_nuclide) // Perform collision physics for inelastic scattering const auto& rx {nuc->reactions_[i]}; inelastic_scatter(*nuc, *rx, p); - p.event_mt_ = rx->mt_; + p.event_mt() = rx->mt_; } // Set event component - p.event_ = TallyEvent::SCATTER; + p.event() = TallyEvent::SCATTER; // Sample new outgoing angle for isotropic-in-lab scattering - const auto& mat {model::materials[p.material_]}; + const auto& mat {model::materials[p.material()]}; if (!mat->p0_.empty()) { int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide]; if (mat->p0_[i_nuc_mat]) { @@ -719,7 +724,7 @@ void scatter(Particle& p, int i_nuclide) p.u().x = mu; p.u().y = std::sqrt(1.0 - mu*mu)*std::cos(phi); p.u().z = std::sqrt(1.0 - mu*mu)*std::sin(phi); - p.mu_ = u_old.dot(p.u()); + p.mu() = u_old.dot(p.u()); } } } @@ -730,7 +735,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, // get pointer to nuclide const auto& nuc {data::nuclides[i_nuclide]}; - double vel = std::sqrt(p.E_); + double vel = std::sqrt(p.E()); double awr = nuc->awr_; // Neutron velocity in LAB @@ -738,9 +743,9 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, // Sample velocity of target nucleus Direction v_t {}; - if (!p.neutron_xs_[i_nuclide].use_ptable) { - v_t = sample_target_velocity(*nuc, p.E_, p.u(), v_n, - p.neutron_xs_[i_nuclide].elastic, kT, p.current_seed()); + if (!p.neutron_xs(i_nuclide).use_ptable) { + v_t = sample_target_velocity(*nuc, p.E(), p.u(), v_n, + p.neutron_xs(i_nuclide).elastic, kT, p.current_seed()); } // Velocity of center-of-mass @@ -758,7 +763,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, auto& d = rx.products_[0].distribution_[0]; auto d_ = dynamic_cast(d.get()); if (!d_->angle().empty()) { - mu_cm = d_->angle().sample(p.E_, p.current_seed()); + mu_cm = d_->angle().sample(p.E(), p.current_seed()); } else { mu_cm = 2.0*prn(p.current_seed()) - 1.0; } @@ -774,12 +779,12 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, // Transform back to LAB frame v_n += v_cm; - p.E_ = v_n.dot(v_n); - vel = std::sqrt(p.E_); + p.E() = v_n.dot(v_n); + vel = std::sqrt(p.E()); // compute cosine of scattering angle in LAB frame by taking dot product of // neutron's pre- and post-collision angle - p.mu_ = p.u().dot(v_n) / vel; + p.mu() = p.u().dot(v_n) / vel; // Set energy and direction of particle in LAB frame p.u() = v_n / vel; @@ -787,22 +792,24 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, // Because of floating-point roundoff, it may be possible for mu_lab to be // outside of the range [-1,1). In these cases, we just set mu_lab to exactly // -1 or 1 - if (std::abs(p.mu_) > 1.0) p.mu_ = std::copysign(1.0, p.mu_); + if (std::abs(p.mu()) > 1.0) + p.mu() = std::copysign(1.0, p.mu()); } void sab_scatter(int i_nuclide, int i_sab, Particle& p) { // Determine temperature index - const auto& micro {p.neutron_xs_[i_nuclide]}; + const auto& micro {p.neutron_xs(i_nuclide)}; int i_temp = micro.index_temp_sab; // Sample energy and angle double E_out; - data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, p.E_, &E_out, &p.mu_, p.current_seed()); + data::thermal_scatt[i_sab]->data_[i_temp].sample( + micro, p.E(), &E_out, &p.mu(), p.current_seed()); // Set energy to outgoing, change direction of particle - p.E_ = E_out; - p.u() = rotate_angle(p.u(), p.mu_, nullptr, p.current_seed()); + p.E() = E_out; + p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); } Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u, @@ -1101,7 +1108,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) { // copy energy of neutron - double E_in = p.E_; + double E_in = p.E(); // sample outgoing energy and scattering cosine double E; @@ -1128,8 +1135,8 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu); // Set outgoing energy and scattering angle - p.E_ = E; - p.mu_ = mu; + p.E() = E; + p.mu() = mu; // change direction of particle p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed()); @@ -1139,19 +1146,19 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p) if (std::floor(yield) == yield) { // If yield is integral, create exactly that many secondary particles for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { - p.create_secondary(p.wgt_, p.u(), p.E_, Particle::Type::neutron); + p.create_secondary(p.wgt(), p.u(), p.E(), Particle::Type::neutron); } } else { // Otherwise, change weight of particle based on yield - p.wgt_ *= yield; + p.wgt() *= yield; } } void sample_secondary_photons(Particle& p, int i_nuclide) { // Sample the number of photons produced - double y_t = p.neutron_xs_[i_nuclide].photon_prod / - p.neutron_xs_[i_nuclide].total; + double y_t = + p.neutron_xs(i_nuclide).photon_prod / p.neutron_xs(i_nuclide).total; int y = static_cast(y_t); if (prn(p.current_seed()) <= y_t - y) ++y; @@ -1166,7 +1173,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide) auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; double E; double mu; - rx->products_[i_product].sample(p.E_, E, mu, p.current_seed()); + rx->products_[i_product].sample(p.E(), E, mu, p.current_seed()); // Sample the new direction Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed()); @@ -1178,9 +1185,9 @@ void sample_secondary_photons(Particle& p, int i_nuclide) // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. double wgt; if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) { - wgt = simulation::keff * p.wgt_; + wgt = simulation::keff * p.wgt(); } else { - wgt = p.wgt_; + wgt = p.wgt(); } // Create the secondary photon diff --git a/src/physics_common.cpp b/src/physics_common.cpp index 32e00056d..2a63b188d 100644 --- a/src/physics_common.cpp +++ b/src/physics_common.cpp @@ -11,14 +11,14 @@ namespace openmc { void russian_roulette(Particle& p) { - if (p.wgt_ < settings::weight_cutoff) { - if (prn(p.current_seed()) < p.wgt_ / settings::weight_survive) { - p.wgt_ = settings::weight_survive; - p.wgt_last_ = p.wgt_; + if (p.wgt() < settings::weight_cutoff) { + if (prn(p.current_seed()) < p.wgt() / settings::weight_survive) { + p.wgt() = settings::weight_survive; + p.wgt_last() = p.wgt(); } else { - p.wgt_ = 0.; - p.wgt_last_ = 0.; - p.alive_ = false; + p.wgt() = 0.; + p.wgt_last() = 0.; + p.alive() = false; } } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index b2dd4078e..bd4966dd8 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -25,14 +25,14 @@ void collision_mg(Particle& p) { // Add to the collision counter for the particle - p.n_collision_++; + p.n_collision()++; // Sample the reaction type sample_reaction(p); // Display information about collision - if ((settings::verbosity >= 10) || p.trace_) { - write_message(fmt::format(" Energy Group = {}", p.g_), 1); + if ((settings::verbosity >= 10) || p.trace()) { + write_message(fmt::format(" Energy Group = {}", p.g()), 1); } } @@ -44,7 +44,7 @@ sample_reaction(Particle& p) // change when sampling fission sites. The following block handles all // absorption (including fission) - if (model::materials[p.material_]->fissionable_) { + if (model::materials[p.material()]->fissionable_) { if (settings::run_mode == RunMode::EIGENVALUE || (settings::run_mode == RunMode::FIXED_SOURCE && settings::create_fission_neutrons)) { @@ -54,12 +54,13 @@ sample_reaction(Particle& p) // If survival biasing is being used, the following subroutine adjusts the // weight of the particle. Otherwise, it checks to see if absorption occurs. - if (p.macro_xs_.absorption > 0.) { + if (p.macro_xs().absorption > 0.) { absorption(p); } else { - p.wgt_absorb_ = 0.; + p.wgt_absorb() = 0.; } - if (!p.alive_) return; + if (!p.alive()) + return; // Sample a scattering event to determine the energy of the exiting neutron scatter(p); @@ -67,24 +68,25 @@ sample_reaction(Particle& p) // Play Russian roulette if survival biasing is turned on if (settings::survival_biasing) { russian_roulette(p); - if (!p.alive_) return; + if (!p.alive()) + return; } } void scatter(Particle& p) { - data::mg.macro_xs_[p.material_].sample_scatter(p.g_last_, p.g_, p.mu_, - p.wgt_, p.current_seed()); + data::mg.macro_xs_[p.material()].sample_scatter( + p.g_last(), p.g(), p.mu(), p.wgt(), p.current_seed()); // Rotate the angle - p.u() = rotate_angle(p.u(), p.mu_, nullptr, p.current_seed()); + p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed()); // Update energy value for downstream compatability (in tallying) - p.E_ = data::mg.energy_bin_avg_[p.g_]; + p.E() = data::mg.energy_bin_avg_[p.g()]; // Set event component - p.event_ = TallyEvent::SCATTER; + p.event() = TallyEvent::SCATTER; } void @@ -95,8 +97,8 @@ create_fission_sites(Particle& p) double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p.wgt_ / simulation::keff * weight * - p.macro_xs_.nu_fission / p.macro_xs_.total; + double nu_t = p.wgt() / simulation::keff * weight * p.macro_xs().nu_fission / + p.macro_xs().total; // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -113,9 +115,9 @@ create_fission_sites(Particle& p) double nu_d[MAX_DELAYED_GROUPS] = {0.}; // Clear out particle's nu fission bank - p.nu_bank_.clear(); + p.nu_bank().clear(); - p.fission_ = true; + p.fission() = true; int skipped = 0; // Determine whether to place fission sites into the shared fission bank @@ -128,8 +130,8 @@ create_fission_sites(Particle& p) site.r = p.r(); site.particle = Particle::Type::neutron; site.wgt = 1. / weight; - site.parent_id = p.id_; - site.progeny_id = p.n_progeny_++; + site.parent_id = p.id(); + site.progeny_id = p.n_progeny()++; // Sample the cosine of the angle, assuming fission neutrons are emitted // isotropically @@ -144,8 +146,8 @@ create_fission_sites(Particle& p) // Sample secondary energy distribution for the fission reaction int dg; int gout; - data::mg.macro_xs_[p.material_].sample_fission_energy(p.g_, dg, gout, - p.current_seed()); + data::mg.macro_xs_[p.material()].sample_fission_energy( + p.g(), dg, gout, p.current_seed()); // Store the energy and delayed groups on the fission bank site.E = gout; @@ -164,20 +166,20 @@ create_fission_sites(Particle& p) break; } } else { - p.secondary_bank_.push_back(site); + p.secondary_bank().push_back(site); } // Set the delayed group on the particle as well - p.delayed_group_ = dg + 1; + p.delayed_group() = dg + 1; // Increment the number of neutrons born delayed - if (p.delayed_group_ > 0) { + if (p.delayed_group() > 0) { nu_d[dg]++; } // Write fission particles to nuBank - p.nu_bank_.emplace_back(); - Particle::NuBank* nu_bank_entry = &p.nu_bank_.back(); + p.nu_bank().emplace_back(); + Particle::NuBank* nu_bank_entry = &p.nu_bank().back(); nu_bank_entry->wgt = site.wgt; nu_bank_entry->E = site.E; nu_bank_entry->delayed_group = site.delayed_group; @@ -186,7 +188,7 @@ create_fission_sites(Particle& p) // If shared fission bank was full, and no fissions could be added, // set the particle fission flag to false. if (nu == skipped) { - p.fission_ = false; + p.fission() = false; return; } @@ -195,10 +197,10 @@ create_fission_sites(Particle& p) nu -= skipped; // Store the total weight banked for analog fission tallies - p.n_bank_ = nu; - p.wgt_bank_ = nu / weight; + p.n_bank() = nu; + p.wgt_bank() = nu / weight; for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { - p.n_delayed_bank_[d] = nu_d[d]; + p.n_delayed_bank(d) = nu_d[d]; } } @@ -207,23 +209,22 @@ absorption(Particle& p) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p.wgt_absorb_ = p.wgt_ * p.macro_xs_.absorption / p.macro_xs_.total; + p.wgt_absorb() = p.wgt() * p.macro_xs().absorption / p.macro_xs().total; // Adjust weight of particle by the probability of absorption - p.wgt_ -= p.wgt_absorb_; - p.wgt_last_ = p.wgt_; + p.wgt() -= p.wgt_absorb(); + p.wgt_last() = p.wgt(); // Score implicit absorpion estimate of keff - p.keff_tally_absorption_ += p.wgt_absorb_ * p.macro_xs_.nu_fission / - p.macro_xs_.absorption; + p.keff_tally_absorption() += + p.wgt_absorb() * p.macro_xs().nu_fission / p.macro_xs().absorption; } else { - if (p.macro_xs_.absorption > prn(p.current_seed()) * p.macro_xs_.total) { - p.keff_tally_absorption_ += p.wgt_ * p.macro_xs_.nu_fission / - p.macro_xs_.absorption; - p.alive_ = false; - p.event_ = TallyEvent::ABSORB; + if (p.macro_xs().absorption > prn(p.current_seed()) * p.macro_xs().total) { + p.keff_tally_absorption() += + p.wgt() * p.macro_xs().nu_fission / p.macro_xs().absorption; + p.alive() = false; + p.event() = TallyEvent::ABSORB; } - } } diff --git a/src/plot.cpp b/src/plot.cpp index c7739c2c6..9333a6fd3 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -42,19 +42,19 @@ IdData::IdData(size_t h_res, size_t v_res) void IdData::set_value(size_t y, size_t x, const Particle& p, int level) { // set cell data - if (p.n_coord_ <= level) { + if (p.n_coord() <= level) { data_(y, x, 0) = NOT_FOUND; } else { - data_(y, x, 0) = model::cells.at(p.coord_.at(level).cell)->id_; + data_(y, x, 0) = model::cells.at(p.coord(level).cell)->id_; } // set material data - Cell* c = model::cells.at(p.coord_.at(p.n_coord_ - 1).cell).get(); - if (p.material_ == MATERIAL_VOID) { + Cell* c = model::cells.at(p.coord(p.n_coord() - 1).cell).get(); + if (p.material() == MATERIAL_VOID) { data_(y, x, 1) = MATERIAL_VOID; return; } else if (c->type_ == Fill::MATERIAL) { - Material* m = model::materials.at(p.material_).get(); + Material* m = model::materials.at(p.material()).get(); data_(y, x, 1) = m->id_; } } @@ -69,10 +69,10 @@ PropertyData::PropertyData(size_t h_res, size_t v_res) void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level) { - Cell* c = model::cells.at(p.coord_.at(p.n_coord_ - 1).cell).get(); - data_(y,x,0) = (p.sqrtkT_ * p.sqrtkT_) / K_BOLTZMANN; - if (c->type_ != Fill::UNIVERSE && p.material_ != MATERIAL_VOID) { - Material* m = model::materials.at(p.material_).get(); + Cell* c = model::cells.at(p.coord(p.n_coord() - 1).cell).get(); + data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN; + if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) { + Material* m = model::materials.at(p.material()).get(); data_(y,x,1) = m->density_gpcc_; } } diff --git a/src/simulation.cpp b/src/simulation.cpp index 21630ad25..d5fdfccf4 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -464,51 +464,52 @@ void initialize_history(Particle& p, int64_t index_source) auto site = sample_external_source(&seed); p.from_source(&site); } - p.current_work_ = index_source; + p.current_work() = index_source; // set identifier for particle - p.id_ = simulation::work_index[mpi::rank] + index_source; + p.id() = simulation::work_index[mpi::rank] + index_source; // set progeny count to zero - p.n_progeny_ = 0; + p.n_progeny() = 0; // Reset particle event counter - p.n_event_ = 0; + p.n_event() = 0; // set random number seed - int64_t particle_seed = (simulation::total_gen + overall_generation() - 1) - * settings::n_particles + p.id_; - init_particle_seeds(particle_seed, p.seeds_); + int64_t particle_seed = + (simulation::total_gen + overall_generation() - 1) * settings::n_particles + + p.id(); + init_particle_seeds(particle_seed, p.seeds()); // set particle trace - p.trace_ = false; + p.trace() = false; if (simulation::current_batch == settings::trace_batch && simulation::current_gen == settings::trace_gen && - p.id_ == settings::trace_particle) p.trace_ = true; + p.id() == settings::trace_particle) + p.trace() = true; // Set particle track. - p.write_track_ = false; + p.write_track() = false; if (settings::write_all_tracks) { - p.write_track_ = true; + p.write_track() = true; } else if (settings::track_identifiers.size() > 0) { for (const auto& t : settings::track_identifiers) { if (simulation::current_batch == t[0] && - simulation::current_gen == t[1] && - p.id_ == t[2]) { - p.write_track_ = true; + simulation::current_gen == t[1] && p.id() == t[2]) { + p.write_track() = true; break; } } } // Display message if high verbosity or trace is on - if (settings::verbosity >= 9 || p.trace_) { - write_message("Simulating Particle {}", p.id_); + if (settings::verbosity >= 9 || p.trace()) { + write_message("Simulating Particle {}", p.id()); } // Add paricle's starting weight to count for normalizing tallies later #pragma omp atomic - simulation::total_weight += p.wgt_; + simulation::total_weight += p.wgt(); initialize_history_partial(p); } @@ -517,21 +518,22 @@ void initialize_history_partial(Particle& p) { // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { - for (auto& micro : p.neutron_xs_) micro.last_E = 0.0; + p.invalidate_neutron_xs(); } // Prepare to write out particle track. - if (p.write_track_) add_particle_track(p); + if (p.write_track()) + add_particle_track(p); // Every particle starts with no accumulated flux derivative. if (!model::active_tallies.empty()) { - p.flux_derivs_.resize(model::tally_derivs.size(), 0.0); - std::fill(p.flux_derivs_.begin(), p.flux_derivs_.end(), 0.0); + p.flux_derivs().resize(model::tally_derivs.size(), 0.0); + std::fill(p.flux_derivs().begin(), p.flux_derivs().end(), 0.0); } // Allocate space for tally filter matches - p.filter_matches_.resize(model::tally_filters.size()); + p.filter_matches().resize(model::tally_filters.size()); } int overall_generation() @@ -678,13 +680,13 @@ void transport_history_based_single_particle(Particle& p) while (true) { p.event_calculate_xs(); p.event_advance(); - if (p.collision_distance_ > p.boundary_.distance) { + if (p.collision_distance() > p.boundary().distance) { p.event_cross_surface(); } else { p.event_collide(); } p.event_revive_from_secondary(); - if (!p.alive_) + if (!p.alive()) break; } p.event_death(); diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index fecfa3f1c..1a950c8d6 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -109,17 +109,17 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, // perturbated variable. const auto& deriv {model::tally_derivs[tally.deriv_]}; - const auto flux_deriv = p.flux_derivs_[tally.deriv_]; + const auto flux_deriv = p.flux_derivs()[tally.deriv_]; // Handle special cases where we know that d_c/d_p must be zero. if (score_bin == SCORE_FLUX) { score *= flux_deriv; return; - } else if (p.material_ == MATERIAL_VOID) { + } else if (p.material() == MATERIAL_VOID) { score *= flux_deriv; return; } - const Material& material {*model::materials[p.material_]}; + const Material& material {*model::materials[p.material()]}; if (material.id_ != deriv.diff_material) { score *= flux_deriv; return; @@ -179,7 +179,7 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, switch (tally.estimator_) { case TallyEstimator::ANALOG: - if (p.event_nuclide_ != deriv.diff_nuclide) { + if (p.event_nuclide() != deriv.diff_nuclide) { score *= flux_deriv; return; } @@ -210,12 +210,11 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (i_nuclide == -1 && p.macro_xs_.total > 0.0) { - score *= flux_deriv - + p.neutron_xs_[deriv.diff_nuclide].total - / p.macro_xs_.total; - } else if (i_nuclide == deriv.diff_nuclide - && p.neutron_xs_[i_nuclide].total) { + if (i_nuclide == -1 && p.macro_xs().total > 0.0) { + score *= flux_deriv + + p.neutron_xs(deriv.diff_nuclide).total / p.macro_xs().total; + } else if (i_nuclide == deriv.diff_nuclide && + p.neutron_xs(i_nuclide).total) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -223,13 +222,12 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, break; case SCORE_SCATTER: - if (i_nuclide == -1 && (p.macro_xs_.total - - p.macro_xs_.absorption) > 0.0) { - score *= flux_deriv - + (p.neutron_xs_[deriv.diff_nuclide].total - - p.neutron_xs_[deriv.diff_nuclide].absorption) - / (p.macro_xs_.total - - p.macro_xs_.absorption); + if (i_nuclide == -1 && + (p.macro_xs().total - p.macro_xs().absorption) > 0.0) { + score *= + flux_deriv + (p.neutron_xs(deriv.diff_nuclide).total - + p.neutron_xs(deriv.diff_nuclide).absorption) / + (p.macro_xs().total - p.macro_xs().absorption); } else if (i_nuclide == deriv.diff_nuclide) { score *= flux_deriv + 1. / atom_density; } else { @@ -238,12 +236,11 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, break; case SCORE_ABSORPTION: - if (i_nuclide == -1 && p.macro_xs_.absorption > 0.0) { - score *= flux_deriv - + p.neutron_xs_[deriv.diff_nuclide].absorption - / p.macro_xs_.absorption; - } else if (i_nuclide == deriv.diff_nuclide - && p.neutron_xs_[i_nuclide].absorption) { + if (i_nuclide == -1 && p.macro_xs().absorption > 0.0) { + score *= flux_deriv + p.neutron_xs(deriv.diff_nuclide).absorption / + p.macro_xs().absorption; + } else if (i_nuclide == deriv.diff_nuclide && + p.neutron_xs(i_nuclide).absorption) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -251,12 +248,11 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, break; case SCORE_FISSION: - if (i_nuclide == -1 && p.macro_xs_.fission > 0.0) { - score *= flux_deriv - + p.neutron_xs_[deriv.diff_nuclide].fission - / p.macro_xs_.fission; - } else if (i_nuclide == deriv.diff_nuclide - && p.neutron_xs_[i_nuclide].fission) { + if (i_nuclide == -1 && p.macro_xs().fission > 0.0) { + score *= flux_deriv + p.neutron_xs(deriv.diff_nuclide).fission / + p.macro_xs().fission; + } else if (i_nuclide == deriv.diff_nuclide && + p.neutron_xs(i_nuclide).fission) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -264,12 +260,11 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, break; case SCORE_NU_FISSION: - if (i_nuclide == -1 && p.macro_xs_.nu_fission > 0.0) { - score *= flux_deriv - + p.neutron_xs_[deriv.diff_nuclide].nu_fission - / p.macro_xs_.nu_fission; - } else if (i_nuclide == deriv.diff_nuclide - && p.neutron_xs_[i_nuclide].nu_fission) { + if (i_nuclide == -1 && p.macro_xs().nu_fission > 0.0) { + score *= flux_deriv + p.neutron_xs(deriv.diff_nuclide).nu_fission / + p.macro_xs().nu_fission; + } else if (i_nuclide == deriv.diff_nuclide && + p.neutron_xs(i_nuclide).nu_fission) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -309,10 +304,11 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, // Find the index of the event nuclide. int i; for (i = 0; i < material.nuclide_.size(); ++i) - if (material.nuclide_[i] == p.event_nuclide_) break; + if (material.nuclide_[i] == p.event_nuclide()) + break; - const auto& nuc {*data::nuclides[p.event_nuclide_]}; - if (!multipole_in_range(nuc, p.E_last_)) { + const auto& nuc {*data::nuclides[p.event_nuclide()]}; + if (!multipole_in_range(nuc, p.E_last())) { score *= flux_deriv; break; } @@ -320,63 +316,65 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (p.neutron_xs_[p.event_nuclide_].total) { + if (p.neutron_xs(p.event_nuclide()).total) { double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i) - / p.macro_xs_.total; + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= flux_deriv + (dsig_s + dsig_a) * + material.atom_density_(i) / + p.macro_xs().total; } else { score *= flux_deriv; } break; case SCORE_SCATTER: - if (p.neutron_xs_[p.event_nuclide_].total - - p.neutron_xs_[p.event_nuclide_].absorption) { + if (p.neutron_xs(p.event_nuclide()).total - + p.neutron_xs(p.event_nuclide()).absorption) { double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv + dsig_s * material.atom_density_(i) - / (p.macro_xs_.total - p.macro_xs_.absorption); + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= + flux_deriv + dsig_s * material.atom_density_(i) / + (p.macro_xs().total - p.macro_xs().absorption); } else { score *= flux_deriv; } break; case SCORE_ABSORPTION: - if (p.neutron_xs_[p.event_nuclide_].absorption) { + if (p.neutron_xs(p.event_nuclide()).absorption) { double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv + dsig_a * material.atom_density_(i) - / p.macro_xs_.absorption; + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= flux_deriv + dsig_a * material.atom_density_(i) / + p.macro_xs().absorption; } else { score *= flux_deriv; } break; case SCORE_FISSION: - if (p.neutron_xs_[p.event_nuclide_].fission) { + if (p.neutron_xs(p.event_nuclide()).fission) { double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv + dsig_f * material.atom_density_(i) - / p.macro_xs_.fission; + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= flux_deriv + + dsig_f * material.atom_density_(i) / p.macro_xs().fission; } else { score *= flux_deriv; } break; case SCORE_NU_FISSION: - if (p.neutron_xs_[p.event_nuclide_].fission) { - double nu = p.neutron_xs_[p.event_nuclide_].nu_fission - / p.neutron_xs_[p.event_nuclide_].fission; + if (p.neutron_xs(p.event_nuclide()).fission) { + double nu = p.neutron_xs(p.event_nuclide()).nu_fission / + p.neutron_xs(p.event_nuclide()).fission; double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv + nu * dsig_f * material.atom_density_(i) - / p.macro_xs_.nu_fission; + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= flux_deriv + nu * dsig_f * material.atom_density_(i) / + p.macro_xs().nu_fission; } else { score *= flux_deriv; } @@ -392,7 +390,7 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, case TallyEstimator::COLLISION: if (i_nuclide != -1) { const auto& nuc {data::nuclides[i_nuclide]}; - if (!multipole_in_range(*nuc, p.E_last_)) { + if (!multipole_in_range(*nuc, p.E_last())) { score *= flux_deriv; return; } @@ -401,141 +399,136 @@ apply_derivative_to_score(const Particle& p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (i_nuclide == -1 && p.macro_xs_.total > 0.0) { + if (i_nuclide == -1 && p.macro_xs().total > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(nuc, p.E_last_) - && p.neutron_xs_[i_nuc].total) { + if (multipole_in_range(nuc, p.E_last()) && + p.neutron_xs(i_nuc).total) { double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / p.macro_xs_.total; - } else if (p.neutron_xs_[i_nuclide].total) { + score *= flux_deriv + cum_dsig / p.macro_xs().total; + } else if (p.neutron_xs(i_nuclide).total) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv - + (dsig_s + dsig_a) / p.neutron_xs_[i_nuclide].total; + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= + flux_deriv + (dsig_s + dsig_a) / p.neutron_xs(i_nuclide).total; } else { score *= flux_deriv; } break; case SCORE_SCATTER: - if (i_nuclide == -1 && (p.macro_xs_.total - - p.macro_xs_.absorption)) { + if (i_nuclide == -1 && (p.macro_xs().total - p.macro_xs().absorption)) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(nuc, p.E_last_) - && (p.neutron_xs_[i_nuc].total - - p.neutron_xs_[i_nuc].absorption)) { + if (multipole_in_range(nuc, p.E_last()) && + (p.neutron_xs(i_nuc).total - p.neutron_xs(i_nuc).absorption)) { double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); cum_dsig += dsig_s * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / (p.macro_xs_.total - - p.macro_xs_.absorption); - } else if (p.neutron_xs_[i_nuclide].total - - p.neutron_xs_[i_nuclide].absorption) { + score *= flux_deriv + + cum_dsig / (p.macro_xs().total - p.macro_xs().absorption); + } else if (p.neutron_xs(i_nuclide).total - + p.neutron_xs(i_nuclide).absorption) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv + dsig_s / (p.neutron_xs_[i_nuclide].total - - p.neutron_xs_[i_nuclide].absorption); + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= flux_deriv + dsig_s / (p.neutron_xs(i_nuclide).total - + p.neutron_xs(i_nuclide).absorption); } else { score *= flux_deriv; } break; case SCORE_ABSORPTION: - if (i_nuclide == -1 && p.macro_xs_.absorption > 0.0) { + if (i_nuclide == -1 && p.macro_xs().absorption > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(nuc, p.E_last_) - && p.neutron_xs_[i_nuc].absorption) { + if (multipole_in_range(nuc, p.E_last()) && + p.neutron_xs(i_nuc).absorption) { double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); cum_dsig += dsig_a * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / p.macro_xs_.absorption; - } else if (p.neutron_xs_[i_nuclide].absorption) { + score *= flux_deriv + cum_dsig / p.macro_xs().absorption; + } else if (p.neutron_xs(i_nuclide).absorption) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv - + dsig_a / p.neutron_xs_[i_nuclide].absorption; + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= flux_deriv + dsig_a / p.neutron_xs(i_nuclide).absorption; } else { score *= flux_deriv; } break; case SCORE_FISSION: - if (i_nuclide == -1 && p.macro_xs_.fission > 0.0) { + if (i_nuclide == -1 && p.macro_xs().fission > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(nuc, p.E_last_) - && p.neutron_xs_[i_nuc].fission) { + if (multipole_in_range(nuc, p.E_last()) && + p.neutron_xs(i_nuc).fission) { double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); cum_dsig += dsig_f * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / p.macro_xs_.fission; - } else if (p.neutron_xs_[i_nuclide].fission) { + score *= flux_deriv + cum_dsig / p.macro_xs().fission; + } else if (p.neutron_xs(i_nuclide).fission) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv - + dsig_f / p.neutron_xs_[i_nuclide].fission; + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= flux_deriv + dsig_f / p.neutron_xs(i_nuclide).fission; } else { score *= flux_deriv; } break; case SCORE_NU_FISSION: - if (i_nuclide == -1 && p.macro_xs_.nu_fission > 0.0) { + if (i_nuclide == -1 && p.macro_xs().nu_fission > 0.0) { double cum_dsig = 0; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto i_nuc = material.nuclide_[i]; const auto& nuc {*data::nuclides[i_nuc]}; - if (multipole_in_range(nuc, p.E_last_) - && p.neutron_xs_[i_nuc].fission) { - double nu = p.neutron_xs_[i_nuc].nu_fission - / p.neutron_xs_[i_nuc].fission; + if (multipole_in_range(nuc, p.E_last()) && + p.neutron_xs(i_nuc).fission) { + double nu = + p.neutron_xs(i_nuc).nu_fission / p.neutron_xs(i_nuc).fission; double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); cum_dsig += nu * dsig_f * material.atom_density_(i); } } - score *= flux_deriv + cum_dsig / p.macro_xs_.nu_fission; - } else if (p.neutron_xs_[i_nuclide].fission) { + score *= flux_deriv + cum_dsig / p.macro_xs().nu_fission; + } else if (p.neutron_xs(i_nuclide).fission) { const auto& nuc {*data::nuclides[i_nuclide]}; double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); - score *= flux_deriv - + dsig_f / p.neutron_xs_[i_nuclide].fission; + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); + score *= flux_deriv + dsig_f / p.neutron_xs(i_nuclide).fission; } else { score *= flux_deriv; } @@ -558,13 +551,14 @@ void score_track_derivative(Particle& p, double distance) { // A void material cannot be perturbed so it will not affect flux derivatives. - if (p.material_ == MATERIAL_VOID) return; + if (p.material() == MATERIAL_VOID) + return; - const Material& material {*model::materials[p.material_]}; + const Material& material {*model::materials[p.material()]}; for (auto idx = 0; idx < model::tally_derivs.size(); idx++) { const auto& deriv = model::tally_derivs[idx]; - auto& flux_deriv = p.flux_derivs_[idx]; + auto& flux_deriv = p.flux_derivs()[idx]; if (deriv.diff_material != material.id_) continue; switch (deriv.variable) { @@ -573,28 +567,26 @@ score_track_derivative(Particle& p, double distance) // phi is proportional to e^(-Sigma_tot * dist) // (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist // (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist - flux_deriv -= distance * p.macro_xs_.total - / material.density_gpcc_; + flux_deriv -= distance * p.macro_xs().total / material.density_gpcc_; break; case DerivativeVariable::NUCLIDE_DENSITY: // phi is proportional to e^(-Sigma_tot * dist) // (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist // (1 / phi) * (d_phi / d_N) = - sigma_tot * dist - flux_deriv -= distance - * p.neutron_xs_[deriv.diff_nuclide].total; + flux_deriv -= distance * p.neutron_xs(deriv.diff_nuclide).total; break; case DerivativeVariable::TEMPERATURE: for (auto i = 0; i < material.nuclide_.size(); ++i) { const auto& nuc {*data::nuclides[material.nuclide_[i]]}; - if (multipole_in_range(nuc, p.E_last_)) { + if (multipole_in_range(nuc, p.E_last())) { // phi is proportional to e^(-Sigma_tot * dist) // (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist // (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_, p.sqrtkT_); + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E(), p.sqrtkT()); flux_deriv -= distance * (dsig_s + dsig_a) * material.atom_density_(i); } @@ -607,13 +599,14 @@ score_track_derivative(Particle& p, double distance) void score_collision_derivative(Particle& p) { // A void material cannot be perturbed so it will not affect flux derivatives. - if (p.material_ == MATERIAL_VOID) return; + if (p.material() == MATERIAL_VOID) + return; - const Material& material {*model::materials[p.material_]}; + const Material& material {*model::materials[p.material()]}; for (auto idx = 0; idx < model::tally_derivs.size(); idx++) { const auto& deriv = model::tally_derivs[idx]; - auto& flux_deriv = p.flux_derivs_[idx]; + auto& flux_deriv = p.flux_derivs()[idx]; if (deriv.diff_material != material.id_) continue; @@ -627,7 +620,8 @@ void score_collision_derivative(Particle& p) break; case DerivativeVariable::NUCLIDE_DENSITY: - if (p.event_nuclide_ != deriv.diff_nuclide) continue; + if (p.event_nuclide() != deriv.diff_nuclide) + continue; // Find the index in this material for the diff_nuclide. int i; for (i = 0; i < material.nuclide_.size(); ++i) @@ -649,14 +643,14 @@ void score_collision_derivative(Particle& p) // Loop over the material's nuclides until we find the event nuclide. for (auto i_nuc : material.nuclide_) { const auto& nuc {*data::nuclides[i_nuc]}; - if (i_nuc == p.event_nuclide_ && multipole_in_range(nuc, p.E_last_)) { + if (i_nuc == p.event_nuclide() && multipole_in_range(nuc, p.E_last())) { // phi is proportional to Sigma_s // (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s // (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s - const auto& micro_xs {p.neutron_xs_[i_nuc]}; + const auto& micro_xs {p.neutron_xs(i_nuc)}; double dsig_s, dsig_a, dsig_f; - std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p.E_last_, p.sqrtkT_); + std::tie(dsig_s, dsig_a, dsig_f) = + nuc.multipole_->evaluate_deriv(p.E_last(), p.sqrtkT()); flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption); // Note that this is an approximation! The real scattering cross // section is diff --git a/src/tallies/filter_azimuthal.cpp b/src/tallies/filter_azimuthal.cpp index 4c64587ee..7bf5343b7 100644 --- a/src/tallies/filter_azimuthal.cpp +++ b/src/tallies/filter_azimuthal.cpp @@ -54,7 +54,7 @@ void AzimuthalFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - Direction u = (estimator == TallyEstimator::TRACKLENGTH) ? p.u() : p.u_last_; + Direction u = (estimator == TallyEstimator::TRACKLENGTH) ? p.u() : p.u_last(); double phi = std::atan2(u.y, u.x); if (phi >= bins_.front() && phi <= bins_.back()) { diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp index f6b804ab9..861a11aa4 100644 --- a/src/tallies/filter_cell.cpp +++ b/src/tallies/filter_cell.cpp @@ -49,8 +49,8 @@ void CellFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - for (int i = 0; i < p.n_coord_; i++) { - auto search = map_.find(p.coord_[i].cell); + for (int i = 0; i < p.n_coord(); i++) { + auto search = map_.find(p.coord(i).cell); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_cell_instance.cpp b/src/tallies/filter_cell_instance.cpp index d44cdfcf6..82af5bb9f 100644 --- a/src/tallies/filter_cell_instance.cpp +++ b/src/tallies/filter_cell_instance.cpp @@ -69,8 +69,8 @@ void CellInstanceFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - gsl::index index_cell = p.coord_[p.n_coord_ - 1].cell; - gsl::index instance = p.cell_instance_; + gsl::index index_cell = p.coord(p.n_coord() - 1).cell; + gsl::index instance = p.cell_instance(); auto search = map_.find({index_cell, instance}); if (search != map_.end()) { int index_bin = search->second; diff --git a/src/tallies/filter_cellborn.cpp b/src/tallies/filter_cellborn.cpp index 572b1503c..ac5b8be9d 100644 --- a/src/tallies/filter_cellborn.cpp +++ b/src/tallies/filter_cellborn.cpp @@ -8,7 +8,7 @@ void CellbornFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - auto search = map_.find(p.cell_born_); + auto search = map_.find(p.cell_born()); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_cellfrom.cpp b/src/tallies/filter_cellfrom.cpp index 39fde3976..9d6af07dd 100644 --- a/src/tallies/filter_cellfrom.cpp +++ b/src/tallies/filter_cellfrom.cpp @@ -8,8 +8,8 @@ void CellFromFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - for (int i = 0; i < p.n_coord_last_; i++) { - auto search = map_.find(p.cell_last_[i]); + for (int i = 0; i < p.n_coord_last(); i++) { + auto search = map_.find(p.cell_last(i)); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp index d3a2ffe5f..ea6fef944 100644 --- a/src/tallies/filter_distribcell.cpp +++ b/src/tallies/filter_distribcell.cpp @@ -43,8 +43,8 @@ DistribcellFilter::get_all_bins(const Particle& p, TallyEstimator estimator, { int offset = 0; auto distribcell_index = model::cells[cell_]->distribcell_index_; - for (int i = 0; i < p.n_coord_; i++) { - auto& c {*model::cells[p.coord_[i].cell]}; + for (int i = 0; i < p.n_coord(); i++) { + auto& c {*model::cells[p.coord(i).cell]}; if (c.type_ == Fill::UNIVERSE) { offset += c.offset_[distribcell_index]; } else if (c.type_ == Fill::LATTICE) { @@ -56,7 +56,7 @@ DistribcellFilter::get_all_bins(const Particle& p, TallyEstimator estimator, offset += lat.offset(distribcell_index, i_xyz); } } - if (cell_ == p.coord_[i].cell) { + if (cell_ == p.coord(i).cell) { match.bins_.push_back(offset); match.weights_.push_back(1.0); return; diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp index 097e8429e..9169c4d85 100644 --- a/src/tallies/filter_energy.cpp +++ b/src/tallies/filter_energy.cpp @@ -61,17 +61,17 @@ void EnergyFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - if (p.g_ != F90_NONE && matches_transport_groups_) { + if (p.g() != F90_NONE && matches_transport_groups_) { if (estimator == TallyEstimator::TRACKLENGTH) { - match.bins_.push_back(data::mg.num_energy_groups_ - p.g_ - 1); + match.bins_.push_back(data::mg.num_energy_groups_ - p.g() - 1); } else { - match.bins_.push_back(data::mg.num_energy_groups_ - p.g_last_ - 1); + match.bins_.push_back(data::mg.num_energy_groups_ - p.g_last() - 1); } match.weights_.push_back(1.0); } else { // Get the pre-collision energy of the particle. - auto E = p.E_last_; + auto E = p.E_last(); // Bin the energy. if (E >= bins_.front() && E <= bins_.back()) { @@ -103,13 +103,13 @@ void EnergyoutFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - if (p.g_ != F90_NONE && matches_transport_groups_) { - match.bins_.push_back(data::mg.num_energy_groups_ - p.g_ - 1); + if (p.g() != F90_NONE && matches_transport_groups_) { + match.bins_.push_back(data::mg.num_energy_groups_ - p.g() - 1); match.weights_.push_back(1.0); } else { - if (p.E_ >= bins_.front() && p.E_ <= bins_.back()) { - auto bin = lower_bound_index(bins_.begin(), bins_.end(), p.E_); + if (p.E() >= bins_.front() && p.E() <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), p.E()); match.bins_.push_back(bin); match.weights_.push_back(1.0); } diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index 498fd52e8..bea35c8e0 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -56,12 +56,12 @@ void EnergyFunctionFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - if (p.E_last_ >= energy_.front() && p.E_last_ <= energy_.back()) { + if (p.E_last() >= energy_.front() && p.E_last() <= energy_.back()) { // Search for the incoming energy bin. - auto i = lower_bound_index(energy_.begin(), energy_.end(), p.E_last_); + auto i = lower_bound_index(energy_.begin(), energy_.end(), p.E_last()); // Compute the interpolation factor between the nearest bins. - double f = (p.E_last_ - energy_[i]) / (energy_[i+1] - energy_[i]); + double f = (p.E_last() - energy_[i]) / (energy_[i + 1] - energy_[i]); // Interpolate on the lin-lin grid. match.bins_.push_back(0); diff --git a/src/tallies/filter_legendre.cpp b/src/tallies/filter_legendre.cpp index 57912983c..314f98a29 100644 --- a/src/tallies/filter_legendre.cpp +++ b/src/tallies/filter_legendre.cpp @@ -28,7 +28,7 @@ LegendreFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { std::vector wgt(n_bins_); - calc_pn_c(order_, p.mu_, wgt.data()); + calc_pn_c(order_, p.mu(), wgt.data()); for (int i = 0; i < n_bins_; i++) { match.bins_.push_back(i); match.weights_.push_back(wgt[i]); diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp index df50ff697..0dec223c5 100644 --- a/src/tallies/filter_material.cpp +++ b/src/tallies/filter_material.cpp @@ -48,7 +48,7 @@ void MaterialFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - auto search = map_.find(p.material_); + auto search = map_.find(p.material()); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_mu.cpp b/src/tallies/filter_mu.cpp index d1f5569cd..b157ff2db 100644 --- a/src/tallies/filter_mu.cpp +++ b/src/tallies/filter_mu.cpp @@ -52,8 +52,8 @@ void MuFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - if (p.mu_ >= bins_.front() && p.mu_ <= bins_.back()) { - auto bin = lower_bound_index(bins_.begin(), bins_.end(), p.mu_); + if (p.mu() >= bins_.front() && p.mu() <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), p.mu()); match.bins_.push_back(bin); match.weights_.push_back(1.0); } diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index c2b045a67..eba56421b 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -38,7 +38,7 @@ ParticleFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { for (auto i = 0; i < particles_.size(); i++) { - if (particles_[i] == p.type_) { + if (particles_[i] == p.type()) { match.bins_.push_back(i); match.weights_.push_back(1.0); } diff --git a/src/tallies/filter_polar.cpp b/src/tallies/filter_polar.cpp index 6d4d899cf..51ff496cd 100644 --- a/src/tallies/filter_polar.cpp +++ b/src/tallies/filter_polar.cpp @@ -53,7 +53,8 @@ void PolarFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - double z = (estimator == TallyEstimator::TRACKLENGTH) ? p.u().z : p.u_last_.z; + double z = + (estimator == TallyEstimator::TRACKLENGTH) ? p.u().z : p.u_last().z; double theta = std::acos(z); if (theta >= bins_.front() && theta <= bins_.back()) { diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp index 44ec5f160..0dbc1e59b 100644 --- a/src/tallies/filter_sph_harm.cpp +++ b/src/tallies/filter_sph_harm.cpp @@ -51,7 +51,7 @@ SphericalHarmonicsFilter::get_all_bins(const Particle& p, TallyEstimator estimat // Determine cosine term for scatter expansion if necessary std::vector wgt(order_ + 1); if (cosine_ == SphericalHarmonicsCosine::scatter) { - calc_pn_c(order_, p.mu_, wgt.data()); + calc_pn_c(order_, p.mu(), wgt.data()); } else { for (int i = 0; i < order_ + 1; i++) { wgt[i] = 1; @@ -60,7 +60,7 @@ SphericalHarmonicsFilter::get_all_bins(const Particle& p, TallyEstimator estimat // Find the Rn,m values std::vector rn(n_bins_); - calc_rn(order_, p.u_last_, rn.data()); + calc_rn(order_, p.u_last(), rn.data()); int j = 0; for (int n = 0; n < order_ + 1; n++) { diff --git a/src/tallies/filter_surface.cpp b/src/tallies/filter_surface.cpp index 141e1e847..d65530d6f 100644 --- a/src/tallies/filter_surface.cpp +++ b/src/tallies/filter_surface.cpp @@ -50,10 +50,10 @@ void SurfaceFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - auto search = map_.find(std::abs(p.surface_)-1); + auto search = map_.find(std::abs(p.surface()) - 1); if (search != map_.end()) { match.bins_.push_back(search->second); - if (p.surface_ < 0) { + if (p.surface() < 0) { match.weights_.push_back(-1.0); } else { match.weights_.push_back(1.0); diff --git a/src/tallies/filter_universe.cpp b/src/tallies/filter_universe.cpp index d1b0d057e..bcd20d795 100644 --- a/src/tallies/filter_universe.cpp +++ b/src/tallies/filter_universe.cpp @@ -48,8 +48,8 @@ void UniverseFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - for (int i = 0; i < p.n_coord_; i++) { - auto search = map_.find(p.coord_[i].universe); + for (int i = 0; i < p.n_coord(); i++) { + auto search = map_.find(p.coord(i).universe); if (search != map_.end()) { match.bins_.push_back(search->second); match.weights_.push_back(1.0); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 91f10f559..c643d014d 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -27,7 +27,7 @@ namespace openmc { //============================================================================== FilterBinIter::FilterBinIter(const Tally& tally, Particle& p) - : filter_matches_{p.filter_matches_}, tally_{tally} + : filter_matches_ {p.filter_matches()}, tally_ {tally} { // Find all valid bins in each relevant filter if they have not already been // found for this event. @@ -181,11 +181,11 @@ double get_nuc_fission_q(const Nuclide& nuc, const Particle& p, int score_bin) { if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) { - return (*nuc.fission_q_prompt_)(p.E_last_); + return (*nuc.fission_q_prompt_)(p.E_last()); } } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) { - return (*nuc.fission_q_recov_)(p.E_last_); + return (*nuc.fission_q_recov_)(p.E_last()); } } return 0.0; @@ -200,42 +200,43 @@ double score_fission_q(const Particle& p, int score_bin, const Tally& tally, double flux, int i_nuclide, double atom_density) { if (tally.estimator_ == TallyEstimator::ANALOG) { - const Nuclide& nuc {*data::nuclides[p.event_nuclide_]}; + const Nuclide& nuc {*data::nuclides[p.event_nuclide()]}; if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { - return p.wgt_absorb_ * get_nuc_fission_q(nuc, p, score_bin) - * p.neutron_xs_[p.event_nuclide_].fission * flux - / p.neutron_xs_[p.event_nuclide_].absorption; + if (p.neutron_xs(p.event_nuclide()).absorption > 0) { + return p.wgt_absorb() * get_nuc_fission_q(nuc, p, score_bin) * + p.neutron_xs(p.event_nuclide()).fission * flux / + p.neutron_xs(p.event_nuclide()).absorption; } } else { // Skip any non-absorption events - if (p.event_ == TallyEvent::SCATTER) return 0.0; + if (p.event() == TallyEvent::SCATTER) + return 0.0; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { - return p.wgt_last_ * get_nuc_fission_q(nuc, p, score_bin) - * p.neutron_xs_[p.event_nuclide_].fission * flux - / p.neutron_xs_[p.event_nuclide_].absorption; + if (p.neutron_xs(p.event_nuclide()).absorption > 0) { + return p.wgt_last() * get_nuc_fission_q(nuc, p, score_bin) * + p.neutron_xs(p.event_nuclide()).fission * flux / + p.neutron_xs(p.event_nuclide()).absorption; } } } else { if (i_nuclide >= 0) { const Nuclide& nuc {*data::nuclides[i_nuclide]}; - return get_nuc_fission_q(nuc, p, score_bin) * atom_density * flux - * p.neutron_xs_[i_nuclide].fission; - } else if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + return get_nuc_fission_q(nuc, p, score_bin) * atom_density * flux * + p.neutron_xs(i_nuclide).fission; + } else if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; double score {0.0}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); const Nuclide& nuc {*data::nuclides[j_nuclide]}; - score += get_nuc_fission_q(nuc, p, score_bin) * atom_density - * p.neutron_xs_[j_nuclide].fission; + score += get_nuc_fission_q(nuc, p, score_bin) * atom_density * + p.neutron_xs(j_nuclide).fission; } return score * flux; } @@ -251,24 +252,26 @@ double get_nuclide_neutron_heating(const Particle& p, const Nuclide& nuc, size_t mt = nuc.reaction_index_[rxn_index]; if (mt == C_NONE) return 0.0; - auto i_temp = p.neutron_xs_[i_nuclide].index_temp; + auto i_temp = p.neutron_xs(i_nuclide).index_temp; if (i_temp < 0) return 0.0; // Can be true due to multipole const auto& rxn {*nuc.reactions_[mt]}; const auto& xs {rxn.xs_[i_temp]}; - auto i_grid = p.neutron_xs_[i_nuclide].index_grid; + auto i_grid = p.neutron_xs(i_nuclide).index_grid; if (i_grid < xs.threshold) return 0.0; // Determine total kerma - auto f = p.neutron_xs_[i_nuclide].interp_factor; + auto f = p.neutron_xs(i_nuclide).interp_factor; double kerma = (1.0 - f) * xs.value[i_grid-xs.threshold] + f * xs.value[i_grid-xs.threshold+1]; if (settings::run_mode == RunMode::EIGENVALUE) { // Determine kerma for fission as (EFR + EB)*sigma_f - double kerma_fission = nuc.fragments_ ? - ((*nuc.fragments_)(p.E_last_) + (*nuc.betas_)(p.E_last_)) - *p.neutron_xs_[i_nuclide].fission : 0.0; + double kerma_fission = + nuc.fragments_ + ? ((*nuc.fragments_)(p.E_last()) + (*nuc.betas_)(p.E_last())) * + p.neutron_xs(i_nuclide).fission + : 0.0; // Determine non-fission kerma as difference double kerma_non_fission = kerma - kerma_fission; @@ -294,14 +297,14 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux, const Nuclide& nuc {*data::nuclides[i_nuclide]}; heating_xs = get_nuclide_neutron_heating(p, nuc, rxn_bin, i_nuclide); if (tally.estimator_ == TallyEstimator::ANALOG) { - heating_xs /= p.neutron_xs_[i_nuclide].total; + heating_xs /= p.neutron_xs(i_nuclide).total; } else { heating_xs *= atom_density; } } else { - if (p.material_ != MATERIAL_VOID) { + if (p.material() != MATERIAL_VOID) { heating_xs = 0.0; - const Material& material {*model::materials[p.material_]}; + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i< material.nuclide_.size(); ++i) { int j_nuclide = material.nuclide_[i]; double atom_density {material.atom_density_(i)}; @@ -309,7 +312,7 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux, heating_xs += atom_density * get_nuclide_neutron_heating(p, nuc, rxn_bin, j_nuclide); } if (tally.estimator_ == TallyEstimator::ANALOG) { - heating_xs /= p.macro_xs_.total; + heating_xs /= p.macro_xs().total; } } } @@ -320,9 +323,9 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux, // reaction-wise heating cross section if (settings::survival_biasing) { // Account for the fact that some weight has been absorbed - score *= p.wgt_last_ + p.wgt_absorb_; + score *= p.wgt_last() + p.wgt_absorb(); } else { - score *= p.wgt_last_; + score *= p.wgt_last(); } } return score; @@ -338,8 +341,8 @@ score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) { auto& tally {*model::tallies[i_tally]}; auto i_eout_filt = tally.filters()[tally.energyout_filter_]; - auto i_bin = p.filter_matches_[i_eout_filt].i_bin_; - auto bin_energyout = p.filter_matches_[i_eout_filt].bins_[i_bin]; + auto i_bin = p.filter_matches(i_eout_filt).i_bin_; + auto bin_energyout = p.filter_matches(i_eout_filt).bins_[i_bin]; const EnergyoutFilter& eo_filt {*dynamic_cast(model::tally_filters[i_eout_filt].get())}; @@ -350,8 +353,8 @@ score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) // rate. Otherwise, the sum of all nu-fission rates would be ~1.0. // loop over number of particles banked - for (auto i = 0; i < p.n_bank_; ++i) { - const auto& bank = p.nu_bank_[i]; + for (auto i = 0; i < p.n_bank(); ++i) { + const auto& bank = p.nu_bank(i); // get the delayed group auto g = bank.delayed_group; @@ -375,7 +378,7 @@ score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) g_out = eo_filt.n_bins() - g_out - 1; // change outgoing energy bin - p.filter_matches_[i_eout_filt].bins_[i_bin] = g_out; + p.filter_matches(i_eout_filt).bins_[i_bin] = g_out; } else { @@ -392,7 +395,7 @@ score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) } else { auto i_match = lower_bound_index(eo_filt.bins().begin(), eo_filt.bins().end(), E_out); - p.filter_matches_[i_eout_filt].bins_[i_bin] = i_match; + p.filter_matches(i_eout_filt).bins_[i_bin] = i_match; } } @@ -406,7 +409,7 @@ score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) double filter_weight = 1.0; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {p.filter_matches_[i_filt]}; + auto& match {p.filter_matches(i_filt)}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(j); filter_weight *= match.weights_[i_bin]; @@ -435,13 +438,13 @@ score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {p.filter_matches_[i_filt]}; + auto& match {p.filter_matches(i_filt)}; auto i_bin = match.i_bin_; filter_weight *= match.weights_[i_bin]; } - score_fission_delayed_dg(i_tally, d_bin, score*filter_weight, - i_score, p.filter_matches_); + score_fission_delayed_dg(i_tally, d_bin, score * filter_weight, + i_score, p.filter_matches()); } } @@ -453,7 +456,7 @@ score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {p.filter_matches_[i_filt]}; + auto& match {p.filter_matches(i_filt)}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(j); filter_weight *= match.weights_[i_bin]; @@ -467,7 +470,7 @@ score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) } // Reset outgoing energy bin and score index - p.filter_matches_[i_eout_filt].bins_[i_bin] = bin_energyout; + p.filter_matches(i_eout_filt).bins_[i_bin] = bin_energyout; } double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin) { @@ -477,7 +480,7 @@ double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin) { auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) return 0.0; const auto& rxn {*nuc.reactions_[m]}; - const auto& micro {p.neutron_xs_[i_nuclide]}; + const auto& micro {p.neutron_xs(i_nuclide)}; // In the URR, the (n,gamma) cross section is sampled randomly from // probability tables. Make sure we use the sampled value (which is equal to @@ -504,10 +507,13 @@ double get_nuclide_xs(const Particle& p, int i_nuclide, int score_bin) { if (settings::run_mode == RunMode::EIGENVALUE && score_bin == HEATING_LOCAL) { // Determine kerma for fission as (EFR + EGP + EGD + EB)*sigma_f - double kerma_fission = nuc.fragments_ ? - ((*nuc.fragments_)(p.E_last_) + (*nuc.betas_)(p.E_last_) + - (*nuc.prompt_photons_)(p.E_last_) + (*nuc.delayed_photons_)(p.E_last_)) - * micro.fission : 0.0; + double kerma_fission = + nuc.fragments_ + ? ((*nuc.fragments_)(p.E_last()) + (*nuc.betas_)(p.E_last()) + + (*nuc.prompt_photons_)(p.E_last()) + + (*nuc.delayed_photons_)(p.E_last())) * + micro.fission + : 0.0; // Determine non-fission kerma as difference double kerma_non_fission = value - kerma_fission; @@ -539,7 +545,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, Tally& tally {*model::tallies[i_tally]}; // Get the pre-collision energy of the particle. - auto E = p.E_last_; + auto E = p.E_last(); using Type = Particle::Type; @@ -557,13 +563,13 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p.wgt_last_ + p.wgt_absorb_; + score = p.wgt_last() + p.wgt_absorb(); } else { - score = p.wgt_last_; + score = p.wgt_last(); } - if (p.type_ == Type::neutron || p.type_ == Type::photon) { - score *= flux / p.macro_xs_.total; + if (p.type() == Type::neutron || p.type() == Type::photon) { + score *= flux / p.macro_xs().total; } else { score = 0.; } @@ -580,27 +586,28 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = (p.wgt_last_ + p.wgt_absorb_) * flux; + score = (p.wgt_last() + p.wgt_absorb()) * flux; } else { - score = p.wgt_last_ * flux; + score = p.wgt_last() * flux; } } else { if (i_nuclide >= 0) { - if (p.type_ == Type::neutron) { - score = p.neutron_xs_[i_nuclide].total * atom_density * flux; - } else if (p.type_ == Type::photon) { - score = p.photon_xs_[i_nuclide].total * atom_density * flux; + if (p.type() == Type::neutron) { + score = p.neutron_xs(i_nuclide).total * atom_density * flux; + } else if (p.type() == Type::photon) { + score = p.photon_xs(i_nuclide).total * atom_density * flux; } } else { - score = p.macro_xs_.total * flux; + score = p.macro_xs().total * flux; } } break; case SCORE_INVERSE_VELOCITY: - if (p.type_ != Type::neutron) continue; + if (p.type() != Type::neutron) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // All events score to an inverse velocity bin. We actually use a @@ -609,11 +616,11 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p.wgt_last_ + p.wgt_absorb_; + score = p.wgt_last() + p.wgt_absorb(); } else { - score = p.wgt_last_; + score = p.wgt_last(); } - score *= flux / p.macro_xs_.total; + score *= flux / p.macro_xs().total; } else { score = flux; } @@ -623,28 +630,30 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_SCATTER: - if (p.type_ != Type::neutron && p.type_ != Type::photon) continue; + if (p.type() != Type::neutron && p.type() != Type::photon) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // Skip any event where the particle didn't scatter - if (p.event_ != TallyEvent::SCATTER) continue; + if (p.event() != TallyEvent::SCATTER) + continue; // Since only scattering events make it here, again we can use the // weight entering the collision as the estimator for the reaction rate - score = p.wgt_last_ * flux; + score = p.wgt_last() * flux; } else { if (i_nuclide >= 0) { - if (p.type_ == Type::neutron) { - const auto& micro = p.neutron_xs_[i_nuclide]; + if (p.type() == Type::neutron) { + const auto& micro = p.neutron_xs(i_nuclide); score = (micro.total - micro.absorption) * atom_density * flux; } else { - const auto& micro = p.photon_xs_[i_nuclide]; + const auto& micro = p.photon_xs(i_nuclide); score = (micro.coherent + micro.incoherent) * atom_density * flux; } } else { - if (p.type_ == Type::neutron) { - score = (p.macro_xs_.total - p.macro_xs_.absorption) * flux; + if (p.type() == Type::neutron) { + score = (p.macro_xs().total - p.macro_xs().absorption) * flux; } else { - score = (p.macro_xs_.coherent + p.macro_xs_.incoherent) * flux; + score = (p.macro_xs().coherent + p.macro_xs().incoherent) * flux; } } } @@ -652,56 +661,63 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_NU_SCATTER: - if (p.type_ != Type::neutron) continue; + if (p.type() != Type::neutron) + continue; // Only analog estimators are available. // Skip any event where the particle didn't scatter - if (p.event_ != TallyEvent::SCATTER) continue; + if (p.event() != TallyEvent::SCATTER) + continue; // For scattering production, we need to use the pre-collision weight // times the yield as the estimate for the number of neutrons exiting a // reaction with neutrons in the exit channel - if (p.event_mt_ == ELASTIC || p.event_mt_ == N_LEVEL || - (p.event_mt_ >= N_N1 && p.event_mt_ <= N_NC)) { + if (p.event_mt() == ELASTIC || p.event_mt() == N_LEVEL || + (p.event_mt() >= N_N1 && p.event_mt() <= N_NC)) { // Don't waste time on very common reactions we know have // multiplicities of one. - score = p.wgt_last_ * flux; + score = p.wgt_last() * flux; } else { // Get yield and apply to score - auto m = data::nuclides[p.event_nuclide_]->reaction_index_[p.event_mt_]; - const auto& rxn {*data::nuclides[p.event_nuclide_]->reactions_[m]}; - score = p.wgt_last_ * flux * (*rxn.products_[0].yield_)(E); + auto m = + data::nuclides[p.event_nuclide()]->reaction_index_[p.event_mt()]; + const auto& rxn {*data::nuclides[p.event_nuclide()]->reactions_[m]}; + score = p.wgt_last() * flux * (*rxn.products_[0].yield_)(E); } break; case SCORE_ABSORPTION: - if (p.type_ != Type::neutron && p.type_ != Type::photon) continue; + if (p.type() != Type::neutron && p.type() != Type::photon) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p.wgt_absorb_ * flux; + score = p.wgt_absorb() * flux; } else { // Skip any event where the particle wasn't absorbed - if (p.event_ == TallyEvent::SCATTER) continue; + if (p.event() == TallyEvent::SCATTER) + continue; // All fission and absorption events will contribute here, so we // can just use the particle's weight entering the collision - score = p.wgt_last_ * flux; + score = p.wgt_last() * flux; } } else { if (i_nuclide >= 0) { - if (p.type_ == Type::neutron) { - score = p.neutron_xs_[i_nuclide].absorption * atom_density * flux; + if (p.type() == Type::neutron) { + score = p.neutron_xs(i_nuclide).absorption * atom_density * flux; } else { - const auto& xs = p.photon_xs_[i_nuclide]; + const auto& xs = p.photon_xs(i_nuclide); score = (xs.total - xs.coherent - xs.incoherent) * atom_density * flux; } } else { - if (p.type_ == Type::neutron) { - score = p.macro_xs_.absorption * flux; + if (p.type() == Type::neutron) { + score = p.macro_xs().absorption * flux; } else { - score = (p.macro_xs_.photoelectric + p.macro_xs_.pair_production) * flux; + score = + (p.macro_xs().photoelectric + p.macro_xs().pair_production) * + flux; } } } @@ -709,43 +725,44 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_FISSION: - if (p.macro_xs_.absorption == 0) continue; + if (p.macro_xs().absorption == 0) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { - score = p.wgt_absorb_ - * p.neutron_xs_[p.event_nuclide_].fission - / p.neutron_xs_[p.event_nuclide_].absorption * flux; + if (p.neutron_xs(p.event_nuclide()).absorption > 0) { + score = p.wgt_absorb() * p.neutron_xs(p.event_nuclide()).fission / + p.neutron_xs(p.event_nuclide()).absorption * flux; } else { score = 0.; } } else { // Skip any non-absorption events - if (p.event_ == TallyEvent::SCATTER) continue; + if (p.event() == TallyEvent::SCATTER) + continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - score = p.wgt_last_ - * p.neutron_xs_[p.event_nuclide_].fission - / p.neutron_xs_[p.event_nuclide_].absorption * flux; + score = p.wgt_last() * p.neutron_xs(p.event_nuclide()).fission / + p.neutron_xs(p.event_nuclide()).absorption * flux; } } else { if (i_nuclide >= 0) { - score = p.neutron_xs_[i_nuclide].fission * atom_density * flux; + score = p.neutron_xs(i_nuclide).fission * atom_density * flux; } else { - score = p.macro_xs_.fission * flux; + score = p.macro_xs().fission * flux; } } break; case SCORE_NU_FISSION: - if (p.macro_xs_.absorption == 0) continue; + if (p.macro_xs().absorption == 0) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p.fission_) { + if (settings::survival_biasing || p.fission()) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -757,38 +774,39 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { - score = p.wgt_absorb_ - * p.neutron_xs_[p.event_nuclide_].nu_fission - / p.neutron_xs_[p.event_nuclide_].absorption * flux; + if (p.neutron_xs(p.event_nuclide()).absorption > 0) { + score = p.wgt_absorb() * + p.neutron_xs(p.event_nuclide()).nu_fission / + p.neutron_xs(p.event_nuclide()).absorption * flux; } else { score = 0.; } } else { // Skip any non-fission events - if (!p.fission_) continue; + if (!p.fission()) + continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - score = simulation::keff * p.wgt_bank_ * flux; + score = simulation::keff * p.wgt_bank() * flux; } } else { if (i_nuclide >= 0) { - score = p.neutron_xs_[i_nuclide].nu_fission * atom_density - * flux; + score = p.neutron_xs(i_nuclide).nu_fission * atom_density * flux; } else { - score = p.macro_xs_.nu_fission * flux; + score = p.macro_xs().nu_fission * flux; } } break; case SCORE_PROMPT_NU_FISSION: - if (p.macro_xs_.absorption == 0) continue; + if (p.macro_xs().absorption == 0) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p.fission_) { + if (settings::survival_biasing || p.fission()) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -800,46 +818,46 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { - score = p.wgt_absorb_ - * p.neutron_xs_[p.event_nuclide_].fission - * data::nuclides[p.event_nuclide_] - ->nu(E, ReactionProduct::EmissionMode::prompt) - / p.neutron_xs_[p.event_nuclide_].absorption * flux; + if (p.neutron_xs(p.event_nuclide()).absorption > 0) { + score = p.wgt_absorb() * p.neutron_xs(p.event_nuclide()).fission * + data::nuclides[p.event_nuclide()]->nu( + E, ReactionProduct::EmissionMode::prompt) / + p.neutron_xs(p.event_nuclide()).absorption * flux; } else { score = 0.; } } else { // Skip any non-fission events - if (!p.fission_) continue; + if (!p.fission()) + continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - auto n_delayed = std::accumulate(p.n_delayed_bank_, - p.n_delayed_bank_+MAX_DELAYED_GROUPS, 0); - auto prompt_frac = 1. - n_delayed / static_cast(p.n_bank_); - score = simulation::keff * p.wgt_bank_ * prompt_frac * flux; + auto n_delayed = std::accumulate( + p.n_delayed_bank(), p.n_delayed_bank() + MAX_DELAYED_GROUPS, 0); + auto prompt_frac = 1. - n_delayed / static_cast(p.n_bank()); + score = simulation::keff * p.wgt_bank() * prompt_frac * flux; } } else { if (i_nuclide >= 0) { - score = p.neutron_xs_[i_nuclide].fission - * data::nuclides[i_nuclide] - ->nu(E, ReactionProduct::EmissionMode::prompt) - * atom_density * flux; + score = p.neutron_xs(i_nuclide).fission * + data::nuclides[i_nuclide]->nu( + E, ReactionProduct::EmissionMode::prompt) * + atom_density * flux; } else { score = 0.; // Add up contributions from each nuclide in the material. - if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += p.neutron_xs_[j_nuclide].fission - * data::nuclides[j_nuclide] - ->nu(E, ReactionProduct::EmissionMode::prompt) - * atom_density * flux; + score += p.neutron_xs(j_nuclide).fission * + data::nuclides[j_nuclide]->nu( + E, ReactionProduct::EmissionMode::prompt) * + atom_density * flux; } } } @@ -848,9 +866,10 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_DELAYED_NU_FISSION: - if (p.macro_xs_.absorption == 0) continue; + if (p.macro_xs().absorption == 0) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p.fission_) { + if (settings::survival_biasing || p.fission()) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -862,8 +881,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - if (p.neutron_xs_[p.event_nuclide_].absorption > 0 - && data::nuclides[p.event_nuclide_]->fissionable_) { + if (p.neutron_xs(p.event_nuclide()).absorption > 0 && + data::nuclides[p.event_nuclide()]->fissionable_) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; const DelayedGroupFilter& filt @@ -872,29 +891,29 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto dg = filt.groups()[d_bin]; - auto yield = data::nuclides[p.event_nuclide_] - ->nu(E, ReactionProduct::EmissionMode::delayed, dg); - score = p.wgt_absorb_ * yield - * p.neutron_xs_[p.event_nuclide_].fission - / p.neutron_xs_[p.event_nuclide_].absorption * flux; - score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p.filter_matches_); + auto yield = data::nuclides[p.event_nuclide()]->nu( + E, ReactionProduct::EmissionMode::delayed, dg); + score = p.wgt_absorb() * yield * + p.neutron_xs(p.event_nuclide()).fission / + p.neutron_xs(p.event_nuclide()).absorption * flux; + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { // If the delayed group filter is not present, compute the score // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs - score = p.wgt_absorb_ - * p.neutron_xs_[p.event_nuclide_].fission - * data::nuclides[p.event_nuclide_] - ->nu(E, ReactionProduct::EmissionMode::delayed) - / p.neutron_xs_[p.event_nuclide_].absorption *flux; + score = p.wgt_absorb() * p.neutron_xs(p.event_nuclide()).fission * + data::nuclides[p.event_nuclide()]->nu( + E, ReactionProduct::EmissionMode::delayed) / + p.neutron_xs(p.event_nuclide()).absorption * flux; } } } else { // Skip any non-fission events - if (!p.fission_) continue; + if (!p.fission()) + continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -910,17 +929,18 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin]; - score = simulation::keff * p.wgt_bank_ / p.n_bank_ - * p.n_delayed_bank_[d-1] * flux; - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); + score = simulation::keff * p.wgt_bank() / p.n_bank() * + p.n_delayed_bank(d - 1) * flux; + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { // Add the contribution from all delayed groups - auto n_delayed = std::accumulate(p.n_delayed_bank_, - p.n_delayed_bank_+MAX_DELAYED_GROUPS, 0); - score = simulation::keff * p.wgt_bank_ / p.n_bank_ * n_delayed - * flux; + auto n_delayed = std::accumulate( + p.n_delayed_bank(), p.n_delayed_bank() + MAX_DELAYED_GROUPS, 0); + score = + simulation::keff * p.wgt_bank() / p.n_bank() * n_delayed * flux; } } } else { @@ -935,18 +955,19 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, auto d = filt.groups()[d_bin]; auto yield = data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = p.neutron_xs_[i_nuclide].fission * yield - * atom_density * flux; - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); + score = + p.neutron_xs(i_nuclide).fission * yield * atom_density * flux; + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { // If the delayed group filter is not present, compute the score // by multiplying the delayed-nu-fission macro xs by the flux - score = p.neutron_xs_[i_nuclide].fission - * data::nuclides[i_nuclide] - ->nu(E, ReactionProduct::EmissionMode::delayed) - * atom_density * flux; + score = p.neutron_xs(i_nuclide).fission * + data::nuclides[i_nuclide]->nu( + E, ReactionProduct::EmissionMode::delayed) * + atom_density * flux; } } else { // Need to add up contributions for each nuclide @@ -955,8 +976,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, const DelayedGroupFilter& filt {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; - if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -965,25 +986,25 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, auto d = filt.groups()[d_bin]; auto yield = data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = p.neutron_xs_[j_nuclide].fission * yield - * atom_density * flux; - score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p.filter_matches_); + score = p.neutron_xs(j_nuclide).fission * yield * + atom_density * flux; + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } } } continue; } else { score = 0.; - if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += p.neutron_xs_[j_nuclide].fission - * data::nuclides[j_nuclide] - ->nu(E, ReactionProduct::EmissionMode::delayed) - * atom_density * flux; + score += p.neutron_xs(j_nuclide).fission * + data::nuclides[j_nuclide]->nu( + E, ReactionProduct::EmissionMode::delayed) * + atom_density * flux; } } } @@ -993,15 +1014,16 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_DECAY_RATE: - if (p.macro_xs_.absorption == 0) continue; + if (p.macro_xs().absorption == 0) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - const auto& nuc {*data::nuclides[p.event_nuclide_]}; - if (p.neutron_xs_[p.event_nuclide_].absorption > 0 - && nuc.fissionable_) { + const auto& nuc {*data::nuclides[p.event_nuclide()]}; + if (p.neutron_xs(p.event_nuclide()).absorption > 0 && + nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; @@ -1014,12 +1036,12 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p.wgt_absorb_ * yield - * p.neutron_xs_[p.event_nuclide_].fission - / p.neutron_xs_[p.event_nuclide_].absorption - * rate * flux; - score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p.filter_matches_); + score = p.wgt_absorb() * yield * + p.neutron_xs(p.event_nuclide()).fission / + p.neutron_xs(p.event_nuclide()).absorption * rate * + flux; + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { @@ -1037,15 +1059,16 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += rate * p.wgt_absorb_ - * p.neutron_xs_[p.event_nuclide_].fission * yield - / p.neutron_xs_[p.event_nuclide_].absorption * flux; + score += rate * p.wgt_absorb() * + p.neutron_xs(p.event_nuclide()).fission * yield / + p.neutron_xs(p.event_nuclide()).absorption * flux; } } } } else { // Skip any non-fission events - if (!p.fission_) continue; + if (!p.fission()) + continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -1054,11 +1077,11 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // ones are delayed. If a delayed neutron is encountered, add its // contribution to the fission bank to the score. score = 0.; - for (auto i = 0; i < p.n_bank_; ++i) { - const auto& bank = p.nu_bank_[i]; + for (auto i = 0; i < p.n_bank(); ++i) { + const auto& bank = p.nu_bank(i); auto g = bank.delayed_group; if (g != 0) { - const auto& nuc {*data::nuclides[p.event_nuclide_]}; + const auto& nuc {*data::nuclides[p.event_nuclide()]}; const auto& rxn {*nuc.fission_rx_[0]}; auto rate = rxn.products_[g].decay_rate_; score += simulation::keff * bank.wgt * rate * flux; @@ -1071,8 +1094,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin]; if (d == g) - score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p.filter_matches_); + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } score = 0.; } @@ -1095,9 +1118,10 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p.neutron_xs_[i_nuclide].fission * yield * flux - * atom_density * rate; - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); + score = p.neutron_xs(i_nuclide).fission * yield * flux * + atom_density * rate; + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { @@ -1111,8 +1135,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += p.neutron_xs_[i_nuclide].fission * flux - * yield * atom_density * rate; + score += p.neutron_xs(i_nuclide).fission * flux * yield * + atom_density * rate; } } } else { @@ -1121,8 +1145,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, const DelayedGroupFilter& filt {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; - if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1135,10 +1159,10 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p.neutron_xs_[j_nuclide].fission * yield - * flux * atom_density * rate; - score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p.filter_matches_); + score = p.neutron_xs(j_nuclide).fission * yield * flux * + atom_density * rate; + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } } } @@ -1146,8 +1170,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, continue; } else { score = 0.; - if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1163,8 +1187,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += p.neutron_xs_[j_nuclide].fission - * yield * atom_density * flux * rate; + score += p.neutron_xs(j_nuclide).fission * yield * + atom_density * flux * rate; } } } @@ -1176,7 +1200,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_KAPPA_FISSION: - if (p.macro_xs_.absorption == 0.) continue; + if (p.macro_xs().absorption == 0.) + continue; score = 0.; // Kappa-fission values are determined from the Q-value listed for the // fission cross section. @@ -1185,27 +1210,28 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - const auto& nuc {*data::nuclides[p.event_nuclide_]}; - if (p.neutron_xs_[p.event_nuclide_].absorption > 0 - && nuc.fissionable_) { + const auto& nuc {*data::nuclides[p.event_nuclide()]}; + if (p.neutron_xs(p.event_nuclide()).absorption > 0 && + nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = p.wgt_absorb_ * rxn.q_value_ - * p.neutron_xs_[p.event_nuclide_].fission - / p.neutron_xs_[p.event_nuclide_].absorption * flux; + score = p.wgt_absorb() * rxn.q_value_ * + p.neutron_xs(p.event_nuclide()).fission / + p.neutron_xs(p.event_nuclide()).absorption * flux; } } else { // Skip any non-absorption events - if (p.event_ == TallyEvent::SCATTER) continue; + if (p.event() == TallyEvent::SCATTER) + continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - const auto& nuc {*data::nuclides[p.event_nuclide_]}; - if (p.neutron_xs_[p.event_nuclide_].absorption > 0 - && nuc.fissionable_) { + const auto& nuc {*data::nuclides[p.event_nuclide()]}; + if (p.neutron_xs(p.event_nuclide()).absorption > 0 && + nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = p.wgt_last_ * rxn.q_value_ - * p.neutron_xs_[p.event_nuclide_].fission - / p.neutron_xs_[p.event_nuclide_].absorption * flux; + score = p.wgt_last() * rxn.q_value_ * + p.neutron_xs(p.event_nuclide()).fission / + p.neutron_xs(p.event_nuclide()).absorption * flux; } } } else { @@ -1213,19 +1239,19 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, const auto& nuc {*data::nuclides[i_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = rxn.q_value_ * p.neutron_xs_[i_nuclide].fission - * atom_density * flux; + score = rxn.q_value_ * p.neutron_xs(i_nuclide).fission * + atom_density * flux; } - } else if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + } else if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); const auto& nuc {*data::nuclides[j_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score += rxn.q_value_ * p.neutron_xs_[j_nuclide].fission - * atom_density * flux; + score += rxn.q_value_ * p.neutron_xs(j_nuclide).fission * + atom_density * flux; } } } @@ -1241,28 +1267,29 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case ELASTIC: - if (p.type_ != Type::neutron) continue; + if (p.type() != Type::neutron) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // Check if event MT matches - if (p.event_mt_ != ELASTIC) continue; - score = p.wgt_last_ * flux; + if (p.event_mt() != ELASTIC) + continue; + score = p.wgt_last() * flux; } else { if (i_nuclide >= 0) { - if (p.neutron_xs_[i_nuclide].elastic == CACHE_INVALID) + if (p.neutron_xs(i_nuclide).elastic == CACHE_INVALID) data::nuclides[i_nuclide]->calculate_elastic_xs(p); - score = p.neutron_xs_[i_nuclide].elastic * atom_density * flux; + score = p.neutron_xs(i_nuclide).elastic * atom_density * flux; } else { score = 0.; - if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - if (p.neutron_xs_[j_nuclide].elastic == CACHE_INVALID) + if (p.neutron_xs(j_nuclide).elastic == CACHE_INVALID) data::nuclides[j_nuclide]->calculate_elastic_xs(p); - score += p.neutron_xs_[j_nuclide].elastic * atom_density - * flux; + score += p.neutron_xs(j_nuclide).elastic * atom_density * flux; } } } @@ -1271,7 +1298,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_FISS_Q_PROMPT: case SCORE_FISS_Q_RECOV: - if (p.macro_xs_.absorption == 0.) continue; + if (p.macro_xs().absorption == 0.) + continue; score = score_fission_q(p, score_bin, tally, flux, i_nuclide, atom_density); break; @@ -1287,12 +1315,14 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // which looks up cross sections if (!simulation::need_depletion_rx) goto default_case; - if (p.type_ != Type::neutron) continue; + if (p.type() != Type::neutron) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // Check if the event MT matches - if (p.event_mt_ != score_bin) continue; - score = p.wgt_last_ * flux; + if (p.event_mt() != score_bin) + continue; + score = p.wgt_last() * flux; } else { int m; switch (score_bin) { @@ -1304,17 +1334,16 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case N_4N: m = 5; break; } if (i_nuclide >= 0) { - score = p.neutron_xs_[i_nuclide].reaction[m] * atom_density - * flux; + score = p.neutron_xs(i_nuclide).reaction[m] * atom_density * flux; } else { score = 0.; - if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += p.neutron_xs_[j_nuclide].reaction[m] - * atom_density * flux; + score += + p.neutron_xs(j_nuclide).reaction[m] * atom_density * flux; } } } @@ -1326,22 +1355,25 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case INCOHERENT: case PHOTOELECTRIC: case PAIR_PROD: - if (p.type_ != Type::photon) continue; + if (p.type() != Type::photon) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { if (score_bin == PHOTOELECTRIC) { // Photoelectric events are assigned an MT value corresponding to the // shell cross section. Also, photons below the energy cutoff are // assumed to have been absorbed via photoelectric absorption - if ((p.event_mt_ < 534 || p.event_mt_ > 572) && - p.event_mt_ != REACTION_NONE) continue; + if ((p.event_mt() < 534 || p.event_mt() > 572) && + p.event_mt() != REACTION_NONE) + continue; } else { - if (p.event_mt_ != score_bin) continue; + if (p.event_mt() != score_bin) + continue; } - score = p.wgt_last_ * flux; + score = p.wgt_last() * flux; } else { if (i_nuclide >= 0) { - const auto& micro = p.photon_xs_[i_nuclide]; + const auto& micro = p.photon_xs(i_nuclide); double xs = (score_bin == COHERENT) ? micro.coherent : (score_bin == INCOHERENT) ? micro.incoherent : @@ -1349,11 +1381,11 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, micro.pair_production; score = xs * atom_density * flux; } else { - double xs = - (score_bin == COHERENT) ? p.macro_xs_.coherent : - (score_bin == INCOHERENT) ? p.macro_xs_.incoherent : - (score_bin == PHOTOELECTRIC) ? p.macro_xs_.photoelectric : - p.macro_xs_.pair_production; + double xs = (score_bin == COHERENT) ? p.macro_xs().coherent + : (score_bin == INCOHERENT) ? p.macro_xs().incoherent + : (score_bin == PHOTOELECTRIC) + ? p.macro_xs().photoelectric + : p.macro_xs().pair_production; score = xs * flux; } } @@ -1361,22 +1393,22 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case HEATING: - if (p.type_ == Type::neutron) { + if (p.type() == Type::neutron) { score = score_neutron_heating(p, tally, flux, HEATING, i_nuclide, atom_density); } else { // The energy deposited is the difference between the pre-collision and // post-collision energy... - score = E - p.E_; + score = E - p.E(); // ...less the energy of any secondary particles since they will be // transported individually later - const auto& bank = p.secondary_bank_; - for (auto it = bank.end() - p.n_bank_second_; it < bank.end(); ++it) { + const auto& bank = p.secondary_bank(); + for (auto it = bank.end() - p.n_bank_second(); it < bank.end(); ++it) { score -= it->E; } - score *= p.wgt_last_; + score *= p.wgt_last(); } break; @@ -1385,14 +1417,16 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // The default block is really only meant for redundant neutron reactions // (e.g. 444, 901) - if (p.type_ != Type::neutron) continue; + if (p.type() != Type::neutron) + continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // Any other score is assumed to be a MT number. Thus, we just need // to check if it matches the MT number of the event - if (p.event_mt_ != score_bin) continue; - score = p.wgt_last_*flux; + if (p.event_mt() != score_bin) + continue; + score = p.wgt_last() * flux; } else { // Any other cross section has to be calculated on-the-fly if (score_bin < 2) fatal_error("Invalid score type on tally " @@ -1400,8 +1434,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, score = 0.; if (i_nuclide >= 0) { score = get_nuclide_xs(p, i_nuclide, score_bin) * atom_density * flux; - } else if (p.material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p.material_]}; + } else if (p.material() != MATERIAL_VOID) { + const Material& material {*model::materials[p.material()]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1443,40 +1477,40 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, // Then we either are alive and had a scatter (and so g changed), // or are dead and g did not change - if (p.alive_) { - p_u = p.u_last_; - p_g = p.g_last_; + if (p.alive()) { + p_u = p.u_last(); + p_g = p.g_last(); } else { p_u = p.u_local(); - p_g = p.g_; + p_g = p.g(); } - } else if (p.event_ == TallyEvent::SCATTER) { + } else if (p.event() == TallyEvent::SCATTER) { // Then the energy group has been changed by the scattering routine // meaning gin is now in p % last_g - p_u = p.u_last_; - p_g = p.g_last_; + p_u = p.u_last(); + p_g = p.g_last(); } else { // No scatter, no change in g. p_u = p.u_local(); - p_g = p.g_; + p_g = p.g(); } } else { // No actual collision so g has not changed. p_u = p.u_local(); - p_g = p.g_; + p_g = p.g(); } // For shorthand, assign pointers to the material and nuclide xs set auto& nuc_xs = data::mg.nuclides_[i_nuclide]; - auto& macro_xs = data::mg.macro_xs_[p.material_]; + auto& macro_xs = data::mg.macro_xs_[p.material()]; // Find the temperature and angle indices of interest macro_xs.set_angle_index(p_u); if (i_nuclide >= 0) { - nuc_xs.set_temperature_index(p.sqrtkT_); + nuc_xs.set_temperature_index(p.sqrtkT()); nuc_xs.set_angle_index(p_u); } @@ -1497,11 +1531,11 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p.wgt_last_ + p.wgt_absorb_; + score = p.wgt_last() + p.wgt_absorb(); } else { - score = p.wgt_last_; + score = p.wgt_last(); } - score *= flux / p.macro_xs_.total; + score *= flux / p.macro_xs().total; } else { score = flux; } @@ -1515,9 +1549,9 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p.wgt_last_ + p.wgt_absorb_; + score = p.wgt_last() + p.wgt_absorb(); } else { - score = p.wgt_last_; + score = p.wgt_last(); } //TODO: should flux be multiplied in above instead of below? if (i_nuclide >= 0) { @@ -1528,7 +1562,7 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, if (i_nuclide >= 0) { score = atom_density * flux * nuc_xs.get_xs(MgxsType::TOTAL, p_g); } else { - score = p.macro_xs_.total * flux; + score = p.macro_xs().total * flux; } } break; @@ -1543,9 +1577,9 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p.wgt_last_ + p.wgt_absorb_; + score = p.wgt_last() + p.wgt_absorb(); } else { - score = p.wgt_last_; + score = p.wgt_last(); } if (i_nuclide >= 0) { score *= flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) @@ -1567,23 +1601,26 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_SCATTER: if (tally.estimator_ == TallyEstimator::ANALOG) { // Skip any event where the particle didn't scatter - if (p.event_ != TallyEvent::SCATTER) continue; + if (p.event() != TallyEvent::SCATTER) + continue; // Since only scattering events make it here, again we can use the // weight entering the collision as the estimator for the reaction rate - score = p.wgt_last_ * flux; + score = p.wgt_last() * flux; if (i_nuclide >= 0) { - score *= atom_density * nuc_xs.get_xs( - MgxsType::SCATTER_FMU, p.g_last_, &p.g_, &p.mu_, nullptr) - / macro_xs.get_xs( - MgxsType::SCATTER_FMU, p.g_last_, &p.g_, &p.mu_, nullptr); + score *= atom_density * + nuc_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last(), &p.g(), + &p.mu(), nullptr) / + macro_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last(), &p.g(), + &p.mu(), nullptr); } } else { if (i_nuclide >= 0) { - score = atom_density * flux * nuc_xs.get_xs( - MgxsType::SCATTER, p_g, nullptr, &p.mu_, nullptr); + score = + atom_density * flux * + nuc_xs.get_xs(MgxsType::SCATTER, p_g, nullptr, &p.mu(), nullptr); } else { score = flux * macro_xs.get_xs( - MgxsType::SCATTER, p_g, nullptr, &p.mu_, nullptr); + MgxsType::SCATTER, p_g, nullptr, &p.mu(), nullptr); } } break; @@ -1592,20 +1629,21 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_NU_SCATTER: if (tally.estimator_ == TallyEstimator::ANALOG) { // Skip any event where the particle didn't scatter - if (p.event_ != TallyEvent::SCATTER) continue; + if (p.event() != TallyEvent::SCATTER) + continue; // For scattering production, we need to use the pre-collision weight // times the multiplicity as the estimate for the number of neutrons // exiting a reaction with neutrons in the exit channel - score = p.wgt_ * flux; + score = p.wgt() * flux; // Since we transport based on material data, the angle selected // was not selected from the f(mu) for the nuclide. Therefore // adjust the score by the actual probability for that nuclide. if (i_nuclide >= 0) { - score *= atom_density - * nuc_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last_, &p.g_, - &p.mu_, nullptr) - / macro_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last_, &p.g_, - &p.mu_, nullptr); + score *= atom_density * + nuc_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last(), &p.g(), + &p.mu(), nullptr) / + macro_xs.get_xs(MgxsType::NU_SCATTER_FMU, p.g_last(), &p.g(), + &p.mu(), nullptr); } } else { if (i_nuclide >= 0) { @@ -1622,13 +1660,14 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p.wgt_absorb_ * flux; + score = p.wgt_absorb() * flux; } else { // Skip any event where the particle wasn't absorbed - if (p.event_ == TallyEvent::SCATTER) continue; + if (p.event() == TallyEvent::SCATTER) + continue; // All fission and absorption events will contribute here, so we // can just use the particle's weight entering the collision - score = p.wgt_last_ * flux; + score = p.wgt_last() * flux; } if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g) @@ -1639,7 +1678,7 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, score = atom_density * flux * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g); } else { - score = p.macro_xs_.absorption * flux; + score = p.macro_xs().absorption * flux; } } break; @@ -1651,14 +1690,15 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - score = p.wgt_absorb_ * flux; + score = p.wgt_absorb() * flux; } else { // Skip any non-absorption events - if (p.event_ == TallyEvent::SCATTER) continue; + if (p.event() == TallyEvent::SCATTER) + continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - score = p.wgt_last_ * flux; + score = p.wgt_last() * flux; } if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) @@ -1679,7 +1719,7 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_NU_FISSION: if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p.fission_) { + if (settings::survival_biasing || p.fission()) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -1691,7 +1731,7 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - score = p.wgt_absorb_ * flux; + score = p.wgt_absorb() * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g) / macro_xs.get_xs(MgxsType::ABSORPTION, p_g); @@ -1701,13 +1741,14 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, } } else { // Skip any non-fission events - if (!p.fission_) continue; + if (!p.fission()) + continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - score = simulation::keff * p.wgt_bank_ * flux; + score = simulation::keff * p.wgt_bank() * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / macro_xs.get_xs(MgxsType::FISSION, p_g); @@ -1726,7 +1767,7 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_PROMPT_NU_FISSION: if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p.fission_) { + if (settings::survival_biasing || p.fission()) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -1738,7 +1779,7 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - score = p.wgt_absorb_ * flux; + score = p.wgt_absorb() * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g) @@ -1749,16 +1790,17 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, } } else { // Skip any non-fission events - if (!p.fission_) continue; + if (!p.fission()) + continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - auto n_delayed = std::accumulate(p.n_delayed_bank_, - p.n_delayed_bank_+MAX_DELAYED_GROUPS, 0); - auto prompt_frac = 1. - n_delayed / static_cast(p.n_bank_); - score = simulation::keff * p.wgt_bank_ * prompt_frac * flux; + auto n_delayed = std::accumulate( + p.n_delayed_bank(), p.n_delayed_bank() + MAX_DELAYED_GROUPS, 0); + auto prompt_frac = 1. - n_delayed / static_cast(p.n_bank()); + score = simulation::keff * p.wgt_bank() * prompt_frac * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / macro_xs.get_xs(MgxsType::FISSION, p_g); @@ -1777,7 +1819,7 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_DELAYED_NU_FISSION: if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p.fission_) { + if (settings::survival_biasing || p.fission()) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -1799,7 +1841,7 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin] - 1; - score = p.wgt_absorb_ * flux; + score = p.wgt_absorb() * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) @@ -1809,14 +1851,15 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, nullptr, nullptr, &d) / abs_xs; } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { // If the delayed group filter is not present, compute the score // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs - score = p.wgt_absorb_ * flux; + score = p.wgt_absorb() * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g) / abs_xs; @@ -1828,7 +1871,8 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, } } else { // Skip any non-fission events - if (!p.fission_) continue; + if (!p.fission()) + continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -1844,21 +1888,22 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin]; - score = simulation::keff * p.wgt_bank_ / p.n_bank_ - * p.n_delayed_bank_[d-1] * flux; + score = simulation::keff * p.wgt_bank() / p.n_bank() * + p.n_delayed_bank(d - 1) * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / macro_xs.get_xs(MgxsType::FISSION, p_g); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { // Add the contribution from all delayed groups - auto n_delayed = std::accumulate(p.n_delayed_bank_, - p.n_delayed_bank_+MAX_DELAYED_GROUPS, 0); - score = simulation::keff * p.wgt_bank_ / p.n_bank_ * n_delayed - * flux; + auto n_delayed = std::accumulate( + p.n_delayed_bank(), p.n_delayed_bank() + MAX_DELAYED_GROUPS, 0); + score = + simulation::keff * p.wgt_bank() / p.n_bank() * n_delayed * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / macro_xs.get_xs(MgxsType::FISSION, p_g); @@ -1881,7 +1926,8 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, score = flux * macro_xs.get_xs( MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { @@ -1912,7 +1958,7 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin] - 1; - score = p.wgt_absorb_ * flux; + score = p.wgt_absorb() * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) @@ -1924,7 +1970,8 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) / abs_xs; } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { @@ -1935,24 +1982,27 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, score = 0.; for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) { if (i_nuclide >= 0) { - score += p.wgt_absorb_ * flux - * nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, - nullptr, &d) - * nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d) / abs_xs; + score += p.wgt_absorb() * flux * + nuc_xs.get_xs( + MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * + nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, + nullptr, nullptr, &d) / + abs_xs; } else { - score += p.wgt_absorb_ * flux - * macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, - nullptr, &d) - * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, - nullptr, nullptr, &d) / abs_xs; + score += p.wgt_absorb() * flux * + macro_xs.get_xs( + MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * + macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, + nullptr, nullptr, &d) / + abs_xs; } } } } } else { // Skip any non-fission events - if (!p.fission_) continue; + if (!p.fission()) + continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -1961,8 +2011,8 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, // ones are delayed. If a delayed neutron is encountered, add its // contribution to the fission bank to the score. score = 0.; - for (auto i = 0; i < p.n_bank_; ++i) { - const auto& bank = p.nu_bank_[i]; + for (auto i = 0; i < p.n_bank(); ++i) { + const auto& bank = p.nu_bank(i); auto d = bank.delayed_group - 1; if (d != -1) { if (i_nuclide >= 0) { @@ -1983,8 +2033,8 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto dg = filt.groups()[d_bin]; if (dg == d + 1) - score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p.filter_matches_); + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } score = 0.; } @@ -2014,7 +2064,8 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); + score_fission_delayed_dg( + i_tally, d_bin, score, score_index, p.filter_matches()); } continue; } else { @@ -2041,14 +2092,15 @@ score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - score = p.wgt_absorb_ * flux; + score = p.wgt_absorb() * flux; } else { // Skip any non-absorption events - if (p.event_ == TallyEvent::SCATTER) continue; + if (p.event() == TallyEvent::SCATTER) + continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - score = p.wgt_last_ * flux; + score = p.wgt_last() * flux; } if (i_nuclide >= 0) { score *= atom_density @@ -2092,7 +2144,7 @@ score_all_nuclides(Particle& p, int i_tally, double flux, int filter_index, double filter_weight) { const Tally& tally {*model::tallies[i_tally]}; - const Material& material {*model::materials[p.material_]}; + const Material& material {*model::materials[p.material()]}; // Score all individual nuclide reaction rates. for (auto i = 0; i < material.nuclide_.size(); ++i) { @@ -2129,8 +2181,9 @@ void score_analog_tally_ce(Particle& p) // Note that the heating score does NOT use the flux and will be non-zero for // electrons/positrons. double flux = - (p.type_ == Particle::Type::neutron || p.type_ == Particle::Type::photon) ? - 1.0 : 0.0; + (p.type() == Particle::Type::neutron || p.type() == Particle::Type::photon) + ? 1.0 + : 0.0; for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2139,7 +2192,7 @@ void score_analog_tally_ce(Particle& p) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p.filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches()); if (filter_iter == end) continue; // Loop over filter bins. @@ -2155,7 +2208,7 @@ void score_analog_tally_ce(Particle& p) // Tally this event in the present nuclide bin if that bin represents // the event nuclide or the total material. Note that the atomic // density argument for score_general is not used for analog tallies. - if (i_nuclide == p.event_nuclide_ || i_nuclide == -1) + if (i_nuclide == p.event_nuclide() || i_nuclide == -1) score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, filter_weight, i_nuclide, -1.0, flux); } @@ -2164,7 +2217,7 @@ void score_analog_tally_ce(Particle& p) // In the case that the user has requested to tally all nuclides, we // can take advantage of the fact that we know exactly how nuclide // bins correspond to nuclide indices. First, tally the nuclide. - auto i = p.event_nuclide_; + auto i = p.event_nuclide(); score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, filter_weight, -1, -1.0, flux); @@ -2183,7 +2236,7 @@ void score_analog_tally_ce(Particle& p) } // Reset all the filter matches for the next tally event. - for (auto& match : p.filter_matches_) + for (auto& match : p.filter_matches()) match.bins_present_ = false; } @@ -2196,7 +2249,7 @@ void score_analog_tally_mg(Particle& p) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p.filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches()); if (filter_iter == end) continue; // Loop over filter bins. @@ -2210,9 +2263,10 @@ void score_analog_tally_mg(Particle& p) double atom_density = 0.; if (i_nuclide >= 0) { - auto j = model::materials[p.material_]->mat_nuclide_index_[i_nuclide]; + auto j = + model::materials[p.material()]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p.material_]->atom_density_(j); + atom_density = model::materials[p.material()]->atom_density_(j); } score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index, @@ -2228,7 +2282,7 @@ void score_analog_tally_mg(Particle& p) } // Reset all the filter matches for the next tally event. - for (auto& match : p.filter_matches_) + for (auto& match : p.filter_matches()) match.bins_present_ = false; } @@ -2236,7 +2290,7 @@ void score_tracklength_tally(Particle& p, double distance) { // Determine the tracklength estimate of the flux - double flux = p.wgt_ * distance; + double flux = p.wgt() * distance; for (auto i_tally : model::active_tracklength_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2245,7 +2299,7 @@ score_tracklength_tally(Particle& p, double distance) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p.filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches()); if (filter_iter == end) continue; // Loop over filter bins. @@ -2255,7 +2309,7 @@ score_tracklength_tally(Particle& p, double distance) // Loop over nuclide bins. if (tally.all_nuclides_) { - if (p.material_ != MATERIAL_VOID) + if (p.material() != MATERIAL_VOID) score_all_nuclides(p, i_tally, flux*filter_weight, filter_index, filter_weight); @@ -2265,10 +2319,11 @@ score_tracklength_tally(Particle& p, double distance) double atom_density = 0.; if (i_nuclide >= 0) { - if (p.material_ != MATERIAL_VOID) { - auto j = model::materials[p.material_]->mat_nuclide_index_[i_nuclide]; + if (p.material() != MATERIAL_VOID) { + auto j = + model::materials[p.material()]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p.material_]->atom_density_(j); + atom_density = model::materials[p.material()]->atom_density_(j); } } @@ -2293,7 +2348,7 @@ score_tracklength_tally(Particle& p, double distance) } // Reset all the filter matches for the next tally event. - for (auto& match : p.filter_matches_) + for (auto& match : p.filter_matches()) match.bins_present_ = false; } @@ -2301,11 +2356,12 @@ void score_collision_tally(Particle& p) { // Determine the collision estimate of the flux double flux = 0.0; - if (p.type_ == Particle::Type::neutron || p.type_ == Particle::Type::photon) { + if (p.type() == Particle::Type::neutron || + p.type() == Particle::Type::photon) { if (!settings::survival_biasing) { - flux = p.wgt_last_ / p.macro_xs_.total; + flux = p.wgt_last() / p.macro_xs().total; } else { - flux = (p.wgt_last_ + p.wgt_absorb_) / p.macro_xs_.total; + flux = (p.wgt_last() + p.wgt_absorb()) / p.macro_xs().total; } } @@ -2316,7 +2372,7 @@ void score_collision_tally(Particle& p) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p.filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches()); if (filter_iter == end) continue; // Loop over filter bins. @@ -2335,9 +2391,10 @@ void score_collision_tally(Particle& p) double atom_density = 0.; if (i_nuclide >= 0) { - auto j = model::materials[p.material_]->mat_nuclide_index_[i_nuclide]; + auto j = + model::materials[p.material()]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p.material_]->atom_density_(j); + atom_density = model::materials[p.material()]->atom_density_(j); } //TODO: consider replacing this "if" with pointers or templates @@ -2360,7 +2417,7 @@ void score_collision_tally(Particle& p) } // Reset all the filter matches for the next tally event. - for (auto& match : p.filter_matches_) + for (auto& match : p.filter_matches()) match.bins_present_ = false; } @@ -2368,7 +2425,7 @@ void score_surface_tally(Particle& p, const std::vector& tallies) { // No collision, so no weight change when survival biasing - double flux = p.wgt_; + double flux = p.wgt(); for (auto i_tally : tallies) { auto& tally {*model::tallies[i_tally]}; @@ -2377,7 +2434,7 @@ score_surface_tally(Particle& p, const std::vector& tallies) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p.filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches()); if (filter_iter == end) continue; // Loop over filter bins. @@ -2404,7 +2461,7 @@ score_surface_tally(Particle& p, const std::vector& tallies) } // Reset all the filter matches for the next tally event. - for (auto& match : p.filter_matches_) + for (auto& match : p.filter_matches()) match.bins_present_ = false; } diff --git a/src/track_output.cpp b/src/track_output.cpp index 314e43351..0a3e06f00 100644 --- a/src/track_output.cpp +++ b/src/track_output.cpp @@ -25,35 +25,35 @@ namespace openmc { void add_particle_track(Particle& p) { - p.tracks_.emplace_back(); + p.tracks().emplace_back(); } void write_particle_track(Particle& p) { - p.tracks_.back().push_back(p.r()); + p.tracks().back().push_back(p.r()); } void finalize_particle_track(Particle& p) { - std::string filename = fmt::format("{}track_{}_{}_{}.h5", - settings::path_output, simulation::current_batch, simulation::current_gen, - p.id_); + std::string filename = + fmt::format("{}track_{}_{}_{}.h5", settings::path_output, + simulation::current_batch, simulation::current_gen, p.id()); // Determine number of coordinates for each particle std::vector n_coords; - for (auto& coords : p.tracks_) { + for (auto& coords : p.tracks()) { n_coords.push_back(coords.size()); } - #pragma omp critical (FinalizeParticleTrack) +#pragma omp critical (FinalizeParticleTrack) { hid_t file_id = file_open(filename, 'w'); write_attribute(file_id, "filetype", "track"); write_attribute(file_id, "version", VERSION_TRACK); - write_attribute(file_id, "n_particles", p.tracks_.size()); + write_attribute(file_id, "n_particles", p.tracks().size()); write_attribute(file_id, "n_coords", n_coords); - for (auto i = 1; i <= p.tracks_.size(); ++i) { - const auto& t {p.tracks_[i-1]}; + for (auto i = 1; i <= p.tracks().size(); ++i) { + const auto& t {p.tracks()[i - 1]}; size_t n = t.size(); xt::xtensor data({n,3}); for (int j = 0; j < n; ++j) { @@ -68,7 +68,7 @@ void finalize_particle_track(Particle& p) } // Clear particle tracks - p.tracks_.clear(); + p.tracks().clear(); } } // namespace openmc diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 15096544d..ecb610ec4 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -129,7 +129,7 @@ std::vector VolumeCalculation::execute() const int64_t id = iterations * n_samples_ + i; uint64_t seed = init_seed(id, STREAM_VOLUME); - p.n_coord_ = 1; + p.n_coord() = 1; Position xi {prn(&seed), prn(&seed), prn(&seed)}; p.r() = lower_left_ + xi*(upper_right_ - lower_left_); p.u() = {0.5, 0.5, 0.5}; @@ -139,28 +139,33 @@ std::vector VolumeCalculation::execute() const continue; if (domain_type_ == TallyDomain::MATERIAL) { - if (p.material_ != MATERIAL_VOID) { + if (p.material() != MATERIAL_VOID) { for (int i_domain = 0; i_domain < n; i_domain++) { - if (model::materials[p.material_]->id_ == domain_ids_[i_domain]) { - this->check_hit(p.material_, indices[i_domain], hits[i_domain]); + if (model::materials[p.material()]->id_ == + domain_ids_[i_domain]) { + this->check_hit( + p.material(), indices[i_domain], hits[i_domain]); break; } } } } else if (domain_type_ == TallyDomain::CELL) { - for (int level = 0; level < p.n_coord_; ++level) { + for (int level = 0; level < p.n_coord(); ++level) { for (int i_domain=0; i_domain < n; i_domain++) { - if (model::cells[p.coord_[level].cell]->id_ == domain_ids_[i_domain]) { - this->check_hit(p.material_, indices[i_domain], hits[i_domain]); + if (model::cells[p.coord(level).cell]->id_ == + domain_ids_[i_domain]) { + this->check_hit( + p.material(), indices[i_domain], hits[i_domain]); break; } } } } else if (domain_type_ == TallyDomain::UNIVERSE) { - for (int level = 0; level < p.n_coord_; ++level) { + for (int level = 0; level < p.n_coord(); ++level) { for (int i_domain = 0; i_domain < n; ++i_domain) { - if (model::universes[p.coord_[level].universe]->id_ == domain_ids_[i_domain]) { - check_hit(p.material_, indices[i_domain], hits[i_domain]); + if (model::universes[p.coord(level).universe]->id_ == + domain_ids_[i_domain]) { + check_hit(p.material(), indices[i_domain], hits[i_domain]); break; } }