Use 0-based indices for materials consistently

This commit is contained in:
Paul Romano 2019-02-25 08:54:19 -06:00
parent 5c6633c2c3
commit 4a89075890
17 changed files with 82 additions and 107 deletions

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

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)

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

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

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;

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

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

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

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

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

View file

@ -425,8 +425,7 @@ int sample_nuclide(const Particle* p)
double cutoff = prn() * simulation::material_xs.total;
// Get pointers to nuclide/density arrays
// TODO: off-by-one
const auto& mat {model::materials[p->material - 1]};
const auto& mat {model::materials[p->material]};
int n = mat->nuclide_.size();
double prob = 0.0;
@ -451,7 +450,7 @@ int sample_element(Particle* p)
double cutoff = prn() * simulation::material_xs.total;
// Get pointers to elements, densities
const auto& mat {model::materials[p->material - 1]};
const auto& mat {model::materials[p->material]};
int n = mat->nuclide_.size();
int i = 0;
@ -675,7 +674,7 @@ void scatter(Particle* p, int i_nuclide)
// Sample new outgoing angle for isotropic-in-lab scattering
// TODO: off-by-one
const auto& mat {model::materials[p->material - 1]};
const auto& mat {model::materials[p->material]};
if (!mat->p0_.empty()) {
int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
if (mat->p0_[i_nuc_mat]) {

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;

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

View file

@ -131,8 +131,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
score *= flux_deriv;
return;
}
//TODO: off-by-one
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
if (material.id_ != deriv.diff_material) {
score *= flux_deriv;
return;
@ -574,8 +573,7 @@ score_track_derivative(const Particle* p, double distance)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material == MATERIAL_VOID) return;
//TODO: off-by-one
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto& deriv : model::tally_derivs) {
if (deriv.diff_material != material.id_) continue;
@ -622,8 +620,7 @@ void score_collision_derivative(const Particle* p)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material == MATERIAL_VOID) return;
//TODO: off-by-one
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto& deriv : model::tally_derivs) {
if (deriv.diff_material != material.id_) continue;

View file

@ -42,7 +42,7 @@ void
MaterialFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
auto search = map_.find(p->material - 1);
auto search = map_.find(p->material);
if (search != map_.end()) {
//TODO: off-by-one
match.bins_.push_back(search->second + 1);

View file

@ -600,7 +600,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
score = 0.;
// Add up contributions from each nuclide in the material.
if (p->material != MATERIAL_VOID) {
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto j_nuclide = material.nuclide_[i];
auto atom_density = material.atom_density_(i);
@ -723,7 +723,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<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);
@ -743,7 +743,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
} else {
score = 0.;
if (p->material != MATERIAL_VOID) {
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto j_nuclide = material.nuclide_[i];
auto atom_density = material.atom_density_(i);
@ -890,7 +890,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
{*dynamic_cast<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);
@ -915,7 +915,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
} else {
score = 0.;
if (p->material != MATERIAL_VOID) {
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto j_nuclide = material.nuclide_[i];
auto atom_density = material.atom_density_(i);
@ -986,7 +986,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
}
} else {
if (p->material != MATERIAL_VOID) {
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto j_nuclide = material.nuclide_[i];
auto atom_density = material.atom_density_(i);
@ -1022,7 +1022,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
} else {
score = 0.;
if (p->material != MATERIAL_VOID) {
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto j_nuclide = material.nuclide_[i];
auto atom_density = material.atom_density_(i);
@ -1096,7 +1096,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
* atom_density * flux;
} else {
if (p->material != MATERIAL_VOID) {
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto j_nuclide = material.nuclide_[i];
auto atom_density = material.atom_density_(i);
@ -1144,7 +1144,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
} else {
score = 0.;
if (p->material != MATERIAL_VOID) {
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto j_nuclide = material.nuclide_[i];
auto atom_density = material.atom_density_(i);
@ -1186,7 +1186,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
}
} else {
if (p->material != MATERIAL_VOID) {
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto j_nuclide = material.nuclide_[i];
auto atom_density = material.atom_density_(i);
@ -1274,7 +1274,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
// To significantly reduce de-referencing, point matxs to the macroscopic
// Mgxs for the material of interest
data::macro_xs[p->material - 1].set_angle_index(p_uvw);
data::macro_xs[p->material].set_angle_index(p_uvw);
// Do same for nucxs, point it to the microscopic nuclide data of interest
if (i_nuclide >= 0) {
@ -1937,7 +1937,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux,
int filter_index)
{
const Tally& tally {*model::tallies[i_tally]};
const Material& material {*model::materials[p->material-1]};
const Material& material {*model::materials[p->material]};
// Score all individual nuclide reaction rates.
for (auto i = 0; i < material.nuclide_.size(); ++i) {
@ -2051,12 +2051,12 @@ void score_analog_tally_mg(const Particle* p)
double atom_density = 0.;
if (i_nuclide >= 0) {
//TODO: off-by-one
auto j = model::materials[p->material-1]
auto j = model::materials[p->material]
->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) continue;
//atom_density = material_atom_density(p->material, j);
//TODO: off-by-one
atom_density = model::materials[p->material-1]->atom_density_(j);
atom_density = model::materials[p->material]->atom_density_(j);
}
score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index,
@ -2110,12 +2110,12 @@ score_tracklength_tally(const Particle* p, double distance)
if (i_nuclide >= 0) {
if (p->material != MATERIAL_VOID) {
//TODO: off-by-one
auto j = model::materials[p->material-1]
auto j = model::materials[p->material]
->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) continue;
//atom_density = material_atom_density(p->material, j);
//TODO: off-by-one
atom_density = model::materials[p->material-1]->atom_density_(j);
atom_density = model::materials[p->material]->atom_density_(j);
}
}
@ -2180,12 +2180,12 @@ void score_collision_tally(const Particle* p)
double atom_density = 0.;
if (i_nuclide >= 0) {
//TODO: off-by-one
auto j = model::materials[p->material-1]
auto j = model::materials[p->material]
->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) continue;
//atom_density = material_atom_density(p->material, j);
//TODO: off-by-one
atom_density = model::materials[p->material-1]->atom_density_(j);
atom_density = model::materials[p->material]->atom_density_(j);
}
//TODO: consider replacing this "if" with pointers or templates

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