From 5c6633c2c36b9138ee750440435883b1b5dfae29 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Feb 2019 08:22:16 -0600 Subject: [PATCH 1/9] Make sure temperatures/volumes are read correctly in Material.from_xml --- openmc/material.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index d7c93ebae..f8f122308 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -930,7 +930,10 @@ class Material(IDManagerMixin): mat_id = int(elem.get('id')) mat = cls(mat_id) mat.name = elem.get('name') - mat.temperature = elem.get('temperature') + if 'temperature' in elem.attrib: + mat.temperature = float(elem.get('temperature')) + if 'volume' in elem.attrib: + mat.volume = float(elem.get('volume')) mat.depletable = bool(elem.get('depletable')) # Get each nuclide From 4a890758905cc091b1613ee3988fadfc4587abce Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Feb 2019 08:54:19 -0600 Subject: [PATCH 2/9] Use 0-based indices for materials consistently --- openmc/capi/cell.py | 7 ++--- openmc/capi/filter.py | 7 ++--- openmc/capi/material.py | 2 +- src/bremsstrahlung.cpp | 4 +-- src/cell.cpp | 5 ++-- src/geometry.cpp | 10 ++----- src/geometry_aux.cpp | 15 +++++------ src/material.cpp | 46 ++++++++++++++++----------------- src/mgxs_interface.cpp | 4 +-- src/particle.cpp | 12 ++++----- src/physics.cpp | 7 +++-- src/physics_mg.cpp | 6 ++--- src/plot.cpp | 4 +-- src/tallies/derivative.cpp | 9 +++---- src/tallies/filter_material.cpp | 2 +- src/tallies/tally_scoring.cpp | 36 +++++++++++++------------- src/volume_calc.cpp | 13 ++++------ 17 files changed, 82 insertions(+), 107 deletions(-) diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 894e4c79e..7ad5c75d7 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -107,12 +107,9 @@ class Cell(_FortranObjectWithID): if fill_type.value == 1: if n.value > 1: - #TODO: off-by-one - return [Material(index=i+1 if i >= 0 else i) - for i in indices[:n.value]] + return [Material(index=i) for i in indices[:n.value]] else: - #TODO: off-by-one - index = indices[0] + 1 if indices[0] >= 0 else indices[0] + index = indices[0] return Material(index=index) else: raise NotImplementedError diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 0f5631519..0b85e026e 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -233,16 +233,13 @@ class MaterialFilter(Filter): materials = POINTER(c_int32)() n = c_int32() _dll.openmc_material_filter_get_bins(self._index, materials, n) - #TODO: fix this off-by-one when materials become 0-indexed - return [Material(index=materials[i]+1) for i in range(n.value)] + return [Material(index=materials[i]) for i in range(n.value)] @bins.setter def bins(self, materials): # Get material indices as int32_t[] n = len(materials) - #TODO: fix this off-by-one when materials become 0-indexed - bins = (c_int32*n)(*(m._index-1 for m in materials)) - + bins = (c_int32*n)(*(m._index for m in materials)) _dll.openmc_material_filter_set_bins(self._index, n, bins) diff --git a/openmc/capi/material.py b/openmc/capi/material.py index 4bf19a1e0..981185338 100644 --- a/openmc/capi/material.py +++ b/openmc/capi/material.py @@ -227,7 +227,7 @@ class _MaterialMapping(Mapping): def __iter__(self): for i in range(len(self)): - yield Material(index=i + 1).id + yield Material(index=i).id def __len__(self): return _dll.n_materials() diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index f0ece1a53..2b44154da 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -36,9 +36,9 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) // Get bremsstrahlung data for this material and particle type BremsstrahlungData* mat; if (p.type == static_cast(ParticleType::positron)) { - mat = &model::materials[p.material -1]->ttb_->positron; + mat = &model::materials[p.material]->ttb_->positron; } else { - mat = &model::materials[p.material -1]->ttb_->electron; + mat = &model::materials[p.material]->ttb_->electron; } double e = std::log(p.E); diff --git a/src/cell.cpp b/src/cell.cpp index 30ca4b7cf..2608d8f6f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -727,9 +727,8 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n, int i_mat = indices[i]; if (i_mat == MATERIAL_VOID) { c.material_.push_back(MATERIAL_VOID); - } else if (i_mat >= 1 && i_mat <= model::materials.size()) { - //TODO: off-by-one - c.material_.push_back(i_mat - 1); + } else if (i_mat >= 0 && i_mat < model::materials.size()) { + c.material_.push_back(i_mat); } else { set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; diff --git a/src/geometry.cpp b/src/geometry.cpp index 38dba164b..fb081c313 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -147,16 +147,10 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) // Set the material and temperature. p->last_material = p->material; - int32_t mat; if (c.material_.size() > 1) { - mat = c.material_[p->cell_instance]; + p->material = c.material_[p->cell_instance]; } else { - mat = c.material_[0]; - } - if (mat == MATERIAL_VOID) { - p->material = MATERIAL_VOID; - } else { - p->material = mat + 1; + p->material = c.material_[0]; } p->last_sqrtkT = p->sqrtkT; if (c.sqrtkT_.size() > 1) { diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index d38b5f8f5..c118c04ee 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -85,18 +85,17 @@ adjust_indices() } } else { c->type_ = FILL_MATERIAL; - for (auto it = c->material_.begin(); it != c->material_.end(); it++) { - int32_t mid = *it; - if (mid != MATERIAL_VOID) { - auto search = model::material_map.find(mid); - if (search != model::material_map.end()) { - *it = search->second; - } else { + for (auto& mat_id : c->material_) { + if (mat_id != MATERIAL_VOID) { + auto search = model::material_map.find(mat_id); + if (search == model::material_map.end()) { std::stringstream err_msg; - err_msg << "Could not find material " << mid + err_msg << "Could not find material " << mat_id << " specified on cell " << c->id_; fatal_error(err_msg); } + // Change from ID to index + mat_id = search->second; } } } diff --git a/src/material.cpp b/src/material.cpp index 3153d585a..8236a4ef0 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -932,7 +932,7 @@ openmc_get_material_index(int32_t id, int32_t* index) set_errmsg("No material exists with ID=" + std::to_string(id) + "."); return OPENMC_E_INVALID_ID; } else { - *index = it->second + 1; + *index = it->second; return 0; } } @@ -941,8 +941,8 @@ extern "C" int openmc_material_add_nuclide(int32_t index, const char* name, double density) { int err = 0; - if (index >= 1 && index <= model::materials.size()) { - Material* m = model::materials[index - 1]; + if (index >= 0 && index < model::materials.size()) { + Material* m = model::materials[index]; // Check if nuclide is already in material for (int i = 0; i < m->nuclide_.size(); ++i) { @@ -987,8 +987,8 @@ openmc_material_add_nuclide(int32_t index, const char* name, double density) extern "C" int openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n) { - if (index >= 1 && index <= model::materials.size()) { - auto& mat = model::materials[index - 1]; + if (index >= 0 && index < model::materials.size()) { + auto& mat = model::materials[index]; if (!mat->nuclide_.empty()) { *nuclides = mat->nuclide_.data(); *densities = mat->atom_density_.data(); @@ -1007,8 +1007,8 @@ openmc_material_get_densities(int32_t index, int** nuclides, double** densities, extern "C" int openmc_material_get_fissionable(int32_t index, bool* fissionable) { - if (index >= 1 && index <= model::materials.size()) { - *fissionable = model::materials[index - 1]->fissionable_; + if (index >= 0 && index < model::materials.size()) { + *fissionable = model::materials[index]->fissionable_; return 0; } else { set_errmsg("Index in materials array is out of bounds."); @@ -1019,8 +1019,8 @@ openmc_material_get_fissionable(int32_t index, bool* fissionable) extern "C" int openmc_material_get_id(int32_t index, int32_t* id) { - if (index >= 1 && index <= model::materials.size()) { - *id = model::materials[index - 1]->id_; + if (index >= 0 && index < model::materials.size()) { + *id = model::materials[index]->id_; return 0; } else { set_errmsg("Index in materials array is out of bounds."); @@ -1031,8 +1031,8 @@ openmc_material_get_id(int32_t index, int32_t* id) extern "C" int openmc_material_get_volume(int32_t index, double* volume) { - if (index >= 1 && index <= model::materials.size()) { - Material* m = model::materials[index - 1]; + if (index >= 0 && index < model::materials.size()) { + Material* m = model::materials[index]; if (m->volume_ >= 0.0) { *volume = m->volume_; return 0; @@ -1051,8 +1051,8 @@ openmc_material_get_volume(int32_t index, double* volume) extern "C" int openmc_material_set_density(int32_t index, double density, const char* units) { - if (index >= 1 && index <= model::materials.size()) { - return model::materials[index - 1]->set_density(density, units); + if (index >= 0 && index < model::materials.size()) { + return model::materials[index]->set_density(density, units); } else { set_errmsg("Index in materials array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; @@ -1062,9 +1062,8 @@ openmc_material_set_density(int32_t index, double density, const char* units) extern "C" int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density) { - if (index >= 1 && index <= model::materials.size()) { - // TODO: off-by-one - auto& mat {model::materials[index - 1]}; + if (index >= 0 && index < model::materials.size()) { + auto& mat {model::materials[index]}; if (n != mat->nuclide_.size()) { mat->nuclide_.resize(n); mat->atom_density_ = xt::zeros({n}); @@ -1099,9 +1098,9 @@ openmc_material_set_densities(int32_t index, int n, const char** name, const dou extern "C" int openmc_material_set_id(int32_t index, int32_t id) { - if (index >= 1 && index <= model::materials.size()) { - model::materials[index - 1]->id_ = id; - model::material_map[id] = index - 1; + if (index >= 0 && index < model::materials.size()) { + model::materials[index]->id_ = id; + model::material_map[id] = index; return 0; } else { set_errmsg("Index in materials array is out of bounds."); @@ -1112,8 +1111,8 @@ openmc_material_set_id(int32_t index, int32_t id) extern "C" int openmc_material_set_volume(int32_t index, double volume) { - if (index >= 1 && index <= model::materials.size()) { - Material* m = model::materials[index - 1]; + if (index >= 0 && index < model::materials.size()) { + Material* m = model::materials[index]; if (volume >= 0.0) { m->volume_ = volume; return 0; @@ -1130,9 +1129,8 @@ openmc_material_set_volume(int32_t index, double volume) extern "C" int openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end) { - // TODO: off-by-one - if (index_start) *index_start = model::materials.size() + 1; - if (index_end) *index_end = model::materials.size() + n; + if (index_start) *index_start = model::materials.size(); + if (index_end) *index_end = model::materials.size() + n - 1; for (int32_t i = 0; i < n; i++) { model::materials.push_back(new Material()); } diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 402799c50..581e6b178 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -251,7 +251,7 @@ void calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs) { - data::macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs, + data::macro_xs[i_mat].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs, nu_fiss_xs); } @@ -302,7 +302,7 @@ get_macro_xs(int index, int xstype, int gin, const int* gout, } else { dg_c_p = dg; } - return data::macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p); + return data::macro_xs[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p); } //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index bcd981b26..ccce1502f 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -84,8 +84,8 @@ Particle::initialize() // clear attributes surface = 0; cell_born = C_NONE; - material = 0; - last_material = 0; + material = C_NONE; + last_material = C_NONE; last_sqrtkT = 0; wgt = 1.0; last_wgt = 1.0; @@ -201,7 +201,7 @@ Particle::transport() // 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 - 1]->calculate_xs(*this); + model::materials[material]->calculate_xs(*this); } } else { // Get the MG data @@ -346,7 +346,7 @@ Particle::transport() // Set last material to none since cross sections will need to be // re-evaluated - last_material = F90_NONE; + last_material = C_NONE; // Set all uvws to base level -- right now, after a collision, only the // base level uvws are changed @@ -584,9 +584,7 @@ Particle::cross_surface() // set new cell value coord[0].cell = i_cell; cell_instance = 0; - // TODO: off-by-one - int mat = model::cells[i_cell]->material_[0]; - material = (mat == MATERIAL_VOID) ? mat : mat + 1; + material = model::cells[i_cell]->material_[0]; sqrtkT = model::cells[i_cell]->sqrtkT_[0]; return; } diff --git a/src/physics.cpp b/src/physics.cpp index 4d0fbf27b..8b65fe852 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -425,8 +425,7 @@ int sample_nuclide(const Particle* p) double cutoff = prn() * simulation::material_xs.total; // Get pointers to nuclide/density arrays - // TODO: off-by-one - const auto& mat {model::materials[p->material - 1]}; + const auto& mat {model::materials[p->material]}; int n = mat->nuclide_.size(); double prob = 0.0; @@ -451,7 +450,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 - 1]}; + const auto& mat {model::materials[p->material]}; int n = mat->nuclide_.size(); int i = 0; @@ -675,7 +674,7 @@ void scatter(Particle* p, int i_nuclide) // Sample new outgoing angle for isotropic-in-lab scattering // TODO: off-by-one - const auto& mat {model::materials[p->material - 1]}; + 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]) { diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index ab5232c2a..87dd37ef8 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -46,7 +46,7 @@ sample_reaction(Particle* p) // change when sampling fission sites. The following block handles all // absorption (including fission) - if (model::materials[p->material - 1]->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, @@ -84,7 +84,7 @@ scatter(Particle* p) // TODO: Remove when no longer needed int gin = p->last_g - 1; int gout = p->g - 1; - int i_mat = p->material - 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 @@ -176,7 +176,7 @@ 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 - 1].sample_fission_energy(p->g - 1, dg, gout); + data::macro_xs[p->material].sample_fission_energy(p->g - 1, dg, gout); bank_array[i].E = static_cast(gout + 1); bank_array[i].delayed_group = dg + 1; diff --git a/src/plot.cpp b/src/plot.cpp index ec486ae27..6f8ddd5a6 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -676,8 +676,8 @@ void position_rgb(Particle p, Plot pl, RGBColor& rgb, int& id) rgb = WHITE; id = -1; } else { - rgb = pl.colors_[p.material - 1]; - id = model::materials[p.material - 1]->id_; + rgb = pl.colors_[p.material]; + id = model::materials[p.material]->id_; } } else if (PlotColorBy::cells == pl.color_by_) { // Assign color based on cell diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index b5b85fd6f..1944b39dd 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -131,8 +131,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, score *= flux_deriv; return; } - //TODO: off-by-one - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; if (material.id_ != deriv.diff_material) { score *= flux_deriv; return; @@ -574,8 +573,7 @@ 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; - //TODO: off-by-one - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto& deriv : model::tally_derivs) { if (deriv.diff_material != material.id_) continue; @@ -622,8 +620,7 @@ 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; - //TODO: off-by-one - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto& deriv : model::tally_derivs) { if (deriv.diff_material != material.id_) continue; diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp index 05842b056..0dee964de 100644 --- a/src/tallies/filter_material.cpp +++ b/src/tallies/filter_material.cpp @@ -42,7 +42,7 @@ void MaterialFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - auto search = map_.find(p->material - 1); + auto search = map_.find(p->material); if (search != map_.end()) { //TODO: off-by-one match.bins_.push_back(search->second + 1); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 9edac5ec2..cdd80adba 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -600,7 +600,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, score = 0.; // Add up contributions from each nuclide in the material. if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -723,7 +723,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -743,7 +743,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } else { score = 0.; if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -890,7 +890,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -915,7 +915,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } else { score = 0.; if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -986,7 +986,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } } else { if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1022,7 +1022,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } else { score = 0.; if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1096,7 +1096,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, * atom_density * flux; } else { if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1144,7 +1144,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } else { score = 0.; if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1186,7 +1186,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } } else { if (p->material != MATERIAL_VOID) { - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1274,7 +1274,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // To significantly reduce de-referencing, point matxs to the macroscopic // Mgxs for the material of interest - data::macro_xs[p->material - 1].set_angle_index(p_uvw); + data::macro_xs[p->material].set_angle_index(p_uvw); // Do same for nucxs, point it to the microscopic nuclide data of interest if (i_nuclide >= 0) { @@ -1937,7 +1937,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux, int filter_index) { const Tally& tally {*model::tallies[i_tally]}; - const Material& material {*model::materials[p->material-1]}; + const Material& material {*model::materials[p->material]}; // Score all individual nuclide reaction rates. for (auto i = 0; i < material.nuclide_.size(); ++i) { @@ -2051,12 +2051,12 @@ void score_analog_tally_mg(const Particle* p) double atom_density = 0.; if (i_nuclide >= 0) { //TODO: off-by-one - auto j = model::materials[p->material-1] + auto j = model::materials[p->material] ->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; //atom_density = material_atom_density(p->material, j); //TODO: off-by-one - atom_density = model::materials[p->material-1]->atom_density_(j); + atom_density = model::materials[p->material]->atom_density_(j); } score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index, @@ -2110,12 +2110,12 @@ score_tracklength_tally(const Particle* p, double distance) if (i_nuclide >= 0) { if (p->material != MATERIAL_VOID) { //TODO: off-by-one - auto j = model::materials[p->material-1] + auto j = model::materials[p->material] ->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; //atom_density = material_atom_density(p->material, j); //TODO: off-by-one - atom_density = model::materials[p->material-1]->atom_density_(j); + atom_density = model::materials[p->material]->atom_density_(j); } } @@ -2180,12 +2180,12 @@ void score_collision_tally(const Particle* p) double atom_density = 0.; if (i_nuclide >= 0) { //TODO: off-by-one - auto j = model::materials[p->material-1] + auto j = model::materials[p->material] ->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; //atom_density = material_atom_density(p->material, j); //TODO: off-by-one - atom_density = model::materials[p->material-1]->atom_density_(j); + atom_density = model::materials[p->material]->atom_density_(j); } //TODO: consider replacing this "if" with pointers or templates diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index e3a9bb479..0f18e245f 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -117,14 +117,11 @@ std::vector VolumeCalculation::execute() const // If this location is not in the geometry at all, move on to next block if (!find_cell(&p, false)) continue; - // TODO: off-by-one - int i_material = p.material == MATERIAL_VOID ? p.material : p.material - 1; - if (domain_type_ == FILTER_MATERIAL) { - if (i_material != MATERIAL_VOID) { + if (p.material != MATERIAL_VOID) { for (int i_domain = 0; i_domain < n; i_domain++) { - if (model::materials[i_material]->id_ == domain_ids_[i_domain]) { - this->check_hit(i_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; } } @@ -133,7 +130,7 @@ std::vector VolumeCalculation::execute() const 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(i_material, indices[i_domain], hits[i_domain]); + this->check_hit(p.material, indices[i_domain], hits[i_domain]); break; } } @@ -142,7 +139,7 @@ std::vector VolumeCalculation::execute() const 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(i_material, indices[i_domain], hits[i_domain]); + check_hit(p.material, indices[i_domain], hits[i_domain]); break; } } From 6ce85a9f0ce67a858c8a5234ad3e88ccab905a2c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Feb 2019 09:55:28 -0600 Subject: [PATCH 3/9] Fix off-by-one for index_grid --- src/nuclide.cpp | 13 +++++-------- src/physics.cpp | 9 ++++----- src/reaction.cpp | 2 ++ src/tallies/tally_scoring.cpp | 16 +++++++--------- 4 files changed, 18 insertions(+), 22 deletions(-) diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 592bcceb5..a279abf2b 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -282,8 +282,7 @@ void Nuclide::create_derived() reaction_index_[rx->mt_] = i; for (int t = 0; t < kTs_.size(); ++t) { - // TODO: off-by-one - int j = rx->xs_[t].threshold - 1; + int j = rx->xs_[t].threshold; int n = rx->xs_[t].value.size(); auto xs = xt::adapt(rx->xs_[t].value); @@ -467,7 +466,7 @@ void Nuclide::calculate_elastic_xs() const // Get temperature index, grid index, and interpolation factor auto& micro = simulation::micro_xs[i_nuclide_]; int i_temp = micro.index_temp; - int i_grid = micro.index_grid - 1; + int i_grid = micro.index_grid; double f = micro.interp_factor; if (i_temp >= 0) { @@ -549,7 +548,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, // set to -1 to force a segfault in case a developer messes up and tries // to use it with multipole. micro_xs.index_temp = -1; - micro_xs.index_grid = 0; + micro_xs.index_grid = -1; micro_xs.interp_factor = 0.0; } else { @@ -613,8 +612,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, (grid.energy[i_grid + 1]- grid.energy[i_grid]); micro_xs.index_temp = i_temp; - // TODO: off-by-one - micro_xs.index_grid = i_grid + 1; + micro_xs.index_grid = i_grid; micro_xs.interp_factor = f; // Calculate microscopic nuclide total cross section @@ -665,8 +663,7 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, continue; } - // TODO: off-by-one - int threshold = rx->xs_[i_temp].threshold - 1; + int threshold = rx->xs_[i_temp].threshold; if (i_grid >= threshold) { micro_xs.reaction[j] = (1.0 - f)*rx_xs[i_grid - threshold] + f*rx_xs[i_grid - threshold + 1]; diff --git a/src/physics.cpp b/src/physics.cpp index 8b65fe852..b93c2c232 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -597,7 +597,7 @@ void scatter(Particle* p, int i_nuclide) const auto& nuc {data::nuclides[i_nuclide]}; const auto& micro {simulation::micro_xs[i_nuclide]}; int i_temp = micro.index_temp; - int i_grid = micro.index_grid - 1; + int i_grid = micro.index_grid; double f = micro.interp_factor; // For tallying purposes, this routine might be called directly. In that @@ -655,12 +655,11 @@ void scatter(Particle* p, int i_nuclide) // if energy is below threshold for this reaction, skip it const auto& xs {nuc->reactions_[i]->xs_[i_temp]}; - int threshold = xs.threshold - 1; - if (i_grid < threshold) continue; + if (i_grid < xs.threshold) continue; // add to cumulative probability - prob += (1.0 - f)*xs.value[i_grid - threshold] + - f*xs.value[i_grid - threshold + 1]; + prob += (1.0 - f)*xs.value[i_grid - xs.threshold] + + f*xs.value[i_grid - xs.threshold + 1]; } // Perform collision physics for inelastic scattering diff --git a/src/reaction.cpp b/src/reaction.cpp index 66b841ef4..0fcac2e1d 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -42,6 +42,8 @@ Reaction::Reaction(hid_t group, const std::vector& temperatures) // Get threshold index TemperatureXS xs; read_attribute(dset, "threshold_idx", xs.threshold); + // TODO: change HDF5 format so that threshold_idx is 0-based + --xs.threshold; // Read cross section values read_dataset(dset, xs.value); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index cdd80adba..a0c2883b0 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1175,13 +1175,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index, const auto& rxn {*nuc.reactions_[m]}; auto i_temp = simulation::micro_xs[i_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = simulation::micro_xs[i_nuclide].index_grid - 1; + auto i_grid = simulation::micro_xs[i_nuclide].index_grid; auto f = simulation::micro_xs[i_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; - auto threshold = xs.threshold - 1; if (i_grid >= xs.threshold) { - score = ((1.0 - f) * xs.value[i_grid-threshold] - + f * xs.value[i_grid-threshold+1]) * atom_density * flux; + score = ((1.0 - f) * xs.value[i_grid-xs.threshold] + + f * xs.value[i_grid-xs.threshold+1]) * atom_density * flux; } } } else { @@ -1196,13 +1195,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index, const auto& rxn {*nuc.reactions_[m]}; auto i_temp = simulation::micro_xs[j_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = simulation::micro_xs[j_nuclide].index_grid - 1; + auto i_grid = simulation::micro_xs[j_nuclide].index_grid; auto f = simulation::micro_xs[j_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; - auto threshold = xs.threshold - 1; - if (i_grid >= threshold) { - score += ((1.0 - f) * xs.value[i_grid-threshold] - + f * xs.value[i_grid-threshold+1]) * atom_density + if (i_grid >= xs.threshold) { + score += ((1.0 - f) * xs.value[i_grid-xs.threshold] + + f * xs.value[i_grid-xs.threshold+1]) * atom_density * flux; } } From a2e058c8cf6cb125e8c42eb322ca2d15c3b9e858 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Feb 2019 10:07:04 -0600 Subject: [PATCH 4/9] Fix off-by-one on event_nuclide and diff_nuclide --- src/output.cpp | 3 +- src/physics.cpp | 10 ++-- src/state_point.cpp | 2 +- src/tallies/derivative.cpp | 61 ++++++++++------------ src/tallies/tally_scoring.cpp | 96 +++++++++++++++++------------------ 5 files changed, 80 insertions(+), 92 deletions(-) diff --git a/src/output.cpp b/src/output.cpp index 35cd9f9aa..97436dc24 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -657,8 +657,7 @@ write_tallies() case DIFF_NUCLIDE_DENSITY: tallies_out << " Nuclide density derivative Material " << std::to_string(deriv.diff_material) << " Nuclide " - // TODO: off-by-one - << data::nuclides[deriv.diff_nuclide-1]->name_ << "\n"; + << data::nuclides[deriv.diff_nuclide]->name_ << "\n"; break; case DIFF_TEMPERATURE: tallies_out << " Temperature derivative Material " diff --git a/src/physics.cpp b/src/physics.cpp index b93c2c232..23ce8f560 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -63,10 +63,10 @@ void collision(Particle* p) std::stringstream msg; if (static_cast(p->type) == ParticleType::neutron) { msg << " " << reaction_name(p->event_MT) << " with " << - data::nuclides[p->event_nuclide-1]->name_ << ". Energy = " << p->E << " eV."; + data::nuclides[p->event_nuclide]->name_ << ". Energy = " << p->E << " eV."; } else { msg << " " << reaction_name(p->event_MT) << " with " << - data::elements[p->event_nuclide-1].name_ << ". Energy = " << p->E << " eV."; + data::elements[p->event_nuclide].name_ << ". Energy = " << p->E << " eV."; } write_message(msg, 1); } @@ -78,8 +78,7 @@ void sample_neutron_reaction(Particle* p) int i_nuclide = sample_nuclide(p); // Save which nuclide particle had collision with - // TODO: off-by-one - p->event_nuclide = i_nuclide + 1; + 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 @@ -228,8 +227,7 @@ void sample_photon_reaction(Particle* p) // Sample element within material int i_element = sample_element(p); - // TODO: off-by-one - p->event_nuclide = i_element + 1; + p->event_nuclide = i_element; const auto& micro {simulation::micro_photon_xs[i_element]}; const auto& element {data::elements[i_element]}; diff --git a/src/state_point.cpp b/src/state_point.cpp index 2a3f9d2bb..3730c3216 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -125,7 +125,7 @@ openmc_statepoint_write(const char* filename, bool* write_source) write_dataset(deriv_group, "independent variable", "nuclide_density"); //TODO: off-by-one write_dataset(deriv_group, "nuclide", - data::nuclides[deriv.diff_nuclide-1]->name_); + data::nuclides[deriv.diff_nuclide]->name_); } else if (deriv.variable == DIFF_TEMPERATURE) { write_dataset(deriv_group, "independent variable", "temperature"); } else { diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 1944b39dd..1098adc55 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -51,8 +51,7 @@ TallyDerivative::TallyDerivative(pugi::xml_node node) for (auto i = 0; i < data::nuclides.size(); ++i) { if (data::nuclides[i]->name_ == nuclide_name) { found = true; - //TODO: off-by-one - diff_nuclide = i + 1; + diff_nuclide = i; } } if (!found) { @@ -188,7 +187,6 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, // where i is the perturbed nuclide. case DIFF_NUCLIDE_DENSITY: - //TODO: off-by-one throughout on diff_nuclide switch (tally.estimator_) { case ESTIMATOR_ANALOG: @@ -208,7 +206,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, // Find the index of the perturbed nuclide. int i; for (i = 0; i < material.nuclide_.size(); ++i) - if (material.nuclide_[i] == deriv.diff_nuclide - 1) break; + if (material.nuclide_[i] == deriv.diff_nuclide) break; score *= flux_deriv + 1. / material.atom_density_(i); } break; @@ -225,9 +223,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, case SCORE_TOTAL: if (i_nuclide == -1 && simulation::material_xs.total > 0.0) { score *= flux_deriv - + simulation::micro_xs[deriv.diff_nuclide-1].total + + simulation::micro_xs[deriv.diff_nuclide].total / simulation::material_xs.total; - } else if (i_nuclide == deriv.diff_nuclide-1 + } else if (i_nuclide == deriv.diff_nuclide && simulation::micro_xs[i_nuclide].total) { score *= flux_deriv + 1. / atom_density; } else { @@ -239,11 +237,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, if (i_nuclide == -1 && (simulation::material_xs.total - simulation::material_xs.absorption) > 0.0) { score *= flux_deriv - + (simulation::micro_xs[deriv.diff_nuclide-1].total - - simulation::micro_xs[deriv.diff_nuclide-1].absorption) + + (simulation::micro_xs[deriv.diff_nuclide].total + - simulation::micro_xs[deriv.diff_nuclide].absorption) / (simulation::material_xs.total - simulation::material_xs.absorption); - } else if (i_nuclide == deriv.diff_nuclide-1) { + } else if (i_nuclide == deriv.diff_nuclide) { score *= flux_deriv + 1. / atom_density; } else { score *= flux_deriv; @@ -253,9 +251,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, case SCORE_ABSORPTION: if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) { score *= flux_deriv - + simulation::micro_xs[deriv.diff_nuclide-1].absorption + + simulation::micro_xs[deriv.diff_nuclide].absorption / simulation::material_xs.absorption; - } else if (i_nuclide == deriv.diff_nuclide-1 + } else if (i_nuclide == deriv.diff_nuclide && simulation::micro_xs[i_nuclide].absorption) { score *= flux_deriv + 1. / atom_density; } else { @@ -266,9 +264,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, case SCORE_FISSION: if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) { score *= flux_deriv - + simulation::micro_xs[deriv.diff_nuclide-1].fission + + simulation::micro_xs[deriv.diff_nuclide].fission / simulation::material_xs.fission; - } else if (i_nuclide == deriv.diff_nuclide-1 + } else if (i_nuclide == deriv.diff_nuclide && simulation::micro_xs[i_nuclide].fission) { score *= flux_deriv + 1. / atom_density; } else { @@ -279,9 +277,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, case SCORE_NU_FISSION: if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) { score *= flux_deriv - + simulation::micro_xs[deriv.diff_nuclide-1].nu_fission + + simulation::micro_xs[deriv.diff_nuclide].nu_fission / simulation::material_xs.nu_fission; - } else if (i_nuclide == deriv.diff_nuclide-1 + } else if (i_nuclide == deriv.diff_nuclide && simulation::micro_xs[i_nuclide].nu_fission) { score *= flux_deriv + 1. / atom_density; } else { @@ -322,9 +320,9 @@ 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-1) break; + if (material.nuclide_[i] == p->event_nuclide) break; - const auto& nuc {*data::nuclides[p->event_nuclide-1]}; + const auto& nuc {*data::nuclides[p->event_nuclide]}; if (!multipole_in_range(&nuc, p->last_E)) { score *= flux_deriv; break; @@ -333,7 +331,7 @@ 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-1].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); @@ -345,8 +343,8 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_SCATTER: - if (simulation::micro_xs[p->event_nuclide-1].total - - simulation::micro_xs[p->event_nuclide-1].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); @@ -359,7 +357,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_ABSORPTION: - if (simulation::micro_xs[p->event_nuclide-1].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); @@ -371,7 +369,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, break; case SCORE_FISSION: - if (simulation::micro_xs[p->event_nuclide-1].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); @@ -383,9 +381,9 @@ 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-1].fission) { - double nu = simulation::micro_xs[p->event_nuclide-1].nu_fission - / simulation::micro_xs[p->event_nuclide-1].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); @@ -592,9 +590,8 @@ score_track_derivative(const Particle* p, double distance) // phi is proportional to e^(-Sigma_tot * dist) // (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist // (1 / phi) * (d_phi / d_N) = - sigma_tot * dist - //TODO: off-by-one deriv.flux_deriv -= distance - * simulation::micro_xs[deriv.diff_nuclide-1].total; + * simulation::micro_xs[deriv.diff_nuclide].total; break; case DIFF_TEMPERATURE: @@ -635,17 +632,16 @@ void score_collision_derivative(const Particle* p) break; case DIFF_NUCLIDE_DENSITY: - //TODO: off-by-one throughout on diff_nuclide 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) - if (material.nuclide_[i] == deriv.diff_nuclide - 1) break; + if (material.nuclide_[i] == deriv.diff_nuclide) break; // Make sure we found the nuclide. - if (material.nuclide_[i] != deriv.diff_nuclide - 1) { + if (material.nuclide_[i] != deriv.diff_nuclide) { std::stringstream err_msg; err_msg << "Could not find nuclide " - << data::nuclides[deriv.diff_nuclide-1]->name_ << " in material " + << data::nuclides[deriv.diff_nuclide]->name_ << " in material " << material.id_ << " for tally derivative " << deriv.id; fatal_error(err_msg); } @@ -661,8 +657,7 @@ void score_collision_derivative(const Particle* p) for (auto i_nuc : material.nuclide_) { const auto& nuc {*data::nuclides[i_nuc]}; //TODO: off-by-one - if (i_nuc == p->event_nuclide - 1 - && 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 diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index a0c2883b0..96aa1abaa 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -337,9 +337,6 @@ score_general_ce(const Particle* p, int i_tally, int start_index, double score; - //TODO: off-by-one throughout on p->event_nuclide - //TODO: off-by-one throughout on p->material - switch (score_bin) { @@ -445,8 +442,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, } else { // Get yield and apply to score auto m = - data::nuclides[p->event_nuclide-1]->reaction_index_[p->event_MT]; - const auto& rxn {*data::nuclides[p->event_nuclide-1]->reactions_[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); } break; @@ -483,10 +480,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + if (simulation::micro_xs[p->event_nuclide].absorption > 0) { score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide-1].fission - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + * simulation::micro_xs[p->event_nuclide].fission + / simulation::micro_xs[p->event_nuclide].absorption * flux; } else { score = 0.; } @@ -497,8 +494,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // weight entering the collision as the estimate for the fission // reaction rate score = p->last_wgt - * simulation::micro_xs[p->event_nuclide-1].fission - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + * simulation::micro_xs[p->event_nuclide].fission + / simulation::micro_xs[p->event_nuclide].absorption * flux; } } else { if (i_nuclide >= 0) { @@ -525,10 +522,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + if (simulation::micro_xs[p->event_nuclide].absorption > 0) { score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide-1].nu_fission - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + * simulation::micro_xs[p->event_nuclide].nu_fission + / simulation::micro_xs[p->event_nuclide].absorption * flux; } else { score = 0.; } @@ -568,12 +565,12 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + if (simulation::micro_xs[p->event_nuclide].absorption > 0) { score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide-1].fission - * data::nuclides[p->event_nuclide-1] + * simulation::micro_xs[p->event_nuclide].fission + * data::nuclides[p->event_nuclide] ->nu(E, ReactionProduct::EmissionMode::prompt) - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + / simulation::micro_xs[p->event_nuclide].absorption * flux; } else { score = 0.; } @@ -630,8 +627,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - if (simulation::micro_xs[p->event_nuclide-1].absorption > 0 - && data::nuclides[p->event_nuclide-1]->fissionable_) { + if (simulation::micro_xs[p->event_nuclide].absorption > 0 + && data::nuclides[p->event_nuclide]->fissionable_) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; const DelayedGroupFilter& filt @@ -640,11 +637,11 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { auto d = filt.groups_[d_bin]; - auto yield = data::nuclides[p->event_nuclide-1] + auto yield = data::nuclides[p->event_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = p->absorb_wgt * yield - * simulation::micro_xs[p->event_nuclide-1].fission - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + * simulation::micro_xs[p->event_nuclide].fission + / simulation::micro_xs[p->event_nuclide].absorption * flux; score_fission_delayed_dg(i_tally, d_bin+1, score, score_index); } @@ -654,10 +651,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs score = p->absorb_wgt - * simulation::micro_xs[p->event_nuclide-1].fission - * data::nuclides[p->event_nuclide-1] + * simulation::micro_xs[p->event_nuclide].fission + * data::nuclides[p->event_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed) - / simulation::micro_xs[p->event_nuclide-1].absorption *flux; + / simulation::micro_xs[p->event_nuclide].absorption *flux; } } } else { @@ -766,8 +763,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - const auto& nuc {*data::nuclides[p->event_nuclide-1]}; - if (simulation::micro_xs[p->event_nuclide-1].absorption > 0 + const auto& nuc {*data::nuclides[p->event_nuclide]}; + if (simulation::micro_xs[p->event_nuclide].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; if (tally.delayedgroup_filter_ != C_NONE) { @@ -782,8 +779,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; score = p->absorb_wgt * yield - * simulation::micro_xs[p->event_nuclide-1].fission - / simulation::micro_xs[p->event_nuclide-1].absorption + * simulation::micro_xs[p->event_nuclide].fission + / simulation::micro_xs[p->event_nuclide].absorption * rate * flux; score_fission_delayed_dg(i_tally, d_bin+1, score, score_index); @@ -805,8 +802,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; score += rate * p->absorb_wgt - * simulation::micro_xs[p->event_nuclide-1].fission * yield - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + * simulation::micro_xs[p->event_nuclide].fission * yield + / simulation::micro_xs[p->event_nuclide].absorption * flux; } } } @@ -826,7 +823,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, const auto& bank = simulation::fission_bank[i_bank]; auto g = bank.delayed_group; if (g != 0) { - const auto& nuc {*data::nuclides[p->event_nuclide-1]}; + const auto& nuc {*data::nuclides[p->event_nuclide]}; const auto& rxn {*nuc.fission_rx_[0]}; auto rate = rxn.products_[g].decay_rate_; score += simulation::keff * bank.wgt * rate * flux; @@ -953,13 +950,13 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - const auto& nuc {*data::nuclides[p->event_nuclide-1]}; - if (simulation::micro_xs[p->event_nuclide-1].absorption > 0 + const auto& nuc {*data::nuclides[p->event_nuclide]}; + if (simulation::micro_xs[p->event_nuclide].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; score = p->absorb_wgt * rxn.q_value_ - * simulation::micro_xs[p->event_nuclide-1].fission - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + * simulation::micro_xs[p->event_nuclide].fission + / simulation::micro_xs[p->event_nuclide].absorption * flux; } } else { // Skip any non-absorption events @@ -967,13 +964,13 @@ 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 - const auto& nuc {*data::nuclides[p->event_nuclide-1]}; - if (simulation::micro_xs[p->event_nuclide-1].absorption > 0 + const auto& nuc {*data::nuclides[p->event_nuclide]}; + if (simulation::micro_xs[p->event_nuclide].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; score = p->last_wgt * rxn.q_value_ - * simulation::micro_xs[p->event_nuclide-1].fission - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + * simulation::micro_xs[p->event_nuclide].fission + / simulation::micro_xs[p->event_nuclide].absorption * flux; } } } else { @@ -1046,8 +1043,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - const auto& nuc {*data::nuclides[p->event_nuclide-1]}; - if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + const auto& nuc {*data::nuclides[p->event_nuclide]}; + if (simulation::micro_xs[p->event_nuclide].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) @@ -1057,8 +1054,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, q_value = (*nuc.fission_q_recov_)(p->last_E); } score = p->absorb_wgt * q_value - * simulation::micro_xs[p->event_nuclide-1].fission - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + * simulation::micro_xs[p->event_nuclide].fission + / simulation::micro_xs[p->event_nuclide].absorption * flux; } } else { // Skip any non-absorption events @@ -1066,8 +1063,8 @@ 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 - const auto& nuc {*data::nuclides[p->event_nuclide-1]}; - if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) { + const auto& nuc {*data::nuclides[p->event_nuclide]}; + if (simulation::micro_xs[p->event_nuclide].absorption > 0) { double q_value = 0.; if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) @@ -1077,8 +1074,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index, q_value = (*nuc.fission_q_recov_)(p->last_E); } score = p->last_wgt * q_value - * simulation::micro_xs[p->event_nuclide-1].fission - / simulation::micro_xs[p->event_nuclide-1].absorption * flux; + * simulation::micro_xs[p->event_nuclide].fission + / simulation::micro_xs[p->event_nuclide].absorption * flux; } } } else { @@ -1992,8 +1989,7 @@ void score_analog_tally_ce(const Particle* p) // the event nuclide or the total material. Note that the i_nuclide // and flux arguments for score_general are not used for analog // tallies. - //TODO: off-by-one - if (i_nuclide == p->event_nuclide-1 || i_nuclide == -1) + if (i_nuclide == p->event_nuclide || i_nuclide == -1) score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, -1, -1., filter_weight); } From 392fbd7b98382d764bc9987d52dfb6cedf197a5c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Feb 2019 10:30:18 -0600 Subject: [PATCH 5/9] Use 0-based indices for tallies --- openmc/capi/tally.py | 2 +- src/tallies/tally.cpp | 104 +++++++++++++++++------------------------- 2 files changed, 44 insertions(+), 62 deletions(-) diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 0102436e0..b8ce0bee5 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -381,7 +381,7 @@ class _TallyMapping(Mapping): def __iter__(self): for i in range(len(self)): - yield Tally(index=i + 1).id + yield Tally(index=i).id def __len__(self): return _dll.tallies_size() diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index f47285545..4d55350a0 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -1039,9 +1039,8 @@ free_memory_tally() extern "C" int openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end) { - // TODO: off-by-one - if (index_start) *index_start = model::tallies.size() + 1; - if (index_end) *index_end = model::tallies.size() + n; + if (index_start) *index_start = model::tallies.size(); + if (index_end) *index_end = model::tallies.size() + n - 1; for (int i = 0; i < n; ++i) { model::tallies.push_back(std::make_unique()); } @@ -1057,8 +1056,7 @@ openmc_get_tally_index(int32_t id, int32_t* index) return OPENMC_E_INVALID_ID; } - // TODO: off-by-one - *index = it->second + 1; + *index = it->second; return 0; } @@ -1075,26 +1073,24 @@ openmc_get_tally_next_id(int32_t* id) extern "C" int openmc_tally_get_estimator(int32_t index, int* estimator) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - // TODO: off-by-one - *estimator = model::tallies[index-1]->estimator_; + *estimator = model::tallies[index]->estimator_; return 0; } extern "C" int openmc_tally_set_estimator(int32_t index, const char* estimator) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - // TODO: off-by-one - auto& t {model::tallies[index-1]}; + auto& t {model::tallies[index]}; std::string est = estimator; if (est == "analog") { @@ -1113,20 +1109,19 @@ openmc_tally_set_estimator(int32_t index, const char* estimator) extern "C" int openmc_tally_get_id(int32_t index, int32_t* id) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - // TODO: off-by-one - *id = model::tallies[index-1]->id_; + *id = model::tallies[index]->id_; return 0; } extern "C" int openmc_tally_set_id(int32_t index, int32_t id) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } @@ -1137,21 +1132,19 @@ openmc_tally_set_id(int32_t index, int32_t id) return OPENMC_E_INVALID_ID; } - // TODO: off-by-one - model::tallies[index-1]->id_ = id; - model::tally_map[id] = index - 1; + model::tallies[index]->id_ = id; + model::tally_map[id] = index; return 0; } extern "C" int openmc_tally_get_type(int32_t index, int32_t* type) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - //TODO: off-by-one - *type = model::tallies[index-1]->type_; + *type = model::tallies[index]->type_; return 0; } @@ -1159,16 +1152,16 @@ openmc_tally_get_type(int32_t index, int32_t* type) extern "C" int openmc_tally_set_type(int32_t index, const char* type) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } if (strcmp(type, "volume") == 0) { - model::tallies[index-1]->type_ = TALLY_VOLUME; + model::tallies[index]->type_ = TALLY_VOLUME; } else if (strcmp(type, "mesh-surface") == 0) { - model::tallies[index-1]->type_ = TALLY_MESH_SURFACE; + model::tallies[index]->type_ = TALLY_MESH_SURFACE; } else if (strcmp(type, "surface") == 0) { - model::tallies[index-1]->type_ = TALLY_SURFACE; + model::tallies[index]->type_ = TALLY_SURFACE; } else { std::stringstream errmsg; errmsg << "Unknown tally type: " << type; @@ -1182,12 +1175,11 @@ openmc_tally_set_type(int32_t index, const char* type) extern "C" int openmc_tally_get_active(int32_t index, bool* active) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - //TODO: off-by-one - *active = model::tallies[index-1]->active_; + *active = model::tallies[index]->active_; return 0; } @@ -1195,12 +1187,11 @@ openmc_tally_get_active(int32_t index, bool* active) extern "C" int openmc_tally_set_active(int32_t index, bool active) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - //TODO: off-by-one - model::tallies[index-1]->active_ = active; + model::tallies[index]->active_ = active; return 0; } @@ -1208,29 +1199,27 @@ openmc_tally_set_active(int32_t index, bool active) extern "C" int openmc_tally_get_scores(int32_t index, int** scores, int* n) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - //TODO: off-by-one - *scores = model::tallies[index-1]->scores_.data(); - *n = model::tallies[index-1]->scores_.size(); + *scores = model::tallies[index]->scores_.data(); + *n = model::tallies[index]->scores_.size(); return 0; } extern "C" int openmc_tally_set_scores(int32_t index, int n, const char** scores) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } std::vector scores_str(scores, scores+n); try { - //TODO: off-by-one - model::tallies[index-1]->set_scores(scores_str); + model::tallies[index]->set_scores(scores_str); } catch (const std::invalid_argument& ex) { set_errmsg(ex.what()); return OPENMC_E_INVALID_ARGUMENT; @@ -1243,14 +1232,13 @@ extern "C" int openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n) { // Make sure the index fits in the array bounds. - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - //TODO: off-by-one - *n = model::tallies[index-1]->nuclides_.size(); - *nuclides = model::tallies[index-1]->nuclides_.data(); + *n = model::tallies[index]->nuclides_.size(); + *nuclides = model::tallies[index]->nuclides_.data(); return 0; } @@ -1259,7 +1247,7 @@ extern "C" int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides) { // Make sure the index fits in the array bounds. - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } @@ -1279,8 +1267,7 @@ openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides) } } - //TODO: off-by-one - model::tallies[index-1]->nuclides_ = nucs; + model::tallies[index]->nuclides_ = nucs; return 0; } @@ -1288,14 +1275,13 @@ openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides) extern "C" int openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n) { - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - //TODO: off-by-one - *indices = model::tallies[index-1]->filters().data(); - *n = model::tallies[index-1]->filters().size(); + *indices = model::tallies[index]->filters().data(); + *n = model::tallies[index]->filters().size(); return 0; } @@ -1303,15 +1289,14 @@ extern "C" int openmc_tally_set_filters(int32_t index, int n, const int32_t* indices) { // Make sure the index fits in the array bounds. - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } // Set the filters. try { - //TODO: off-by-one - model::tallies[index-1]->set_filters(indices, n); + model::tallies[index]->set_filters(indices, n); } catch (const std::out_of_range& ex) { set_errmsg(ex.what()); return OPENMC_E_OUT_OF_BOUNDS; @@ -1325,13 +1310,12 @@ extern "C" int openmc_tally_reset(int32_t index) { // Make sure the index fits in the array bounds. - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - // TODO: off-by-one - model::tallies[index-1]->reset(); + model::tallies[index]->reset(); return 0; } @@ -1339,13 +1323,12 @@ extern "C" int openmc_tally_get_n_realizations(int32_t index, int32_t* n) { // Make sure the index fits in the array bounds. - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - // TODO: off-by-one - *n = model::tallies[index - 1]->n_realizations_; + *n = model::tallies[index]->n_realizations_; return 0; } @@ -1355,13 +1338,12 @@ extern "C" int openmc_tally_results(int32_t index, double** results, size_t* shape) { // Make sure the index fits in the array bounds. - if (index < 1 || index > model::tallies.size()) { + if (index < 0 || index >= model::tallies.size()) { set_errmsg("Index in tallies array is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } - // TODO: off-by-one - const auto& t {model::tallies[index - 1]}; + const auto& t {model::tallies[index]}; // TODO: Change to zero when xtensor is updated if (t->results_.size() == 1) { set_errmsg("Tally results have not been allocated yet."); From eeaa38c2ab78541576db588bea02504729e85a85 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Feb 2019 10:55:03 -0600 Subject: [PATCH 6/9] Use 0-based indices for filters --- openmc/capi/filter.py | 2 +- openmc/capi/tally.py | 4 ++-- src/tallies/filter.cpp | 19 +++++++------------ src/tallies/filter_cell.cpp | 3 +-- src/tallies/filter_energy.cpp | 4 ++-- src/tallies/filter_legendre.cpp | 4 ++-- src/tallies/filter_material.cpp | 4 ++-- src/tallies/filter_mesh.cpp | 4 ++-- src/tallies/filter_sph_harm.cpp | 2 +- src/tallies/filter_sptl_legendre.cpp | 2 +- src/tallies/filter_zernike.cpp | 2 +- 11 files changed, 22 insertions(+), 28 deletions(-) diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 0b85e026e..05ec85d01 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -402,7 +402,7 @@ class _FilterMapping(Mapping): def __iter__(self): for i in range(len(self)): - yield _get_filter(i + 1).id + yield _get_filter(i).id def __len__(self): return _dll.tally_filters_size() diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index b8ce0bee5..eec3a2b66 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -251,13 +251,13 @@ class Tally(_FortranObjectWithID): filt_idx = POINTER(c_int32)() n = c_int() _dll.openmc_tally_get_filters(self._index, filt_idx, n) - return [_get_filter(filt_idx[i]+1) for i in range(n.value)] + return [_get_filter(filt_idx[i]) for i in range(n.value)] @filters.setter def filters(self, filters): # Get filter indices as int32_t[] n = len(filters) - indices = (c_int32*n)(*(f._index-1 for f in filters)) + indices = (c_int32*n)(*(f._index for f in filters)) _dll.openmc_tally_set_filters(self._index, n, indices) diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 0351fd4dc..e9aae53a7 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -114,7 +114,7 @@ extern "C" size_t tally_filters_size() int verify_filter(int32_t index) { - if (index < 1 || index > model::tally_filters.size()) { + if (index < 0 || index >= model::tally_filters.size()) { set_errmsg("Filter index is out of bounds."); return OPENMC_E_OUT_OF_BOUNDS; } @@ -126,8 +126,7 @@ openmc_filter_get_id(int32_t index, int32_t* id) { if (int err = verify_filter(index)) return err; - // TODO: off-by-one - *id = model::tally_filters[index-1]->id_; + *id = model::tally_filters[index]->id_; return 0; } @@ -141,9 +140,8 @@ openmc_filter_set_id(int32_t index, int32_t id) return OPENMC_E_INVALID_ID; } - // TODO: off-by-one - model::tally_filters[index-1]->id_ = id; - model::filter_map[id] = index - 1; + model::tally_filters[index]->id_ = id; + model::filter_map[id] = index; return 0; } @@ -152,8 +150,7 @@ openmc_filter_get_type(int32_t index, char* type) { if (int err = verify_filter(index)) return err; - // TODO: off-by-one - std::strcpy(type, model::tally_filters[index-1]->type().c_str()); + std::strcpy(type, model::tally_filters[index]->type().c_str()); return 0; } @@ -166,8 +163,7 @@ openmc_get_filter_index(int32_t id, int32_t* index) return OPENMC_E_INVALID_ID; } - // TODO: off-by-one - *index = it->second + 1; + *index = it->second; return 0; } @@ -184,9 +180,8 @@ openmc_get_filter_next_id(int32_t* id) extern "C" int openmc_new_filter(const char* type, int32_t* index) { - allocate_filter(type); - // TODO: off-by-one *index = model::tally_filters.size(); + allocate_filter(type); return 0; } diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp index 80b6678a5..4e25e8541 100644 --- a/src/tallies/filter_cell.cpp +++ b/src/tallies/filter_cell.cpp @@ -76,8 +76,7 @@ openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n) { if (int err = verify_filter(index)) return err; - // TODO: off-by-one - const auto& filt = model::tally_filters[index-1].get(); + const auto& filt = model::tally_filters[index].get(); if (filt->type() != "cell") { set_errmsg("Tried to get cells from a non-cell filter."); return OPENMC_E_INVALID_TYPE; diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp index bd1de4b46..90cd7b3bf 100644 --- a/src/tallies/filter_energy.cpp +++ b/src/tallies/filter_energy.cpp @@ -123,7 +123,7 @@ openmc_energy_filter_get_bins(int32_t index, double** energies, int32_t* n) if (int err = verify_filter(index)) return err; // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. @@ -145,7 +145,7 @@ openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies) if (int err = verify_filter(index)) return err; // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. diff --git a/src/tallies/filter_legendre.cpp b/src/tallies/filter_legendre.cpp index 55f225b77..a5f583a15 100644 --- a/src/tallies/filter_legendre.cpp +++ b/src/tallies/filter_legendre.cpp @@ -52,7 +52,7 @@ openmc_legendre_filter_get_order(int32_t index, int* order) if (int err = verify_filter(index)) return err; // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. @@ -73,7 +73,7 @@ openmc_legendre_filter_set_order(int32_t index, int order) if (int err = verify_filter(index)) return err; // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp index 0dee964de..488f962fa 100644 --- a/src/tallies/filter_material.cpp +++ b/src/tallies/filter_material.cpp @@ -77,7 +77,7 @@ openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n) if (int err = verify_filter(index)) return err; // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. @@ -99,7 +99,7 @@ openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins) if (int err = verify_filter(index)) return err; // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index c17055a9e..324e38cc7 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -90,7 +90,7 @@ openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh) if (int err = verify_filter(index)) return err; // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. @@ -111,7 +111,7 @@ openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh) if (int err = verify_filter(index)) return err; // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp index c772d13ac..d75be0f9d 100644 --- a/src/tallies/filter_sph_harm.cpp +++ b/src/tallies/filter_sph_harm.cpp @@ -100,7 +100,7 @@ check_sphharm_filter(int32_t index) } // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. diff --git a/src/tallies/filter_sptl_legendre.cpp b/src/tallies/filter_sptl_legendre.cpp index 7757d1f69..b5d75f4d2 100644 --- a/src/tallies/filter_sptl_legendre.cpp +++ b/src/tallies/filter_sptl_legendre.cpp @@ -107,7 +107,7 @@ check_sptl_legendre_filter(int32_t index) } // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. diff --git a/src/tallies/filter_zernike.cpp b/src/tallies/filter_zernike.cpp index 467c64d65..3ea449238 100644 --- a/src/tallies/filter_zernike.cpp +++ b/src/tallies/filter_zernike.cpp @@ -132,7 +132,7 @@ check_zernike_filter(int32_t index) } // Get a pointer to the filter and downcast. - const auto& filt_base = model::tally_filters[index-1].get(); + const auto& filt_base = model::tally_filters[index].get(); auto* filt = dynamic_cast(filt_base); // Check the filter type. From 2729d986bf29456435c992fc98bf254184a33e0d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Feb 2019 21:09:17 -0600 Subject: [PATCH 7/9] Fix off-by-one for filter matches --- src/eigenvalue.cpp | 3 +- src/mesh.cpp | 29 ++++----- src/output.cpp | 8 +-- src/physics.cpp | 1 - src/state_point.cpp | 1 - src/tallies/derivative.cpp | 1 - src/tallies/filter_azimuthal.cpp | 6 +- src/tallies/filter_cell.cpp | 6 +- src/tallies/filter_cellborn.cpp | 6 +- src/tallies/filter_cellfrom.cpp | 6 +- src/tallies/filter_delayedgroup.cpp | 6 +- src/tallies/filter_distribcell.cpp | 6 +- src/tallies/filter_energy.cpp | 20 ++---- src/tallies/filter_energyfunc.cpp | 3 +- src/tallies/filter_legendre.cpp | 6 +- src/tallies/filter_material.cpp | 6 +- src/tallies/filter_mesh.cpp | 3 +- src/tallies/filter_meshsurface.cpp | 5 +- src/tallies/filter_mu.cpp | 6 +- src/tallies/filter_particle.cpp | 3 +- src/tallies/filter_polar.cpp | 6 +- src/tallies/filter_sph_harm.cpp | 9 ++- src/tallies/filter_sptl_legendre.cpp | 6 +- src/tallies/filter_surface.cpp | 6 +- src/tallies/filter_universe.cpp | 6 +- src/tallies/filter_zernike.cpp | 14 ++-- src/tallies/tally_scoring.cpp | 96 +++++++++++----------------- 27 files changed, 102 insertions(+), 172 deletions(-) diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index d0e493608..ccf70ee2e 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -608,8 +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 - // TODO: off by one - int mesh_bin = m->get_bin({p->coord[0].xyz}) - 1; + int mesh_bin = m->get_bin({p->coord[0].xyz}); if (mesh_bin < 0) { p->write_restart(); fatal_error("Source site outside UFS mesh!"); diff --git a/src/mesh.cpp b/src/mesh.cpp index b02a49ff4..ab263cd77 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -168,11 +168,11 @@ int RegularMesh::get_bin_from_indices(const int* ijk) const { switch (n_dimension_) { case 1: - return ijk[0]; + return ijk[0] - 1; case 2: - return (ijk[1] - 1)*shape_[0] + ijk[0]; + return (ijk[1] - 1)*shape_[0] + ijk[0] - 1; case 3: - return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0]; + return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0] - 1; default: throw std::runtime_error{"Invalid number of mesh dimensions"}; } @@ -193,14 +193,14 @@ void RegularMesh::get_indices(Position r, int* ijk, bool* in_mesh) const void RegularMesh::get_indices_from_bin(int bin, int* ijk) const { if (n_dimension_ == 1) { - ijk[0] = bin; + ijk[0] = bin + 1; } else if (n_dimension_ == 2) { - ijk[0] = (bin - 1) % shape_[0] + 1; - ijk[1] = (bin - 1) / shape_[0] + 1; + ijk[0] = bin % shape_[0] + 1; + ijk[1] = bin / shape_[0] + 1; } else if (n_dimension_ == 3) { - ijk[0] = (bin - 1) % shape_[0] + 1; - ijk[1] = ((bin - 1) % (shape_[0] * shape_[1])) / shape_[0] + 1; - ijk[2] = (bin - 1) / (shape_[0] * shape_[1]) + 1; + ijk[0] = bin % shape_[0] + 1; + ijk[1] = (bin % (shape_[0] * shape_[1])) / shape_[0] + 1; + ijk[2] = bin / (shape_[0] * shape_[1]) + 1; } } @@ -586,7 +586,7 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_)) { int i_surf = 4*i + 3; int i_mesh = get_bin_from_indices(ijk0.data()); - int i_bin = 4*n*(i_mesh - 1) + i_surf; + int i_bin = 4*n*i_mesh + i_surf - 1; bins.push_back(i_bin); } @@ -600,7 +600,7 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_)) { int i_surf = 4*i + 2; int i_mesh = get_bin_from_indices(ijk0.data()); - int i_bin = 4*n*(i_mesh - 1) + i_surf; + int i_bin = 4*n*i_mesh + i_surf - 1; bins.push_back(i_bin); } @@ -612,7 +612,7 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_) ){ int i_surf = 4*i + 1; int i_mesh = get_bin_from_indices(ijk0.data()); - int i_bin = 4*n*(i_mesh - 1) + i_surf; + int i_bin = 4*n*i_mesh + i_surf - 1; bins.push_back(i_bin); } @@ -626,7 +626,7 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_)) { int i_surf = 4*i + 4; int i_mesh = get_bin_from_indices(ijk0.data()); - int i_bin = 4*n*(i_mesh - 1) + i_surf; + int i_bin = 4*n*i_mesh + i_surf - 1; bins.push_back(i_bin); } @@ -670,8 +670,7 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, for (int64_t i = 0; i < n; ++i) { // determine scoring bin for entropy mesh - // TODO: off-by-one - int mesh_bin = get_bin({bank[i].xyz}) - 1; + int mesh_bin = get_bin({bank[i].xyz}); // if outside mesh, skip particle if (mesh_bin < 0) { diff --git a/src/output.cpp b/src/output.cpp index 97436dc24..095a17526 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -679,13 +679,12 @@ write_tallies() // prevents redundant output. int indent = 0; for (auto i = 0; i < tally.filters().size(); ++i) { - if ((filter_index-1) % tally.strides(i) == 0) { + if (filter_index % tally.strides(i) == 0) { auto i_filt = tally.filters(i); const auto& filt {*model::tally_filters[i_filt]}; auto& match {simulation::filter_matches[i_filt]}; tallies_out << std::string(indent+1, ' ') - // TODO: off-by-one - << filt.text_label(match.i_bin_+1) << "\n"; + << filt.text_label(match.i_bin_) << "\n"; } indent += 2; } @@ -712,9 +711,8 @@ write_tallies() std::string score_name = score > 0 ? reaction_name(score) : score_names.at(score); double mean, stdev; - //TODO: off-by-one std::tie(mean, stdev) = mean_stdev( - &tally.results_(filter_index-1, score_index, 0), tally.n_realizations_); + &tally.results_(filter_index, score_index, 0), tally.n_realizations_); tallies_out << std::string(indent+1, ' ') << std::left << std::setw(36) << score_name << " " << mean << " +/- " << t_value * stdev << "\n"; diff --git a/src/physics.cpp b/src/physics.cpp index 23ce8f560..5d6ca2dee 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -670,7 +670,6 @@ void scatter(Particle* p, int i_nuclide) p->event = EVENT_SCATTER; // Sample new outgoing angle for isotropic-in-lab scattering - // TODO: off-by-one const auto& mat {model::materials[p->material]}; if (!mat->p0_.empty()) { int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide]; diff --git a/src/state_point.cpp b/src/state_point.cpp index 3730c3216..cca622950 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -123,7 +123,6 @@ openmc_statepoint_write(const char* filename, bool* write_source) write_dataset(deriv_group, "independent variable", "density"); } else if (deriv.variable == DIFF_NUCLIDE_DENSITY) { write_dataset(deriv_group, "independent variable", "nuclide_density"); - //TODO: off-by-one write_dataset(deriv_group, "nuclide", data::nuclides[deriv.diff_nuclide]->name_); } else if (deriv.variable == DIFF_TEMPERATURE) { diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 1098adc55..2ad251b73 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -656,7 +656,6 @@ 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]}; - //TODO: off-by-one 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 diff --git a/src/tallies/filter_azimuthal.cpp b/src/tallies/filter_azimuthal.cpp index dc6528482..f4e265dd8 100644 --- a/src/tallies/filter_azimuthal.cpp +++ b/src/tallies/filter_azimuthal.cpp @@ -48,8 +48,7 @@ AzimuthalFilter::get_all_bins(const Particle* p, int estimator, } if (phi >= bins_.front() && phi <= bins_.back()) { - //TODO: off-by-one - auto bin = lower_bound_index(bins_.begin(), bins_.end(), phi) + 1; + auto bin = lower_bound_index(bins_.begin(), bins_.end(), phi); match.bins_.push_back(bin); match.weights_.push_back(1.0); } @@ -66,8 +65,7 @@ std::string AzimuthalFilter::text_label(int bin) const { std::stringstream out; - //TODO: off-by-one - out << "Azimuthal Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + out << "Azimuthal Angle [" << bins_[bin] << ", " << bins_[bin+1] << ")"; return out.str(); } diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp index 4e25e8541..1b2b05416 100644 --- a/src/tallies/filter_cell.cpp +++ b/src/tallies/filter_cell.cpp @@ -44,8 +44,7 @@ CellFilter::get_all_bins(const Particle* p, int estimator, for (int i = 0; i < p->n_coord; i++) { auto search = map_.find(p->coord[i].cell); if (search != map_.end()) { - //TODO: off-by-one - match.bins_.push_back(search->second + 1); + match.bins_.push_back(search->second); match.weights_.push_back(1.0); } } @@ -63,8 +62,7 @@ CellFilter::to_statepoint(hid_t filter_group) const std::string CellFilter::text_label(int bin) const { - //TODO: off-by-one - return "Cell " + std::to_string(model::cells[cells_[bin-1]]->id_); + return "Cell " + std::to_string(model::cells[cells_[bin]]->id_); } //============================================================================== diff --git a/src/tallies/filter_cellborn.cpp b/src/tallies/filter_cellborn.cpp index b173b2e39..005d3844b 100644 --- a/src/tallies/filter_cellborn.cpp +++ b/src/tallies/filter_cellborn.cpp @@ -10,8 +10,7 @@ CellbornFilter::get_all_bins(const Particle* p, int estimator, { auto search = map_.find(p->cell_born); if (search != map_.end()) { - //TODO: off-by-one - match.bins_.push_back(search->second + 1); + match.bins_.push_back(search->second); match.weights_.push_back(1.0); } } @@ -19,8 +18,7 @@ CellbornFilter::get_all_bins(const Particle* p, int estimator, std::string CellbornFilter::text_label(int bin) const { - //TODO: off-by-one - return "Birth Cell " + std::to_string(model::cells[cells_[bin-1]]->id_); + return "Birth Cell " + std::to_string(model::cells[cells_[bin]]->id_); } } // namespace openmc diff --git a/src/tallies/filter_cellfrom.cpp b/src/tallies/filter_cellfrom.cpp index 1262ebecd..9a43f7c43 100644 --- a/src/tallies/filter_cellfrom.cpp +++ b/src/tallies/filter_cellfrom.cpp @@ -11,8 +11,7 @@ CellFromFilter::get_all_bins(const Particle* p, int estimator, for (int i = 0; i < p->last_n_coord; i++) { auto search = map_.find(p->last_cell[i]); if (search != map_.end()) { - //TODO: off-by-one - match.bins_.push_back(search->second + 1); + match.bins_.push_back(search->second); match.weights_.push_back(1.0); } } @@ -21,8 +20,7 @@ CellFromFilter::get_all_bins(const Particle* p, int estimator, std::string CellFromFilter::text_label(int bin) const { - //TODO: off-by-one - return "Cell from " + std::to_string(model::cells[cells_[bin-1]]->id_); + return "Cell from " + std::to_string(model::cells[cells_[bin]]->id_); } } // namespace openmc diff --git a/src/tallies/filter_delayedgroup.cpp b/src/tallies/filter_delayedgroup.cpp index f1876fae9..4f448358e 100644 --- a/src/tallies/filter_delayedgroup.cpp +++ b/src/tallies/filter_delayedgroup.cpp @@ -29,8 +29,7 @@ void DelayedGroupFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { - //TODO: off-by-one - match.bins_.push_back(1); + match.bins_.push_back(0); match.weights_.push_back(1.0); } @@ -44,8 +43,7 @@ DelayedGroupFilter::to_statepoint(hid_t filter_group) const std::string DelayedGroupFilter::text_label(int bin) const { - //TODO: off-by-one - return "Delayed Group " + std::to_string(groups_[bin-1]); + return "Delayed Group " + std::to_string(groups_[bin]); } } // namespace openmc diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp index ee4f08dd4..ab7ea67cf 100644 --- a/src/tallies/filter_distribcell.cpp +++ b/src/tallies/filter_distribcell.cpp @@ -54,8 +54,7 @@ DistribcellFilter::get_all_bins(const Particle* p, int estimator, } } if (cell_ == p->coord[i].cell) { - //TODO: off-by-one - match.bins_.push_back(offset + 1); + match.bins_.push_back(offset); match.weights_.push_back(1.0); return; } @@ -73,8 +72,7 @@ std::string DistribcellFilter::text_label(int bin) const { auto map = model::cells[cell_]->distribcell_index_; - //TODO: off-by-one - auto path = distribcell_path(cell_, map, bin-1); + auto path = distribcell_path(cell_, map, bin); return "Distributed Cell " + path; } diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp index 90cd7b3bf..4cbde1cc4 100644 --- a/src/tallies/filter_energy.cpp +++ b/src/tallies/filter_energy.cpp @@ -43,11 +43,9 @@ const { if (p->g != F90_NONE && matches_transport_groups_) { if (estimator == ESTIMATOR_TRACKLENGTH) { - //TODO: off-by-one - match.bins_.push_back(data::num_energy_groups - p->g + 1); + match.bins_.push_back(data::num_energy_groups - p->g); } else { - //TODO: off-by-one - match.bins_.push_back(data::num_energy_groups - p->last_g + 1); + match.bins_.push_back(data::num_energy_groups - p->last_g); } match.weights_.push_back(1.0); @@ -57,8 +55,7 @@ const // Bin the energy. if (E >= bins_.front() && E <= bins_.back()) { - //TODO: off-by-one - auto bin = lower_bound_index(bins_.begin(), bins_.end(), E) + 1; + auto bin = lower_bound_index(bins_.begin(), bins_.end(), E); match.bins_.push_back(bin); match.weights_.push_back(1.0); } @@ -76,8 +73,7 @@ std::string EnergyFilter::text_label(int bin) const { std::stringstream out; - //TODO: off-by-one - out << "Incoming Energy [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + out << "Incoming Energy [" << bins_[bin] << ", " << bins_[bin+1] << ")"; return out.str(); } @@ -90,13 +86,12 @@ 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 + 1); + 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()) { - //TODO: off-by-one - auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->E) + 1; + auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->E); match.bins_.push_back(bin); match.weights_.push_back(1.0); } @@ -107,8 +102,7 @@ std::string EnergyoutFilter::text_label(int bin) const { std::stringstream out; - //TODO: off-by-one - out << "Outgoing Energy [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + out << "Outgoing Energy [" << bins_[bin] << ", " << bins_[bin+1] << ")"; return out.str(); } diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp index b47e7b7b6..3c616aabf 100644 --- a/src/tallies/filter_energyfunc.cpp +++ b/src/tallies/filter_energyfunc.cpp @@ -41,8 +41,7 @@ EnergyFunctionFilter::get_all_bins(const Particle* p, int estimator, double f = (p->last_E - energy_[i]) / (energy_[i+1] - energy_[i]); // Interpolate on the lin-lin grid. - //TODO: off-by-one - match.bins_.push_back(1); + match.bins_.push_back(0); match.weights_.push_back((1-f) * y_[i] + f * y_[i+1]); } } diff --git a/src/tallies/filter_legendre.cpp b/src/tallies/filter_legendre.cpp index a5f583a15..1e799932a 100644 --- a/src/tallies/filter_legendre.cpp +++ b/src/tallies/filter_legendre.cpp @@ -21,8 +21,7 @@ LegendreFilter::get_all_bins(const Particle* p, int estimator, double wgt[n_bins_]; calc_pn_c(order_, p->mu, wgt); for (int i = 0; i < n_bins_; i++) { - //TODO: off-by-one - match.bins_.push_back(i + 1); + match.bins_.push_back(i); match.weights_.push_back(wgt[i]); } } @@ -37,8 +36,7 @@ LegendreFilter::to_statepoint(hid_t filter_group) const std::string LegendreFilter::text_label(int bin) const { - //TODO: off-by-one - return "Legendre expansion, P" + std::to_string(bin - 1); + return "Legendre expansion, P" + std::to_string(bin); } //============================================================================== diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp index 488f962fa..a4da2deb1 100644 --- a/src/tallies/filter_material.cpp +++ b/src/tallies/filter_material.cpp @@ -44,8 +44,7 @@ MaterialFilter::get_all_bins(const Particle* p, int estimator, { auto search = map_.find(p->material); if (search != map_.end()) { - //TODO: off-by-one - match.bins_.push_back(search->second + 1); + match.bins_.push_back(search->second); match.weights_.push_back(1.0); } } @@ -62,8 +61,7 @@ MaterialFilter::to_statepoint(hid_t filter_group) const std::string MaterialFilter::text_label(int bin) const { - //TODO: off-by-one - return "Material " + std::to_string(model::materials[materials_[bin-1]]->id_); + return "Material " + std::to_string(model::materials[materials_[bin]]->id_); } //============================================================================== diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index 324e38cc7..177dc407e 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -25,7 +25,7 @@ MeshFilter::from_xml(pugi::xml_node node) set_mesh(search->second); } else{ std::stringstream err_msg; - err_msg << "Could not find cell " << id << " specified on tally filter."; + err_msg << "Could not find mesh " << id << " specified on tally filter."; fatal_error(err_msg); } } @@ -37,7 +37,6 @@ const if (estimator != ESTIMATOR_TRACKLENGTH) { auto bin = model::meshes[mesh_]->get_bin(p->coord[0].xyz); if (bin >= 0) { - //TODO: off-by-one match.bins_.push_back(bin); match.weights_.push_back(1.0); } diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp index a8854c7ec..57f10b09c 100644 --- a/src/tallies/filter_meshsurface.cpp +++ b/src/tallies/filter_meshsurface.cpp @@ -22,9 +22,8 @@ MeshSurfaceFilter::text_label(int bin) const int n_dim = mesh.n_dimension_; // Get flattend mesh index and surface index. - //TODO: off-by-one - int i_mesh = (bin - 1) / (4 * n_dim) + 1; - int i_surf = ((bin - 1) % (4 * n_dim)) + 1; + int i_mesh = bin / (4 * n_dim); + int i_surf = (bin % (4 * n_dim)) + 1; // Get mesh index part of label. std::string out = MeshFilter::text_label(i_mesh); diff --git a/src/tallies/filter_mu.cpp b/src/tallies/filter_mu.cpp index 521deefe3..78a648e23 100644 --- a/src/tallies/filter_mu.cpp +++ b/src/tallies/filter_mu.cpp @@ -39,8 +39,7 @@ MuFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) const { if (p->mu >= bins_.front() && p->mu <= bins_.back()) { - //TODO: off-by-one - auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->mu) + 1; + auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->mu); match.bins_.push_back(bin); match.weights_.push_back(1.0); } @@ -57,8 +56,7 @@ std::string MuFilter::text_label(int bin) const { std::stringstream out; - //TODO: off-by-one - out << "Change-in-Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + out << "Change-in-Angle [" << bins_[bin] << ", " << bins_[bin+1] << ")"; return out.str(); } diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index 403ef9500..967febc06 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -18,8 +18,7 @@ ParticleFilter::get_all_bins(const Particle* p, int estimator, { for (auto i = 0; i < particles_.size(); i++) { if (particles_[i] == p->type) { - //TODO: off-by-one - match.bins_.push_back(i + 1); + match.bins_.push_back(i); match.weights_.push_back(1.0); } } diff --git a/src/tallies/filter_polar.cpp b/src/tallies/filter_polar.cpp index 50fc47999..5216511a9 100644 --- a/src/tallies/filter_polar.cpp +++ b/src/tallies/filter_polar.cpp @@ -47,8 +47,7 @@ const } if (theta >= bins_.front() && theta <= bins_.back()) { - //TODO: off-by-one - auto bin = lower_bound_index(bins_.begin(), bins_.end(), theta) + 1; + auto bin = lower_bound_index(bins_.begin(), bins_.end(), theta); match.bins_.push_back(bin); match.weights_.push_back(1.0); } @@ -65,8 +64,7 @@ std::string PolarFilter::text_label(int bin) const { std::stringstream out; - //TODO: off-by-one - out << "Polar Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + out << "Polar Angle [" << bins_[bin] << ", " << bins_[bin+1] << ")"; return out.str(); } diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp index d75be0f9d..253bf7ff3 100644 --- a/src/tallies/filter_sph_harm.cpp +++ b/src/tallies/filter_sph_harm.cpp @@ -54,8 +54,8 @@ SphericalHarmonicsFilter::get_all_bins(const Particle* p, int estimator, // Append the matching (bin,weight) for each moment for (int i = 0; i < num_nm; i++) { match.weights_.push_back(wgt[n] * rn[j]); - //TODO: off-by-one - match.bins_.push_back(++j); + match.bins_.push_back(j); + ++j; } } } @@ -77,9 +77,8 @@ SphericalHarmonicsFilter::text_label(int bin) const { std::stringstream out; for (int n = 0; n < order_ + 1; n++) { - if (bin <= (n + 1) * (n + 1)) { - //TODO: off-by-one - int m = (bin - n*n - 1) - n; + if (bin < (n + 1) * (n + 1)) { + int m = (bin - n*n) - n; out << "Spherical harmonic expansion, Y" << n << "," << m; return out.str(); } diff --git a/src/tallies/filter_sptl_legendre.cpp b/src/tallies/filter_sptl_legendre.cpp index b5d75f4d2..dae62d1d2 100644 --- a/src/tallies/filter_sptl_legendre.cpp +++ b/src/tallies/filter_sptl_legendre.cpp @@ -53,8 +53,7 @@ SpatialLegendreFilter::get_all_bins(const Particle* p, int estimator, double wgt[order_ + 1]; calc_pn_c(order_, x_norm, wgt); for (int i = 0; i < order_ + 1; i++) { - //TODO: off-by-one - match.bins_.push_back(i + 1); + match.bins_.push_back(i); match.weights_.push_back(wgt[i]); } } @@ -88,8 +87,7 @@ SpatialLegendreFilter::text_label(int bin) const } else { out << "z"; } - //TODO: off-by-one - out << " axis, P" << std::to_string(bin - 1); + out << " axis, P" << std::to_string(bin); return out.str(); } diff --git a/src/tallies/filter_surface.cpp b/src/tallies/filter_surface.cpp index 5b175ca70..774571d75 100644 --- a/src/tallies/filter_surface.cpp +++ b/src/tallies/filter_surface.cpp @@ -43,8 +43,7 @@ SurfaceFilter::get_all_bins(const Particle* p, int estimator, { auto search = map_.find(std::abs(p->surface)-1); if (search != map_.end()) { - //TODO: off-by-one - match.bins_.push_back(search->second + 1); + match.bins_.push_back(search->second); if (p->surface < 0) { match.weights_.push_back(-1.0); } else { @@ -65,8 +64,7 @@ SurfaceFilter::to_statepoint(hid_t filter_group) const std::string SurfaceFilter::text_label(int bin) const { - //TODO: off-by-one - return "Surface " + std::to_string(model::surfaces[surfaces_[bin-1]]->id_); + return "Surface " + std::to_string(model::surfaces[surfaces_[bin]]->id_); } } // namespace openmc diff --git a/src/tallies/filter_universe.cpp b/src/tallies/filter_universe.cpp index ecd12b237..5a6d5e1eb 100644 --- a/src/tallies/filter_universe.cpp +++ b/src/tallies/filter_universe.cpp @@ -44,8 +44,7 @@ UniverseFilter::get_all_bins(const Particle* p, int estimator, for (int i = 0; i < p->n_coord; i++) { auto search = map_.find(p->coord[i].universe); if (search != map_.end()) { - //TODO: off-by-one - match.bins_.push_back(search->second + 1); + match.bins_.push_back(search->second); match.weights_.push_back(1.0); } } @@ -63,8 +62,7 @@ UniverseFilter::to_statepoint(hid_t filter_group) const std::string UniverseFilter::text_label(int bin) const { - //TODO: off-by-one - return "Universe " + std::to_string(model::universes[universes_[bin-1]]->id_); + return "Universe " + std::to_string(model::universes[universes_[bin]]->id_); } } // namespace openmc diff --git a/src/tallies/filter_zernike.cpp b/src/tallies/filter_zernike.cpp index 3ea449238..a005f7ddc 100644 --- a/src/tallies/filter_zernike.cpp +++ b/src/tallies/filter_zernike.cpp @@ -39,8 +39,7 @@ ZernikeFilter::get_all_bins(const Particle* p, int estimator, double zn[n_bins_]; calc_zn(order_, r, theta, zn); for (int i = 0; i < n_bins_; i++) { - //TODO: off-by-one - match.bins_.push_back(i+1); + match.bins_.push_back(i); match.weights_.push_back(zn[i]); } } @@ -62,9 +61,8 @@ ZernikeFilter::text_label(int bin) const std::stringstream out; for (int n = 0; n < order_+1; n++) { int last = (n + 1) * (n + 2) / 2; - //TODO: off-by-one - if (bin <= last) { - int first = last - n; + if (bin < last) { + int first = last - (n + 1); int m = -n + (bin - first) * 2; out << "Zernike expansion, Z" << n << "," << m; return out.str(); @@ -97,8 +95,7 @@ ZernikeRadialFilter::get_all_bins(const Particle* p, int estimator, double zn[n_bins_]; calc_zn_rad(order_, r, zn); for (int i = 0; i < n_bins_; i++) { - //TODO: off-by-one - match.bins_.push_back(i+1); + match.bins_.push_back(i); match.weights_.push_back(zn[i]); } } @@ -107,8 +104,7 @@ ZernikeRadialFilter::get_all_bins(const Particle* p, int estimator, std::string ZernikeRadialFilter::text_label(int bin) const { - //TODO: off-by-one - return "Zernike expansion, Z" + std::to_string(2*(bin-1)) + ",0"; + return "Zernike expansion, Z" + std::to_string(2*bin) + ",0"; } void diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 96aa1abaa..5882d241a 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -49,7 +49,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, const Particle* p) } // Compute the initial index and weight. - compute_index_weight(); + this->compute_index_weight(); } FilterBinIter::FilterBinIter(const Tally& tally, bool end) @@ -67,8 +67,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, bool end) match.bins_.clear(); match.weights_.clear(); for (auto i = 0; i < model::tally_filters[i_filt]->n_bins_; ++i) { - // TODO: off-by-one - match.bins_.push_back(i+1); + match.bins_.push_back(i); match.weights_.push_back(1.0); } match.bins_present_ = true; @@ -83,7 +82,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, bool end) } // Compute the initial index and weight. - compute_index_weight(); + this->compute_index_weight(); } FilterBinIter& @@ -124,14 +123,13 @@ FilterBinIter::operator++() void FilterBinIter::compute_index_weight() { - index_ = 1; + index_ = 0; weight_ = 1.; for (auto i = 0; i < tally_.filters().size(); ++i) { auto i_filt = tally_.filters(i); auto& match {simulation::filter_matches[i_filt]}; auto i_bin = match.i_bin_; - //TODO: off-by-one - index_ += (match.bins_[i_bin] - 1) * tally_.strides(i); + index_ += match.bins_[i_bin] * tally_.strides(i); weight_ *= match.weights_[i_bin]; } } @@ -154,19 +152,17 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index) dg_match.bins_[i_bin] = d_bin; // Determine the filter scoring index - auto filter_index = 1; + auto filter_index = 0; for (auto i = 0; i < tally.filters().size(); ++i) { auto i_filt = tally.filters(i); auto& match {simulation::filter_matches[i_filt]}; auto i_bin = match.i_bin_; - //TODO: off-by-one - filter_index += (match.bins_[i_bin] - 1) * tally.strides(i); + filter_index += match.bins_[i_bin] * tally.strides(i); } // Update the tally result - // TODO: off-by-one #pragma omp atomic - tally.results_(filter_index-1, score_index, RESULT_VALUE) += score; + tally.results_(filter_index, score_index, RESULT_VALUE) += score; // Reset the original delayed group bin dg_match.bins_[i_bin] = original_bin; @@ -217,7 +213,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) // modify the value so that g_out = 1 corresponds to the highest energy // bin - g_out = eo_filt.n_bins_ - g_out + 1; + g_out = eo_filt.n_bins_ - g_out; // change outgoing energy bin simulation::filter_matches[i_eout_filt].bins_[i_bin] = g_out; @@ -235,9 +231,8 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) if (E_out < eo_filt.bins_.front() || E_out > eo_filt.bins_.back()) { continue; } else { - //TODO: off-by-one auto i_match = lower_bound_index(eo_filt.bins_.begin(), - eo_filt.bins_.end(), E_out) + 1; + eo_filt.bins_.end(), E_out); simulation::filter_matches[i_eout_filt].bins_[i_bin] = i_match; } @@ -249,17 +244,17 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) // Find the filter scoring index for this filter combination //TODO: should this include a weight? - int filter_index = 1; + int filter_index = 0; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); auto& match {simulation::filter_matches[i_filt]}; auto i_bin = match.i_bin_; - filter_index += (match.bins_[i_bin] - 1) * tally.strides(j); + filter_index += match.bins_[i_bin] * tally.strides(j); } // Update tally results #pragma omp atomic - tally.results_(filter_index-1, i_score, RESULT_VALUE) += score; + tally.results_(filter_index, i_score, RESULT_VALUE) += score; } else if (score_bin == SCORE_DELAYED_NU_FISSION && g != 0) { @@ -276,17 +271,15 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) for (auto d_bin = 0; d_bin < dg_filt.n_bins_; ++d_bin) { if (dg_filt.groups_[d_bin] == g) { // Find the filter index and weight for this filter combination - int filter_index = 1; double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); auto& match {simulation::filter_matches[i_filt]}; auto i_bin = match.i_bin_; - filter_index += (match.bins_[i_bin] - 1) * tally.strides(j); filter_weight *= match.weights_[i_bin]; } - score_fission_delayed_dg(i_tally, d_bin+1, score*filter_weight, + score_fission_delayed_dg(i_tally, d_bin, score*filter_weight, i_score); } } @@ -295,19 +288,19 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) } else { // Find the filter index and weight for this filter combination - int filter_index = 1; + int filter_index = 0; double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); auto& match {simulation::filter_matches[i_filt]}; auto i_bin = match.i_bin_; - filter_index += (match.bins_[i_bin] - 1) * tally.strides(j); + filter_index += match.bins_[i_bin] * tally.strides(j); filter_weight *= match.weights_[i_bin]; } // Update tally results #pragma omp atomic - tally.results_(filter_index-1, i_score, RESULT_VALUE) += score*filter_weight; + tally.results_(filter_index, i_score, RESULT_VALUE) += score*filter_weight; } } } @@ -642,7 +635,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, score = p->absorb_wgt * yield * simulation::micro_xs[p->event_nuclide].fission / simulation::micro_xs[p->event_nuclide].absorption * flux; - score_fission_delayed_dg(i_tally, d_bin+1, score, + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; @@ -677,7 +670,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto d = filt.groups_[d_bin]; score = simulation::keff * p->wgt_bank / p->n_bank * p->n_delayed_bank[d-1] * flux; - score_fission_delayed_dg(i_tally, d_bin+1, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; } else { @@ -702,7 +695,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = simulation::micro_xs[i_nuclide].fission * yield * atom_density * flux; - score_fission_delayed_dg(i_tally, d_bin+1, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; } else { @@ -731,7 +724,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = simulation::micro_xs[j_nuclide].fission * yield * atom_density * flux; - score_fission_delayed_dg(i_tally, d_bin+1, score, + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } } @@ -782,7 +775,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, * simulation::micro_xs[p->event_nuclide].fission / simulation::micro_xs[p->event_nuclide].absorption * rate * flux; - score_fission_delayed_dg(i_tally, d_bin+1, score, + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; @@ -836,7 +829,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { auto d = filt.groups_[d_bin]; if (d == g) - score_fission_delayed_dg(i_tally, d_bin+1, score, + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } score = 0.; @@ -862,7 +855,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto rate = rxn.products_[d].decay_rate_; score = simulation::micro_xs[i_nuclide].fission * yield * flux * atom_density * rate; - score_fission_delayed_dg(i_tally, d_bin+1, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; } else { @@ -902,7 +895,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, auto rate = rxn.products_[d].decay_rate_; score = simulation::micro_xs[j_nuclide].fission * yield * flux * atom_density * rate; - score_fission_delayed_dg(i_tally, d_bin+1, score, + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } } @@ -1214,7 +1207,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index, // Update tally results #pragma omp atomic - tally.results_(filter_index-1, score_index, RESULT_VALUE) += score; + tally.results_(filter_index, score_index, RESULT_VALUE) += score; } } @@ -1229,8 +1222,6 @@ score_general_mg(const Particle* p, int i_tally, int start_index, { auto& tally {*model::tallies[i_tally]}; - //TODO: off-by-one throughout on p->material - // Set the direction and group to use with get_xs const double* p_uvw; int p_g; @@ -1626,7 +1617,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, p_g, nullptr, nullptr, &d) / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); } - score_fission_delayed_dg(i_tally, d_bin+1, score, + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; @@ -1671,7 +1662,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); } - score_fission_delayed_dg(i_tally, d_bin+1, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; } else { @@ -1705,7 +1696,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin+1, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; } else { @@ -1752,7 +1743,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, p_g, nullptr, nullptr, &d) / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); } - score_fission_delayed_dg(i_tally, d_bin+1, score, + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; @@ -1817,7 +1808,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) { auto d = filt.groups_[d_bin]; if (d == g) - score_fission_delayed_dg(i_tally, d_bin+1, score, + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } score = 0.; @@ -1848,7 +1839,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin+1, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index); } continue; } else { @@ -1921,7 +1912,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // Update tally results #pragma omp atomic - tally.results_(filter_index-1, score_index, RESULT_VALUE) += score; + tally.results_(filter_index, score_index, RESULT_VALUE) += score; } } @@ -2044,12 +2035,8 @@ void score_analog_tally_mg(const Particle* p) double atom_density = 0.; if (i_nuclide >= 0) { - //TODO: off-by-one - auto j = model::materials[p->material] - ->mat_nuclide_index_[i_nuclide]; + auto j = model::materials[p->material]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - //atom_density = material_atom_density(p->material, j); - //TODO: off-by-one atom_density = model::materials[p->material]->atom_density_(j); } @@ -2103,12 +2090,8 @@ score_tracklength_tally(const Particle* p, double distance) double atom_density = 0.; if (i_nuclide >= 0) { if (p->material != MATERIAL_VOID) { - //TODO: off-by-one - auto j = model::materials[p->material] - ->mat_nuclide_index_[i_nuclide]; + auto j = model::materials[p->material]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - //atom_density = material_atom_density(p->material, j); - //TODO: off-by-one atom_density = model::materials[p->material]->atom_density_(j); } } @@ -2173,12 +2156,8 @@ void score_collision_tally(const Particle* p) double atom_density = 0.; if (i_nuclide >= 0) { - //TODO: off-by-one - auto j = model::materials[p->material] - ->mat_nuclide_index_[i_nuclide]; + auto j = model::materials[p->material]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - //atom_density = material_atom_density(p->material, j); - //TODO: off-by-one atom_density = model::materials[p->material]->atom_density_(j); } @@ -2233,9 +2212,8 @@ score_surface_tally(const Particle* p, const std::vector& tallies) double score = flux * filter_weight; for (auto score_index = 0; score_index < tally.scores_.size(); ++score_index) { - //TODO: off-by-one #pragma omp atomic - tally.results_(filter_index-1, score_index, RESULT_VALUE) += score; + tally.results_(filter_index, score_index, RESULT_VALUE) += score; } } From 964664893758db912e1e779a06cd1381c93cff7f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Feb 2019 22:48:57 -0600 Subject: [PATCH 8/9] Fix off-by-one for get_nuclide_xs --- src/mgxs_interface.cpp | 2 +- src/tallies/tally_scoring.cpp | 72 +++++++++++++++++------------------ 2 files changed, 37 insertions(+), 37 deletions(-) diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 581e6b178..d2913a17a 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -277,7 +277,7 @@ get_nuclide_xs(int index, int xstype, int gin, const int* gout, } else { dg_c_p = dg; } - return data::nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p); + return data::nuclides_MG[index].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p); } //============================================================================== diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 5882d241a..5288ea597 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1311,12 +1311,12 @@ score_general_mg(const Particle* p, int i_tally, int start_index, //TODO: should flux be multiplied in above instead of below? if (i_nuclide >= 0) { score *= flux * atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_TOTAL, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_TOTAL, p_g) / get_macro_xs(p->material, MG_GET_XS_TOTAL, p_g); } } else { if (i_nuclide >= 0) { - score = get_nuclide_xs(i_nuclide+1, MG_GET_XS_TOTAL, p_g) + score = get_nuclide_xs(i_nuclide, MG_GET_XS_TOTAL, p_g) * atom_density * flux; } else { score = simulation::material_xs.total * flux; @@ -1340,7 +1340,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } if (i_nuclide >= 0) { score *= flux - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_INVERSE_VELOCITY, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_INVERSE_VELOCITY, p_g) / get_macro_xs(p->material, MG_GET_XS_TOTAL, p_g); } else { score *= flux @@ -1350,7 +1350,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } else { if (i_nuclide >= 0) { score = flux - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_INVERSE_VELOCITY, p_g); + * get_nuclide_xs(i_nuclide, MG_GET_XS_INVERSE_VELOCITY, p_g); } else { score = flux * get_macro_xs(p->material, MG_GET_XS_INVERSE_VELOCITY, p_g); @@ -1368,7 +1368,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = p->last_wgt * flux; if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_SCATTER_FMU_MULT, + * get_nuclide_xs(i_nuclide, MG_GET_XS_SCATTER_FMU_MULT, p->last_g, &p->g, &p->mu, nullptr) / get_macro_xs(p->material, MG_GET_XS_SCATTER_FMU_MULT, p->last_g, &p->g, &p->mu, nullptr); @@ -1376,7 +1376,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } else { if (i_nuclide >= 0) { score = atom_density * flux * get_nuclide_xs( - i_nuclide+1, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu, nullptr); + i_nuclide, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu, nullptr); } else { score = flux * get_macro_xs( p->material, MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu, nullptr); @@ -1398,7 +1398,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, // adjust the score by the actual probability for that nuclide. if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_SCATTER_FMU, + * get_nuclide_xs(i_nuclide, MG_GET_XS_SCATTER_FMU, p->last_g, &p->g, &p->mu, nullptr) / get_macro_xs(p->material, MG_GET_XS_SCATTER_FMU, p->last_g, &p->g, &p->mu, nullptr); @@ -1406,7 +1406,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } else { if (i_nuclide >= 0) { score = atom_density * flux * get_nuclide_xs( - i_nuclide+1, MG_GET_XS_SCATTER, p_g); + i_nuclide, MG_GET_XS_SCATTER, p_g); } else { score = flux * get_macro_xs(p->material, MG_GET_XS_SCATTER, p_g); } @@ -1429,13 +1429,13 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_ABSORPTION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_ABSORPTION, p_g) / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); } } else { if (i_nuclide >= 0) { score = atom_density * flux - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_ABSORPTION, p_g); + * get_nuclide_xs(i_nuclide, MG_GET_XS_ABSORPTION, p_g); } else { score = simulation::material_xs.absorption * flux; } @@ -1460,7 +1460,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); } else { score *= @@ -1469,7 +1469,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } } else { if (i_nuclide >= 0) { - score = get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + score = get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) * atom_density * flux; } else { score = get_macro_xs(p->material, MG_GET_XS_FISSION, p_g) * flux; @@ -1495,7 +1495,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = p->absorb_wgt * flux; if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_NU_FISSION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_NU_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); } else { score *= @@ -1513,13 +1513,13 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = simulation::keff * p->wgt_bank * flux; if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); } } } else { if (i_nuclide >= 0) { - score = get_nuclide_xs(i_nuclide+1, MG_GET_XS_NU_FISSION, p_g) + score = get_nuclide_xs(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * atom_density * flux; } else { score = get_macro_xs(p->material, MG_GET_XS_NU_FISSION, p_g) * flux; @@ -1545,7 +1545,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = p->absorb_wgt * flux; if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_PROMPT_NU_FISSION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); } else { score *= @@ -1566,13 +1566,13 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = simulation::keff * p->wgt_bank * prompt_frac * flux; if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); } } } else { if (i_nuclide >= 0) { - score = get_nuclide_xs(i_nuclide+1, MG_GET_XS_PROMPT_NU_FISSION, p_g) + score = get_nuclide_xs(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * atom_density * flux; } else { score = get_macro_xs(p->material, MG_GET_XS_PROMPT_NU_FISSION, p_g) @@ -1608,7 +1608,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = p->absorb_wgt * flux; if (i_nuclide >= 0) { score *= - get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + 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 { @@ -1628,7 +1628,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = p->absorb_wgt * flux; if (i_nuclide >= 0) { score *= - get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g) + get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); } else { score *= @@ -1659,7 +1659,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * p->n_delayed_bank[d-1] * flux; if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); } score_fission_delayed_dg(i_tally, d_bin, score, score_index); @@ -1673,7 +1673,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * flux; if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); } } @@ -1689,7 +1689,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, auto d = filt.groups_[d_bin]; if (i_nuclide >= 0) { score = flux * atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } else { score = flux @@ -1702,7 +1702,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } else { if (i_nuclide >= 0) { score = flux * atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, p_g); + * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g); } else { score = flux * get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION, p_g); @@ -1730,9 +1730,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = p->absorb_wgt * flux; if (i_nuclide >= 0) { score *= - get_nuclide_xs(i_nuclide+1, MG_GET_XS_DECAY_RATE, + get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + * 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 { @@ -1756,9 +1756,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, for (auto d = 0; d < data::num_delayed_groups; ++d) { if (i_nuclide >= 0) { score += p->absorb_wgt * flux - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DECAY_RATE, + * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + * 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 { @@ -1790,9 +1790,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, if (g != 0) { if (i_nuclide >= 0) { score += simulation::keff * atom_density * bank.wgt * flux - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DECAY_RATE, p_g, + * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &g) - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_FISSION, p_g); } else { score += simulation::keff * bank.wgt * flux @@ -1828,9 +1828,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, auto d = filt.groups_[d_bin]; if (i_nuclide >= 0) { score += atom_density * flux - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DECAY_RATE, + * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } else { score += flux @@ -1847,9 +1847,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index, for (auto d = 0; d < data::num_delayed_groups; ++d) { if (i_nuclide >= 0) { score += atom_density * flux - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DECAY_RATE, + * get_nuclide_xs(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d) - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_DELAYED_NU_FISSION, + * get_nuclide_xs(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } else { score += flux @@ -1881,7 +1881,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } if (i_nuclide >= 0) { score *= atom_density - * get_nuclide_xs(i_nuclide+1, MG_GET_XS_KAPPA_FISSION, p_g) + * get_nuclide_xs(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) / get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g); } else { score *= @@ -1890,7 +1890,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, } } else { if (i_nuclide >= 0) { - score = get_nuclide_xs(i_nuclide+1, MG_GET_XS_KAPPA_FISSION, p_g) + score = get_nuclide_xs(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * atom_density * flux; } else { score = get_macro_xs(p->material, MG_GET_XS_KAPPA_FISSION, p_g) From 6ac86503efe915f57cf0801648680f717bc4b97d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 26 Feb 2019 08:51:48 -0600 Subject: [PATCH 9/9] Add one @wbinventor publication in list --- docs/source/publications.rst | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/docs/source/publications.rst b/docs/source/publications.rst index 012560a28..990a8fc21 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -196,9 +196,14 @@ Miscellaneous Reaction Rate Kernel Density Estimators in OpenMC," *Trans. Am. Nucl. Soc.*, **109**, 683-686 (2013). ------------------------------------- -Multi-group Cross Section Generation ------------------------------------- +----------------------------------- +Multigroup Cross Section Generation +----------------------------------- + +- William Boyd, Benoit Forget, and Kord Smith, "`A single-step framework to + generate spatially self-shielded multi-group cross sections from Monte Carlo + transport simulations `_," + *Ann. Nucl. Energy*, **125**, 261-271 (2019). - Changho Lee and Yeon Sang Jung, "Verification of the Cross Section Library Generated Using OpenMC and MC\ :sup:`2`-3 for PROTEUS," *Proc. PHYSOR*, Cancun,