crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1179 from paulromano/off-by-one

Browse source 
Resolve most off-by-ones in C++ code
This commit is contained in:
Sterling Harper 2019-02-28 13:54:30 -05:00 • committed by GitHub
commit 241651bb28
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
47 changed files with 394 additions and 521 deletions
Download patch file Download diff file Expand all files Collapse all files

11
docs/source/publications.rst
View file

@ -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 <https://doi.org/10.1016/j.anucene.2018.11.017>`_,"
*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,

7
openmc/capi/cell.py
View file

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

9
openmc/capi/filter.py
View file

@ -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)
@ -405,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()

2
openmc/capi/material.py
View file

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

6
openmc/capi/tally.py
View file

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

5
openmc/material.py
View file

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

4
src/bremsstrahlung.cpp
View file

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

5
src/cell.cpp
View file

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

3
src/eigenvalue.cpp
View file

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

10
src/geometry.cpp
View file

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

15
src/geometry_aux.cpp
View file

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

46
src/material.cpp
View file

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

29
src/mesh.cpp
View file

@ -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<int>& 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<int>& 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<int>& 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<int>& 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<double> 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) {

6
src/mgxs_interface.cpp
View file

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

13
src/nuclide.cpp
View file

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

11
src/output.cpp
View file

@ -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 "
@ -680,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;
}
@ -713,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";

12
src/particle.cpp
View file

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

27
src/physics.cpp
View file

@ -63,10 +63,10 @@ void collision(Particle* p)
std::stringstream msg;
if (static_cast<ParticleType>(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]};
@ -425,8 +423,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 +448,7 @@ int sample_element(Particle* p)
double cutoff = prn() * simulation::material_xs.total;
// Get pointers to elements, densities
const auto& mat {model::materials[p->material - 1]};
const auto& mat {model::materials[p->material]};
int n = mat->nuclide_.size();
int i = 0;
@ -598,7 +595,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
@ -656,12 +653,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
@ -674,8 +670,7 @@ 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 - 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]) {

6
src/physics_mg.cpp
View file

@ -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<double>(gout + 1);
bank_array[i].delayed_group = dg + 1;

4
src/plot.cpp
View file

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

2
src/reaction.cpp
View file

@ -42,6 +42,8 @@ Reaction::Reaction(hid_t group, const std::vector<int>& 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);

3
src/state_point.cpp
View file

@ -123,9 +123,8 @@ 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-1]->name_);
data::nuclides[deriv.diff_nuclide]->name_);
} else if (deriv.variable == DIFF_TEMPERATURE) {
write_dataset(deriv_group, "independent variable", "temperature");
} else {

71
src/tallies/derivative.cpp
View file

@ -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) {
@ -131,8 +130,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;
@ -189,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:
@ -209,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;
@ -226,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 {
@ -240,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;
@ -254,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 {
@ -267,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 {
@ -280,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 {
@ -323,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;
@ -334,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);
@ -346,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);
@ -360,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);
@ -372,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);
@ -384,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);
@ -574,8 +571,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;
@ -594,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:
@ -622,8 +617,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;
@ -638,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);
}
@ -663,9 +656,7 @@ 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 - 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

19
src/tallies/filter.cpp
View file

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

6
src/tallies/filter_azimuthal.cpp
View file

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

9
src/tallies/filter_cell.cpp
View file

@ -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_);
}
//==============================================================================
@ -76,8 +74,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;

6
src/tallies/filter_cellborn.cpp
View file

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

6
src/tallies/filter_cellfrom.cpp
View file

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

6
src/tallies/filter_delayedgroup.cpp
View file

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

6
src/tallies/filter_distribcell.cpp
View file

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

24
src/tallies/filter_energy.cpp
View file

@ -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();
}
@ -123,7 +117,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<EnergyFilter*>(filt_base);
// Check the filter type.
@ -145,7 +139,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<EnergyFilter*>(filt_base);
// Check the filter type.

3
src/tallies/filter_energyfunc.cpp
View file

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

10
src/tallies/filter_legendre.cpp
View file

@ -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);
}
//==============================================================================
@ -52,7 +50,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<LegendreFilter*>(filt_base);
// Check the filter type.
@ -73,7 +71,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<LegendreFilter*>(filt_base);
// Check the filter type.

12
src/tallies/filter_material.cpp
View file

@ -42,10 +42,9 @@ 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);
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_);
}
//==============================================================================
@ -77,7 +75,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<MaterialFilter*>(filt_base);
// Check the filter type.
@ -99,7 +97,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<MaterialFilter*>(filt_base);
// Check the filter type.

7
src/tallies/filter_mesh.cpp
View file

@ -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);
}
@ -90,7 +89,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<MeshFilter*>(filt_base);
// Check the filter type.
@ -111,7 +110,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<MeshFilter*>(filt_base);
// Check the filter type.

5
src/tallies/filter_meshsurface.cpp
View file

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

6
src/tallies/filter_mu.cpp
View file

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

3
src/tallies/filter_particle.cpp
View file

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

6
src/tallies/filter_polar.cpp
View file

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

11
src/tallies/filter_sph_harm.cpp
View file

@ -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();
}
@ -100,7 +99,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<SphericalHarmonicsFilter*>(filt_base);
// Check the filter type.

8
src/tallies/filter_sptl_legendre.cpp
View file

@ -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();
}
@ -107,7 +105,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<SpatialLegendreFilter*>(filt_base);
// Check the filter type.

6
src/tallies/filter_surface.cpp
View file

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

6
src/tallies/filter_universe.cpp
View file

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

16
src/tallies/filter_zernike.cpp
View file

@ -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
@ -132,7 +128,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<ZernikeFilter*>(filt_base);
// Check the filter type.

104
src/tallies/tally.cpp
View file

@ -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<Tally>());
}
@ -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<std::string> 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.");

310
src/tallies/tally_scoring.cpp
View file

@ -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;
}
}
}
@ -337,9 +330,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 +435,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 +473,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 +487,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 +515,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 +558,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.;
}
@ -600,7 +590,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);
@ -630,8 +620,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// delayed-nu-fission
if (simulation::micro_xs[p->event_nuclide-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,12 +630,12 @@ 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;
score_fission_delayed_dg(i_tally, d_bin+1, score,
* simulation::micro_xs[p->event_nuclide].fission
/ simulation::micro_xs[p->event_nuclide].absorption * flux;
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
continue;
@ -654,10 +644,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 {
@ -680,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 {
@ -705,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 {
@ -723,7 +713,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
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);
@ -734,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);
}
}
@ -743,7 +733,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);
@ -766,8 +756,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// No fission events occur if survival biasing is on -- need to
// calculate fraction of absorptions that would have resulted in
// delayed-nu-fission
const auto& nuc {*data::nuclides[p->event_nuclide-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,10 +772,10 @@ 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_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
continue;
@ -805,8 +795,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 +816,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;
@ -839,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.;
@ -865,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 {
@ -890,7 +880,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<DelayedGroupFilter*>(
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);
@ -905,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);
}
}
@ -915,7 +905,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);
@ -953,13 +943,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 +957,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 {
@ -986,7 +976,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 +1012,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);
@ -1046,8 +1036,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 +1047,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 +1056,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 +1067,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 {
@ -1096,7 +1086,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 +1134,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);
@ -1175,18 +1165,17 @@ 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 {
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);
@ -1196,13 +1185,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;
}
}
@ -1219,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;
}
}
@ -1234,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;
@ -1274,7 +1260,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) {
@ -1325,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;
@ -1354,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
@ -1364,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);
@ -1382,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);
@ -1390,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);
@ -1412,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);
@ -1420,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);
}
@ -1443,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;
}
@ -1474,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 *=
@ -1483,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;
@ -1509,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 *=
@ -1527,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;
@ -1559,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 *=
@ -1580,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)
@ -1622,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 {
@ -1631,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;
@ -1642,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 *=
@ -1673,10 +1659,10 @@ 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+1, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
}
continue;
} else {
@ -1687,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);
}
}
@ -1703,20 +1689,20 @@ 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
* 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 {
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);
@ -1744,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 {
@ -1757,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;
@ -1770,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 {
@ -1804,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
@ -1822,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.;
@ -1842,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
@ -1853,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 {
@ -1861,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
@ -1895,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 *=
@ -1904,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)
@ -1926,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;
}
}
@ -1937,7 +1923,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux,
int filter_index)
{
const Tally& tally {*model::tallies[i_tally]};
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
// Score all individual nuclide reaction rates.
for (auto i = 0; i < material.nuclide_.size(); ++i) {
@ -1994,8 +1980,7 @@ void score_analog_tally_ce(const Particle* p)
// the event nuclide or the total material. Note that the i_nuclide
// and flux arguments for score_general are not used for analog
// tallies.
//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);
}
@ -2050,13 +2035,9 @@ 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]
->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-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,
@ -2109,13 +2090,9 @@ 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-1]
->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-1]->atom_density_(j);
atom_density = model::materials[p->material]->atom_density_(j);
}
}
@ -2179,13 +2156,9 @@ 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]
->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-1]->atom_density_(j);
atom_density = model::materials[p->material]->atom_density_(j);
}
//TODO: consider replacing this "if" with pointers or templates
@ -2239,9 +2212,8 @@ score_surface_tally(const Particle* p, const std::vector<int>& 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;
}
}

13
src/volume_calc.cpp
View file

@ -117,14 +117,11 @@ std::vector<VolumeCalculation::Result> 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::Result> 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::Result> 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;
}
}