diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 6f2d58832..84e628fe2 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -35,7 +35,7 @@ constexpr double REL_MAX_LOST_PARTICLES {1.0e-6}; //! Particle types enum class ParticleType { - neutron = 1, photon = 2, electron = 3, positron = 4 + neutron, photon, electron, positron }; extern "C" { diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index 1ba1073a5..ce699a3f4 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -30,12 +30,10 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) { if (p.material == MATERIAL_VOID) return; - // TODO: off-by-one - int photon = static_cast(ParticleType::photon) - 1; + int photon = static_cast(ParticleType::photon); if (p.E < settings::energy_cutoff[photon]) return; // Get bremsstrahlung data for this material and particle type - // TODO: off-by-one BremsstrahlungData* mat; if (p.type == static_cast(ParticleType::positron)) { mat = &model::materials[p.material -1]->ttb_->positron; diff --git a/src/constants.F90 b/src/constants.F90 index 93a40abdb..0b3927f03 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -122,10 +122,10 @@ module constants ! Particle type integer, parameter :: & - NEUTRON = 1, & - PHOTON = 2, & - ELECTRON = 3, & - POSITRON = 4 + NEUTRON = 0, & + PHOTON = 1, & + ELECTRON = 2, & + POSITRON = 3 ! Angular distribution type integer, parameter :: & diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index a73dbbcb6..87d5a2f50 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -218,8 +218,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, // Determine if minimum/maximum energy for this nuclide is greater/less // than the previous if (data::nuclides[i_nuclide]->grid_.size() >= 1) { - // TODO: off-by-one - int neutron = static_cast(ParticleType::neutron) - 1; + int neutron = static_cast(ParticleType::neutron); data::energy_min[neutron] = std::max(data::energy_min[neutron], data::nuclides[i_nuclide]->grid_[0].energy.front()); data::energy_max[neutron] = std::min(data::energy_max[neutron], @@ -253,8 +252,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, // the previous const auto& elem {data::elements.back()}; if (elem.energy_.size() >= 1) { - // TODO: off-by-one - int photon = static_cast(ParticleType::photon) - 1; + int photon = static_cast(ParticleType::photon); int n = elem.energy_.size(); data::energy_min[photon] = std::max(data::energy_min[photon], std::exp(elem.energy_(1))); @@ -313,8 +311,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, for (auto& nuc : data::nuclides) { nuc->init_grid(); } - // TODO: off-by-one - int neutron = static_cast(ParticleType::neutron) - 1; + int neutron = static_cast(ParticleType::neutron); simulation::log_spacing = std::log(data::energy_max[neutron] / data::energy_min[neutron]) / settings::n_log_bins; @@ -322,8 +319,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, // Determine if minimum/maximum energy for bremsstrahlung is greater/less // than the current minimum/maximum if (data::ttb_e_grid.size() >= 1) { - // TODO: off-by-one - int photon = static_cast(ParticleType::photon) - 1; + int photon = static_cast(ParticleType::photon); int n_e = data::ttb_e_grid.size(); data::energy_min[photon] = std::max(data::energy_min[photon], data::ttb_e_grid(1)); data::energy_max[photon] = std::min(data::energy_max[photon], data::ttb_e_grid(n_e - 1)); @@ -339,7 +335,7 @@ read_ce_cross_sections(const std::vector>& nuc_temps, // grid has not been allocated if (nuc->grid_.size() > 0) { double max_E = nuc->grid_[0].energy.back(); - int neutron = static_cast(ParticleType::neutron) - 1; + int neutron = static_cast(ParticleType::neutron); if (max_E == data::energy_max[neutron]) { write_message("Maximum neutron transport energy: " + std::to_string(data::energy_max[neutron]) + " eV for " + diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ca6ccca34..102a95dfb 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1456,10 +1456,8 @@ contains call read_mg_cross_sections_header_c(file_id) ! Get the minimum and maximum energies - energy_min(NEUTRON) = energy_bins(num_energy_groups + 1) - energy_max(NEUTRON) = energy_bins(1) - call set_particle_energy_bounds(NEUTRON, energy_min(NEUTRON), & - energy_max(NEUTRON)) + call set_particle_energy_bounds(NEUTRON, & + energy_bins(num_energy_groups + 1), energy_bins(1)) ! Close MGXS HDF5 file call file_close(file_id) diff --git a/src/material.cpp b/src/material.cpp index 8f2cbab5a..51b273ca9 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -659,8 +659,7 @@ void Material::calculate_xs(const Particle& p) const void Material::calculate_neutron_xs(const Particle& p) const { - // TODO: off-by-one - int neutron = static_cast(ParticleType::neutron) - 1; + int neutron = static_cast(ParticleType::neutron); // Find energy index on energy grid int i_grid = std::log(p.E/data::energy_min[neutron])/simulation::log_spacing; diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 6c689977b..db7343172 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -392,7 +392,7 @@ void Nuclide::create_derived() void Nuclide::init_grid() { - int neutron = static_cast(ParticleType::neutron) - 1; + int neutron = static_cast(ParticleType::neutron); double E_min = data::energy_min[neutron]; double E_max = data::energy_max[neutron]; int M = settings::n_log_bins; @@ -939,8 +939,8 @@ extern "C" int openmc_load_nuclide(const char* name) extern "C" void set_particle_energy_bounds(int particle, double E_min, double E_max) { - data::energy_min[particle - 1] = E_min; - data::energy_max[particle - 1] = E_max; + data::energy_min[particle] = E_min; + data::energy_max[particle] = E_max; } extern "C" int nuclides_size() { return data::nuclide_map.size(); } diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index af3dab9a6..78d436bdd 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -150,11 +150,6 @@ module nuclide_header type(MaterialMacroXS), bind(C) :: material_xs ! Cache for current material !$omp threadprivate(micro_xs, material_xs) - ! Minimum/maximum energies - real(8) :: energy_min(2) = [ZERO, ZERO] - real(8) :: energy_max(2) = [INFINITY, INFINITY] - - interface function library_present_c(type, name) result(b) bind(C, name='library_present') import C_INT, C_CHAR, C_BOOL diff --git a/src/particle.cpp b/src/particle.cpp index c8a753801..e5a38a983 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -203,11 +203,6 @@ Particle::write_restart() const void reset_coord(LocalCoord* c) { c->reset(); } void particle_clear(Particle* p) { p->clear(); } -void particle_create_secondary(Particle* p, const double* uvw, double E, - int type, bool run_CE) -{ - p->create_secondary(uvw, E, type, run_CE); -} void particle_initialize(Particle* p) { p->initialize(); } void particle_from_source(Particle* p, const Bank* src) { diff --git a/src/particle_header.F90 b/src/particle_header.F90 index db2317f4b..268fbb285 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -115,15 +115,6 @@ module particle_header type(Particle), intent(inout) :: p end subroutine particle_clear - subroutine particle_create_secondary(p, uvw, E, type, run_CE) bind(C) - import Particle, C_DOUBLE, C_INT, C_BOOL - type(Particle), intent(inout) :: p - real(C_DOUBLE), intent(in) :: uvw(3) - real(C_DOUBLE), value :: E - integer(C_INT), value :: type - logical(C_BOOL), value :: run_CE - end subroutine particle_create_secondary - subroutine particle_initialize(p) bind(C) import Particle type(Particle), intent(inout) :: p diff --git a/src/photon.cpp b/src/photon.cpp index 52d3118dc..0dff1d47f 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -223,7 +223,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) } // Truncate the bremsstrahlung data at the cutoff energy - int photon = static_cast(ParticleType::photon) - 1; + int photon = static_cast(ParticleType::photon); const auto& E {electron_energy}; double cutoff = settings::energy_cutoff[photon]; if (cutoff > E(0)) { @@ -782,8 +782,7 @@ extern "C" void photon_from_hdf5(hid_t group) // the previous const auto& element {data::elements.back()}; if (element.energy_.size() >= 1) { - // TODO: off-by-one - int photon = static_cast(ParticleType::photon) - 1; + int photon = static_cast(ParticleType::photon); int n = element.energy_.size(); data::energy_min[photon] = std::max(data::energy_min[photon], std::exp(element.energy_(1))); diff --git a/src/physics.cpp b/src/physics.cpp index 4823dc0ab..4d0fbf27b 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -52,7 +52,7 @@ void collision(Particle* p) } // Kill particle if energy falls below cutoff - if (p->E < settings::energy_cutoff[p->type - 1]) { + if (p->E < settings::energy_cutoff[p->type]) { p->alive = false; p->wgt = 0.0; p->last_wgt = 0.0; @@ -219,7 +219,7 @@ void sample_photon_reaction(Particle* p) // Kill photon if below energy cutoff -- an extra check is made here because // photons with energy below the cutoff may have been produced by neutrons // reactions or atomic relaxation - int photon = static_cast(ParticleType::photon) - 1; + int photon = static_cast(ParticleType::photon); if (p->E < settings::energy_cutoff[photon]) { p->E = 0.0; p->alive = false; @@ -1034,8 +1034,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank rx->products_[group].sample(E_in, site->E, mu); // resample if energy is greater than maximum neutron energy - // TODO: off-by-one - constexpr int neutron = static_cast(ParticleType::neutron) - 1; + constexpr int neutron = static_cast(ParticleType::neutron); if (site->E < data::energy_max[neutron]) break; // check for large number of resamples @@ -1060,8 +1059,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank rx->products_[0].sample(E_in, site->E, mu); // resample if energy is greater than maximum neutron energy - // TODO: off-by-one - constexpr int neutron = static_cast(ParticleType::neutron) - 1; + constexpr int neutron = static_cast(ParticleType::neutron); if (site->E < data::energy_max[neutron]) break; // check for large number of resamples diff --git a/src/reaction.cpp b/src/reaction.cpp index 918d61fea..4e418e7fe 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -302,16 +302,7 @@ int reaction_product_emission_mode(Reaction* rx, int product) int reaction_product_particle(Reaction* rx, int product) { - switch (rx->products_[product - 1].particle_) { - case ParticleType::neutron: - return 1; - case ParticleType::photon: - return 2; - case ParticleType::electron: - return 3; - case ParticleType::positron: - return 4; - } + return static_cast(rx->products_[product - 1].particle_); } void reaction_product_sample(Reaction* rx, int product, double E_in, double* E_out, double* mu) diff --git a/src/source.cpp b/src/source.cpp index bfa119823..3997c3cc0 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -174,14 +174,11 @@ Bank SourceDistribution::sample() const // Determine material auto c = model::cells[cell_index - 1]; int32_t mat_index = c->material_[instance]; - auto m = model::materials[mat_index]; if (mat_index == MATERIAL_VOID) { found = false; } else { - bool fissionable; - openmc_material_get_fissionable(mat_index + 1, &fissionable); - if (!fissionable) found = false; + if (!model::materials[mat_index]->fissionable_) found = false; } } } @@ -212,10 +209,10 @@ Bank SourceDistribution::sample() const auto energy_ptr = dynamic_cast(energy_.get()); if (energy_ptr) { auto energies = xt::adapt(energy_ptr->x()); - if (xt::any(energies > data::energy_max[p-1])) { + if (xt::any(energies > data::energy_max[p])) { fatal_error("Source energy above range of energies of at least " "one cross section table"); - } else if (xt::any(energies < data::energy_min[p-1])) { + } else if (xt::any(energies < data::energy_min[p])) { fatal_error("Source energy below range of energies of at least " "one cross section table"); } @@ -226,7 +223,7 @@ Bank SourceDistribution::sample() const site.E = energy_->sample(); // Resample if energy falls outside minimum or maximum particle energy - if (site.E < data::energy_max[p-1] && site.E > data::energy_min[p-1]) break; + if (site.E < data::energy_max[p] && site.E > data::energy_min[p]) break; } // Set delayed group diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index 27558c59b..b8bc32286 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -8,6 +8,7 @@ void ParticleFilter::from_xml(pugi::xml_node node) { particles_ = get_node_array(node, "bins"); + for (auto& p : particles_) --p; n_bins_ = particles_.size(); } @@ -34,8 +35,7 @@ ParticleFilter::to_statepoint(hid_t filter_group) const std::string ParticleFilter::text_label(int bin) const { - //TODO: off-by-one - return "Particle " + std::to_string(particles_[bin-1]); + return "Particle " + std::to_string(particles_[bin]); } //============================================================================== @@ -43,6 +43,6 @@ ParticleFilter::text_label(int bin) const //============================================================================== extern "C" int particle_filter_particles(ParticleFilter* filt, int i) -{return filt->particles_[i-1];} +{return filt->particles_[i];} } // namespace openmc