From e57916896934ed0530c41f06af11df66db9a4240 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 17 Jun 2020 15:49:20 -0500 Subject: [PATCH] Calculate min/max energies during simulation_init --- include/openmc/simulation.h | 3 ++ src/cross_sections.cpp | 51 ---------------------------- src/simulation.cpp | 68 +++++++++++++++++++++++++++++++++++++ 3 files changed, 71 insertions(+), 51 deletions(-) diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index a3032ee01..96f2c632e 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -57,6 +57,9 @@ void allocate_banks(); //! Determine number of particles to transport per process void calculate_work(); +//! Determine energy limits for incident neutron/photon data +void determine_energy_limits(); + //! Initialize a batch void initialize_batch(); diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index 3654d85e3..0b19a51ac 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -226,16 +226,6 @@ read_ce_cross_sections(const std::vector>& nuc_temps, close_group(group); file_close(file_id); - // Determine if minimum/maximum energy for this nuclide is greater/less - // than the previous - if (data::nuclides[i_nuclide]->grid_.size() >= 1) { - int neutron = static_cast(Particle::Type::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], - data::nuclides[i_nuclide]->grid_[0].energy.back()); - } - // Add name and alias to dictionary already_read.insert(name); @@ -257,18 +247,6 @@ read_ce_cross_sections(const std::vector>& nuc_temps, hid_t group = open_group(file_id, element.c_str()); data::elements.emplace_back(group); - // Determine if minimum/maximum energy for this element is greater/less than - // the previous - const auto& elem {data::elements.back()}; - if (elem.energy_.size() >= 1) { - int photon = static_cast(Particle::Type::photon); - int n = elem.energy_.size(); - data::energy_min[photon] = std::max(data::energy_min[photon], - std::exp(elem.energy_(1))); - data::energy_max[photon] = std::min(data::energy_max[photon], - std::exp(elem.energy_(n - 1))); - } - close_group(group); file_close(file_id); @@ -317,39 +295,10 @@ read_ce_cross_sections(const std::vector>& nuc_temps, } // materials if (settings::photon_transport && settings::electron_treatment == ElectronTreatment::TTB) { - // Determine if minimum/maximum energy for bremsstrahlung is greater/less - // than the current minimum/maximum - if (data::ttb_e_grid.size() >= 1) { - int photon = static_cast(Particle::Type::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)); - } - // Take logarithm of energies since they are log-log interpolated data::ttb_e_grid = xt::log(data::ttb_e_grid); } - // Show which nuclide results in lowest energy for neutron transport - for (const auto& nuc : data::nuclides) { - // If a nuclide is present in a material that's not used in the model, its - // grid has not been allocated - if (nuc->grid_.size() > 0) { - double max_E = nuc->grid_[0].energy.back(); - int neutron = static_cast(Particle::Type::neutron); - if (max_E == data::energy_max[neutron]) { - write_message("Maximum neutron transport energy: " + - std::to_string(data::energy_max[neutron]) + " eV for " + - nuc->name_, 7); - if (mpi::master && data::energy_max[neutron] < 20.0e6) { - warning("Maximum neutron energy is below 20 MeV. This may bias " - " the results."); - } - break; - } - } - } - // Show minimum/maximum temperature write_message("Minimum neutron data temperature: " + std::to_string(data::temperature_min) + " K", 4); diff --git a/src/simulation.cpp b/src/simulation.cpp index fa497e5e0..37b9bd882 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -34,6 +34,7 @@ #endif #include +#include #include @@ -68,6 +69,11 @@ int openmc_simulation_init() // Skip if simulation has already been initialized if (simulation::initialized) return 0; + // Determine limits for incident neutron/photon data + if (settings::run_CE) { + determine_energy_limits(); + } + // Set up logarithmic grid for nuclides for (auto& nuc : data::nuclides) { nuc->init_grid(); @@ -563,6 +569,68 @@ void calculate_work() } } +void determine_energy_limits() +{ + // Determine minimum/maximum energy for incident neutron/photon data + data::energy_max = {INFTY, INFTY}; + data::energy_min = {0.0, 0.0}; + for (const auto& nuc : data::nuclides) { + if (nuc->grid_.size() >= 1) { + int neutron = static_cast(Particle::Type::neutron); + data::energy_min[neutron] = std::max(data::energy_min[neutron], + nuc->grid_[0].energy.front()); + data::energy_max[neutron] = std::min(data::energy_max[neutron], + nuc->grid_[0].energy.back()); + } + } + + if (settings::photon_transport) { + for (const auto& elem : data::elements) { + if (elem.energy_.size() >= 1) { + int photon = static_cast(Particle::Type::photon); + int n = elem.energy_.size(); + data::energy_min[photon] = std::max(data::energy_min[photon], + std::exp(elem.energy_(1))); + data::energy_max[photon] = std::min(data::energy_max[photon], + std::exp(elem.energy_(n - 1))); + } + } + + if (settings::electron_treatment == ElectronTreatment::TTB) { + // Determine if minimum/maximum energy for bremsstrahlung is greater/less + // than the current minimum/maximum + if (data::ttb_e_grid.size() >= 1) { + int photon = static_cast(Particle::Type::photon); + int n_e = data::ttb_e_grid.size(); + data::energy_min[photon] = std::max(data::energy_min[photon], + std::exp(data::ttb_e_grid(1))); + data::energy_max[photon] = std::min(data::energy_max[photon], + std::exp(data::ttb_e_grid(n_e - 1))); + } + } + } + + // Show which nuclide results in lowest energy for neutron transport + for (const auto& nuc : data::nuclides) { + // If a nuclide is present in a material that's not used in the model, its + // grid has not been allocated + if (nuc->grid_.size() > 0) { + double max_E = nuc->grid_[0].energy.back(); + int neutron = static_cast(Particle::Type::neutron); + if (max_E == data::energy_max[neutron]) { + write_message("Maximum neutron transport energy: " + + std::to_string(data::energy_max[neutron]) + " eV for " + + nuc->name_, 7); + if (mpi::master && data::energy_max[neutron] < 20.0e6) { + warning("Maximum neutron energy is below 20 MeV. This may bias " + "the results."); + } + break; + } + } + } +} + #ifdef OPENMC_MPI void broadcast_results() { // Broadcast tally results so that each process has access to results