Move the EnergyFunctionFilter to C++

This commit is contained in:
Sterling Harper 2018-10-17 19:59:45 -04:00
parent 9b5e402a11
commit 960c64bc18
3 changed files with 93 additions and 110 deletions

View file

@ -0,0 +1,89 @@
#ifndef OPENMC_TALLY_FILTER_ENERGYFUNC_H
#define OPENMC_TALLY_FILTER_ENERGYFUNC_H
#include <iomanip> // for setprecision
#include <ios> // for scientific
#include <sstream>
#include <vector>
#include "openmc/error.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/tallies/tally_filter.h"
namespace openmc {
class EnergyFunctionFilter : public TallyFilter
{
public:
virtual std::string type() const override {return "energyfunction";}
EnergyFunctionFilter()
: TallyFilter {}
{
n_bins_ = 1;
}
virtual ~EnergyFunctionFilter() override = default;
virtual void
from_xml(pugi::xml_node node) override
{
if (!settings::run_CE)
fatal_error("EnergyFunction filters are only supported for "
"continuous-energy transport calculations");
if (!check_for_node(node, "energy"))
fatal_error("Energy grid not specified for EnergyFunction filter.");
energy_ = get_node_array<double>(node, "energy");
if (!check_for_node(node, "y"))
fatal_error("y values not specified for EnergyFunction filter.");
y_ = get_node_array<double>(node, "y");
}
virtual void
get_all_bins(Particle* p, int estimator, TallyFilterMatch& match)
const override
{
if (p->last_E >= energy_.front() && p->last_E <= energy_.back()) {
// Search for the incoming energy bin.
auto i = lower_bound_index(energy_.begin(), energy_.end(), p->last_E);
// Compute the interpolation factor between the nearest bins.
double f = (p->last_E - energy_[i]) / (energy_[i+1] - energy_[i]);
// Interpolate on the lin-lin grid.
match.bins.push_back(1);
match.weights.push_back((1-f) * y_[i] + f * y_[i+1]);
}
}
virtual void
to_statepoint(hid_t filter_group) const override
{
TallyFilter::to_statepoint(filter_group);
write_dataset(filter_group, "energy", energy_);
write_dataset(filter_group, "y", y_);
}
virtual std::string
text_label(int bin) const override
{
std::stringstream out;
out << std::scientific << std::setprecision(1)
<< "Energy Function f"
<< "([ " << energy_.front() << ", ..., " << energy_.back() << "]) = "
<< "[" << y_.front() << ", ..., " << y_.back() << "]";
return out.str();
}
std::vector<double> energy_;
std::vector<double> y_;
};
} // namespace openmc
#endif // OPENMC_TALLY_FILTER_ENERGYFUNC_H

View file

@ -8,6 +8,7 @@
#include "openmc/tallies/tally_filter_cellborn.h"
#include "openmc/tallies/tally_filter_cellfrom.h"
#include "openmc/tallies/tally_filter_distribcell.h"
#include "openmc/tallies/tally_filter_energyfunc.h"
#include "openmc/tallies/tally_filter_legendre.h"
#include "openmc/tallies/tally_filter_material.h"
#include "openmc/tallies/tally_filter_mesh.h"
@ -98,6 +99,8 @@ extern "C" {
tally_filters.push_back(new CellFromFilter());
} else if (type_ == "distribcell") {
tally_filters.push_back(new DistribcellFilter());
} else if (type_ == "energyfunction") {
tally_filters.push_back(new EnergyFunctionFilter());
} else if (type_ == "legendre") {
tally_filters.push_back(new LegendreFilter());
} else if (type_ == "material") {

View file

@ -1,124 +1,15 @@
module tally_filter_energyfunc
use, intrinsic :: ISO_C_BINDING
use algorithm, only: binary_search
use constants
use error, only: fatal_error
use hdf5_interface
use particle_header, only: Particle
use settings, only: run_CE
use string, only: to_str
use tally_filter_header
use xml_interface
implicit none
private
!===============================================================================
! EnergyFunctionFilter multiplies tally scores by an arbitrary function of
! incident energy described by a piecewise linear-linear interpolation.
!===============================================================================
type, public, extends(TallyFilter) :: EnergyFunctionFilter
real(8), allocatable :: energy(:)
real(8), allocatable :: y(:)
contains
procedure :: from_xml
procedure :: get_all_bins => get_all_bins_energyfunction
procedure :: to_statepoint => to_statepoint_energyfunction
procedure :: text_label => text_label_energyfunction
type, extends(CppTallyFilter) :: EnergyFunctionFilter
end type EnergyFunctionFilter
contains
subroutine from_xml(this, node)
class(EnergyFunctionFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
integer :: n
this % n_bins = 1
! Make sure this is continuous-energy mode.
if (.not. run_CE) then
call fatal_error("EnergyFunction filters are only supported for &
&continuous-energy transport calculations")
end if
! Allocate and store energy grid.
if (.not. check_for_node(node, "energy")) then
call fatal_error("Energy grid not specified for EnergyFunction &
&filter.")
end if
n = node_word_count(node, "energy")
allocate(this % energy(n))
call get_node_array(node, "energy", this % energy)
! Allocate and store interpolant values.
if (.not. check_for_node(node, "y")) then
call fatal_error("y values not specified for EnergyFunction &
&filter.")
end if
n = node_word_count(node, "y")
allocate(this % y(n))
call get_node_array(node, "y", this % y)
end subroutine from_xml
subroutine get_all_bins_energyfunction(this, p, estimator, match)
class(EnergyFunctionFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
integer :: n, indx
real(8) :: E, f, weight
select type(this)
type is (EnergyFunctionFilter)
n = size(this % energy)
! Get pre-collision energy of particle
E = p % last_E
! Search to find incoming energy bin.
indx = binary_search(this % energy, n, E)
! Compute an interpolation factor between nearest bins.
f = (E - this % energy(indx)) &
/ (this % energy(indx+1) - this % energy(indx))
! Interpolate on the lin-lin grid.
call match % bins_push_back(1)
weight = (ONE - f) * this % y(indx) + f * this % y(indx+1)
call match % weights_push_back(weight)
end select
end subroutine get_all_bins_energyfunction
subroutine to_statepoint_energyfunction(this, filter_group)
class(EnergyFunctionFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
select type(this)
type is (EnergyFunctionFilter)
call write_dataset(filter_group, "type", "energyfunction")
call write_dataset(filter_group, "energy", this % energy)
call write_dataset(filter_group, "y", this % y)
end select
end subroutine to_statepoint_energyfunction
function text_label_energyfunction(this, bin) result(label)
class(EnergyFunctionFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
select type(this)
type is (EnergyFunctionFilter)
write(label, FMT="(A, ES8.1, A, ES8.1, A, ES8.1, A, ES8.1, A)") &
"Energy Function f([", this % energy(1), ", ..., ", &
this % energy(size(this % energy)), "]) = [", this % y(1), &
", ..., ", this % y(size(this % y)), "]"
end select
end function text_label_energyfunction
end module tally_filter_energyfunc