From b4ed267d4bbd8c02b6c627029874718aab736733 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 27 Feb 2019 08:23:43 -0600 Subject: [PATCH] Add trailing underscore on Particle data members --- include/openmc/particle.h | 79 ++-- src/bremsstrahlung.cpp | 14 +- src/cross_sections.cpp | 2 +- src/dagmc.cpp | 2 +- src/eigenvalue.cpp | 2 +- src/geometry.cpp | 228 ++++++------ src/material.cpp | 20 +- src/mesh.cpp | 12 +- src/output.cpp | 46 +-- src/particle.cpp | 324 ++++++++-------- src/particle_restart.cpp | 32 +- src/physics.cpp | 176 ++++----- src/physics_common.cpp | 14 +- src/physics_mg.cpp | 66 ++-- src/plot.cpp | 44 +-- src/simulation.cpp | 16 +- src/tallies/derivative.cpp | 88 ++--- src/tallies/filter_azimuthal.cpp | 4 +- src/tallies/filter_cell.cpp | 4 +- src/tallies/filter_cellborn.cpp | 2 +- src/tallies/filter_cellfrom.cpp | 4 +- src/tallies/filter_distribcell.cpp | 14 +- src/tallies/filter_energy.cpp | 16 +- src/tallies/filter_energyfunc.cpp | 6 +- src/tallies/filter_legendre.cpp | 2 +- src/tallies/filter_material.cpp | 2 +- src/tallies/filter_mesh.cpp | 2 +- src/tallies/filter_mu.cpp | 4 +- src/tallies/filter_particle.cpp | 2 +- src/tallies/filter_polar.cpp | 4 +- src/tallies/filter_sph_harm.cpp | 4 +- src/tallies/filter_sptl_legendre.cpp | 6 +- src/tallies/filter_surface.cpp | 4 +- src/tallies/filter_universe.cpp | 4 +- src/tallies/filter_zernike.cpp | 8 +- src/tallies/tally_scoring.cpp | 534 +++++++++++++-------------- src/track_output.cpp | 4 +- src/volume_calc.cpp | 28 +- 38 files changed, 912 insertions(+), 911 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 32d6c843b..c747b74a7 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -57,70 +57,71 @@ struct LocalCoord { //! State of a particle being transported through geometry //============================================================================ -struct Particle { - int64_t id; //!< Unique ID - int type; //!< Particle type (n, p, e, etc.) +class Particle { +public: + int64_t id_; //!< Unique ID + int type_; //!< Particle type (n, p, e, etc.) - int n_coord; //!< number of current coordinate levels - int cell_instance; //!< offset for distributed properties - LocalCoord coord[MAX_COORD]; //!< coordinates for all levels + int n_coord_; //!< number of current coordinate levels + int cell_instance_; //!< offset for distributed properties + LocalCoord coord_[MAX_COORD]; //!< coordinates for all levels // Particle coordinates before crossing a surface - int last_n_coord; //!< number of current coordinates - int last_cell[MAX_COORD]; //!< coordinates for all levels + int last_n_coord_; //!< number of current coordinates + int last_cell_[MAX_COORD]; //!< coordinates for all levels // Energy data - double E; //!< post-collision energy in eV - double last_E; //!< pre-collision energy in eV - int g; //!< post-collision energy group (MG only) - int last_g; //!< pre-collision energy group (MG only) + double E_; //!< post-collision energy in eV + double last_E_; //!< pre-collision energy in eV + int g_; //!< post-collision energy group (MG only) + int last_g_; //!< pre-collision energy group (MG only) // Other physical data - double wgt; //!< particle weight - double mu; //!< angle of scatter - bool alive; //!< is particle alive? + double wgt_; //!< particle weight + double mu_; //!< angle of scatter + bool alive_; //!< is particle alive? // Other physical data - double last_xyz_current[3]; //!< coordinates of the last collision or + double last_xyz_current_[3]; //!< coordinates of the last collision or //!< reflective/periodic surface crossing for //!< current tallies - double last_xyz[3]; //!< previous coordinates - double last_uvw[3]; //!< previous direction coordinates - double last_wgt; //!< pre-collision particle weight - double absorb_wgt; //!< weight absorbed for survival biasing + double last_xyz_[3]; //!< previous coordinates + double last_uvw_[3]; //!< previous direction coordinates + double last_wgt_; //!< pre-collision particle weight + double absorb_wgt_; //!< weight absorbed for survival biasing // What event took place - bool fission; //!< did particle cause implicit fission - int event; //!< scatter, absorption - int event_nuclide; //!< index in nuclides array - int event_MT; //!< reaction MT - int delayed_group; //!< delayed group + bool fission_; //!< did particle cause implicit fission + int event_; //!< scatter, absorption + int event_nuclide_; //!< index in nuclides array + int event_mt_; //!< reaction MT + int delayed_group_; //!< delayed group // Post-collision physical data - int n_bank; //!< number of fission sites banked - double wgt_bank; //!< weight of fission sites banked - int n_delayed_bank[MAX_DELAYED_GROUPS]; //!< number of delayed fission + int n_bank_; //!< number of fission sites banked + double wgt_bank_; //!< weight of fission sites banked + int n_delayed_bank_[MAX_DELAYED_GROUPS]; //!< number of delayed fission //!< sites banked // Indices for various arrays - int surface; //!< index for surface particle is on - int cell_born; //!< index for cell particle was born in - int material; //!< index for current material - int last_material; //!< index for last material + int surface_; //!< index for surface particle is on + int cell_born_; //!< index for cell particle was born in + int material_; //!< index for current material + int last_material_; //!< index for last material // Temperature of current cell - double sqrtkT; //!< sqrt(k_Boltzmann * temperature) in eV - double last_sqrtkT; //!< last temperature + double sqrtkT_; //!< sqrt(k_Boltzmann * temperature) in eV + double last_sqrtkT_; //!< last temperature // Statistical data - int n_collision; //!< number of collisions + int n_collision_; //!< number of collisions - // Track output - bool write_track {false}; +// Track output + bool write_track_ {false}; // Secondary particles created - int64_t n_secondary {}; - Bank secondary_bank[MAX_SECONDARY]; + int64_t n_secondary_ {}; + Bank secondary_bank_[MAX_SECONDARY]; //! resets all coordinate levels for the particle void clear(); diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index 2b44154da..1b22c24ae 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -28,20 +28,20 @@ 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(ParticleType::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 == static_cast(ParticleType::positron)) { - mat = &model::materials[p.material]->ttb_->positron; + if (p.type_ == static_cast(ParticleType::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 @@ -109,7 +109,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) if (w > settings::energy_cutoff[photon]) { // Create secondary photon int photon_ = static_cast(ParticleType::photon); - p.create_secondary(p.coord[0].uvw, w, photon_, true); + p.create_secondary(p.coord_[0].uvw, w, photon_, true); *E_lost += w; } } diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index 7cbb0789c..4027a6d38 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -312,7 +312,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, } } // thermal_tables_ - // Finish setting up materials (normalizing densities, etc.) + // Finish setting up.material_s (normalizing densities, etc.) mat->finalize(); } // materials diff --git a/src/dagmc.cpp b/src/dagmc.cpp index a9d116dfa..e87752261 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -221,7 +221,7 @@ void load_dagmc_geometry() c->material_.push_back(MATERIAL_VOID); } else { if (using_uwuw) { - // lookup material in uwuw if the were present + // lookup.material_ in uwuw if the were present std::string uwuw_mat = DMD.volume_material_property_data_eh[vol_handle]; if (uwuw.material_library.count(uwuw_mat) != 0) { // Note: material numbers are set by UWUW diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index ccf70ee2e..1cf159889 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -608,7 +608,7 @@ double ufs_get_weight(const Particle* p) auto& m = model::meshes[settings::index_ufs_mesh]; // Determine indices on ufs mesh for current location - int mesh_bin = m->get_bin({p->coord[0].xyz}); + int mesh_bin = m->get_bin({p->coord_[0].xyz}); if (mesh_bin < 0) { p->write_restart(); fatal_error("Source site outside UFS mesh!"); diff --git a/src/geometry.cpp b/src/geometry.cpp index fb081c313..3598f8d5b 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -34,21 +34,21 @@ std::vector overlap_check_count; extern "C" bool check_cell_overlap(Particle* p) { - 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]; int n = univ.cells_.size(); // 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].xyz, p->coord[j].uvw, p->surface)) { - if (index_cell != p->coord[j].cell) { + if (c.contains(p->coord_[j].xyz, p->coord_[j].uvw, p->surface_)) { + if (index_cell != p->coord_[j].cell) { std::stringstream err_msg; err_msg << "Overlapping cells detected: " << c.id_ << ", " - << model::cells[p->coord[j].cell]->id_ << " on universe " + << model::cells[p->coord_[j].cell]->id_ << " on universe " << univ.id_; fatal_error(err_msg); } @@ -74,36 +74,36 @@ 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->coord[p->n_coord-1].xyz}; - Direction u {p->coord[p->n_coord-1].uvw}; - auto surf = p->surface; + Position r {p->coord_[p->n_coord_-1].xyz}; + Direction u {p->coord_[p->n_coord_-1].uvw}; + 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; } } } else { - int i_universe = p->coord[p->n_coord-1].universe; + int i_universe = p->coord_[p->n_coord_-1].universe; const auto& cells {model::universes[i_universe]->cells_}; for (auto it = cells.cbegin(); it != cells.cend(); it++) { 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->coord[p->n_coord-1].xyz}; - Direction u {p->coord[p->n_coord-1].uvw}; - auto surf = p->surface; + Position r {p->coord_[p->n_coord_-1].xyz}; + Direction u {p->coord_[p->n_coord_-1].uvw}; + 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; } @@ -126,37 +126,37 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) // Find the distribcell instance number. if (c.material_.size() > 1 || c.sqrtkT_.size() > 1) { int offset = 0; - for (int i = 0; i < p->n_coord; i++) { - Cell& c_i {*model::cells[p->coord[i].cell]}; + for (int i = 0; i < p->n_coord_; i++) { + Cell& 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) { - Lattice& lat {*model::lattices[p->coord[i+1].lattice-1]}; - int i_xyz[3] {p->coord[i+1].lattice_x, - p->coord[i+1].lattice_y, - p->coord[i+1].lattice_z}; + Lattice& lat {*model::lattices[p->coord_[i+1].lattice-1]}; + int i_xyz[3] {p->coord_[i+1].lattice_x, + p->coord_[i+1].lattice_y, + p->coord_[i+1].lattice_z}; if (lat.are_valid_indices(i_xyz)) { offset += lat.offset(c.distribcell_index_, i_xyz); } } } - p->cell_instance = offset; + p->cell_instance_ = offset; } else { - p->cell_instance = 0; + p->cell_instance_ = 0; } // Set the material and temperature. - p->last_material = p->material; + p->last_material_ = 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->last_sqrtkT = p->sqrtkT; + p->last_sqrtkT_ = 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; @@ -166,44 +166,44 @@ 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. - p->coord[p->n_coord].universe = c.fill_; + p->coord_[p->n_coord_].universe = c.fill_; // Set the position and direction. for (int i = 0; i < 3; i++) { - p->coord[p->n_coord].xyz[i] = p->coord[p->n_coord-1].xyz[i]; - p->coord[p->n_coord].uvw[i] = p->coord[p->n_coord-1].uvw[i]; + p->coord_[p->n_coord_].xyz[i] = p->coord_[p->n_coord_-1].xyz[i]; + p->coord_[p->n_coord_].uvw[i] = p->coord_[p->n_coord_-1].uvw[i]; } // Apply translation. - p->coord[p->n_coord].xyz[0] -= c.translation_.x; - p->coord[p->n_coord].xyz[1] -= c.translation_.y; - p->coord[p->n_coord].xyz[2] -= c.translation_.z; + p->coord_[p->n_coord_].xyz[0] -= c.translation_.x; + p->coord_[p->n_coord_].xyz[1] -= c.translation_.y; + p->coord_[p->n_coord_].xyz[2] -= c.translation_.z; // Apply rotation. if (!c.rotation_.empty()) { - auto x = p->coord[p->n_coord].xyz[0]; - auto y = p->coord[p->n_coord].xyz[1]; - auto z = p->coord[p->n_coord].xyz[2]; - p->coord[p->n_coord].xyz[0] = x*c.rotation_[3] + y*c.rotation_[4] + auto x = p->coord_[p->n_coord_].xyz[0]; + auto y = p->coord_[p->n_coord_].xyz[1]; + auto z = p->coord_[p->n_coord_].xyz[2]; + p->coord_[p->n_coord_].xyz[0] = x*c.rotation_[3] + y*c.rotation_[4] + z*c.rotation_[5]; - p->coord[p->n_coord].xyz[1] = x*c.rotation_[6] + y*c.rotation_[7] + p->coord_[p->n_coord_].xyz[1] = x*c.rotation_[6] + y*c.rotation_[7] + z*c.rotation_[8]; - p->coord[p->n_coord].xyz[2] = x*c.rotation_[9] + y*c.rotation_[10] + p->coord_[p->n_coord_].xyz[2] = x*c.rotation_[9] + y*c.rotation_[10] + z*c.rotation_[11]; - auto u = p->coord[p->n_coord].uvw[0]; - auto v = p->coord[p->n_coord].uvw[1]; - auto w = p->coord[p->n_coord].uvw[2]; - p->coord[p->n_coord].uvw[0] = u*c.rotation_[3] + v*c.rotation_[4] + auto u = p->coord_[p->n_coord_].uvw[0]; + auto v = p->coord_[p->n_coord_].uvw[1]; + auto w = p->coord_[p->n_coord_].uvw[2]; + p->coord_[p->n_coord_].uvw[0] = u*c.rotation_[3] + v*c.rotation_[4] + w*c.rotation_[5]; - p->coord[p->n_coord].uvw[1] = u*c.rotation_[6] + v*c.rotation_[7] + p->coord_[p->n_coord_].uvw[1] = u*c.rotation_[6] + v*c.rotation_[7] + w*c.rotation_[8]; - p->coord[p->n_coord].uvw[2] = u*c.rotation_[9] + v*c.rotation_[10] + p->coord_[p->n_coord_].uvw[2] = u*c.rotation_[9] + v*c.rotation_[10] + w*c.rotation_[11]; - p->coord[p->n_coord].rotated = true; + p->coord_[p->n_coord_].rotated = true; } // Update the coordinate level and recurse. - ++p->n_coord; + ++p->n_coord_; return find_cell_inner(p, nullptr); } else if (c.type_ == FILL_LATTICE) { @@ -213,35 +213,35 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) Lattice& lat {*model::lattices[c.fill_]}; // Determine lattice indices. - Position r {p->coord[p->n_coord-1].xyz}; - Direction u {p->coord[p->n_coord-1].uvw}; + Position r {p->coord_[p->n_coord_-1].xyz}; + Direction u {p->coord_[p->n_coord_-1].uvw}; r += TINY_BIT * u; auto i_xyz = lat.get_indices(r); // Store lower level coordinates. - r = lat.get_local_position(p->coord[p->n_coord-1].xyz, i_xyz); - p->coord[p->n_coord].xyz[0] = r.x; - p->coord[p->n_coord].xyz[1] = r.y; - p->coord[p->n_coord].xyz[2] = r.z; - p->coord[p->n_coord].uvw[0] = u.x; - p->coord[p->n_coord].uvw[1] = u.y; - p->coord[p->n_coord].uvw[2] = u.z; + r = lat.get_local_position(p->coord_[p->n_coord_-1].xyz, i_xyz); + p->coord_[p->n_coord_].xyz[0] = r.x; + p->coord_[p->n_coord_].xyz[1] = r.y; + p->coord_[p->n_coord_].xyz[2] = r.z; + p->coord_[p->n_coord_].uvw[0] = u.x; + p->coord_[p->n_coord_].uvw[1] = u.y; + p->coord_[p->n_coord_].uvw[2] = u.z; // Set lattice indices. - p->coord[p->n_coord].lattice = c.fill_ + 1; - p->coord[p->n_coord].lattice_x = i_xyz[0]; - p->coord[p->n_coord].lattice_y = i_xyz[1]; - p->coord[p->n_coord].lattice_z = i_xyz[2]; + p->coord_[p->n_coord_].lattice = c.fill_ + 1; + p->coord_[p->n_coord_].lattice_x = i_xyz[0]; + p->coord_[p->n_coord_].lattice_y = i_xyz[1]; + p->coord_[p->n_coord_].lattice_z = i_xyz[2]; // Set the lower coordinate level universe. if (lat.are_valid_indices(i_xyz)) { - p->coord[p->n_coord].universe = lat[i_xyz]; + p->coord_[p->n_coord_].universe = lat[i_xyz]; } else { if (lat.outer_ != NO_OUTER_UNIVERSE) { - p->coord[p->n_coord].universe = lat.outer_; + p->coord_[p->n_coord_].universe = lat.outer_; } else { std::stringstream err_msg; - err_msg << "Particle " << p->id << " is outside lattice " + err_msg << "Particle " << p->id_ << " is outside lattice " << lat.id_ << " but the lattice has no defined outer " "universe."; warning(err_msg); @@ -250,7 +250,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) } // Update the coordinate level and recurse. - ++p->n_coord; + ++p->n_coord_; return find_cell_inner(p, nullptr); } } @@ -264,22 +264,22 @@ extern "C" bool find_cell(Particle* p, bool use_neighbor_lists) { // Determine universe (if not yet set, use root universe). - 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 < MAX_COORD; i++) { - p->coord[i].reset(); + for (int i = p->n_coord_; i < MAX_COORD; i++) { + p->coord_[i].reset(); } if (use_neighbor_lists) { // 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 @@ -291,7 +291,7 @@ find_cell(Particle* p, bool use_neighbor_lists) // cells in this universe, and update the neighbor list if we find a new // neighboring cell. found = find_cell_inner(p, nullptr); - if (found) c.neighbors_.push_back(p->coord[coord_lvl].cell); + if (found) c.neighbors_.push_back(p->coord_[coord_lvl].cell); return found; } else { @@ -305,55 +305,55 @@ find_cell(Particle* p, bool use_neighbor_lists) extern "C" void cross_lattice(Particle* p, int lattice_translation[3]) { - Lattice& lat {*model::lattices[p->coord[p->n_coord-1].lattice-1]}; + Lattice& lat {*model::lattices[p->coord_[p->n_coord_-1].lattice-1]}; if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; msg << " Crossing lattice " << lat.id_ << ". Current position (" - << p->coord[p->n_coord-1].lattice_x << "," - << p->coord[p->n_coord-1].lattice_y << "," - << p->coord[p->n_coord-1].lattice_z << ")"; + << p->coord_[p->n_coord_-1].lattice_x << "," + << p->coord_[p->n_coord_-1].lattice_y << "," + << p->coord_[p->n_coord_-1].lattice_z << ")"; write_message(msg, 1); } // Set the lattice indices. - p->coord[p->n_coord-1].lattice_x += lattice_translation[0]; - p->coord[p->n_coord-1].lattice_y += lattice_translation[1]; - p->coord[p->n_coord-1].lattice_z += lattice_translation[2]; - std::array i_xyz {p->coord[p->n_coord-1].lattice_x, - p->coord[p->n_coord-1].lattice_y, - p->coord[p->n_coord-1].lattice_z}; + p->coord_[p->n_coord_-1].lattice_x += lattice_translation[0]; + p->coord_[p->n_coord_-1].lattice_y += lattice_translation[1]; + p->coord_[p->n_coord_-1].lattice_z += lattice_translation[2]; + std::array i_xyz {p->coord_[p->n_coord_-1].lattice_x, + p->coord_[p->n_coord_-1].lattice_y, + p->coord_[p->n_coord_-1].lattice_z}; // Set the new coordinate position. - auto r = lat.get_local_position(p->coord[p->n_coord-2].xyz, i_xyz); - p->coord[p->n_coord-1].xyz[0] = r.x; - p->coord[p->n_coord-1].xyz[1] = r.y; - p->coord[p->n_coord-1].xyz[2] = r.z; + auto r = lat.get_local_position(p->coord_[p->n_coord_-2].xyz, i_xyz); + p->coord_[p->n_coord_-1].xyz[0] = r.x; + p->coord_[p->n_coord_-1].xyz[1] = r.y; + p->coord_[p->n_coord_-1].xyz[2] = r.z; 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 = find_cell(p, 0); - if (!found && p->alive) { + if (!found && p->alive_) { std::stringstream err_msg; - err_msg << "Could not locate particle " << p->id + err_msg << "Could not locate particle " << p->id_ << " after crossing a lattice boundary"; p->mark_as_lost(err_msg); } } else { // Find cell in next lattice element. - p->coord[p->n_coord-1].universe = lat[i_xyz]; + p->coord_[p->n_coord_-1].universe = lat[i_xyz]; bool found = find_cell(p, 0); 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 = find_cell(p, 0); - if (!found && p->alive) { + if (!found && p->alive_) { std::stringstream err_msg; - err_msg << "Could not locate particle " << p->id + err_msg << "Could not locate particle " << p->id_ << " after crossing a lattice boundary"; p->mark_as_lost(err_msg); } @@ -377,21 +377,21 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, std::array level_lat_trans; // Loop over each coordinate level. - for (int i = 0; i < p->n_coord; i++) { - Position r {p->coord[i].xyz}; - Direction u {p->coord[i].uvw}; - Cell& c {*model::cells[p->coord[i].cell]}; + for (int i = 0; i < p->n_coord_; i++) { + Position r {p->coord_[i].xyz}; + Direction u {p->coord_[i].uvw}; + Cell& c {*model::cells[p->coord_[i].cell]}; // Find the oncoming surface in this cell and the distance to it. - auto surface_distance = c.distance(r, u, p->surface); + auto surface_distance = c.distance(r, u, p->surface_); d_surf = surface_distance.first; level_surf_cross = surface_distance.second; // Find the distance to the next lattice tile crossing. - if (p->coord[i].lattice != F90_NONE) { - Lattice& lat {*model::lattices[p->coord[i].lattice-1]}; - std::array i_xyz {p->coord[i].lattice_x, p->coord[i].lattice_y, - p->coord[i].lattice_z}; + if (p->coord_[i].lattice != F90_NONE) { + Lattice& lat {*model::lattices[p->coord_[i].lattice-1]}; + std::array i_xyz {p->coord_[i].lattice_x, p->coord_[i].lattice_y, + p->coord_[i].lattice_z}; //TODO: refactor so both lattice use the same position argument (which //also means the lat.type attribute can be removed) std::pair> lattice_distance; @@ -400,8 +400,8 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, lattice_distance = lat.distance(r, u, i_xyz); break; case LatticeType::hex: - Position r_hex {p->coord[i-1].xyz[0], p->coord[i-1].xyz[1], - p->coord[i].xyz[2]}; + Position r_hex {p->coord_[i-1].xyz[0], p->coord_[i-1].xyz[1], + p->coord_[i].xyz[2]}; lattice_distance = lat.distance(r_hex, u, i_xyz); break; } @@ -410,7 +410,7 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, if (d_lat < 0) { std::stringstream err_msg; - err_msg << "Particle " << p->id + err_msg << "Particle " << p->id_ << " had a negative distance to a lattice boundary"; p->mark_as_lost(err_msg); } @@ -468,10 +468,10 @@ openmc_find_cell(const double* xyz, int32_t* index, int32_t* instance) Particle p; p.initialize(); - std::copy(xyz, xyz + 3, p.coord[0].xyz); - p.coord[0].uvw[0] = 0.0; - p.coord[0].uvw[1] = 0.0; - p.coord[0].uvw[2] = 1.0; + std::copy(xyz, xyz + 3, p.coord_[0].xyz); + p.coord_[0].uvw[0] = 0.0; + p.coord_[0].uvw[1] = 0.0; + p.coord_[0].uvw[2] = 1.0; if (!find_cell(&p, false)) { std::stringstream msg; @@ -481,8 +481,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 2eed3cc71..5c3f21e5b 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -652,9 +652,9 @@ void Material::calculate_xs(const Particle& p) const simulation::material_xs.fission = 0.0; simulation::material_xs.nu_fission = 0.0; - if (p.type == static_cast(ParticleType::neutron)) { + if (p.type_ == static_cast(ParticleType::neutron)) { this->calculate_neutron_xs(p); - } else if (p.type == static_cast(ParticleType::photon)) { + } else if (p.type_ == static_cast(ParticleType::photon)) { this->calculate_photon_xs(p); } } @@ -664,7 +664,7 @@ void Material::calculate_neutron_xs(const Particle& p) const int neutron = static_cast(ParticleType::neutron); // Find energy index on energy grid - 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); @@ -691,7 +691,7 @@ void Material::calculate_neutron_xs(const 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]->threshold()) i_sab = C_NONE; + if (p.E_ > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE; // Increment position in thermal_tables_ ++j; @@ -709,12 +709,12 @@ void Material::calculate_neutron_xs(const Particle& p) const // Calculate microscopic cross section for this nuclide const auto& micro {simulation::micro_xs[i_nuclide]}; - if (p.E != micro.last_E - || p.sqrtkT != micro.last_sqrtkT + 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, p.E, i_grid, - p.sqrtkT, sab_frac); + data::nuclides[i_nuclide]->calculate_xs(i_sab, p.E_, i_grid, + p.sqrtkT_, sab_frac); } // ====================================================================== @@ -748,8 +748,8 @@ void Material::calculate_photon_xs(const Particle& p) const // Calculate microscopic cross section for this nuclide const auto& micro {simulation::micro_photon_xs[i_element]}; - if (p.E != micro.last_E) { - data::elements[i_element].calculate_xs(p.E); + if (p.E_ != micro.last_E) { + data::elements[i_element].calculate_xs(p.E_); } // ======================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index 068d02918..4f0a746c7 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -385,9 +385,9 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, // just a bit for the purposes of determining if there was an intersection // in case the mesh surfaces coincide with lattice/geometric surfaces which // might produce finite-precision errors. - Position last_r {p->last_xyz}; - Position r {p->coord[0].xyz}; - Direction u {p->coord[0].uvw}; + Position last_r {p->last_xyz_}; + Position r {p->coord_[0].xyz}; + Direction u {p->coord_[0].uvw}; Position r0 = last_r + TINY_BIT*u; Position r1 = r - TINY_BIT*u; @@ -522,9 +522,9 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins // Determine if the track intersects the tally mesh. // Copy the starting and ending coordinates of the particle. - Position r0 {p->last_xyz_current}; - Position r1 {p->coord[0].xyz}; - Direction u {p->coord[0].uvw}; + Position r0 {p->last_xyz_current_}; + Position r1 {p->coord_[0].xyz}; + Direction u {p->coord_[0].uvw}; // Determine indices for starting and ending location. int n = n_dimension_; diff --git a/src/output.cpp b/src/output.cpp index 095a17526..04f4a7d46 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -142,7 +142,7 @@ std::string time_stamp() extern "C" void print_particle(Particle* p) { // Display particle type and ID. - switch (p->type) { + switch (p->type_) { case static_cast(ParticleType::neutron): std::cout << "Neutron "; break; @@ -158,48 +158,48 @@ extern "C" void print_particle(Particle* p) default: std::cout << "Unknown Particle "; } - std::cout << p->id << "\n"; + std::cout << p->id_ << "\n"; // Display particle geometry hierarchy. - for (auto i = 0; i < p->n_coord; i++) { + for (auto i = 0; i < p->n_coord_; i++) { std::cout << " Level " << i << "\n"; - 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]}; std::cout << " Cell = " << c.id_ << "\n"; } - 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]}; std::cout << " Universe = " << u.id_ << "\n"; } - if (p->coord[i].lattice != F90_NONE) { - const Lattice& lat {*model::lattices[p->coord[i].lattice]}; + if (p->coord_[i].lattice != F90_NONE) { + const Lattice& lat {*model::lattices[p->coord_[i].lattice]}; std::cout << " Lattice = " << lat.id_ << "\n"; - std::cout << " Lattice position = (" << p->coord[i].lattice_x - << "," << p->coord[i].lattice_y << "," - << p->coord[i].lattice_z << ")\n"; + std::cout << " Lattice position = (" << p->coord_[i].lattice_x + << "," << p->coord_[i].lattice_y << "," + << p->coord_[i].lattice_z << ")\n"; } - std::cout << " xyz = " << p->coord[i].xyz[0] << " " - << p->coord[i].xyz[1] << " " << p->coord[i].xyz[2] << "\n"; - std::cout << " uvw = " << p->coord[i].uvw[0] << " " - << p->coord[i].uvw[1] << " " << p->coord[i].uvw[2] << "\n"; + std::cout << " xyz = " << p->coord_[i].xyz[0] << " " + << p->coord_[i].xyz[1] << " " << p->coord_[i].xyz[2] << "\n"; + std::cout << " uvw = " << p->coord_[i].uvw[0] << " " + << p->coord_[i].uvw[1] << " " << p->coord_[i].uvw[2] << "\n"; } // Display miscellaneous info. - if (p->surface != ERROR_INT) { - const Surface& surf {*model::surfaces[std::abs(p->surface)-1]}; - std::cout << " Surface = " << std::copysign(surf.id_, p->surface) << "\n"; + if (p->surface_ != ERROR_INT) { + const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]}; + std::cout << " Surface = " << std::copysign(surf.id_, p->surface_) << "\n"; } - std::cout << " Weight = " << p->wgt << "\n"; + std::cout << " Weight = " << p->wgt_ << "\n"; if (settings::run_CE) { - std::cout << " Energy = " << p->E << "\n"; + std::cout << " Energy = " << p->E_ << "\n"; } else { - std::cout << " Energy Group = " << p->g << "\n"; + std::cout << " Energy Group = " << p->g_ << "\n"; } - std::cout << " Delayed Group = " << p->delayed_group << "\n"; + std::cout << " Delayed Group = " << p->delayed_group_ << "\n"; std::cout << "\n"; } diff --git a/src/particle.cpp b/src/particle.cpp index dfb9fca61..ed8403e14 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -51,24 +51,24 @@ void Particle::clear() { // reset any coordinate levels - for (int i=0; i(ParticleType::neutron); - alive = true; + type_ = static_cast(ParticleType::neutron); + alive_ = true; // clear attributes - surface = 0; - cell_born = C_NONE; - material = C_NONE; - last_material = C_NONE; - last_sqrtkT = 0; - wgt = 1.0; - last_wgt = 1.0; - absorb_wgt = 0.0; - n_bank = 0; - wgt_bank = 0.0; - sqrtkT = -1.0; - n_collision = 0; - fission = false; - delayed_group = 0; + surface_ = 0; + cell_born_ = C_NONE; + material_ = C_NONE; + last_material_ = C_NONE; + last_sqrtkT_ = 0; + wgt_ = 1.0; + last_wgt_ = 1.0; + absorb_wgt_ = 0.0; + n_bank_ = 0; + wgt_bank_ = 0.0; + sqrtkT_ = -1.0; + n_collision_ = 0; + fission_ = false; + delayed_group_ = 0; for (int i=0; iparticle; - wgt = src->wgt; - last_wgt = src->wgt; - std::copy(src->xyz, src->xyz + 3, coord[0].xyz); - std::copy(src->uvw, src->uvw + 3, coord[0].uvw); - std::copy(src->xyz, src->xyz + 3, last_xyz_current); - std::copy(src->xyz, src->xyz + 3, last_xyz); - std::copy(src->uvw, src->uvw + 3, last_uvw); + type_ = src->particle; + wgt_ = src->wgt; + last_wgt_ = src->wgt; + std::copy(src->xyz, src->xyz + 3, coord_[0].xyz); + std::copy(src->uvw, src->uvw + 3, coord_[0].uvw); + std::copy(src->xyz, src->xyz + 3, last_xyz_current_); + std::copy(src->xyz, src->xyz + 3, last_xyz_); + std::copy(src->uvw, src->uvw + 3, last_uvw_); if (settings::run_CE) { - E = src->E; - g = 0; + E_ = src->E; + g_ = 0; } else { - g = static_cast(src->E); - last_g = static_cast(src->E); - E = data::energy_bin_avg[g - 1]; + g_ = static_cast(src->E); + last_g_ = static_cast(src->E); + E_ = data::energy_bin_avg[g_ - 1]; } - last_E = E; + last_E_ = E_; } void @@ -138,7 +138,7 @@ Particle::transport() { // Display message if high verbosity or trace is on if (settings::verbosity >= 9 || simulation::trace) { - write_message("Simulating Particle " + std::to_string(id)); + write_message("Simulating Particle " + std::to_string(id_)); } // Initialize number of events to zero @@ -146,7 +146,7 @@ Particle::transport() // Add paricle's starting weight to count for normalizing tallies later #pragma omp atomic - simulation::total_weight += wgt; + simulation::total_weight += wgt_; // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { @@ -156,62 +156,62 @@ Particle::transport() } // Prepare to write out particle track. - if (write_track) add_particle_track(); + if (write_track_) add_particle_track(); // Every particle starts with no accumulated flux derivative. if (!model::active_tallies.empty()) zero_flux_derivs(); while (true) { // Set the random number stream - if (type == static_cast(ParticleType::neutron)) { + if (type_ == static_cast(ParticleType::neutron)) { prn_set_stream(STREAM_TRACKING); } else { prn_set_stream(STREAM_PHOTON); } // Store pre-collision particle properties - last_wgt = wgt; - last_E = E; - std::copy(coord[0].uvw, coord[0].uvw + 3, last_uvw); - std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz); + last_wgt_ = wgt_; + last_E_ = E_; + std::copy(coord_[0].uvw, coord_[0].uvw + 3, last_uvw_); + std::copy(coord_[0].xyz, coord_[0].xyz + 3, last_xyz_); // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles - if (coord[n_coord - 1].cell == C_NONE) { + if (coord_[n_coord_ - 1].cell == C_NONE) { if (!find_cell(this, false)) { this->mark_as_lost("Could not find the cell containing particle " - + std::to_string(id)); + + std::to_string(id_)); return; } // set birth cell attribute - if (cell_born == C_NONE) cell_born = coord[n_coord - 1].cell; + if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell; } // Write particle track. - if (write_track) write_particle_track(*this); + if (write_track_) write_particle_track(*this); if (settings::check_overlaps) check_cell_overlap(this); // Calculate microscopic and macroscopic cross sections - if (material != MATERIAL_VOID) { + if (material_ != MATERIAL_VOID) { if (settings::run_CE) { - if (material != last_material || sqrtkT != last_sqrtkT) { + if (material_ != last_material_ || sqrtkT_ != last_sqrtkT_) { // If the material is the same as the last material and the // temperature hasn't changed, we don't need to lookup cross // sections again. - model::materials[material]->calculate_xs(*this); + model::materials[material_]->calculate_xs(*this); } } else { // Get the MG data - calculate_xs_c(material, g, sqrtkT, coord[n_coord-1].uvw, + calculate_xs_c(material_, g_, sqrtkT_, coord_[n_coord_-1].uvw, simulation::material_xs.total, simulation::material_xs.absorption, simulation::material_xs.nu_fission); // Finally, update the particle group while we have already checked // for if multi-group - last_g = g; + last_g_ = g_; } } else { simulation::material_xs.total = 0.0; @@ -230,8 +230,8 @@ Particle::transport() // Sample a distance to collision double d_collision; - if (type == static_cast(ParticleType::electron) || - type == static_cast(ParticleType::positron)) { + if (type_ == static_cast(ParticleType::electron) || + type_ == static_cast(ParticleType::positron)) { d_collision = 0.0; } else if (simulation::material_xs.total == 0.0) { d_collision = INFINITY; @@ -243,11 +243,11 @@ Particle::transport() double distance = std::min(d_boundary, d_collision); // Advance particle - for (int j = 0; j < n_coord; ++j) { + for (int j = 0; j < n_coord_; ++j) { // TODO: use Position - coord[j].xyz[0] += distance * coord[j].uvw[0]; - coord[j].xyz[1] += distance * coord[j].uvw[1]; - coord[j].xyz[2] += distance * coord[j].uvw[2]; + coord_[j].xyz[0] += distance * coord_[j].uvw[0]; + coord_[j].xyz[1] += distance * coord_[j].uvw[1]; + coord_[j].xyz[2] += distance * coord_[j].uvw[2]; } // Score track-length tallies @@ -257,8 +257,8 @@ Particle::transport() // Score track-length estimate of k-eff if (settings::run_mode == RUN_MODE_EIGENVALUE && - type == static_cast(ParticleType::neutron)) { - global_tally_tracklength += wgt * distance * simulation::material_xs.nu_fission; + type_ == static_cast(ParticleType::neutron)) { + global_tally_tracklength += wgt_ * distance * simulation::material_xs.nu_fission; } // Score flux derivative accumulators for differential tallies. @@ -270,25 +270,25 @@ Particle::transport() // ==================================================================== // PARTICLE CROSSES SURFACE - if (next_level > 0) n_coord = next_level; + if (next_level > 0) n_coord_ = next_level; // Saving previous cell data - for (int j = 0; j < n_coord; ++j) { - last_cell[j] = coord[j].cell; + for (int j = 0; j < n_coord_; ++j) { + last_cell_[j] = coord_[j].cell; } - last_n_coord = n_coord; + last_n_coord_ = n_coord_; if (lattice_translation[0] != 0 || lattice_translation[1] != 0 || lattice_translation[2] != 0) { // Particle crosses lattice boundary - surface = ERROR_INT; + surface_ = ERROR_INT; cross_lattice(this, lattice_translation); - event = EVENT_LATTICE; + event_ = EVENT_LATTICE; } else { // Particle crosses surface - surface = surface_crossed; + surface_ = surface_crossed; this->cross_surface(); - event = EVENT_SURFACE; + event_ = EVENT_SURFACE; } // Score cell to cell partial currents if (!model::active_surface_tallies.empty()) { @@ -300,8 +300,8 @@ Particle::transport() // Score collision estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE && - type == static_cast(ParticleType::neutron)) { - global_tally_collision += wgt * simulation::material_xs.nu_fission + type_ == static_cast(ParticleType::neutron)) { + global_tally_collision += wgt_ * simulation::material_xs.nu_fission / simulation::material_xs.total; } @@ -313,7 +313,7 @@ Particle::transport() score_surface_tally(this, model::active_meshsurf_tallies); // Clear surface component - surface = ERROR_INT; + surface_ = ERROR_INT; if (settings::run_CE) { collision(this); @@ -334,33 +334,33 @@ Particle::transport() } // Reset banked weight during collision - n_bank = 0; - wgt_bank = 0.0; - for (int& v : n_delayed_bank) v = 0; + n_bank_ = 0; + wgt_bank_ = 0.0; + for (int& v : n_delayed_bank_) v = 0; // Reset fission logical - fission = false; + fission_ = false; // Save coordinates for tallying purposes - std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz_current); + std::copy(coord_[0].xyz, coord_[0].xyz + 3, last_xyz_current_); // Set last material to none since cross sections will need to be // re-evaluated - last_material = C_NONE; + last_material_ = C_NONE; // Set all uvws to base level -- right now, after a collision, only the // base level uvws are changed - for (int j = 0; j < n_coord - 1; ++j) { - if (coord[j + 1].rotated) { + for (int j = 0; j < n_coord_ - 1; ++j) { + if (coord_[j + 1].rotated) { // If next level is rotated, apply rotation matrix - const auto& m {model::cells[coord[j].cell]->rotation_}; - Direction u {coord[j].uvw}; - coord[j + 1].uvw[0] = m[3]*u.x + m[4]*u.y + m[5]*u.z; - coord[j + 1].uvw[1] = m[6]*u.x + m[7]*u.y + m[8]*u.z; - coord[j + 1].uvw[2] = m[9]*u.x + m[10]*u.y + m[11]*u.z; + const auto& m {model::cells[coord_[j].cell]->rotation_}; + Direction u {coord_[j].uvw}; + coord_[j + 1].uvw[0] = m[3]*u.x + m[4]*u.y + m[5]*u.z; + coord_[j + 1].uvw[1] = m[6]*u.x + m[7]*u.y + m[8]*u.z; + coord_[j + 1].uvw[2] = m[9]*u.x + m[10]*u.y + m[11]*u.z; } else { // Otherwise, copy this level's direction - std::copy(coord[j].uvw, coord[j].uvw + 3, coord[j + 1].uvw); + std::copy(coord_[j].uvw, coord_[j].uvw + 3, coord_[j + 1].uvw); } } @@ -371,27 +371,27 @@ Particle::transport() // If particle has too many events, display warning and kill it ++n_event; if (n_event == MAX_EVENTS) { - warning("Particle " + std::to_string(id) + + warning("Particle " + std::to_string(id_) + " underwent maximum number of events."); - alive = false; + alive_ = false; } // Check for secondary particles if this particle is dead - if (!alive) { + if (!alive_) { // If no secondary particles, break out of event loop - if (n_secondary == 0) break; + if (n_secondary_ == 0) break; - this->from_source(&secondary_bank[n_secondary - 1]); - --n_secondary; + this->from_source(&secondary_bank_[n_secondary_ - 1]); + --n_secondary_; n_event = 0; // Enter new particle in particle track file - if (write_track) add_particle_track(); + if (write_track_) add_particle_track(); } } // Finish particle track output. - if (write_track) { + if (write_track_) { write_particle_track(*this); finalize_particle_track(*this); } @@ -400,7 +400,7 @@ Particle::transport() void Particle::cross_surface() { - int i_surface = std::abs(surface); + int i_surface = std::abs(surface_); // TODO: off-by-one const auto& surf {model::surfaces[i_surface - 1].get()}; if (settings::verbosity >= 10 || simulation::trace) { @@ -412,7 +412,7 @@ Particle::cross_surface() // PARTICLE LEAKS OUT OF PROBLEM // Kill particle - alive = false; + alive_ = false; // Score any surface current tallies -- note that the particle is moved // forward slightly so that if the mesh boundary is on the surface, it is @@ -423,14 +423,14 @@ Particle::cross_surface() // physically moving the particle forward slightly // TODO: Use Position - coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0]; - coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1]; - coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2]; + coord_[0].xyz[0] += TINY_BIT * coord_[0].uvw[0]; + coord_[0].xyz[1] += TINY_BIT * coord_[0].uvw[1]; + coord_[0].xyz[2] += TINY_BIT * coord_[0].uvw[2]; score_surface_tally(this, model::active_meshsurf_tallies); } // Score to global leakage tally - global_tally_leakage += wgt; + global_tally_leakage += wgt_; // Display message if (settings::verbosity >= 10 || simulation::trace) { @@ -443,8 +443,8 @@ Particle::cross_surface() // PARTICLE REFLECTS FROM SURFACE // Do not handle reflective boundary conditions on lower universes - if (n_coord != 1) { - this->mark_as_lost("Cannot reflect particle " + std::to_string(id) + + if (n_coord_ != 1) { + this->mark_as_lost("Cannot reflect particle " + std::to_string(id_) + " off surface in a lower universe."); return; } @@ -462,32 +462,32 @@ Particle::cross_surface() if (!model::active_meshsurf_tallies.empty()) { - Position r {coord[0].xyz}; - coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0]; - coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1]; - coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2]; + Position r {coord_[0].xyz}; + coord_[0].xyz[0] -= TINY_BIT * coord_[0].uvw[0]; + coord_[0].xyz[1] -= TINY_BIT * coord_[0].uvw[1]; + coord_[0].xyz[2] -= TINY_BIT * coord_[0].uvw[2]; score_surface_tally(this, model::active_meshsurf_tallies); - std::copy(&r.x, &r.x + 3, coord[0].xyz); + std::copy(&r.x, &r.x + 3, coord_[0].xyz); } // Reflect particle off surface - Direction u = surf->reflect(coord[0].xyz, coord[0].uvw); + Direction u = surf->reflect(coord_[0].xyz, coord_[0].uvw); // Make sure new particle direction is normalized double norm = u.norm(); - coord[0].uvw[0] = u.x/norm; - coord[0].uvw[1] = u.y/norm; - coord[0].uvw[2] = u.z/norm; + coord_[0].uvw[0] = u.x/norm; + coord_[0].uvw[1] = u.y/norm; + coord_[0].uvw[2] = u.z/norm; // Reassign particle's cell and surface - coord[0].cell = last_cell[last_n_coord - 1]; - surface = -surface; + coord_[0].cell = last_cell_[last_n_coord_ - 1]; + surface_ = -surface_; // If a reflective surface is coincident with a lattice or universe // boundary, it is necessary to redetermine the particle's coordinates in // the lower universes. - n_coord = 1; + n_coord_ = 1; if (!find_cell(this, true)) { this->mark_as_lost("Couldn't find particle after reflecting from surface " + std::to_string(surf->id_) + "."); @@ -495,9 +495,9 @@ Particle::cross_surface() } // Set previous coordinate going slightly past surface crossing - last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT*coord[0].uvw[0]; - last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT*coord[0].uvw[1]; - last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT*coord[0].uvw[2]; + last_xyz_current_[0] = coord_[0].xyz[0] + TINY_BIT*coord_[0].uvw[0]; + last_xyz_current_[1] = coord_[0].xyz[1] + TINY_BIT*coord_[0].uvw[1]; + last_xyz_current_[2] = coord_[0].xyz[2] + TINY_BIT*coord_[0].uvw[2]; // Diagnostic message if (settings::verbosity >= 10 || simulation::trace) { @@ -510,8 +510,8 @@ Particle::cross_surface() // PERIODIC BOUNDARY // Do not handle periodic boundary conditions on lower universes - if (n_coord != 1) { - this->mark_as_lost("Cannot transfer particle " + std::to_string(id) + + if (n_coord_ != 1) { + this->mark_as_lost("Cannot transfer particle " + std::to_string(id_) + " across surface in a lower universe. Boundary conditions must be " "applied to root universe."); return; @@ -521,12 +521,12 @@ Particle::cross_surface() // particle to change -- artificially move the particle slightly back in // case the surface crossing is coincident with a mesh boundary if (!model::active_meshsurf_tallies.empty()) { - Position r {coord[0].xyz}; - coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0]; - coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1]; - coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2]; + Position r {coord_[0].xyz}; + coord_[0].xyz[0] -= TINY_BIT * coord_[0].uvw[0]; + coord_[0].xyz[1] -= TINY_BIT * coord_[0].uvw[1]; + coord_[0].xyz[2] -= TINY_BIT * coord_[0].uvw[2]; score_surface_tally(this, model::active_meshsurf_tallies); - std::copy(&r.x, &r.x + 3, coord[0].xyz); + std::copy(&r.x, &r.x + 3, coord_[0].xyz); } // Get a pointer to the partner periodic surface @@ -535,20 +535,20 @@ Particle::cross_surface() model::surfaces[surf_p->i_periodic_].get()); // Adjust the particle's location and direction. - Position r {coord[0].xyz}; - Direction u {coord[0].uvw}; + Position r {coord_[0].xyz}; + Direction u {coord_[0].uvw}; bool rotational = other->periodic_translate(surf_p, r, u); - std::copy(&r.x, &r.x + 3, coord[0].xyz); - std::copy(&u.x, &u.x + 3, coord[0].uvw); + std::copy(&r.x, &r.x + 3, coord_[0].xyz); + std::copy(&u.x, &u.x + 3, coord_[0].uvw); // Reassign particle's surface // TODO: off-by-one - surface = rotational ? + surface_ = rotational ? surf_p->i_periodic_ + 1 : - std::copysign(surf_p->i_periodic_ + 1, surface); + std::copysign(surf_p->i_periodic_ + 1, surface_); // Figure out what cell particle is in now - n_coord = 1; + n_coord_ = 1; if (!find_cell(this, true)) { this->mark_as_lost("Couldn't find particle after hitting periodic " @@ -557,9 +557,9 @@ Particle::cross_surface() } // Set previous coordinate going slightly past surface crossing - last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT * coord[0].uvw[0]; - last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT * coord[0].uvw[1]; - last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT * coord[0].uvw[2]; + last_xyz_current_[0] = coord_[0].xyz[0] + TINY_BIT * coord_[0].uvw[0]; + last_xyz_current_[1] = coord_[0].xyz[1] + TINY_BIT * coord_[0].uvw[1]; + last_xyz_current_[2] = coord_[0].xyz[2] + TINY_BIT * coord_[0].uvw[2]; // Diagnostic message if (settings::verbosity >= 10 || simulation::trace) { @@ -574,18 +574,18 @@ Particle::cross_surface() #ifdef DAGMC if (settings::dagmc) { - auto cellp = dynamic_cast(model::cells[last_cell[0]]); + auto cellp = dynamic_cast(model::cells[last_cell_[0]]); // TODO: off-by-one - auto surfp = dynamic_cast(model::surfaces[std::abs(surface) - 1]); + auto surfp = dynamic_cast(model::surfaces[std::abs(surface_) - 1]); int32_t i_cell = next_cell(cellp, surfp) - 1; // save material and temp - last_material = material; - last_sqrtkT = sqrtkT; + last_material_ = material_; + last_sqrtkT_ = sqrtkT_; // set new cell value - coord[0].cell = i_cell; - cell_instance = 0; - material = model::cells[i_cell]->material_[0]; - sqrtkT = model::cells[i_cell]->sqrtkT_[0]; + coord_[0].cell = i_cell; + cell_instance_ = 0; + material_ = model::cells[i_cell]->material_[0]; + sqrtkT_ = model::cells[i_cell]->sqrtkT_[0]; return; } #endif @@ -596,8 +596,8 @@ Particle::cross_surface() // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS // Remove lower coordinate levels and assignment of surface - surface = ERROR_INT; - n_coord = 1; + surface_ = ERROR_INT; + n_coord_ = 1; bool found = find_cell(this, false); if (settings::run_mode != RUN_MODE_PLOTTING && (!found)) { @@ -606,16 +606,16 @@ Particle::cross_surface() // the particle is really traveling tangent to a surface, if we move it // forward a tiny bit it should fix the problem. - n_coord = 1; - coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0]; - coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1]; - coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2]; + n_coord_ = 1; + coord_[0].xyz[0] += TINY_BIT * coord_[0].uvw[0]; + coord_[0].xyz[1] += TINY_BIT * coord_[0].uvw[1]; + coord_[0].xyz[2] += TINY_BIT * coord_[0].uvw[2]; // Couldn't find next cell anywhere! This probably means there is an actual // undefined region in the geometry. if (!find_cell(this, false)) { - this->mark_as_lost("After particle " + std::to_string(id) + + this->mark_as_lost("After particle " + std::to_string(id_) + " crossed surface " + std::to_string(surf->id_) + " it could not be located in any cell and it did not leak."); return; @@ -631,7 +631,7 @@ Particle::mark_as_lost(const char* message) write_restart(); // Increment number of lost particles - alive = false; + alive_ = false; #pragma omp atomic simulation::n_lost_particles += 1; @@ -655,7 +655,7 @@ Particle::write_restart() const // Set up file name std::stringstream filename; filename << settings::path_output << "particle_" << simulation::current_batch - << '_' << id << ".h5"; + << '_' << id_ << ".h5"; #pragma omp critical (WriteParticleRestart) { @@ -686,8 +686,8 @@ Particle::write_restart() const write_dataset(file_id, "run_mode", "particle restart"); break; } - write_dataset(file_id, "id", id); - write_dataset(file_id, "type", type); + write_dataset(file_id, "id", id_); + write_dataset(file_id, "type", type_); int64_t i = simulation::current_work; write_dataset(file_id, "weight", simulation::source_bank[i-1].wgt); diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 19dd50d77..90c37d92f 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -39,29 +39,29 @@ void read_particle_restart(Particle& p, int& previous_run_mode) } else if (mode == "fixed source") { previous_run_mode = RUN_MODE_FIXEDSOURCE; } - read_dataset(file_id, "id", p.id); - read_dataset(file_id, "type", p.type); - read_dataset(file_id, "weight", p.wgt); - read_dataset(file_id, "energy", p.E); + read_dataset(file_id, "id", p.id_); + read_dataset(file_id, "type", p.type_); + read_dataset(file_id, "weight", p.wgt_); + read_dataset(file_id, "energy", p.E_); std::array x; read_dataset(file_id, "xyz", x); - std::copy(x.data(), x.data() + 3, p.coord[0].xyz); + std::copy(x.data(), x.data() + 3, p.coord_[0].xyz); read_dataset(file_id, "uvw", x); - std::copy(x.data(), x.data() + 3, p.coord[0].uvw); + std::copy(x.data(), x.data() + 3, p.coord_[0].uvw); // Set energy group and average energy in multi-group mode if (!settings::run_CE) { - p.g = p.E; - p.E = data::energy_bin_avg[p.g - 1]; + p.g_ = p.E_; + p.E_ = data::energy_bin_avg[p.g_ - 1]; } // Set particle last attributes - p.last_wgt = p.wgt; - std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz_current); - std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz); - std::copy(p.coord[0].uvw, p.coord[0].uvw + 3, p.last_uvw); - p.last_E = p.E; - p.last_g = p.g; + p.last_wgt_ = p.wgt_; + std::copy(p.coord_[0].xyz, p.coord_[0].xyz + 3, p.last_xyz_current_); + std::copy(p.coord_[0].xyz, p.coord_[0].xyz + 3, p.last_xyz_); + std::copy(p.coord_[0].uvw, p.coord_[0].uvw + 3, p.last_uvw_); + p.last_E_ = p.E_; + p.last_g_ = p.g_; // Close hdf5 file file_close(file_id); @@ -94,10 +94,10 @@ void run_particle_restart() int64_t particle_seed; switch (previous_run_mode) { case RUN_MODE_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 RUN_MODE_FIXEDSOURCE: - particle_seed = p.id; + particle_seed = p.id_; break; } set_particle_seed(particle_seed); diff --git a/src/physics.cpp b/src/physics.cpp index 5d6ca2dee..0949532ff 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -33,10 +33,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 (static_cast(p->type)) { + switch (static_cast(p->type_)) { case ParticleType::neutron: sample_neutron_reaction(p); break; @@ -52,21 +52,21 @@ void collision(Particle* p) } // Kill particle if energy falls below cutoff - if (p->E < settings::energy_cutoff[p->type]) { - p->alive = false; - p->wgt = 0.0; - p->last_wgt = 0.0; + if (p->E_ < settings::energy_cutoff[p->type_]) { + p->alive_ = false; + p->wgt_ = 0.0; + p->last_wgt_ = 0.0; } // Display information about collision if (settings::verbosity >= 10 || simulation::trace) { std::stringstream msg; - if (static_cast(p->type) == ParticleType::neutron) { - msg << " " << reaction_name(p->event_MT) << " with " << - data::nuclides[p->event_nuclide]->name_ << ". Energy = " << p->E << " eV."; + if (static_cast(p->type_) == ParticleType::neutron) { + msg << " " << reaction_name(p->event_mt_) << " with " << + data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV."; } else { - msg << " " << reaction_name(p->event_MT) << " with " << - data::elements[p->event_nuclide].name_ << ". Energy = " << p->E << " eV."; + msg << " " << reaction_name(p->event_mt_) << " with " << + data::elements[p->event_nuclide_].name_ << ". Energy = " << p->E_ << " eV."; } write_message(msg, 1); } @@ -78,7 +78,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 @@ -88,14 +88,14 @@ void sample_neutron_reaction(Particle* p) const auto& nuc {data::nuclides[i_nuclide]}; if (nuc->fissionable_) { - Reaction* rx = sample_fission(i_nuclide, p->E); + Reaction* rx = sample_fission(i_nuclide, p->E_); if (settings::run_mode == RUN_MODE_EIGENVALUE) { create_fission_sites(p, i_nuclide, rx, simulation::fission_bank.data(), &simulation::n_bank, simulation::fission_bank.size()); } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE && settings::create_fission_neutrons) { - create_fission_sites(p, i_nuclide, rx, p->secondary_bank, - &p->n_secondary, MAX_SECONDARY); + create_fission_sites(p, i_nuclide, rx, p->secondary_bank_, + &p->n_secondary_, MAX_SECONDARY); } } @@ -112,16 +112,16 @@ void sample_neutron_reaction(Particle* p) if (simulation::micro_xs[i_nuclide].absorption > 0.0) { absorption(p, i_nuclide); } else { - p->absorb_wgt = 0.0; + p->absorb_wgt_ = 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->last_E) { + if (p->E_ != p->last_E_) { prn_set_stream(STREAM_URR_PTABLE); advance_prn_seed(data::nuclides.size()); prn_set_stream(STREAM_TRACKING); @@ -130,7 +130,7 @@ void sample_neutron_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; } } @@ -145,7 +145,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ 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 * simulation::micro_xs[ + double nu_t = p->wgt_ / simulation::keff * weight * simulation::micro_xs[ i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].total; // Sample the number of neutrons produced @@ -177,12 +177,12 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ // group. double nu_d[MAX_DELAYED_GROUPS] = {0.}; - p->fission = true; + p->fission_ = true; for (size_t i = *size_bank; i < std::min(*size_bank + nu, bank_capacity); ++i) { // Bank source neutrons by copying the particle data - bank_array[i].xyz[0] = p->coord[0].xyz[0]; - bank_array[i].xyz[1] = p->coord[0].xyz[1]; - bank_array[i].xyz[2] = p->coord[0].xyz[2]; + bank_array[i].xyz[0] = p->coord_[0].xyz[0]; + bank_array[i].xyz[1] = p->coord_[0].xyz[1]; + bank_array[i].xyz[2] = p->coord_[0].xyz[2]; // Set that the bank particle is a neutron bank_array[i].particle = static_cast(ParticleType::neutron); @@ -191,14 +191,14 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ bank_array[i].wgt = 1. / weight; // Sample delayed group and angle/energy for fission reaction - sample_fission_neutron(i_nuclide, rx, p->E, &bank_array[i]); + sample_fission_neutron(i_nuclide, rx, p->E_, &bank_array[i]); // Set the delayed group on the particle as well - p->delayed_group = bank_array[i].delayed_group; + p->delayed_group_ = bank_array[i].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]++; } } @@ -206,10 +206,10 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_ *size_bank = std::min(*size_bank + nu, bank_capacity); // 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]; } } @@ -219,20 +219,20 @@ 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(ParticleType::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); - p->event_nuclide = i_element; + p->event_nuclide_ = i_element; const auto& micro {simulation::micro_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 @@ -243,8 +243,8 @@ void sample_photon_reaction(Particle* p) prob += micro.coherent; if (prob > cutoff) { double mu = element.rayleigh_scatter(alpha); - rotate_angle_c(p->coord[0].uvw, mu, nullptr); - p->event_MT = COHERENT; + rotate_angle_c(p->coord_[0].uvw, mu, nullptr); + p->event_mt_ = COHERENT; return; } @@ -270,7 +270,7 @@ void sample_photon_reaction(Particle* p) / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu); double phi = 2.0*PI*prn(); double uvw[3]; - std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + std::copy(p->coord_[0].uvw, p->coord_[0].uvw + 3, uvw); rotate_angle_c(uvw, mu_electron, &phi); int electron = static_cast(ParticleType::electron); p->create_secondary(uvw, E_electron, electron, true); @@ -284,9 +284,9 @@ void sample_photon_reaction(Particle* p) } phi += PI; - p->E = alpha_out*MASS_ELECTRON_EV; - rotate_angle_c(p->coord[0].uvw, mu, &phi); - p->event_MT = INCOHERENT; + p->E_ = alpha_out*MASS_ELECTRON_EV; + rotate_angle_c(p->coord_[0].uvw, mu, &phi); + p->event_mt_ = INCOHERENT; return; } @@ -309,7 +309,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 @@ -338,9 +338,9 @@ void sample_photon_reaction(Particle* p) // Allow electrons to fill orbital and produce auger electrons // and fluorescent photons element.atomic_relaxation(shell, *p); - p->event_MT = 533 + shell.index_subshell; - p->alive = false; - p->E = 0.0; + p->event_mt_ = 533 + shell.index_subshell; + p->alive_ = false; + p->E_ = 0.0; return; } } @@ -357,20 +357,20 @@ void sample_photon_reaction(Particle* p) // Create secondary electron double uvw[3]; - std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + std::copy(p->coord_[0].uvw, p->coord_[0].uvw + 3, uvw); rotate_angle_c(uvw, mu_electron, nullptr); int electron = static_cast(ParticleType::electron); p->create_secondary(uvw, E_electron, electron, true); // Create secondary positron - std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + std::copy(p->coord_[0].uvw, p->coord_[0].uvw + 3, uvw); rotate_angle_c(uvw, mu_positron, nullptr); int positron = static_cast(ParticleType::positron); p->create_secondary(uvw, E_positron, positron, true); - p->event_MT = PAIR_PROD; - p->alive = false; - p->E = 0.0; + p->event_mt_ = PAIR_PROD; + p->alive_ = false; + p->E_ = 0.0; } } @@ -383,8 +383,8 @@ void sample_electron_reaction(Particle* p) thick_target_bremsstrahlung(*p, &E_lost); } - p->E = 0.0; - p->alive = false; + p->E_ = 0.0; + p->alive_ = false; } void sample_positron_reaction(Particle* p) @@ -413,8 +413,8 @@ void sample_positron_reaction(Particle* p) uvw[2] = -uvw[2]; p->create_secondary(uvw.data(), MASS_ELECTRON_EV, photon, true); - p->E = 0.0; - p->alive = false; + p->E_ = 0.0; + p->alive_ = false; } int sample_nuclide(const Particle* p) @@ -423,7 +423,7 @@ int sample_nuclide(const Particle* p) double cutoff = prn() * simulation::material_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; @@ -448,7 +448,7 @@ int sample_element(Particle* p) double cutoff = prn() * simulation::material_xs.total; // Get pointers to elements, densities - const auto& mat {model::materials[p->material]}; + const auto& mat {model::materials[p->material_]}; int n = mat->nuclide_.size(); int i = 0; @@ -557,16 +557,16 @@ void absorption(Particle* p, int i_nuclide) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->absorb_wgt = p->wgt * simulation::micro_xs[i_nuclide].absorption / + p->absorb_wgt_ = p->wgt_ * simulation::micro_xs[i_nuclide].absorption / simulation::micro_xs[i_nuclide].total; // Adjust weight of particle by probability of absorption - p->wgt -= p->absorb_wgt; - p->last_wgt = p->wgt; + p->wgt_ -= p->absorb_wgt_; + p->last_wgt_ = p->wgt_; // Score implicit absorption estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE) { - global_tally_absorption += p->absorb_wgt * simulation::micro_xs[ + global_tally_absorption += p->absorb_wgt_ * simulation::micro_xs[ i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption; } } else { @@ -575,13 +575,13 @@ void absorption(Particle* p, int i_nuclide) prn() * simulation::micro_xs[i_nuclide].total) { // Score absorption estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE) { - global_tally_absorption += p->wgt * simulation::micro_xs[ + global_tally_absorption += p->wgt_ * simulation::micro_xs[ i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption; } - p->alive = false; - p->event = EVENT_ABSORB; - p->event_MT = N_DISAPPEAR; + p->alive_ = false; + p->event_ = EVENT_ABSORB; + p->event_mt_ = N_DISAPPEAR; } } } @@ -589,7 +589,7 @@ void absorption(Particle* p, int i_nuclide) void scatter(Particle* p, int i_nuclide) { // copy incoming direction - Direction u_old {p->coord[0].uvw}; + Direction u_old {p->coord_[0].uvw}; // Get pointer to nuclide and grid index/interpolation factor const auto& nuc {data::nuclides[i_nuclide]}; @@ -614,13 +614,13 @@ 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].get(), kT, - &p->E, p->coord[0].uvw, &p->mu, &p->wgt); + &p->E_, p->coord_[0].uvw, &p->mu_, &p->wgt_); - p->event_MT = ELASTIC; + p->event_mt_ = ELASTIC; sampled = true; } @@ -629,9 +629,9 @@ void scatter(Particle* p, int i_nuclide) // ======================================================================= // S(A,B) SCATTERING - sab_scatter(i_nuclide, micro.index_sab, &p->E, p->coord[0].uvw, &p->mu); + sab_scatter(i_nuclide, micro.index_sab, &p->E_, p->coord_[0].uvw, &p->mu_); - p->event_MT = ELASTIC; + p->event_mt_ = ELASTIC; sampled = true; } @@ -663,14 +663,14 @@ void scatter(Particle* p, int i_nuclide) // Perform collision physics for inelastic scattering const auto& rx {nuc->reactions_[i]}; inelastic_scatter(nuc.get(), rx.get(), p); - p->event_MT = rx->mt_; + p->event_mt_ = rx->mt_; } // Set event component - p->event = EVENT_SCATTER; + p->event_ = EVENT_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]) { @@ -682,12 +682,12 @@ void scatter(Particle* p, int i_nuclide) u_new.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); u_new.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); - p->mu = u_old.dot(u_new); + p->mu_ = u_old.dot(u_new); // Change direction of particle - p->coord[0].uvw[0] = u_new.x; - p->coord[0].uvw[1] = u_new.y; - p->coord[0].uvw[2] = u_new.z; + p->coord_[0].uvw[0] = u_new.x; + p->coord_[0].uvw[1] = u_new.y; + p->coord_[0].uvw[2] = u_new.z; } } } @@ -1071,7 +1071,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank 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; @@ -1098,11 +1098,11 @@ 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 - rotate_angle_c(p->coord[0].uvw, mu, nullptr); + rotate_angle_c(p->coord_[0].uvw, mu, nullptr); // evaluate yield double yield = (*rx->products_[0].yield_)(E_in); @@ -1110,18 +1110,18 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) // If yield is integral, create exactly that many secondary particles for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { int neutron = static_cast(ParticleType::neutron); - p->create_secondary(p->coord[0].uvw, p->E, neutron, true); + p->create_secondary(p->coord_[0].uvw, p->E_, neutron, true); } } 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->wgt * simulation::micro_xs[i_nuclide].photon_prod / + double y_t = p->wgt_ * simulation::micro_xs[i_nuclide].photon_prod / simulation::micro_xs[i_nuclide].total; int y = static_cast(y_t); if (prn() <= y_t - y) ++y; @@ -1131,17 +1131,17 @@ void sample_secondary_photons(Particle* p, int i_nuclide) // Sample the reaction and product int i_rx; int i_product; - sample_photon_product(i_nuclide, p->E, &i_rx, &i_product); + sample_photon_product(i_nuclide, p->E_, &i_rx, &i_product); // Sample the outgoing energy and angle auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; double E; double mu; - rx->products_[i_product].sample(p->E, E, mu); + rx->products_[i_product].sample(p->E_, E, mu); // Sample the new direction double uvw[3]; - std::copy(p->coord[0].uvw, p->coord[0].uvw + 3, uvw); + std::copy(p->coord_[0].uvw, p->coord_[0].uvw + 3, uvw); rotate_angle_c(uvw, mu, nullptr); // Create the secondary photon diff --git a/src/physics_common.cpp b/src/physics_common.cpp index 02c77ae7d..42843a515 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->wgt / settings::weight_survive) { - p->wgt = settings::weight_survive; - p->last_wgt = p->wgt; + if (p->wgt_ < settings::weight_cutoff) { + if (prn() < p->wgt_ / settings::weight_survive) { + p->wgt_ = settings::weight_survive; + p->last_wgt_ = p->wgt_; } else { - p->wgt = 0.; - p->last_wgt = 0.; - p->alive = false; + p->wgt_ = 0.; + p->last_wgt_ = 0.; + p->alive_ = false; } } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 87dd37ef8..942f3f475 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -25,7 +25,7 @@ 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); @@ -33,7 +33,7 @@ collision_mg(Particle* p) // Display information about collision if ((settings::verbosity >= 10) || (simulation::trace)) { std::stringstream msg; - msg << " Energy Group = " << p->g; + msg << " Energy Group = " << p->g_; write_message(msg, 1); } } @@ -46,14 +46,14 @@ 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 == RUN_MODE_EIGENVALUE) { create_fission_sites( p, simulation::fission_bank.data(), &simulation::n_bank, simulation::fission_bank.size()); } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) && (settings::create_fission_neutrons)) { - create_fission_sites(p, p->secondary_bank, &(p->n_secondary), + create_fission_sites(p, p->secondary_bank_, &(p->n_secondary_), MAX_SECONDARY); } } @@ -63,9 +63,9 @@ sample_reaction(Particle* p) if (simulation::material_xs.absorption > 0.) { absorption(p); } else { - p->absorb_wgt = 0.; + p->absorb_wgt_ = 0.; } - if (!p->alive) return; + if (!p->alive_) return; // Sample a scattering event to determine the energy of the exiting neutron scatter(p); @@ -73,7 +73,7 @@ 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; } } @@ -82,23 +82,23 @@ scatter(Particle* p) { // Adjust indices for Fortran to C++ indexing // TODO: Remove when no longer needed - int gin = p->last_g - 1; - int gout = p->g - 1; - int i_mat = p->material; - data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt); + int gin = p->last_g_ - 1; + int gout = p->g_ - 1; + int i_mat = p->material_; + data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu_, p->wgt_); // Adjust return value for fortran indexing // TODO: Remove when no longer needed - p->g = gout + 1; + p->g_ = gout + 1; // Rotate the angle - rotate_angle_c(p->coord[0].uvw, p->mu, nullptr); + rotate_angle_c(p->coord_[0].uvw, p->mu_, nullptr); // Update energy value for downstream compatability (in tallying) - p->E = data::energy_bin_avg[gout]; + p->E_ = data::energy_bin_avg[gout]; // Set event component - p->event = EVENT_SCATTER; + p->event_ = EVENT_SCATTER; } void @@ -112,7 +112,7 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, 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 * + double nu_t = p->wgt_ / simulation::keff * weight * simulation::material_xs.nu_fission / simulation::material_xs.total; // Sample the number of neutrons produced @@ -148,13 +148,13 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, // group. double nu_d[MAX_DELAYED_GROUPS] = {0.}; - p->fission = true; + p->fission_ = true; for (size_t i = static_cast(*size_bank); i < static_cast(std::min(*size_bank + nu, bank_array_size)); i++) { // Bank source neutrons by copying the particle data - bank_array[i].xyz[0] = p->coord[0].xyz[0]; - bank_array[i].xyz[1] = p->coord[0].xyz[1]; - bank_array[i].xyz[2] = p->coord[0].xyz[2]; + bank_array[i].xyz[0] = p->coord_[0].xyz[0]; + bank_array[i].xyz[1] = p->coord_[0].xyz[1]; + bank_array[i].xyz[2] = p->coord_[0].xyz[2]; // Set that the bank particle is a neutron bank_array[i].particle = static_cast(ParticleType::neutron); @@ -176,15 +176,15 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, // the energy in the fission bank int dg; int gout; - data::macro_xs[p->material].sample_fission_energy(p->g - 1, dg, gout); + data::macro_xs[p->material_].sample_fission_energy(p->g_ - 1, dg, gout); bank_array[i].E = static_cast(gout + 1); bank_array[i].delayed_group = dg + 1; // 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]++; } } @@ -193,10 +193,10 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank, *size_bank = std::min(*size_bank + nu, bank_array_size); // 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]; } } @@ -205,26 +205,26 @@ absorption(Particle* p) { if (settings::survival_biasing) { // Determine weight absorbed in survival biasing - p->absorb_wgt = p->wgt * + p->absorb_wgt_ = p->wgt_ * simulation::material_xs.absorption / simulation::material_xs.total; // Adjust weight of particle by the probability of absorption - p->wgt -= p->absorb_wgt; - p->last_wgt = p->wgt; + p->wgt_ -= p->absorb_wgt_; + p->last_wgt_ = p->wgt_; // Score implicit absorpion estimate of keff #pragma omp atomic - global_tally_absorption += p->absorb_wgt * + global_tally_absorption += p->absorb_wgt_ * simulation::material_xs.nu_fission / simulation::material_xs.absorption; } else { if (simulation::material_xs.absorption > prn() * simulation::material_xs.total) { #pragma omp atomic - global_tally_absorption += p->wgt * simulation::material_xs.nu_fission / + global_tally_absorption += p->wgt_ * simulation::material_xs.nu_fission / simulation::material_xs.absorption; - p->alive = false; - p->event = EVENT_ABSORB; + p->alive_ = false; + p->event_ = EVENT_ABSORB; } } diff --git a/src/plot.cpp b/src/plot.cpp index e5f5092ea..a1d24c938 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -137,18 +137,18 @@ void create_ppm(Plot pl) { Particle p; p.initialize(); - std::copy(xyz, xyz+3, p.coord[0].xyz); - std::copy(dir, dir+3, p.coord[0].uvw); - p.coord[0].universe = model::root_universe; + std::copy(xyz, xyz+3, p.coord_[0].xyz); + std::copy(dir, dir+3, p.coord_[0].uvw); + p.coord_[0].universe = model::root_universe; #pragma omp for for (int y = 0; y < height; y++) { - p.coord[0].xyz[out_i] = xyz[out_i] - out_pixel * y; + p.coord_[0].xyz[out_i] = xyz[out_i] - out_pixel * y; for (int x = 0; x < width; x++) { // local variables RGBColor rgb; int id; - p.coord[0].xyz[in_i] = xyz[in_i] + in_pixel * x; + p.coord_[0].xyz[in_i] = xyz[in_i] + in_pixel * x; position_rgb(p, pl, rgb, id); data(x,y) = rgb; } @@ -648,11 +648,11 @@ index_meshlines_mesh_(-1) void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id) { - p.n_coord = 1; + p.n_coord_ = 1; bool found_cell = find_cell(&p, 0); - int j = p.n_coord - 1; + int j = p.n_coord_ - 1; if (settings::check_overlaps) {check_cell_overlap(&p);} @@ -666,23 +666,23 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id) } else { if (PlotColorBy::mats == pl.color_by_) { // Assign color based on material - const auto& c = model::cells[p.coord[j].cell]; + const auto& c = model::cells[p.coord_[j].cell]; if (c->type_ == FILL_UNIVERSE) { // If we stopped on a middle universe level, treat as if not found rgb = pl.not_found_; id = -1; - } else if (p.material == MATERIAL_VOID) { + } else if (p.material_ == MATERIAL_VOID) { // By default, color void cells white rgb = WHITE; id = -1; } else { - rgb = pl.colors_[p.material]; - id = model::materials[p.material]->id_; + rgb = pl.colors_[p.material_]; + id = model::materials[p.material_]->id_; } } else if (PlotColorBy::cells == pl.color_by_) { // Assign color based on cell - rgb = pl.colors_[p.coord[j].cell]; - id = model::cells[p.coord[j].cell]->id_; + rgb = pl.colors_[p.coord_[j].cell]; + id = model::cells[p.coord_[j].cell]->id_; } } // endif found_cell } @@ -855,9 +855,9 @@ void create_voxel(Plot pl) double dir[3] = {0.5, 0.5, 0.5}; Particle p; p.initialize(); - std::copy(ll.begin(), ll.begin()+ll.size(), p.coord[0].xyz); - std::copy(dir, dir+3, p.coord[0].uvw); - p.coord[0].universe = model::root_universe; + std::copy(ll.begin(), ll.begin()+ll.size(), p.coord_[0].xyz); + std::copy(dir, dir+3, p.coord_[0].uvw); + p.coord_[0].universe = model::root_universe; // Open binary plot file for writing std::ofstream of; @@ -910,16 +910,16 @@ void create_voxel(Plot pl) // write to plot data data[y][x] = id; // advance particle in x direction - p.coord[0].xyz[0] = p.coord[0].xyz[0] + vox[0]; + p.coord_[0].xyz[0] = p.coord_[0].xyz[0] + vox[0]; } // advance particle in y direction - p.coord[0].xyz[1] = p.coord[0].xyz[1] + vox[1]; - p.coord[0].xyz[0] = ll[0]; + p.coord_[0].xyz[1] = p.coord_[0].xyz[1] + vox[1]; + p.coord_[0].xyz[0] = ll[0]; } // advance particle in z direction - p.coord[0].xyz[2] = p.coord[0].xyz[2] + vox[2]; - p.coord[0].xyz[1] = ll[1]; - p.coord[0].xyz[0] = ll[0]; + p.coord_[0].xyz[2] = p.coord_[0].xyz[2] + vox[2]; + p.coord_[0].xyz[1] = ll[1]; + p.coord_[0].xyz[0] = ll[0]; // Write to HDF5 dataset voxel_write_slice(z, dspace, dset, memspace, &(data[0])); } diff --git a/src/simulation.cpp b/src/simulation.cpp index 50ce1b5c0..b6c9a6b10 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -73,7 +73,7 @@ int openmc_simulation_init() t->init_results(); } - // Set up material nuclide index mapping + // Set up.material_ nuclide index mapping for (auto& mat : model::materials) { mat->init_nuclide_index(); } @@ -474,29 +474,29 @@ void initialize_history(Particle* p, int64_t index_source) p->from_source(&simulation::source_bank[index_source - 1]); // set identifier for particle - p->id = simulation::work_index[mpi::rank] + index_source; + p->id_ = simulation::work_index[mpi::rank] + index_source; // set random number seed int64_t particle_seed = (simulation::total_gen + overall_generation() - 1) - * settings::n_particles + p->id; + * settings::n_particles + p->id_; set_particle_seed(particle_seed); // set particle trace simulation::trace = false; if (simulation::current_batch == settings::trace_batch && simulation::current_gen == settings::trace_gen && - p->id == settings::trace_particle) simulation::trace = true; + p->id_ == settings::trace_particle) simulation::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; + p->id_ == t[2]) { + p->write_track_ = true; break; } } diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 2ad251b73..1ddd7af83 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -126,11 +126,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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; @@ -190,7 +190,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (tally.estimator_) { case ESTIMATOR_ANALOG: - if (p->event_nuclide != deriv.diff_nuclide) { + if (p->event_nuclide_ != deriv.diff_nuclide) { score *= flux_deriv; return; } @@ -320,10 +320,10 @@ 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->last_E)) { + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (!multipole_in_range(&nuc, p->last_E_)) { score *= flux_deriv; break; } @@ -331,10 +331,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, switch (score_bin) { case SCORE_TOTAL: - if (simulation::micro_xs[p->event_nuclide].total) { + if (simulation::micro_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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i) / simulation::material_xs.total; } else { @@ -343,11 +343,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_SCATTER: - if (simulation::micro_xs[p->event_nuclide].total - - simulation::micro_xs[p->event_nuclide].absorption) { + if (simulation::micro_xs[p->event_nuclide_].total + - simulation::micro_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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + dsig_s * material.atom_density_(i) / (simulation::material_xs.total - simulation::material_xs.absorption); @@ -357,10 +357,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_ABSORPTION: - if (simulation::micro_xs[p->event_nuclide].absorption) { + if (simulation::micro_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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + dsig_a * material.atom_density_(i) / simulation::material_xs.absorption; } else { @@ -369,10 +369,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_FISSION: - if (simulation::micro_xs[p->event_nuclide].fission) { + if (simulation::micro_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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + dsig_f * material.atom_density_(i) / simulation::material_xs.fission; } else { @@ -381,12 +381,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_NU_FISSION: - if (simulation::micro_xs[p->event_nuclide].fission) { - double nu = simulation::micro_xs[p->event_nuclide].nu_fission - / simulation::micro_xs[p->event_nuclide].fission; + if (simulation::micro_xs[p->event_nuclide_].fission) { + double nu = simulation::micro_xs[p->event_nuclide_].nu_fission + / simulation::micro_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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + nu * dsig_f * material.atom_density_(i) / simulation::material_xs.nu_fission; } else { @@ -404,7 +404,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, case ESTIMATOR_COLLISION: if (i_nuclide != -1) { const auto& nuc {data::nuclides[i_nuclide]}; - if (!multipole_in_range(nuc.get(), p->last_E)) { + if (!multipole_in_range(nuc.get(), p->last_E_)) { score *= flux_deriv; return; } @@ -418,11 +418,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E) + if (multipole_in_range(&nuc, p->last_E_) && simulation::micro_xs[i_nuc].total) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i); } } @@ -431,7 +431,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + (dsig_s + dsig_a) / simulation::micro_xs[i_nuclide].total; } else { @@ -446,12 +446,12 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E) + if (multipole_in_range(&nuc, p->last_E_) && (simulation::micro_xs[i_nuc].total - simulation::micro_xs[i_nuc].absorption)) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); cum_dsig += dsig_s * material.atom_density_(i); } } @@ -462,7 +462,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total - simulation::micro_xs[i_nuclide].absorption); } else { @@ -476,11 +476,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E) + if (multipole_in_range(&nuc, p->last_E_) && simulation::micro_xs[i_nuc].absorption) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); cum_dsig += dsig_a * material.atom_density_(i); } } @@ -489,7 +489,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + dsig_a / simulation::micro_xs[i_nuclide].absorption; } else { @@ -503,11 +503,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E) + if (multipole_in_range(&nuc, p->last_E_) && simulation::micro_xs[i_nuc].fission) { double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); cum_dsig += dsig_f * material.atom_density_(i); } } @@ -516,7 +516,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + dsig_f / simulation::micro_xs[i_nuclide].fission; } else { @@ -530,13 +530,13 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E) + if (multipole_in_range(&nuc, p->last_E_) && simulation::micro_xs[i_nuc].fission) { double nu = simulation::micro_xs[i_nuc].nu_fission / simulation::micro_xs[i_nuc].fission; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); cum_dsig += nu * dsig_f * material.atom_density_(i); } } @@ -545,7 +545,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, 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->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); score *= flux_deriv + dsig_f / simulation::micro_xs[i_nuclide].fission; } else { @@ -570,8 +570,8 @@ void score_track_derivative(const Particle* p, double distance) { // A void material cannot be perturbed so it will not affect flux derivatives. - if (p->material == MATERIAL_VOID) return; - const Material& material {*model::materials[p->material]}; + if (p->material_ == MATERIAL_VOID) return; + const Material& material {*model::materials[p->material_]}; for (auto& deriv : model::tally_derivs) { if (deriv.diff_material != material.id_) continue; @@ -597,13 +597,13 @@ score_track_derivative(const Particle* p, double distance) case DIFF_TEMPERATURE: for (auto i = 0; i < material.nuclide_.size(); ++i) { const auto& nuc {*data::nuclides[material.nuclide_[i]]}; - if (multipole_in_range(&nuc, p->last_E)) { + if (multipole_in_range(&nuc, p->last_E_)) { // 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); + = nuc.multipole_->evaluate_deriv(p->E_, p->sqrtkT_); deriv.flux_deriv -= distance * (dsig_s + dsig_a) * material.atom_density_(i); } @@ -616,8 +616,8 @@ score_track_derivative(const Particle* p, double distance) void score_collision_derivative(const Particle* p) { // A void material cannot be perturbed so it will not affect flux derivatives. - if (p->material == MATERIAL_VOID) return; - const Material& material {*model::materials[p->material]}; + if (p->material_ == MATERIAL_VOID) return; + const Material& material {*model::materials[p->material_]}; for (auto& deriv : model::tally_derivs) { if (deriv.diff_material != material.id_) continue; @@ -632,7 +632,7 @@ void score_collision_derivative(const Particle* p) break; case DIFF_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) @@ -656,14 +656,14 @@ void score_collision_derivative(const 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->last_E)) { + if (i_nuc == p->event_nuclide_ && multipole_in_range(&nuc, p->last_E_)) { // 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 {simulation::micro_xs[i_nuc]}; double dsig_s, dsig_a, dsig_f; std::tie(dsig_s, dsig_a, dsig_f) - = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT); + = nuc.multipole_->evaluate_deriv(p->last_E_, p->sqrtkT_); deriv.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 f4e265dd8..4b78cabd5 100644 --- a/src/tallies/filter_azimuthal.cpp +++ b/src/tallies/filter_azimuthal.cpp @@ -42,9 +42,9 @@ AzimuthalFilter::get_all_bins(const Particle* p, int estimator, { double phi; if (estimator == ESTIMATOR_TRACKLENGTH) { - phi = std::atan2(p->coord[0].uvw[1], p->coord[0].uvw[0]); + phi = std::atan2(p->coord_[0].uvw[1], p->coord_[0].uvw[0]); } else { - phi = std::atan2(p->last_uvw[1], p->last_uvw[0]); + phi = std::atan2(p->last_uvw_[1], p->last_uvw_[0]); } if (phi >= bins_.front() && phi <= bins_.back()) { diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp index 1b2b05416..8286a6976 100644 --- a/src/tallies/filter_cell.cpp +++ b/src/tallies/filter_cell.cpp @@ -41,8 +41,8 @@ void CellFilter::get_all_bins(const Particle* p, int 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_cellborn.cpp b/src/tallies/filter_cellborn.cpp index 005d3844b..150962297 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, int 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 9a43f7c43..726bd7453 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, int estimator, FilterMatch& match) const { - for (int i = 0; i < p->last_n_coord; i++) { - auto search = map_.find(p->last_cell[i]); + for (int i = 0; i < p->last_n_coord_; i++) { + auto search = map_.find(p->last_cell_[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 ab7ea67cf..ce176cea3 100644 --- a/src/tallies/filter_distribcell.cpp +++ b/src/tallies/filter_distribcell.cpp @@ -40,20 +40,20 @@ DistribcellFilter::get_all_bins(const Particle* p, int 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) { - auto& lat {*model::lattices[p->coord[i+1].lattice-1]}; - int i_xyz[3] {p->coord[i+1].lattice_x, - p->coord[i+1].lattice_y, - p->coord[i+1].lattice_z}; + auto& lat {*model::lattices[p->coord_[i+1].lattice-1]}; + int i_xyz[3] {p->coord_[i+1].lattice_x, + p->coord_[i+1].lattice_y, + p->coord_[i+1].lattice_z}; if (lat.are_valid_indices(i_xyz)) { 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 4cbde1cc4..dd9402e83 100644 --- a/src/tallies/filter_energy.cpp +++ b/src/tallies/filter_energy.cpp @@ -41,17 +41,17 @@ void EnergyFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - if (p->g != F90_NONE && matches_transport_groups_) { + if (p->g_ != F90_NONE && matches_transport_groups_) { if (estimator == ESTIMATOR_TRACKLENGTH) { - match.bins_.push_back(data::num_energy_groups - p->g); + match.bins_.push_back(data::num_energy_groups - p->g_); } else { - match.bins_.push_back(data::num_energy_groups - p->last_g); + match.bins_.push_back(data::num_energy_groups - p->last_g_); } match.weights_.push_back(1.0); } else { // Get the pre-collision energy of the particle. - auto E = p->last_E; + auto E = p->last_E_; // Bin the energy. if (E >= bins_.front() && E <= bins_.back()) { @@ -85,13 +85,13 @@ void EnergyoutFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - if (p->g != F90_NONE && matches_transport_groups_) { - match.bins_.push_back(data::num_energy_groups - p->g); + if (p->g_ != F90_NONE && matches_transport_groups_) { + match.bins_.push_back(data::num_energy_groups - p->g_); 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 3c616aabf..fc51a4a26 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -33,12 +33,12 @@ void EnergyFunctionFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - if (p->last_E >= energy_.front() && p->last_E <= energy_.back()) { + if (p->last_E_ >= energy_.front() && p->last_E_ <= energy_.back()) { // Search for the incoming energy bin. - auto i = lower_bound_index(energy_.begin(), energy_.end(), p->last_E); + auto i = lower_bound_index(energy_.begin(), energy_.end(), p->last_E_); // Compute the interpolation factor between the nearest bins. - double f = (p->last_E - energy_[i]) / (energy_[i+1] - energy_[i]); + double f = (p->last_E_ - 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 1e799932a..1d1d440b2 100644 --- a/src/tallies/filter_legendre.cpp +++ b/src/tallies/filter_legendre.cpp @@ -19,7 +19,7 @@ LegendreFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { double wgt[n_bins_]; - calc_pn_c(order_, p->mu, wgt); + calc_pn_c(order_, p->mu_, wgt); 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 a4da2deb1..fd9e637bd 100644 --- a/src/tallies/filter_material.cpp +++ b/src/tallies/filter_material.cpp @@ -42,7 +42,7 @@ void MaterialFilter::get_all_bins(const Particle* p, int 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_mesh.cpp b/src/tallies/filter_mesh.cpp index 177dc407e..518470bf4 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -35,7 +35,7 @@ MeshFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { if (estimator != ESTIMATOR_TRACKLENGTH) { - auto bin = model::meshes[mesh_]->get_bin(p->coord[0].xyz); + auto bin = model::meshes[mesh_]->get_bin(p->coord_[0].xyz); if (bin >= 0) { match.bins_.push_back(bin); match.weights_.push_back(1.0); diff --git a/src/tallies/filter_mu.cpp b/src/tallies/filter_mu.cpp index 78a648e23..aed8371b8 100644 --- a/src/tallies/filter_mu.cpp +++ b/src/tallies/filter_mu.cpp @@ -38,8 +38,8 @@ void MuFilter::get_all_bins(const Particle* p, int 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 967febc06..98393df0e 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -17,7 +17,7 @@ ParticleFilter::get_all_bins(const Particle* p, int 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 5216511a9..20eb7608b 100644 --- a/src/tallies/filter_polar.cpp +++ b/src/tallies/filter_polar.cpp @@ -41,9 +41,9 @@ const { double theta; if (estimator == ESTIMATOR_TRACKLENGTH) { - theta = std::acos(p->coord[0].uvw[2]); + theta = std::acos(p->coord_[0].uvw[2]); } else { - theta = std::acos(p->last_uvw[2]); + theta = std::acos(p->last_uvw_[2]); } if (theta >= bins_.front() && theta <= bins_.back()) { diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp index 253bf7ff3..3ef2df43c 100644 --- a/src/tallies/filter_sph_harm.cpp +++ b/src/tallies/filter_sph_harm.cpp @@ -37,14 +37,14 @@ SphericalHarmonicsFilter::get_all_bins(const Particle* p, int estimator, // Determine cosine term for scatter expansion if necessary double wgt[order_ + 1]; if (cosine_ == SphericalHarmonicsCosine::scatter) { - calc_pn_c(order_, p->mu, wgt); + calc_pn_c(order_, p->mu_, wgt); } else { for (int i = 0; i < order_ + 1; i++) wgt[i] = 1; } // Find the Rn,m values double rn[n_bins_]; - calc_rn_c(order_, p->last_uvw, rn); + calc_rn_c(order_, p->last_uvw_, rn); int j = 0; for (int n = 0; n < order_ + 1; n++) { diff --git a/src/tallies/filter_sptl_legendre.cpp b/src/tallies/filter_sptl_legendre.cpp index dae62d1d2..14b8fe14e 100644 --- a/src/tallies/filter_sptl_legendre.cpp +++ b/src/tallies/filter_sptl_legendre.cpp @@ -38,11 +38,11 @@ SpatialLegendreFilter::get_all_bins(const Particle* p, int estimator, // Get the coordinate along the axis of interest. double x; if (axis_ == LegendreAxis::x) { - x = p->coord[0].xyz[0]; + x = p->coord_[0].xyz[0]; } else if (axis_ == LegendreAxis::y) { - x = p->coord[0].xyz[1]; + x = p->coord_[0].xyz[1]; } else { - x = p->coord[0].xyz[2]; + x = p->coord_[0].xyz[2]; } if (x >= min_ && x <= max_) { diff --git a/src/tallies/filter_surface.cpp b/src/tallies/filter_surface.cpp index 774571d75..faa8aafbf 100644 --- a/src/tallies/filter_surface.cpp +++ b/src/tallies/filter_surface.cpp @@ -41,10 +41,10 @@ void SurfaceFilter::get_all_bins(const Particle* p, int 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 5a6d5e1eb..50c058c05 100644 --- a/src/tallies/filter_universe.cpp +++ b/src/tallies/filter_universe.cpp @@ -41,8 +41,8 @@ void UniverseFilter::get_all_bins(const Particle* p, int 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/filter_zernike.cpp b/src/tallies/filter_zernike.cpp index a005f7ddc..93ff16109 100644 --- a/src/tallies/filter_zernike.cpp +++ b/src/tallies/filter_zernike.cpp @@ -29,8 +29,8 @@ ZernikeFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { // Determine the normalized (r,theta) coordinates. - double x = p->coord[0].xyz[0] - x_; - double y = p->coord[0].xyz[1] - y_; + double x = p->coord_[0].xyz[0] - x_; + double y = p->coord_[0].xyz[1] - y_; double r = std::sqrt(x*x + y*y) / r_; double theta = std::atan2(y, x); @@ -86,8 +86,8 @@ ZernikeRadialFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { // Determine the normalized radius coordinate. - double x = p->coord[0].xyz[0] - x_; - double y = p->coord[0].xyz[1] - y_; + double x = p->coord_[0].xyz[0] - x_; + double y = p->coord_[0].xyz[1] - y_; double r = std::sqrt(x*x + y*y) / r_; if (r <= 1.0) { diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 5288ea597..c5c02eb7d 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -190,8 +190,8 @@ score_fission_eout(const 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) { - auto i_bank = simulation::n_bank - p->n_bank + i; + for (auto i = 0; i < p->n_bank_; ++i) { + auto i_bank = simulation::n_bank - p->n_bank_ + i; const auto& bank = simulation::fission_bank[i_bank]; // get the delayed group @@ -322,7 +322,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto& tally {*model::tallies[i_tally]}; // Get the pre-collision energy of the particle. - auto E = p->last_E; + auto E = p->last_E_; for (auto i = 0; i < tally.scores_.size(); ++i) { auto score_bin = tally.scores_[i]; @@ -341,13 +341,13 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->last_wgt_ + p->absorb_wgt_; } else { - score = p->last_wgt; + score = p->last_wgt_; } - if (p->type == static_cast(ParticleType::neutron) || - p->type == static_cast(ParticleType::photon)) { + if (p->type_ == static_cast(ParticleType::neutron) || + p->type_ == static_cast(ParticleType::photon)) { score *= flux / simulation::material_xs.total; } else { score = 0.; @@ -365,9 +365,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = (p->last_wgt + p->absorb_wgt) * flux; + score = (p->last_wgt_ + p->absorb_wgt_) * flux; } else { - score = p->last_wgt * flux; + score = p->last_wgt_ * flux; } } else { @@ -388,9 +388,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->last_wgt_ + p->absorb_wgt_; } else { - score = p->last_wgt; + score = p->last_wgt_; } score *= flux / simulation::material_xs.total; } else { @@ -404,10 +404,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_SCATTER: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Skip any event where the particle didn't scatter - if (p->event != EVENT_SCATTER) continue; + if (p->event_ != EVENT_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->last_wgt * flux; + score = p->last_wgt_ * flux; } else { if (i_nuclide >= 0) { score = (simulation::micro_xs[i_nuclide].total @@ -423,21 +423,21 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_NU_SCATTER: // Only analog estimators are available. // Skip any event where the particle didn't scatter - if (p->event != EVENT_SCATTER) continue; + if (p->event_ != EVENT_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->last_wgt * flux; + score = p->last_wgt_ * 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->last_wgt * flux * (*rxn.products_[0].yield_)(E); + data::nuclides[p->event_nuclide_]->reaction_index_[p->event_mt_]; + const auto& rxn {*data::nuclides[p->event_nuclide_]->reactions_[m]}; + score = p->last_wgt_ * flux * (*rxn.products_[0].yield_)(E); } break; @@ -447,13 +447,13 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p->absorb_wgt * flux; + score = p->absorb_wgt_ * flux; } else { // Skip any event where the particle wasn't absorbed - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_SCATTER) continue; // All fission and absorption events will contribute here, so we // can just use the particle's weight entering the collision - score = p->last_wgt * flux; + score = p->last_wgt_ * flux; } } else { if (i_nuclide >= 0) { @@ -473,22 +473,22 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - if (simulation::micro_xs[p->event_nuclide].absorption > 0) { - score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + score = p->absorb_wgt_ + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } else { score = 0.; } } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_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->last_wgt - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->last_wgt_ + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } else { if (i_nuclide >= 0) { @@ -503,7 +503,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_NU_FISSION: if (simulation::material_xs.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_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. @@ -515,22 +515,22 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - if (simulation::micro_xs[p->event_nuclide].absorption > 0) { - score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide].nu_fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + score = p->absorb_wgt_ + * simulation::micro_xs[p->event_nuclide_].nu_fission + / simulation::micro_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) { @@ -546,7 +546,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_PROMPT_NU_FISSION: if (simulation::material_xs.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_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. @@ -558,27 +558,27 @@ score_general_ce(const Particle* p, int i_tally, int start_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 (simulation::micro_xs[p->event_nuclide].absorption > 0) { - score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide].fission - * data::nuclides[p->event_nuclide] + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { + score = p->absorb_wgt_ + * simulation::micro_xs[p->event_nuclide_].fission + * data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::prompt) - / simulation::micro_xs[p->event_nuclide].absorption * flux; + / simulation::micro_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) { @@ -589,8 +589,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } 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); @@ -608,7 +608,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case SCORE_DELAYED_NU_FISSION: if (simulation::material_xs.absorption == 0) continue; if (tally.estimator_ == ESTIMATOR_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. @@ -620,8 +620,8 @@ score_general_ce(const Particle* p, int i_tally, int start_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 (simulation::micro_xs[p->event_nuclide].absorption > 0 - && data::nuclides[p->event_nuclide]->fissionable_) { + if (simulation::micro_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 @@ -630,11 +630,11 @@ score_general_ce(const Particle* p, int i_tally, int start_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]; - auto yield = data::nuclides[p->event_nuclide] + auto yield = data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = p->absorb_wgt * yield - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->absorb_wgt_ * yield + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); } @@ -643,16 +643,16 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // 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->absorb_wgt - * simulation::micro_xs[p->event_nuclide].fission - * data::nuclides[p->event_nuclide] + score = p->absorb_wgt_ + * simulation::micro_xs[p->event_nuclide_].fission + * data::nuclides[p->event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed) - / simulation::micro_xs[p->event_nuclide].absorption *flux; + / simulation::micro_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 @@ -668,16 +668,16 @@ score_general_ce(const Particle* p, int i_tally, int start_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; score_fission_delayed_dg(i_tally, d_bin, score, score_index); } 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 + 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; } } @@ -712,8 +712,8 @@ score_general_ce(const Particle* p, int i_tally, int start_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); @@ -732,8 +732,8 @@ score_general_ce(const Particle* p, int i_tally, int start_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); @@ -756,8 +756,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // 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 (simulation::micro_xs[p->event_nuclide].absorption > 0 + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; if (tally.delayedgroup_filter_ != C_NONE) { @@ -771,9 +771,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p->absorb_wgt * yield - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption + score = p->absorb_wgt_ * yield + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * rate * flux; score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -794,15 +794,15 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += rate * p->absorb_wgt - * simulation::micro_xs[p->event_nuclide].fission * yield - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score += rate * p->absorb_wgt_ + * simulation::micro_xs[p->event_nuclide_].fission * yield + / simulation::micro_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 @@ -811,12 +811,12 @@ score_general_ce(const Particle* p, int i_tally, int start_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) { - auto i_bank = simulation::n_bank - p->n_bank + i; + for (auto i = 0; i < p->n_bank_; ++i) { + auto i_bank = simulation::n_bank - p->n_bank_ + i; const auto& bank = simulation::fission_bank[i_bank]; 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; @@ -879,8 +879,8 @@ score_general_ce(const Particle* p, int i_tally, int start_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); @@ -904,8 +904,8 @@ score_general_ce(const Particle* p, int i_tally, int start_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); @@ -943,27 +943,27 @@ score_general_ce(const Particle* p, int i_tally, int start_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 (simulation::micro_xs[p->event_nuclide].absorption > 0 + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = p->absorb_wgt * rxn.q_value_ - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->absorb_wgt_ * rxn.q_value_ + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_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 (simulation::micro_xs[p->event_nuclide].absorption > 0 + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = p->last_wgt * rxn.q_value_ - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->last_wgt_ * rxn.q_value_ + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } } else { @@ -975,8 +975,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, * atom_density * flux; } } else { - 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); @@ -1002,8 +1002,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case ELASTIC: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Check if event MT matches - if (p->event_MT != ELASTIC) continue; - score = p->last_wgt * flux; + if (p->event_mt_ != ELASTIC) continue; + score = p->last_wgt_ * flux; } else { if (i_nuclide >= 0) { if (simulation::micro_xs[i_nuclide].elastic == CACHE_INVALID) @@ -1011,8 +1011,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, score = simulation::micro_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); @@ -1036,39 +1036,39 @@ score_general_ce(const Particle* p, int i_tally, int start_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 (simulation::micro_xs[p->event_nuclide].absorption > 0) { + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->last_E); + q_value = (*nuc.fission_q_prompt_)(p->last_E_); } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->last_E); + q_value = (*nuc.fission_q_recov_)(p->last_E_); } - score = p->absorb_wgt * q_value - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->absorb_wgt_ * q_value + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_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 (simulation::micro_xs[p->event_nuclide].absorption > 0) { + const auto& nuc {*data::nuclides[p->event_nuclide_]}; + if (simulation::micro_xs[p->event_nuclide_].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->last_E); + q_value = (*nuc.fission_q_prompt_)(p->last_E_); } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->last_E); + q_value = (*nuc.fission_q_recov_)(p->last_E_); } - score = p->last_wgt * q_value - * simulation::micro_xs[p->event_nuclide].fission - / simulation::micro_xs[p->event_nuclide].absorption * flux; + score = p->last_wgt_ * q_value + * simulation::micro_xs[p->event_nuclide_].fission + / simulation::micro_xs[p->event_nuclide_].absorption * flux; } } } else { @@ -1077,16 +1077,16 @@ score_general_ce(const Particle* p, int i_tally, int start_index, double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->last_E); + q_value = (*nuc.fission_q_prompt_)(p->last_E_); } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->last_E); + q_value = (*nuc.fission_q_recov_)(p->last_E_); } score = q_value * simulation::micro_xs[i_nuclide].fission * atom_density * flux; } else { - 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); @@ -1094,10 +1094,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) - q_value = (*nuc.fission_q_prompt_)(p->last_E); + q_value = (*nuc.fission_q_prompt_)(p->last_E_); } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) - q_value = (*nuc.fission_q_recov_)(p->last_E); + q_value = (*nuc.fission_q_recov_)(p->last_E_); } score += q_value * simulation::micro_xs[j_nuclide].fission * atom_density * flux; @@ -1116,8 +1116,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, case N_A: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Check if the event MT matches - if (p->event_MT != score_bin) continue; - score = p->last_wgt * flux; + if (p->event_mt_ != score_bin) continue; + score = p->last_wgt_ * flux; } else { int m; switch (score_bin) { @@ -1133,8 +1133,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, * 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); @@ -1151,8 +1151,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, if (tally.estimator_ == ESTIMATOR_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->last_wgt*flux; + if (p->event_mt_ != score_bin) continue; + score = p->last_wgt_*flux; } else { // Any other cross section has to be calculated on-the-fly if (score_bin < 2) fatal_error("Invalid score type on tally " @@ -1174,8 +1174,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } } } else { - 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); @@ -1232,40 +1232,40 @@ score_general_mg(const Particle* p, int i_tally, int start_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_uvw = p->last_uvw; - p_g = p->last_g; + if (p->alive_) { + p_uvw = p->last_uvw_; + p_g = p->last_g_; } else { - p_uvw = p->coord[p->n_coord-1].uvw; - p_g = p->g; + p_uvw = p->coord_[p->n_coord_-1].uvw; + p_g = p->g_; } - } else if (p->event == EVENT_SCATTER) { + } else if (p->event_ == EVENT_SCATTER) { // Then the energy group has been changed by the scattering routine // meaning gin is now in p % last_g - p_uvw = p->last_uvw; - p_g = p->last_g; + p_uvw = p->last_uvw_; + p_g = p->last_g_; } else { // No scatter, no change in g. - p_uvw = p->coord[p->n_coord-1].uvw; - p_g = p->g; + p_uvw = p->coord_[p->n_coord_-1].uvw; + p_g = p->g_; } } else { // No actual collision so g has not changed. - p_uvw = p->coord[p->n_coord-1].uvw; - p_g = p->g; + p_uvw = p->coord_[p->n_coord_-1].uvw; + p_g = p->g_; } // To significantly reduce de-referencing, point matxs to the macroscopic // Mgxs for the material of interest - data::macro_xs[p->material].set_angle_index(p_uvw); + data::macro_xs[p->material_].set_angle_index(p_uvw); // Do same for nucxs, point it to the microscopic nuclide data of interest if (i_nuclide >= 0) { // And since we haven't calculated this temperature index yet, do so now - data::nuclides_MG[i_nuclide].set_temperature_index(p->sqrtkT); + data::nuclides_MG[i_nuclide].set_temperature_index(p->sqrtkT_); data::nuclides_MG[i_nuclide].set_angle_index(p_uvw); } @@ -1286,9 +1286,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->last_wgt_ + p->absorb_wgt_; } else { - score = p->last_wgt; + score = p->last_wgt_; } score *= flux / simulation::material_xs.total; } else { @@ -1304,15 +1304,15 @@ score_general_mg(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->last_wgt_ + p->absorb_wgt_; } else { - score = p->last_wgt; + score = p->last_wgt_; } //TODO: should flux be multiplied in above instead of below? if (i_nuclide >= 0) { score *= flux * atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_TOTAL, p_g) - / get_macro_xs(p->material, MG_GET_XS_TOTAL, p_g); + / get_macro_xs(p->material_, MG_GET_XS_TOTAL, p_g); } } else { if (i_nuclide >= 0) { @@ -1334,18 +1334,18 @@ score_general_mg(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->last_wgt + p->absorb_wgt; + score = p->last_wgt_ + p->absorb_wgt_; } else { - score = p->last_wgt; + score = p->last_wgt_; } if (i_nuclide >= 0) { score *= flux * get_nuclide_xs(i_nuclide, MG_GET_XS_INVERSE_VELOCITY, p_g) - / get_macro_xs(p->material, MG_GET_XS_TOTAL, p_g); + / get_macro_xs(p->material_, MG_GET_XS_TOTAL, p_g); } else { score *= flux - * get_macro_xs(p->material, MG_GET_XS_INVERSE_VELOCITY, p_g) - / get_macro_xs(p->material, MG_GET_XS_TOTAL, p_g); + * get_macro_xs(p->material_, MG_GET_XS_INVERSE_VELOCITY, p_g) + / get_macro_xs(p->material_, MG_GET_XS_TOTAL, p_g); } } else { if (i_nuclide >= 0) { @@ -1353,7 +1353,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * get_nuclide_xs(i_nuclide, MG_GET_XS_INVERSE_VELOCITY, p_g); } else { score = flux - * get_macro_xs(p->material, MG_GET_XS_INVERSE_VELOCITY, p_g); + * get_macro_xs(p->material_, MG_GET_XS_INVERSE_VELOCITY, p_g); } } break; @@ -1362,24 +1362,24 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_SCATTER: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Skip any event where the particle didn't scatter - if (p->event != EVENT_SCATTER) continue; + if (p->event_ != EVENT_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->last_wgt * flux; + score = p->last_wgt_ * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_SCATTER_FMU_MULT, - p->last_g, &p->g, &p->mu, nullptr) - / get_macro_xs(p->material, MG_GET_XS_SCATTER_FMU_MULT, - p->last_g, &p->g, &p->mu, nullptr); + p->last_g_, &p->g_, &p->mu_, nullptr) + / get_macro_xs(p->material_, MG_GET_XS_SCATTER_FMU_MULT, + p->last_g_, &p->g_, &p->mu_, nullptr); } } else { if (i_nuclide >= 0) { score = atom_density * flux * get_nuclide_xs( - i_nuclide, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu, nullptr); + i_nuclide, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr); } else { score = flux * get_macro_xs( - p->material, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu, nullptr); + p->material_, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr); } } break; @@ -1388,27 +1388,27 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_NU_SCATTER: if (tally.estimator_ == ESTIMATOR_ANALOG) { // Skip any event where the particle didn't scatter - if (p->event != EVENT_SCATTER) continue; + if (p->event_ != EVENT_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 * get_nuclide_xs(i_nuclide, MG_GET_XS_SCATTER_FMU, - p->last_g, &p->g, &p->mu, nullptr) - / get_macro_xs(p->material, MG_GET_XS_SCATTER_FMU, - p->last_g, &p->g, &p->mu, nullptr); + p->last_g_, &p->g_, &p->mu_, nullptr) + / get_macro_xs(p->material_, MG_GET_XS_SCATTER_FMU, + p->last_g_, &p->g_, &p->mu_, nullptr); } } else { if (i_nuclide >= 0) { score = atom_density * flux * get_nuclide_xs( i_nuclide, MG_GET_XS_SCATTER, p_g); } else { - score = flux * get_macro_xs(p->material, MG_GET_XS_SCATTER, p_g); + score = flux * get_macro_xs(p->material_, MG_GET_XS_SCATTER, p_g); } } break; @@ -1419,18 +1419,18 @@ score_general_mg(const Particle* p, int i_tally, int start_index, if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p->absorb_wgt * flux; + score = p->absorb_wgt_ * flux; } else { // Skip any event where the particle wasn't absorbed - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_SCATTER) continue; // All fission and absorption events will contribute here, so we // can just use the particle's weight entering the collision - score = p->last_wgt * flux; + score = p->last_wgt_ * flux; } if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_ABSORPTION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } else { if (i_nuclide >= 0) { @@ -1449,30 +1449,30 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - score = p->absorb_wgt * flux; + score = p->absorb_wgt_ * flux; } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_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->last_wgt * flux; + score = p->last_wgt_ * flux; } if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } else { if (i_nuclide >= 0) { score = get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) * atom_density * flux; } else { - score = get_macro_xs(p->material, MG_GET_XS_FISSION, p_g) * flux; + score = get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g) * flux; } } break; @@ -1480,7 +1480,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_NU_FISSION: if (tally.estimator_ == ESTIMATOR_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. @@ -1492,29 +1492,29 @@ score_general_mg(const Particle* p, int i_tally, int start_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->absorb_wgt * flux; + score = p->absorb_wgt_ * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_NU_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } 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 * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } } } else { @@ -1522,7 +1522,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = get_nuclide_xs(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * atom_density * flux; } else { - score = get_macro_xs(p->material, MG_GET_XS_NU_FISSION, p_g) * flux; + score = get_macro_xs(p->material_, MG_GET_XS_NU_FISSION, p_g) * flux; } } break; @@ -1530,7 +1530,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_PROMPT_NU_FISSION: if (tally.estimator_ == ESTIMATOR_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. @@ -1542,32 +1542,32 @@ score_general_mg(const Particle* p, int i_tally, int start_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->absorb_wgt * flux; + score = p->absorb_wgt_ * flux; if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_PROMPT_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_PROMPT_NU_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } 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 * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } } } else { @@ -1575,7 +1575,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = get_nuclide_xs(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * atom_density * flux; } else { - score = get_macro_xs(p->material, MG_GET_XS_PROMPT_NU_FISSION, p_g) + score = get_macro_xs(p->material_, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux; } } @@ -1584,7 +1584,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, case SCORE_DELAYED_NU_FISSION: if (tally.estimator_ == ESTIMATOR_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. @@ -1596,7 +1596,7 @@ score_general_mg(const Particle* p, int i_tally, int start_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 (get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g) > 0) { + if (get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g) > 0) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; const DelayedGroupFilter& filt @@ -1605,17 +1605,17 @@ score_general_mg(const Particle* p, int i_tally, int start_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 = p->absorb_wgt * flux; + score = p->absorb_wgt_ * flux; if (i_nuclide >= 0) { score *= get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1625,21 +1625,21 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // 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->absorb_wgt * flux; + score = p->absorb_wgt_ * flux; if (i_nuclide >= 0) { score *= get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } } } 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 @@ -1655,26 +1655,26 @@ score_general_mg(const Particle* p, int i_tally, int start_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 * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); } 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 + 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 * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } } } @@ -1693,7 +1693,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, p_g, nullptr, nullptr, &d); } else { score = flux - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1705,7 +1705,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g); } else { score = flux - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, p_g); + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g); } } } @@ -1718,7 +1718,7 @@ score_general_mg(const Particle* p, int i_tally, int start_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 (get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g) > 0) { + if (get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g) > 0) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; const DelayedGroupFilter& filt @@ -1727,21 +1727,21 @@ score_general_mg(const Particle* p, int i_tally, int start_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 = p->absorb_wgt * flux; + score = p->absorb_wgt_ * flux; if (i_nuclide >= 0) { score *= get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, + get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1755,26 +1755,26 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = 0.; for (auto d = 0; d < data::num_delayed_groups; ++d) { if (i_nuclide >= 0) { - score += p->absorb_wgt * flux + score += p->absorb_wgt_ * flux * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { - score += p->absorb_wgt * flux - * get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, + score += p->absorb_wgt_ * flux + * get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } } } } 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 @@ -1783,8 +1783,8 @@ score_general_mg(const Particle* p, int i_tally, int start_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) { - auto i_bank = simulation::n_bank - p->n_bank + i; + for (auto i = 0; i < p->n_bank_; ++i) { + auto i_bank = simulation::n_bank - p->n_bank_ + i; const auto& bank = simulation::fission_bank[i_bank]; auto g = bank.delayed_group; if (g != 0) { @@ -1793,10 +1793,10 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &g) * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_FISSION, p_g); } else { score += simulation::keff * bank.wgt * flux - * get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, p_g, + * get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &g); } if (tally.delayedgroup_filter_ != C_NONE) { @@ -1834,9 +1834,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, p_g, nullptr, nullptr, &d); } else { score += flux - * get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, + * get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1853,9 +1853,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, p_g, nullptr, nullptr, &d); } else { score += flux - * get_macro_xs(p->material, MG_GET_XS_DECAY_RATE, + * get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, + * get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } } @@ -1870,30 +1870,30 @@ score_general_mg(const Particle* p, int i_tally, int start_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->absorb_wgt * flux; + score = p->absorb_wgt_ * flux; } else { // Skip any non-absorption events - if (p->event == EVENT_SCATTER) continue; + if (p->event_ == EVENT_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->last_wgt * flux; + score = p->last_wgt_ * flux; } if (i_nuclide >= 0) { score *= atom_density * get_nuclide_xs(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } else { score *= - get_macro_xs(p->material, MG_GET_XS_KAPPA_FISSION, p_g) - / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); + get_macro_xs(p->material_, MG_GET_XS_KAPPA_FISSION, p_g) + / get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g); } } else { if (i_nuclide >= 0) { score = get_nuclide_xs(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * atom_density * flux; } else { - score = get_macro_xs(p->material, MG_GET_XS_KAPPA_FISSION, p_g) + score = get_macro_xs(p->material_, MG_GET_XS_KAPPA_FISSION, p_g) * flux; } } @@ -1923,7 +1923,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux, int filter_index) { 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) { @@ -1980,7 +1980,7 @@ void score_analog_tally_ce(const Particle* p) // the event nuclide or the total material. Note that the i_nuclide // and flux arguments for score_general are 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, -1, -1., filter_weight); } @@ -1989,7 +1989,7 @@ void score_analog_tally_ce(const 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, -1, -1., filter_weight); @@ -2035,9 +2035,9 @@ void score_analog_tally_mg(const 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, @@ -2061,7 +2061,7 @@ void score_tracklength_tally(const 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]}; @@ -2080,7 +2080,7 @@ score_tracklength_tally(const 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); } else { @@ -2089,10 +2089,10 @@ score_tracklength_tally(const 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); } } @@ -2126,9 +2126,9 @@ void score_collision_tally(const Particle* p) // Determine the collision estimate of the flux double flux; if (!settings::survival_biasing) { - flux = p->last_wgt / simulation::material_xs.total; + flux = p->last_wgt_ / simulation::material_xs.total; } else { - flux = (p->last_wgt + p->absorb_wgt) / simulation::material_xs.total; + flux = (p->last_wgt_ + p->absorb_wgt_) / simulation::material_xs.total; } for (auto i_tally : model::active_collision_tallies) { @@ -2156,9 +2156,9 @@ void score_collision_tally(const 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 @@ -2189,7 +2189,7 @@ void score_surface_tally(const 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]}; diff --git a/src/track_output.cpp b/src/track_output.cpp index df5cef2fc..34cf7fdf6 100644 --- a/src/track_output.cpp +++ b/src/track_output.cpp @@ -40,14 +40,14 @@ void add_particle_track() void write_particle_track(const Particle& p) { - tracks.back().push_back({p.coord[0].xyz}); + tracks.back().push_back({p.coord_[0].xyz}); } void finalize_particle_track(const Particle& p) { std::stringstream filename; filename << settings::path_output << "track_" << simulation::current_batch - << '_' << simulation::current_gen << '_' << p.id << ".h5"; + << '_' << simulation::current_gen << '_' << p.id_ << ".h5"; // Determine number of coordinates for each particle std::vector n_coords; diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 0f18e245f..5ffae4b20 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -105,41 +105,41 @@ std::vector VolumeCalculation::execute() const for (int i = i_start; i < i_end; i++) { set_particle_seed(i); - p.n_coord = 1; + p.n_coord_ = 1; Position xi {prn(), prn(), prn()}; Position r {lower_left_ + xi*(upper_right_ - lower_left_)}; // TODO: assign directly when xyz is Position - std::copy(&r.x, &r.x + 3, p.coord[0].xyz); - p.coord[0].uvw[0] = 0.5; - p.coord[1].uvw[1] = 0.5; - p.coord[2].uvw[2] = 0.5; + std::copy(&r.x, &r.x + 3, p.coord_[0].xyz); + p.coord_[0].uvw[0] = 0.5; + p.coord_[1].uvw[1] = 0.5; + p.coord_[2].uvw[2] = 0.5; // If this location is not in the geometry at all, move on to next block if (!find_cell(&p, false)) continue; if (domain_type_ == FILTER_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_ == FILTER_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_ == FILTER_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; } }