Move most EnergyFilter functionality to C++

This commit is contained in:
Sterling Harper 2018-10-27 15:11:49 -04:00
parent 4ae02e3c1b
commit 8d3f9d3c69
7 changed files with 182 additions and 122 deletions

View file

@ -438,6 +438,7 @@ add_library(libopenmc SHARED
src/tallies/filter_cell.cpp
src/tallies/filter_distribcell.cpp
src/tallies/filter_energyfunc.cpp
src/tallies/filter_energy.cpp
src/tallies/filter_legendre.cpp
src/tallies/filter_material.cpp
src/tallies/filter_mesh.cpp

View file

@ -420,6 +420,7 @@ constexpr int DIFF_NUCLIDE_DENSITY {2};
constexpr int DIFF_TEMPERATURE {3};
constexpr int C_NONE {-1};
constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE
// Interpolation rules
enum class Interpolation {

View file

@ -0,0 +1,55 @@
#ifndef OPENMC_TALLIES_FILTER_ENERGY_H
#define OPENMC_TALLIES_FILTER_ENERGY_H
//#include <cstdint>
//#include <unordered_map>
#include <vector>
#include "openmc/tallies/filter.h"
namespace openmc {
//==============================================================================
//! Bins the incident neutron energy.
//==============================================================================
class EnergyFilter : public Filter
{
public:
~EnergyFilter() = default;
std::string type() const override {return "energy";}
void from_xml(pugi::xml_node node) override;
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
void to_statepoint(hid_t filter_group) const override;
std::string text_label(int bin) const override;
std::vector<double> bins_;
bool matches_transport_groups_ {false};
};
//==============================================================================
//! Bins the outgoing neutron energy.
//!
//! Only scattering events use the get_all_bins functionality. Nu-fission
//! tallies manually iterate over the filter bins.
//==============================================================================
class EnergyoutFilter : public EnergyFilter
{
public:
std::string type() const override {return "energyout";}
void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const override;
std::string text_label(int bin) const override;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_ENERGY_H

View file

@ -16,8 +16,6 @@ namespace openmc {
std::vector<int64_t> overlap_check_count;
constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE
//==============================================================================
extern "C" bool

View file

@ -11,6 +11,7 @@
#include "openmc/tallies/filter_cellfrom.h"
#include "openmc/tallies/filter_distribcell.h"
#include "openmc/tallies/filter_energyfunc.h"
#include "openmc/tallies/filter_energy.h"
#include "openmc/tallies/filter_legendre.h"
#include "openmc/tallies/filter_material.h"
#include "openmc/tallies/filter_mesh.h"
@ -105,6 +106,10 @@ extern "C" {
tally_filters.push_back(new DistribcellFilter());
} else if (type_ == "energyfunction") {
tally_filters.push_back(new EnergyFunctionFilter());
} else if (type_ == "energy") {
tally_filters.push_back(new EnergyFilter());
} else if (type_ == "energyout") {
tally_filters.push_back(new EnergyoutFilter());
} else if (type_ == "legendre") {
tally_filters.push_back(new LegendreFilter());
} else if (type_ == "material") {

View file

@ -0,0 +1,116 @@
#include "openmc/tallies/filter_energy.h"
#include "openmc/constants.h" // For F90_NONE
#include "openmc/mgxs_interface.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/xml_interface.h"
namespace openmc {
//==============================================================================
// EnergyFilter implementation
//==============================================================================
// Used to grab the rev_energy_bins array defined in Fortran.
extern "C" double* rev_energy_bins_ptr();
void
EnergyFilter::from_xml(pugi::xml_node node)
{
bins_ = get_node_array<double>(node, "bins");
n_bins_ = bins_.size() - 1;
// In MG mode, check if the filter bins match the transport bins.
// We can save tallying time if we know that the tally bins match the energy
// group structure. In that case, the matching bin index is simply the group
// (after flipping for the different ordering of the library and tallying
// systems).
if (!settings::run_CE) {
if (n_bins_ == num_energy_groups) {
matches_transport_groups_ = true;
double* rev_energy_bins = rev_energy_bins_ptr();
for (auto i = 0; i < n_bins_ + 1; i++) {
if (rev_energy_bins[i] != bins_[i]) {
matches_transport_groups_ = false;
break;
}
}
}
}
}
void
EnergyFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
const
{
if (p->g != F90_NONE && matches_transport_groups_) {
if (estimator == ESTIMATOR_TRACKLENGTH) {
//TODO: off-by-one
match.bins_.push_back(num_energy_groups - p->g + 1);
} else {
//TODO: off-by-one
match.bins_.push_back(num_energy_groups - p->last_g + 1);
}
match.weights_.push_back(1.0);
} else {
// Get the pre-collision energy of the particle.
auto E = p->last_E;
// 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;
match.bins_.push_back(bin);
match.weights_.push_back(1.0);
}
}
}
void
EnergyFilter::to_statepoint(hid_t filter_group) const
{
Filter::to_statepoint(filter_group);
write_dataset(filter_group, "bins", bins_);
}
std::string
EnergyFilter::text_label(int bin) const
{
std::stringstream out;
out << "Incoming Energy [" << bins_[bin-1] << ", " << bins_[bin] << ")";
return out.str();
}
//==============================================================================
// EnergyoutFilter implementation
//==============================================================================
void
EnergyoutFilter::get_all_bins(const Particle* p, int estimator,
FilterMatch& match) const
{
if (p->g != F90_NONE && matches_transport_groups_) {
match.bins_.push_back(num_energy_groups - p->g + 1);
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;
match.bins_.push_back(bin);
match.weights_.push_back(1.0);
}
}
}
std::string
EnergyoutFilter::text_label(int bin) const
{
std::stringstream out;
out << "Outgoing Energy [" << bins_[bin-1] << ", " << bins_[bin] << ")";
return out.str();
}
}// namespace openmc

View file

@ -10,7 +10,7 @@ module tally_filter_energy
use particle_header, only: Particle
use settings, only: run_CE
use string, only: to_str
use tally_filter_header
use tally_filter_cpp
use xml_interface
implicit none
@ -22,16 +22,13 @@ module tally_filter_energy
! ENERGYFILTER bins the incident neutron energy.
!===============================================================================
type, public, extends(TallyFilter) :: EnergyFilter
type, public, extends(CppTallyFilter) :: EnergyFilter
real(8), allocatable :: bins(:)
! True if transport group number can be used directly to get bin number
logical :: matches_transport_groups = .false.
contains
procedure :: from_xml => from_xml_energy
procedure :: get_all_bins => get_all_bins_energy
procedure :: to_statepoint => to_statepoint_energy
procedure :: text_label => text_label_energy
end type EnergyFilter
!===============================================================================
@ -41,11 +38,6 @@ module tally_filter_energy
!===============================================================================
type, public, extends(EnergyFilter) :: EnergyoutFilter
contains
! Inherit from_xml from EnergyFilter
procedure :: get_all_bins => get_all_bins_energyout
procedure :: to_statepoint => to_statepoint_energyout
procedure :: text_label => text_label_energyout
end type EnergyoutFilter
contains
@ -60,6 +52,8 @@ contains
integer :: n
call this % from_xml_cpp_inner(node)
n = node_word_count(node, "bins")
! Allocate and store bins
@ -81,116 +75,6 @@ contains
end if
end subroutine from_xml_energy
subroutine get_all_bins_energy(this, p, estimator, match)
class(EnergyFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
integer :: n
integer :: bin
real(8) :: E
n = this % n_bins
if (p % g /= NONE .and. this % matches_transport_groups) then
if (estimator == ESTIMATOR_TRACKLENGTH) then
call match % bins_push_back(num_energy_groups - p % g + 1)
call match % weights_push_back(ONE)
else
call match % bins_push_back(num_energy_groups - p % last_g + 1)
call match % weights_push_back(ONE)
end if
else
! Pre-collision energy of particle
E = p % last_E
! Search to find incoming energy bin.
bin = binary_search(this % bins, n + 1, E)
if (bin /= NO_BIN_FOUND) then
call match % bins_push_back(bin)
call match % weights_push_back(ONE)
end if
end if
end subroutine get_all_bins_energy
subroutine to_statepoint_energy(this, filter_group)
class(EnergyFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
call write_dataset(filter_group, "type", "energy")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % bins)
end subroutine to_statepoint_energy
function text_label_energy(this, bin) result(label)
class(EnergyFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
real(8) :: E0, E1
E0 = this % bins(bin)
E1 = this % bins(bin + 1)
label = "Incoming Energy [" // trim(to_str(E0)) // ", " &
// trim(to_str(E1)) // ")"
end function text_label_energy
!===============================================================================
! EnergyoutFilter methods
!===============================================================================
subroutine get_all_bins_energyout(this, p, estimator, match)
class(EnergyoutFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
integer :: n
integer :: bin
n = this % n_bins
if (p % g /= NONE .and. this % matches_transport_groups) then
! Tallies are ordered in increasing groups, group indices
! however are the opposite, so switch
call match % bins_push_back(num_energy_groups - p % g + 1)
call match % weights_push_back(ONE)
else
! Search to find incoming energy bin.
bin = binary_search(this % bins, n + 1, p % E)
if (bin /= NO_BIN_FOUND) then
call match % bins_push_back(bin)
call match % weights_push_back(ONE)
end if
end if
end subroutine get_all_bins_energyout
subroutine to_statepoint_energyout(this, filter_group)
class(EnergyoutFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
call write_dataset(filter_group, "type", "energyout")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % bins)
end subroutine to_statepoint_energyout
function text_label_energyout(this, bin) result(label)
class(EnergyoutFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
real(8) :: E0, E1
E0 = this % bins(bin)
E1 = this % bins(bin + 1)
label = "Outgoing Energy [" // trim(to_str(E0)) // ", " &
// trim(to_str(E1)) // ")"
end function text_label_energyout
!===============================================================================
! C API FUNCTIONS
!===============================================================================