crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Add trailing underscore on Particle data members

Browse source
This commit is contained in:
Paul Romano 2019-02-27 08:23:43 -06:00
parent 38e39c4486
commit b4ed267d4b
38 changed files with 912 additions and 911 deletions
Download patch file Download diff file Expand all files Collapse all files

79
include/openmc/particle.h
View file

@ -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();

14
src/bremsstrahlung.cpp
View file

@ -28,20 +28,20 @@ std::vector<Bremsstrahlung> 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<int>(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<int>(ParticleType::positron)) {
mat = &model::materials[p.material]->ttb_->positron;
if (p.type_ == static_cast<int>(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<int>(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;
}
}

2
src/cross_sections.cpp
View file

@ -312,7 +312,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
}
} // thermal_tables_
// Finish setting up materials (normalizing densities, etc.)
// Finish setting up.material_s (normalizing densities, etc.)
mat->finalize();
} // materials

2
src/dagmc.cpp
View file

@ -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

2
src/eigenvalue.cpp
View file

@ -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!");

228
src/geometry.cpp
View file

@ -34,21 +34,21 @@ std::vector<int64_t> 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<int, 3> 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<int, 3> 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<int, 3> 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<int, 3> 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<int, 3> 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<double, std::array<int, 3>> 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;
}

20
src/material.cpp
View file

@ -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<int>(ParticleType::neutron)) {
if (p.type_ == static_cast<int>(ParticleType::neutron)) {
this->calculate_neutron_xs(p);
} else if (p.type == static_cast<int>(ParticleType::photon)) {
} else if (p.type_ == static_cast<int>(ParticleType::photon)) {
this->calculate_photon_xs(p);
}
}
@ -664,7 +664,7 @@ void Material::calculate_neutron_xs(const Particle& p) const
int neutron = static_cast<int>(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_);
}
// ========================================================================

12
src/mesh.cpp
View file

@ -385,9 +385,9 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector<int>& 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<int>& 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_;

46
src/output.cpp
View file

@ -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<int>(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";
}

324
src/particle.cpp
View file

@ -51,24 +51,24 @@ void
Particle::clear()
{
// reset any coordinate levels
for (int i=0; i<MAX_COORD; ++i) coord[i].reset();
for (int i=0; i<MAX_COORD; ++i) coord_[i].reset();
}
void
Particle::create_secondary(const double* uvw, double E, int type, bool run_CE)
{
if (n_secondary == MAX_SECONDARY) {
if (n_secondary_ == MAX_SECONDARY) {
fatal_error("Too many secondary particles created.");
}
int64_t n = n_secondary;
secondary_bank[n].particle = type;
secondary_bank[n].wgt = wgt;
std::copy(coord[0].xyz, coord[0].xyz + 3, secondary_bank[n].xyz);
std::copy(uvw, uvw + 3, secondary_bank[n].uvw);
secondary_bank[n].E = E;
if (!run_CE) secondary_bank[n].E = g;
n_secondary += 1;
int64_t n = n_secondary_;
secondary_bank_[n].particle = type_;
secondary_bank_[n].wgt = wgt_;
std::copy(coord_[0].xyz, coord_[0].xyz + 3, secondary_bank_[n].xyz);
std::copy(uvw, uvw + 3, secondary_bank_[n].uvw);
secondary_bank_[n].E = E_;
if (!run_CE) secondary_bank_[n].E = g_;
++n_secondary_;
}
void
@ -78,33 +78,33 @@ Particle::initialize()
clear();
// Set particle to neutron that's alive
type = static_cast<int>(ParticleType::neutron);
alive = true;
type_ = static_cast<int>(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; i<MAX_DELAYED_GROUPS; ++i) {
n_delayed_bank[i] = 0;
n_delayed_bank_[i] = 0;
}
g = 0;
g_ = 0;
// Set up base level coordinates
coord[0].universe = C_NONE;
n_coord = 1;
last_n_coord = 1;
coord_[0].universe = C_NONE;
n_coord_ = 1;
last_n_coord_ = 1;
}
void
@ -114,23 +114,23 @@ Particle::from_source(const Bank* src)
initialize();
// copy attributes from source bank site
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);
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<int>(src->E);
last_g = static_cast<int>(src->E);
E = data::energy_bin_avg[g - 1];
g_ = static_cast<int>(src->E);
last_g_ = static_cast<int>(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<int>(ParticleType::neutron)) {
if (type_ == static_cast<int>(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<int>(ParticleType::electron) ||
type == static_cast<int>(ParticleType::positron)) {
if (type_ == static_cast<int>(ParticleType::electron) ||
type_ == static_cast<int>(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<int>(ParticleType::neutron)) {
global_tally_tracklength += wgt * distance * simulation::material_xs.nu_fission;
type_ == static_cast<int>(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<int>(ParticleType::neutron)) {
global_tally_collision += wgt * simulation::material_xs.nu_fission
type_ == static_cast<int>(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<DAGCell*>(model::cells[last_cell[0]]);
auto cellp = dynamic_cast<DAGCell*>(model::cells[last_cell_[0]]);
// TODO: off-by-one
auto surfp = dynamic_cast<DAGSurface*>(model::surfaces[std::abs(surface) - 1]);
auto surfp = dynamic_cast<DAGSurface*>(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);

32
src/particle_restart.cpp
View file

@ -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<double, 3> 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);

176
src/physics.cpp
View file

@ -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<ParticleType>(p->type)) {
switch (static_cast<ParticleType>(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<ParticleType>(p->type) == ParticleType::neutron) {
msg << " " << reaction_name(p->event_MT) << " with " <<
data::nuclides[p->event_nuclide]->name_ << ". Energy = " << p->E << " eV.";
if (static_cast<ParticleType>(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<int>(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<int>(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<int>(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<int>(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<int>(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<int>(std::round(yield)) - 1; ++i) {
int neutron = static_cast<int>(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<int>(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

14
src/physics_common.cpp
View file

@ -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;
}
}
}

66
src/physics_mg.cpp
View file

@ -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_t>(*size_bank);
i < static_cast<size_t>(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<int>(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<double>(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;
}
}

44
src/plot.cpp
View file

@ -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]));
}

16
src/simulation.cpp
View file

@ -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;
}
}

88
src/tallies/derivative.cpp
View file

@ -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

4
src/tallies/filter_azimuthal.cpp
View file

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

4
src/tallies/filter_cell.cpp
View file

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

2
src/tallies/filter_cellborn.cpp
View file

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

4
src/tallies/filter_cellfrom.cpp
View file

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

14
src/tallies/filter_distribcell.cpp
View file

@ -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;

16
src/tallies/filter_energy.cpp
View file

@ -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);
}

6
src/tallies/filter_energyfunc.cpp
View file

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

2
src/tallies/filter_legendre.cpp
View file

@ -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]);

2
src/tallies/filter_material.cpp
View file

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

2
src/tallies/filter_mesh.cpp
View file

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

4
src/tallies/filter_mu.cpp
View file

@ -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);
}

2
src/tallies/filter_particle.cpp
View file

@ -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);
}

4
src/tallies/filter_polar.cpp
View file

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

4
src/tallies/filter_sph_harm.cpp
View file

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

6
src/tallies/filter_sptl_legendre.cpp
View file

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

4
src/tallies/filter_surface.cpp
View file

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

4
src/tallies/filter_universe.cpp
View file

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

8
src/tallies/filter_zernike.cpp
View file

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

534
src/tallies/tally_scoring.cpp
View file

File diff suppressed because it is too large Load diff

4
src/track_output.cpp
View file

@ -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<int> n_coords;

28
src/volume_calc.cpp
View file

@ -105,41 +105,41 @@ std::vector<VolumeCalculation::Result> 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;
}
}