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

Rename other 'last' members in Particle

Browse source
This commit is contained in:
Paul Romano 2019-03-01 15:52:25 -06:00
parent 368f89697d
commit cd0ca81914
12 changed files with 141 additions and 137 deletions
Download patch file Download diff file Expand all files Collapse all files

20
include/openmc/particle.h
View file

@ -84,14 +84,14 @@ public:
LocalCoord coord_[MAX_COORD]; //!< coordinates for all levels
// Particle coordinates before crossing a surface
int last_n_coord_ {1}; //!< number of current coordinates
int last_cell_[MAX_COORD]; //!< coordinates for all levels
int n_coord_last_ {1}; //!< number of current coordinates
int cell_last_[MAX_COORD]; //!< coordinates for all levels
// Energy data
double E_; //!< post-collision energy in eV
double last_E_; //!< pre-collision energy in eV
double E_last_; //!< pre-collision energy in eV
int g_ {0}; //!< post-collision energy group (MG only)
int last_g_; //!< pre-collision energy group (MG only)
int g_last_; //!< pre-collision energy group (MG only)
// Other physical data
double wgt_ {1.0}; //!< particle weight
@ -104,8 +104,8 @@ public:
//!< current tallies
Position r_last_; //!< previous coordinates
Direction u_last_; //!< previous direction coordinates
double last_wgt_ {1.0}; //!< pre-collision particle weight
double absorb_wgt_ {0.0}; //!< weight absorbed for survival biasing
double wgt_last_ {1.0}; //!< pre-collision particle weight
double wgt_absorb_ {0.0}; //!< weight absorbed for survival biasing
// What event took place
bool fission_ {false}; //!< did particle cause implicit fission
@ -124,11 +124,11 @@ public:
int surface_ {0}; //!< index for surface particle is on
int cell_born_ {-1}; //!< index for cell particle was born in
int material_ {-1}; //!< index for current material
int last_material_ {-1}; //!< index for last material
int material_last_ {-1}; //!< index for last material
// Temperature of current cell
double sqrtkT_ {-1.0}; //!< sqrt(k_Boltzmann * temperature) in eV
double last_sqrtkT_ {0.0}; //!< last temperature
double sqrtkT_last_ {0.0}; //!< last temperature
// Statistical data
int n_collision_ {0}; //!< number of collisions
@ -140,15 +140,19 @@ public:
int64_t n_secondary_ {};
Bank secondary_bank_[MAX_SECONDARY];
// Accessors for position in global coordinates
Position& r() { return coord_[0].r; }
const Position& r() const { return coord_[0].r; }
// Accessors for position in local coordinates
Position& r_local() { return coord_[n_coord_ - 1].r; }
const Position& r_local() const { return coord_[n_coord_ - 1].r; }
// Accessors for direction in global coordinates
Direction& u() { return coord_[0].u; }
const Direction& u() const { return coord_[0].u; }
// Accessors for direction in local coordinates
Direction& u_local() { return coord_[n_coord_ - 1].u; }
const Direction& u_local() const { return coord_[n_coord_ - 1].u; }

4
src/geometry.cpp
View file

@ -146,13 +146,13 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
}
// Set the material and temperature.
p->last_material_ = p->material_;
p->material_last_ = p->material_;
if (c.material_.size() > 1) {
p->material_ = c.material_[p->cell_instance_];
} else {
p->material_ = c.material_[0];
}
p->last_sqrtkT_ = p->sqrtkT_;
p->sqrtkT_last_ = p->sqrtkT_;
if (c.sqrtkT_.size() > 1) {
p->sqrtkT_ = c.sqrtkT_[p->cell_instance_];
} else {

28
src/particle.cpp
View file

@ -96,7 +96,7 @@ Particle::from_source(const Bank* src)
// copy attributes from source bank site
type_ = src->particle;
wgt_ = src->wgt;
last_wgt_ = src->wgt;
wgt_last_ = src->wgt;
this->r() = src->r;
this->u() = src->u;
r_last_current_ = src->r;
@ -107,10 +107,10 @@ Particle::from_source(const Bank* src)
g_ = 0;
} else {
g_ = static_cast<int>(src->E);
last_g_ = static_cast<int>(src->E);
g_last_ = static_cast<int>(src->E);
E_ = data::energy_bin_avg[g_ - 1];
}
last_E_ = E_;
E_last_ = E_;
}
void
@ -150,8 +150,8 @@ Particle::transport()
}
// Store pre-collision particle properties
last_wgt_ = wgt_;
last_E_ = E_;
wgt_last_ = wgt_;
E_last_ = E_;
u_last_ = this->u();
r_last_ = this->r();
@ -177,7 +177,7 @@ Particle::transport()
// Calculate microscopic and macroscopic cross sections
if (material_ != MATERIAL_VOID) {
if (settings::run_CE) {
if (material_ != last_material_ || sqrtkT_ != last_sqrtkT_) {
if (material_ != material_last_ || sqrtkT_ != sqrtkT_last_) {
// 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.
@ -191,7 +191,7 @@ Particle::transport()
// Finally, update the particle group while we have already checked
// for if multi-group
last_g_ = g_;
g_last_ = g_;
}
} else {
simulation::material_xs.total = 0.0;
@ -251,9 +251,9 @@ Particle::transport()
// Saving previous cell data
for (int j = 0; j < n_coord_; ++j) {
last_cell_[j] = coord_[j].cell;
cell_last_[j] = coord_[j].cell;
}
last_n_coord_ = n_coord_;
n_coord_last_ = n_coord_;
if (lattice_translation[0] != 0 || lattice_translation[1] != 0 ||
lattice_translation[2] != 0) {
@ -323,7 +323,7 @@ Particle::transport()
// Set last material to none since cross sections will need to be
// re-evaluated
last_material_ = C_NONE;
material_last_ = C_NONE;
// Set all directions to base level -- right now, after a collision, only
// the base level directions are changed
@ -449,7 +449,7 @@ Particle::cross_surface()
this->u() = u / u.norm();
// Reassign particle's cell and surface
coord_[0].cell = last_cell_[last_n_coord_ - 1];
coord_[0].cell = cell_last_[n_coord_last_ - 1];
surface_ = -surface_;
// If a reflective surface is coincident with a lattice or universe
@ -533,13 +533,13 @@ 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[cell_last_[0]]);
// TODO: off-by-one
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_;
material_last_ = material_;
sqrtkT_last_ = sqrtkT_;
// set new cell value
coord_[0].cell = i_cell;
cell_instance_ = 0;

6
src/particle_restart.cpp
View file

@ -55,12 +55,12 @@ void read_particle_restart(Particle& p, int& previous_run_mode)
}
// Set particle last attributes
p.last_wgt_ = p.wgt_;
p.wgt_last_ = p.wgt_;
p.r_last_current_ = p.r();
p.r_last_ = p.r();
p.u_last_ = p.u();
p.last_E_ = p.E_;
p.last_g_ = p.g_;
p.E_last_ = p.E_;
p.g_last_ = p.g_;
// Close hdf5 file
file_close(file_id);

14
src/physics.cpp
View file

@ -56,7 +56,7 @@ void collision(Particle* p)
if (p->E_ < settings::energy_cutoff[type]) {
p->alive_ = false;
p->wgt_ = 0.0;
p->last_wgt_ = 0.0;
p->wgt_last_ = 0.0;
}
// Display information about collision
@ -113,7 +113,7 @@ 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->wgt_absorb_ = 0.0;
}
if (!p->alive_) return;
@ -122,7 +122,7 @@ void sample_neutron_reaction(Particle* p)
scatter(p, i_nuclide);
// Advance URR seed stream 'N' times after energy changes
if (p->E_ != p->last_E_) {
if (p->E_ != p->E_last_) {
prn_set_stream(STREAM_URR_PTABLE);
advance_prn_seed(data::nuclides.size());
prn_set_stream(STREAM_TRACKING);
@ -542,16 +542,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->wgt_absorb_ = 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->wgt_absorb_;
p->wgt_last_ = 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->wgt_absorb_ * simulation::micro_xs[
i_nuclide].nu_fission / simulation::micro_xs[i_nuclide].absorption;
}
} else {

4
src/physics_common.cpp
View file

@ -14,10 +14,10 @@ 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_;
p->wgt_last_ = p->wgt_;
} else {
p->wgt_ = 0.;
p->last_wgt_ = 0.;
p->wgt_last_ = 0.;
p->alive_ = false;
}
}

12
src/physics_mg.cpp
View file

@ -63,7 +63,7 @@ sample_reaction(Particle* p)
if (simulation::material_xs.absorption > 0.) {
absorption(p);
} else {
p->absorb_wgt_ = 0.;
p->wgt_absorb_ = 0.;
}
if (!p->alive_) return;
@ -82,7 +82,7 @@ scatter(Particle* p)
{
// Adjust indices for Fortran to C++ indexing
// TODO: Remove when no longer needed
int gin = p->last_g_ - 1;
int gin = p->g_last_ - 1;
int gout = p->g_ - 1;
int i_mat = p->material_;
data::macro_xs[i_mat].sample_scatter(gin, gout, p->mu_, p->wgt_);
@ -203,16 +203,16 @@ absorption(Particle* p)
{
if (settings::survival_biasing) {
// Determine weight absorbed in survival biasing
p->absorb_wgt_ = p->wgt_ *
p->wgt_absorb_ = 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->wgt_absorb_;
p->wgt_last_ = p->wgt_;
// Score implicit absorpion estimate of keff
#pragma omp atomic
global_tally_absorption += p->absorb_wgt_ *
global_tally_absorption += p->wgt_absorb_ *
simulation::material_xs.nu_fission /
simulation::material_xs.absorption;
} else {

50
src/tallies/derivative.cpp
View file

@ -323,7 +323,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
if (material.nuclide_[i] == p->event_nuclide_) break;
const auto& nuc {*data::nuclides[p->event_nuclide_]};
if (!multipole_in_range(&nuc, p->last_E_)) {
if (!multipole_in_range(&nuc, p->E_last_)) {
score *= flux_deriv;
break;
}
@ -334,7 +334,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
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->E_last_, p->sqrtkT_);
score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i)
/ simulation::material_xs.total;
} else {
@ -347,7 +347,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
- 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->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_s * material.atom_density_(i)
/ (simulation::material_xs.total
- simulation::material_xs.absorption);
@ -360,7 +360,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
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->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_a * material.atom_density_(i)
/ simulation::material_xs.absorption;
} else {
@ -372,7 +372,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
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->E_last_, p->sqrtkT_);
score *= flux_deriv + dsig_f * material.atom_density_(i)
/ simulation::material_xs.fission;
} else {
@ -386,7 +386,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
/ 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->E_last_, 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->E_last_)) {
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->E_last_)
&& 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->E_last_, 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->E_last_, 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->E_last_)
&& (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->E_last_, 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->E_last_, 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->E_last_)
&& 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->E_last_, 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->E_last_, 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->E_last_)
&& 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->E_last_, 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->E_last_, 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->E_last_)
&& 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->E_last_, 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->E_last_, p->sqrtkT_);
score *= flux_deriv
+ dsig_f / simulation::micro_xs[i_nuclide].fission;
} else {
@ -597,7 +597,7 @@ 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->E_last_)) {
// phi is proportional to e^(-Sigma_tot * dist)
// (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist
// (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist
@ -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->E_last_)) {
// phi is proportional to Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
const auto& micro_xs {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->E_last_, 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_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->n_coord_last_; i++) {
auto search = map_.find(p->cell_last_[i]);
if (search != map_.end()) {
match.bins_.push_back(search->second);
match.weights_.push_back(1.0);

4
src/tallies/filter_energy.cpp
View file

@ -45,13 +45,13 @@ const
if (estimator == ESTIMATOR_TRACKLENGTH) {
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->g_last_);
}
match.weights_.push_back(1.0);
} else {
// Get the pre-collision energy of the particle.
auto E = p->last_E_;
auto E = p->E_last_;
// Bin the energy.
if (E >= bins_.front() && E <= bins_.back()) {

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->E_last_ >= energy_.front() && p->E_last_ <= energy_.back()) {
// Search for the incoming energy bin.
auto i = lower_bound_index(energy_.begin(), energy_.end(), p->last_E_);
auto i = lower_bound_index(energy_.begin(), energy_.end(), p->E_last_);
// Compute the interpolation factor between the nearest bins.
double f = (p->last_E_ - energy_[i]) / (energy_[i+1] - energy_[i]);
double f = (p->E_last_ - energy_[i]) / (energy_[i+1] - energy_[i]);
// Interpolate on the lin-lin grid.
match.bins_.push_back(0);

126
src/tallies/tally_scoring.cpp
View file

@ -323,7 +323,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
auto& tally {*model::tallies[i_tally]};
// Get the pre-collision energy of the particle.
auto E = p->last_E_;
auto E = p->E_last_;
for (auto i = 0; i < tally.scores_.size(); ++i) {
auto score_bin = tally.scores_[i];
@ -342,9 +342,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p->last_wgt_ + p->absorb_wgt_;
score = p->wgt_last_ + p->wgt_absorb_;
} else {
score = p->last_wgt_;
score = p->wgt_last_;
}
if (p->type_ == Particle::Type::neutron ||
@ -366,9 +366,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = (p->last_wgt_ + p->absorb_wgt_) * flux;
score = (p->wgt_last_ + p->wgt_absorb_) * flux;
} else {
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
}
} else {
@ -389,9 +389,9 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p->last_wgt_ + p->absorb_wgt_;
score = p->wgt_last_ + p->wgt_absorb_;
} else {
score = p->last_wgt_;
score = p->wgt_last_;
}
score *= flux / simulation::material_xs.total;
} else {
@ -408,7 +408,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
if (p->event_ != EVENT_SCATTER) continue;
// Since only scattering events make it here, again we can use the
// weight entering the collision as the estimator for the reaction rate
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
} else {
if (i_nuclide >= 0) {
score = (simulation::micro_xs[i_nuclide].total
@ -432,13 +432,13 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
(p->event_mt_ >= N_N1 && p->event_mt_ <= N_NC)) {
// Don't waste time on very common reactions we know have
// multiplicities of one.
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
} else {
// Get yield and apply to score
auto m =
data::nuclides[p->event_nuclide_]->reaction_index_[p->event_mt_];
const auto& rxn {*data::nuclides[p->event_nuclide_]->reactions_[m]};
score = p->last_wgt_ * flux * (*rxn.products_[0].yield_)(E);
score = p->wgt_last_ * flux * (*rxn.products_[0].yield_)(E);
}
break;
@ -448,13 +448,13 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
if (settings::survival_biasing) {
// No absorption events actually occur if survival biasing is on --
// just use weight absorbed in survival biasing
score = p->absorb_wgt_ * flux;
score = p->wgt_absorb_ * flux;
} else {
// Skip any event where the particle wasn't absorbed
if (p->event_ == EVENT_SCATTER) continue;
// All fission and absorption events will contribute here, so we
// can just use the particle's weight entering the collision
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
}
} else {
if (i_nuclide >= 0) {
@ -475,7 +475,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// calculate fraction of absorptions that would have resulted in
// fission
if (simulation::micro_xs[p->event_nuclide_].absorption > 0) {
score = p->absorb_wgt_
score = p->wgt_absorb_
* simulation::micro_xs[p->event_nuclide_].fission
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
} else {
@ -487,7 +487,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// All fission events will contribute, so again we can use particle's
// weight entering the collision as the estimate for the fission
// reaction rate
score = p->last_wgt_
score = p->wgt_last_
* simulation::micro_xs[p->event_nuclide_].fission
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
}
@ -517,7 +517,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// calculate fraction of absorptions that would have resulted in
// nu-fission
if (simulation::micro_xs[p->event_nuclide_].absorption > 0) {
score = p->absorb_wgt_
score = p->wgt_absorb_
* simulation::micro_xs[p->event_nuclide_].nu_fission
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
} else {
@ -560,7 +560,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// calculate fraction of absorptions that would have resulted in
// prompt-nu-fission
if (simulation::micro_xs[p->event_nuclide_].absorption > 0) {
score = p->absorb_wgt_
score = p->wgt_absorb_
* simulation::micro_xs[p->event_nuclide_].fission
* data::nuclides[p->event_nuclide_]
->nu(E, ReactionProduct::EmissionMode::prompt)
@ -633,7 +633,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
auto d = filt.groups_[d_bin];
auto yield = data::nuclides[p->event_nuclide_]
->nu(E, ReactionProduct::EmissionMode::delayed, d);
score = p->absorb_wgt_ * yield
score = p->wgt_absorb_ * yield
* simulation::micro_xs[p->event_nuclide_].fission
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
score_fission_delayed_dg(i_tally, d_bin, score,
@ -644,7 +644,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// If the delayed group filter is not present, compute the score
// by multiplying the absorbed weight by the fraction of the
// delayed-nu-fission xs to the absorption xs
score = p->absorb_wgt_
score = p->wgt_absorb_
* simulation::micro_xs[p->event_nuclide_].fission
* data::nuclides[p->event_nuclide_]
->nu(E, ReactionProduct::EmissionMode::delayed)
@ -772,7 +772,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
auto yield
= nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
auto rate = rxn.products_[d].decay_rate_;
score = p->absorb_wgt_ * yield
score = p->wgt_absorb_ * yield
* simulation::micro_xs[p->event_nuclide_].fission
/ simulation::micro_xs[p->event_nuclide_].absorption
* rate * flux;
@ -795,7 +795,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
auto yield
= nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1);
auto rate = rxn.products_[d+1].decay_rate_;
score += rate * p->absorb_wgt_
score += rate * p->wgt_absorb_
* simulation::micro_xs[p->event_nuclide_].fission * yield
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
}
@ -948,7 +948,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
if (simulation::micro_xs[p->event_nuclide_].absorption > 0
&& nuc.fissionable_) {
const auto& rxn {*nuc.fission_rx_[0]};
score = p->absorb_wgt_ * rxn.q_value_
score = p->wgt_absorb_ * rxn.q_value_
* simulation::micro_xs[p->event_nuclide_].fission
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
}
@ -962,7 +962,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
if (simulation::micro_xs[p->event_nuclide_].absorption > 0
&& nuc.fissionable_) {
const auto& rxn {*nuc.fission_rx_[0]};
score = p->last_wgt_ * rxn.q_value_
score = p->wgt_last_ * rxn.q_value_
* simulation::micro_xs[p->event_nuclide_].fission
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
}
@ -1004,7 +1004,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
if (tally.estimator_ == ESTIMATOR_ANALOG) {
// Check if event MT matches
if (p->event_mt_ != ELASTIC) continue;
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
} else {
if (i_nuclide >= 0) {
if (simulation::micro_xs[i_nuclide].elastic == CACHE_INVALID)
@ -1042,12 +1042,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
double q_value = 0.;
if (score_bin == SCORE_FISS_Q_PROMPT) {
if (nuc.fission_q_prompt_)
q_value = (*nuc.fission_q_prompt_)(p->last_E_);
q_value = (*nuc.fission_q_prompt_)(p->E_last_);
} else if (score_bin == SCORE_FISS_Q_RECOV) {
if (nuc.fission_q_recov_)
q_value = (*nuc.fission_q_recov_)(p->last_E_);
q_value = (*nuc.fission_q_recov_)(p->E_last_);
}
score = p->absorb_wgt_ * q_value
score = p->wgt_absorb_ * q_value
* simulation::micro_xs[p->event_nuclide_].fission
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
}
@ -1062,12 +1062,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
double q_value = 0.;
if (score_bin == SCORE_FISS_Q_PROMPT) {
if (nuc.fission_q_prompt_)
q_value = (*nuc.fission_q_prompt_)(p->last_E_);
q_value = (*nuc.fission_q_prompt_)(p->E_last_);
} else if (score_bin == SCORE_FISS_Q_RECOV) {
if (nuc.fission_q_recov_)
q_value = (*nuc.fission_q_recov_)(p->last_E_);
q_value = (*nuc.fission_q_recov_)(p->E_last_);
}
score = p->last_wgt_ * q_value
score = p->wgt_last_ * q_value
* simulation::micro_xs[p->event_nuclide_].fission
/ simulation::micro_xs[p->event_nuclide_].absorption * flux;
}
@ -1078,10 +1078,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
double q_value = 0.;
if (score_bin == SCORE_FISS_Q_PROMPT) {
if (nuc.fission_q_prompt_)
q_value = (*nuc.fission_q_prompt_)(p->last_E_);
q_value = (*nuc.fission_q_prompt_)(p->E_last_);
} else if (score_bin == SCORE_FISS_Q_RECOV) {
if (nuc.fission_q_recov_)
q_value = (*nuc.fission_q_recov_)(p->last_E_);
q_value = (*nuc.fission_q_recov_)(p->E_last_);
}
score = q_value * simulation::micro_xs[i_nuclide].fission
* atom_density * flux;
@ -1095,10 +1095,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
double q_value = 0.;
if (score_bin == SCORE_FISS_Q_PROMPT) {
if (nuc.fission_q_prompt_)
q_value = (*nuc.fission_q_prompt_)(p->last_E_);
q_value = (*nuc.fission_q_prompt_)(p->E_last_);
} else if (score_bin == SCORE_FISS_Q_RECOV) {
if (nuc.fission_q_recov_)
q_value = (*nuc.fission_q_recov_)(p->last_E_);
q_value = (*nuc.fission_q_recov_)(p->E_last_);
}
score += q_value * simulation::micro_xs[j_nuclide].fission
* atom_density * flux;
@ -1118,7 +1118,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
if (tally.estimator_ == ESTIMATOR_ANALOG) {
// Check if the event MT matches
if (p->event_mt_ != score_bin) continue;
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
} else {
int m;
switch (score_bin) {
@ -1153,7 +1153,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// Any other score is assumed to be a MT number. Thus, we just need
// to check if it matches the MT number of the event
if (p->event_mt_ != score_bin) continue;
score = p->last_wgt_*flux;
score = p->wgt_last_*flux;
} else {
// Any other cross section has to be calculated on-the-fly
if (score_bin < 2) fatal_error("Invalid score type on tally "
@ -1235,7 +1235,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// or are dead and g did not change
if (p->alive_) {
p_u = p->u_last_;
p_g = p->last_g_;
p_g = p->g_last_;
} else {
p_u = p->u_local();
p_g = p->g_;
@ -1245,7 +1245,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// Then the energy group has been changed by the scattering routine
// meaning gin is now in p % last_g
p_u = p->u_last_;
p_g = p->last_g_;
p_g = p->g_last_;
} else {
// No scatter, no change in g.
@ -1287,9 +1287,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p->last_wgt_ + p->absorb_wgt_;
score = p->wgt_last_ + p->wgt_absorb_;
} else {
score = p->last_wgt_;
score = p->wgt_last_;
}
score *= flux / simulation::material_xs.total;
} else {
@ -1305,9 +1305,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p->last_wgt_ + p->absorb_wgt_;
score = p->wgt_last_ + p->wgt_absorb_;
} else {
score = p->last_wgt_;
score = p->wgt_last_;
}
//TODO: should flux be multiplied in above instead of below?
if (i_nuclide >= 0) {
@ -1335,9 +1335,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p->last_wgt_ + p->absorb_wgt_;
score = p->wgt_last_ + p->wgt_absorb_;
} else {
score = p->last_wgt_;
score = p->wgt_last_;
}
if (i_nuclide >= 0) {
score *= flux
@ -1366,13 +1366,13 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
if (p->event_ != EVENT_SCATTER) continue;
// Since only scattering events make it here, again we can use the
// weight entering the collision as the estimator for the reaction rate
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
if (i_nuclide >= 0) {
score *= atom_density
* get_nuclide_xs(i_nuclide, MG_GET_XS_SCATTER_FMU_MULT,
p->last_g_, &p->g_, &p->mu_, nullptr)
p->g_last_, &p->g_, &p->mu_, nullptr)
/ get_macro_xs(p->material_, MG_GET_XS_SCATTER_FMU_MULT,
p->last_g_, &p->g_, &p->mu_, nullptr);
p->g_last_, &p->g_, &p->mu_, nullptr);
}
} else {
if (i_nuclide >= 0) {
@ -1400,9 +1400,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
if (i_nuclide >= 0) {
score *= atom_density
* get_nuclide_xs(i_nuclide, MG_GET_XS_SCATTER_FMU,
p->last_g_, &p->g_, &p->mu_, nullptr)
p->g_last_, &p->g_, &p->mu_, nullptr)
/ get_macro_xs(p->material_, MG_GET_XS_SCATTER_FMU,
p->last_g_, &p->g_, &p->mu_, nullptr);
p->g_last_, &p->g_, &p->mu_, nullptr);
}
} else {
if (i_nuclide >= 0) {
@ -1420,13 +1420,13 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
if (settings::survival_biasing) {
// No absorption events actually occur if survival biasing is on --
// just use weight absorbed in survival biasing
score = p->absorb_wgt_ * flux;
score = p->wgt_absorb_ * flux;
} else {
// Skip any event where the particle wasn't absorbed
if (p->event_ == EVENT_SCATTER) continue;
// All fission and absorption events will contribute here, so we
// can just use the particle's weight entering the collision
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
}
if (i_nuclide >= 0) {
score *= atom_density
@ -1450,14 +1450,14 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// fission
score = p->absorb_wgt_ * flux;
score = p->wgt_absorb_ * flux;
} else {
// Skip any non-absorption events
if (p->event_ == EVENT_SCATTER) continue;
// All fission events will contribute, so again we can use particle's
// weight entering the collision as the estimate for the fission
// reaction rate
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
}
if (i_nuclide >= 0) {
score *= atom_density
@ -1493,7 +1493,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// nu-fission
score = p->absorb_wgt_ * flux;
score = p->wgt_absorb_ * flux;
if (i_nuclide >= 0) {
score *= atom_density
* get_nuclide_xs(i_nuclide, MG_GET_XS_NU_FISSION, p_g)
@ -1543,7 +1543,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// prompt-nu-fission
score = p->absorb_wgt_ * flux;
score = p->wgt_absorb_ * flux;
if (i_nuclide >= 0) {
score *= atom_density
* get_nuclide_xs(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g)
@ -1606,7 +1606,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
score = p->absorb_wgt_ * flux;
score = p->wgt_absorb_ * flux;
if (i_nuclide >= 0) {
score *=
get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION,
@ -1626,7 +1626,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// If the delayed group filter is not present, compute the score
// by multiplying the absorbed weight by the fraction of the
// delayed-nu-fission xs to the absorption xs
score = p->absorb_wgt_ * flux;
score = p->wgt_absorb_ * flux;
if (i_nuclide >= 0) {
score *=
get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g)
@ -1728,7 +1728,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// Tally each delayed group bin individually
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
score = p->absorb_wgt_ * flux;
score = p->wgt_absorb_ * flux;
if (i_nuclide >= 0) {
score *=
get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE,
@ -1756,14 +1756,14 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
score = 0.;
for (auto d = 0; d < data::num_delayed_groups; ++d) {
if (i_nuclide >= 0) {
score += p->absorb_wgt_ * flux
score += p->wgt_absorb_ * flux
* get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE,
p_g, nullptr, nullptr, &d)
* get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION,
p_g, nullptr, nullptr, &d)
/ get_macro_xs(p->material_, MG_GET_XS_ABSORPTION, p_g);
} else {
score += p->absorb_wgt_ * flux
score += p->wgt_absorb_ * flux
* get_macro_xs(p->material_, MG_GET_XS_DECAY_RATE,
p_g, nullptr, nullptr, &d)
* get_macro_xs(p->material_, MG_GET_XS_DELAYED_NU_FISSION,
@ -1871,14 +1871,14 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// fission scaled by the Q-value
score = p->absorb_wgt_ * flux;
score = p->wgt_absorb_ * flux;
} else {
// Skip any non-absorption events
if (p->event_ == EVENT_SCATTER) continue;
// All fission events will contribute, so again we can use particle's
// weight entering the collision as the estimate for the fission
// reaction rate
score = p->last_wgt_ * flux;
score = p->wgt_last_ * flux;
}
if (i_nuclide >= 0) {
score *= atom_density
@ -2127,9 +2127,9 @@ void score_collision_tally(const Particle* p)
// Determine the collision estimate of the flux
double flux;
if (!settings::survival_biasing) {
flux = p->last_wgt_ / simulation::material_xs.total;
flux = p->wgt_last_ / simulation::material_xs.total;
} else {
flux = (p->last_wgt_ + p->absorb_wgt_) / simulation::material_xs.total;
flux = (p->wgt_last_ + p->wgt_absorb_) / simulation::material_xs.total;
}
for (auto i_tally : model::active_collision_tallies) {