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

Fix off-by-one for filter matches

Browse source
This commit is contained in:
Paul Romano 2019-02-25 21:09:17 -06:00
parent eeaa38c2ab
commit 2729d986bf
27 changed files with 102 additions and 172 deletions
Download patch file Download diff file Expand all files Collapse all files

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

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

8
src/output.cpp
View file

@ -679,13 +679,12 @@ write_tallies()
// prevents redundant output.
int indent = 0;
for (auto i = 0; i < tally.filters().size(); ++i) {
if ((filter_index-1) % tally.strides(i) == 0) {
if (filter_index % tally.strides(i) == 0) {
auto i_filt = tally.filters(i);
const auto& filt {*model::tally_filters[i_filt]};
auto& match {simulation::filter_matches[i_filt]};
tallies_out << std::string(indent+1, ' ')
// TODO: off-by-one
<< filt.text_label(match.i_bin_+1) << "\n";
<< filt.text_label(match.i_bin_) << "\n";
}
indent += 2;
}
@ -712,9 +711,8 @@ write_tallies()
std::string score_name = score > 0 ? reaction_name(score)
: score_names.at(score);
double mean, stdev;
//TODO: off-by-one
std::tie(mean, stdev) = mean_stdev(
&tally.results_(filter_index-1, score_index, 0), tally.n_realizations_);
&tally.results_(filter_index, score_index, 0), tally.n_realizations_);
tallies_out << std::string(indent+1, ' ') << std::left
<< std::setw(36) << score_name << " " << mean << " +/- "
<< t_value * stdev << "\n";

1
src/physics.cpp
View file

@ -670,7 +670,6 @@ void scatter(Particle* p, int i_nuclide)
p->event = EVENT_SCATTER;
// Sample new outgoing angle for isotropic-in-lab scattering
// TODO: off-by-one
const auto& mat {model::materials[p->material]};
if (!mat->p0_.empty()) {
int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];

1
src/state_point.cpp
View file

@ -123,7 +123,6 @@ openmc_statepoint_write(const char* filename, bool* write_source)
write_dataset(deriv_group, "independent variable", "density");
} else if (deriv.variable == DIFF_NUCLIDE_DENSITY) {
write_dataset(deriv_group, "independent variable", "nuclide_density");
//TODO: off-by-one
write_dataset(deriv_group, "nuclide",
data::nuclides[deriv.diff_nuclide]->name_);
} else if (deriv.variable == DIFF_TEMPERATURE) {

1
src/tallies/derivative.cpp
View file

@ -656,7 +656,6 @@ void score_collision_derivative(const Particle* p)
// Loop over the material's nuclides until we find the event nuclide.
for (auto i_nuc : material.nuclide_) {
const auto& nuc {*data::nuclides[i_nuc]};
//TODO: off-by-one
if (i_nuc == p->event_nuclide && multipole_in_range(&nuc, p->last_E)) {
// phi is proportional to Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s

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

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

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

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

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

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

6
src/tallies/filter_material.cpp
View file

@ -44,8 +44,7 @@ MaterialFilter::get_all_bins(const Particle* p, int estimator,
{
auto search = map_.find(p->material);
if (search != map_.end()) {
//TODO: off-by-one
match.bins_.push_back(search->second + 1);
match.bins_.push_back(search->second);
match.weights_.push_back(1.0);
}
}
@ -62,8 +61,7 @@ MaterialFilter::to_statepoint(hid_t filter_group) const
std::string
MaterialFilter::text_label(int bin) const
{
//TODO: off-by-one
return "Material " + std::to_string(model::materials[materials_[bin-1]]->id_);
return "Material " + std::to_string(model::materials[materials_[bin]]->id_);
}
//==============================================================================

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

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

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

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

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

14
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

96
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;
}
}
}
@ -642,7 +635,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
score = p->absorb_wgt * yield
* simulation::micro_xs[p->event_nuclide].fission
/ simulation::micro_xs[p->event_nuclide].absorption * flux;
score_fission_delayed_dg(i_tally, d_bin+1, score,
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
continue;
@ -677,7 +670,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
auto d = filt.groups_[d_bin];
score = simulation::keff * p->wgt_bank / p->n_bank
* p->n_delayed_bank[d-1] * flux;
score_fission_delayed_dg(i_tally, d_bin+1, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
}
continue;
} else {
@ -702,7 +695,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
->nu(E, ReactionProduct::EmissionMode::delayed, d);
score = simulation::micro_xs[i_nuclide].fission * yield
* atom_density * flux;
score_fission_delayed_dg(i_tally, d_bin+1, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
}
continue;
} else {
@ -731,7 +724,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
->nu(E, ReactionProduct::EmissionMode::delayed, d);
score = simulation::micro_xs[j_nuclide].fission * yield
* atom_density * flux;
score_fission_delayed_dg(i_tally, d_bin+1, score,
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
}
@ -782,7 +775,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
* simulation::micro_xs[p->event_nuclide].fission
/ simulation::micro_xs[p->event_nuclide].absorption
* rate * flux;
score_fission_delayed_dg(i_tally, d_bin+1, score,
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
continue;
@ -836,7 +829,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
if (d == g)
score_fission_delayed_dg(i_tally, d_bin+1, score,
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
score = 0.;
@ -862,7 +855,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
auto rate = rxn.products_[d].decay_rate_;
score = simulation::micro_xs[i_nuclide].fission * yield * flux
* atom_density * rate;
score_fission_delayed_dg(i_tally, d_bin+1, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
}
continue;
} else {
@ -902,7 +895,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
auto rate = rxn.products_[d].decay_rate_;
score = simulation::micro_xs[j_nuclide].fission * yield
* flux * atom_density * rate;
score_fission_delayed_dg(i_tally, d_bin+1, score,
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
}
@ -1214,7 +1207,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
// Update tally results
#pragma omp atomic
tally.results_(filter_index-1, score_index, RESULT_VALUE) += score;
tally.results_(filter_index, score_index, RESULT_VALUE) += score;
}
}
@ -1229,8 +1222,6 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
{
auto& tally {*model::tallies[i_tally]};
//TODO: off-by-one throughout on p->material
// Set the direction and group to use with get_xs
const double* p_uvw;
int p_g;
@ -1626,7 +1617,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
p_g, nullptr, nullptr, &d)
/ get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g);
}
score_fission_delayed_dg(i_tally, d_bin+1, score,
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
continue;
@ -1671,7 +1662,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
* get_nuclide_xs(i_nuclide+1, MG_GET_XS_FISSION, p_g)
/ get_macro_xs(p->material, MG_GET_XS_FISSION, p_g);
}
score_fission_delayed_dg(i_tally, d_bin+1, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
}
continue;
} else {
@ -1705,7 +1696,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
* get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION,
p_g, nullptr, nullptr, &d);
}
score_fission_delayed_dg(i_tally, d_bin+1, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
}
continue;
} else {
@ -1752,7 +1743,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
p_g, nullptr, nullptr, &d)
/ get_macro_xs(p->material, MG_GET_XS_ABSORPTION, p_g);
}
score_fission_delayed_dg(i_tally, d_bin+1, score,
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
continue;
@ -1817,7 +1808,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
for (auto d_bin = 0; d_bin < filt.n_bins_; ++d_bin) {
auto d = filt.groups_[d_bin];
if (d == g)
score_fission_delayed_dg(i_tally, d_bin+1, score,
score_fission_delayed_dg(i_tally, d_bin, score,
score_index);
}
score = 0.;
@ -1848,7 +1839,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
* get_macro_xs(p->material, MG_GET_XS_DELAYED_NU_FISSION,
p_g, nullptr, nullptr, &d);
}
score_fission_delayed_dg(i_tally, d_bin+1, score, score_index);
score_fission_delayed_dg(i_tally, d_bin, score, score_index);
}
continue;
} else {
@ -1921,7 +1912,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// Update tally results
#pragma omp atomic
tally.results_(filter_index-1, score_index, RESULT_VALUE) += score;
tally.results_(filter_index, score_index, RESULT_VALUE) += score;
}
}
@ -2044,12 +2035,8 @@ void score_analog_tally_mg(const Particle* p)
double atom_density = 0.;
if (i_nuclide >= 0) {
//TODO: off-by-one
auto j = model::materials[p->material]
->mat_nuclide_index_[i_nuclide];
auto j = model::materials[p->material]->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) continue;
//atom_density = material_atom_density(p->material, j);
//TODO: off-by-one
atom_density = model::materials[p->material]->atom_density_(j);
}
@ -2103,12 +2090,8 @@ score_tracklength_tally(const Particle* p, double distance)
double atom_density = 0.;
if (i_nuclide >= 0) {
if (p->material != MATERIAL_VOID) {
//TODO: off-by-one
auto j = model::materials[p->material]
->mat_nuclide_index_[i_nuclide];
auto j = model::materials[p->material]->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) continue;
//atom_density = material_atom_density(p->material, j);
//TODO: off-by-one
atom_density = model::materials[p->material]->atom_density_(j);
}
}
@ -2173,12 +2156,8 @@ void score_collision_tally(const Particle* p)
double atom_density = 0.;
if (i_nuclide >= 0) {
//TODO: off-by-one
auto j = model::materials[p->material]
->mat_nuclide_index_[i_nuclide];
auto j = model::materials[p->material]->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) continue;
//atom_density = material_atom_density(p->material, j);
//TODO: off-by-one
atom_density = model::materials[p->material]->atom_density_(j);
}
@ -2233,9 +2212,8 @@ score_surface_tally(const Particle* p, const std::vector<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;
}
}