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

Resolve off-by-one issues for ParticleType

Browse source
This commit is contained in:
Paul Romano 2019-01-31 22:54:32 -06:00
parent 380d23b0e4
commit 4d2284ca7d
15 changed files with 31 additions and 74 deletions
Download patch file Download diff file Expand all files Collapse all files

2
include/openmc/particle.h
View file

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

4
src/bremsstrahlung.cpp
View file

@ -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<int>(ParticleType::photon) - 1;
int photon = static_cast<int>(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<int>(ParticleType::positron)) {
mat = &model::materials[p.material -1]->ttb_->positron;

8
src/constants.F90
View file

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

14
src/cross_sections.cpp
View file

@ -218,8 +218,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& 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<int>(ParticleType::neutron) - 1;
int neutron = static_cast<int>(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<std::vector<double>>& nuc_temps,
// the previous
const auto& elem {data::elements.back()};
if (elem.energy_.size() >= 1) {
// TODO: off-by-one
int photon = static_cast<int>(ParticleType::photon) - 1;
int photon = static_cast<int>(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<std::vector<double>>& nuc_temps,
for (auto& nuc : data::nuclides) {
nuc->init_grid();
}
// TODO: off-by-one
int neutron = static_cast<int>(ParticleType::neutron) - 1;
int neutron = static_cast<int>(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<std::vector<double>>& 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<int>(ParticleType::photon) - 1;
int photon = static_cast<int>(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<std::vector<double>>& nuc_temps,
// grid has not been allocated
if (nuc->grid_.size() > 0) {
double max_E = nuc->grid_[0].energy.back();
int neutron = static_cast<int>(ParticleType::neutron) - 1;
int neutron = static_cast<int>(ParticleType::neutron);
if (max_E == data::energy_max[neutron]) {
write_message("Maximum neutron transport energy: " +
std::to_string(data::energy_max[neutron]) + " eV for " +

6
src/input_xml.F90
View file

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

3
src/material.cpp
View file

@ -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<int>(ParticleType::neutron) - 1;
int neutron = static_cast<int>(ParticleType::neutron);
// Find energy index on energy grid
int i_grid = std::log(p.E/data::energy_min[neutron])/simulation::log_spacing;

6
src/nuclide.cpp
View file

@ -392,7 +392,7 @@ void Nuclide::create_derived()
void Nuclide::init_grid()
{
int neutron = static_cast<int>(ParticleType::neutron) - 1;
int neutron = static_cast<int>(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(); }

5
src/nuclide_header.F90
View file

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

5
src/particle.cpp
View file

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

9
src/particle_header.F90
View file

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

5
src/photon.cpp
View file

@ -223,7 +223,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
}
// Truncate the bremsstrahlung data at the cutoff energy
int photon = static_cast<int>(ParticleType::photon) - 1;
int photon = static_cast<int>(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<int>(ParticleType::photon) - 1;
int photon = static_cast<int>(ParticleType::photon);
int n = element.energy_.size();
data::energy_min[photon] = std::max(data::energy_min[photon],
std::exp(element.energy_(1)));

10
src/physics.cpp
View file

@ -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<int>(ParticleType::photon) - 1;
int photon = static_cast<int>(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<int>(ParticleType::neutron) - 1;
constexpr int neutron = static_cast<int>(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<int>(ParticleType::neutron) - 1;
constexpr int neutron = static_cast<int>(ParticleType::neutron);
if (site->E < data::energy_max[neutron]) break;
// check for large number of resamples

11
src/reaction.cpp
View file

@ -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<int>(rx->products_[product - 1].particle_);
}
void reaction_product_sample(Reaction* rx, int product, double E_in, double* E_out, double* mu)

11
src/source.cpp
View file

@ -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<Discrete*>(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

6
src/tallies/filter_particle.cpp
View file

@ -8,6 +8,7 @@ void
ParticleFilter::from_xml(pugi::xml_node node)
{
particles_ = get_node_array<int>(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