Move SurfaceFilter to C++

This commit is contained in:
Sterling Harper 2018-10-14 13:28:03 -04:00
parent 6c4891d591
commit ea0a91b0e6
5 changed files with 89 additions and 104 deletions

View file

@ -0,0 +1,85 @@
#ifndef OPENMC_TALLY_FILTER_SURFACE_H
#define OPENMC_TALLY_FILTER_SURFACE_H
#include <cstdint>
#include <sstream>
#include <unordered_map>
#include <vector>
#include "openmc/error.h"
#include "openmc/surface.h"
#include "openmc/tallies/tally_filter.h"
namespace openmc {
class SurfaceFilter : public TallyFilter
{
public:
virtual std::string type() const override {return "surface";}
virtual ~SurfaceFilter() override = default;
virtual void
from_xml(pugi::xml_node node) override
{
surfaces_ = get_node_array<int32_t>(node, "bins");
n_bins_ = surfaces_.size();
}
virtual void
initialize() override
{
for (auto& s : surfaces_) {
auto search = surface_map.find(s);
if (search != surface_map.end()) {
s = search->second;
} else {
std::stringstream err_msg;
err_msg << "Could not find surface " << s
<< " specified on tally filter.";
fatal_error(err_msg);
}
}
for (int i = 0; i < surfaces_.size(); i++) {
map_[surfaces_[i]] = i;
}
}
virtual void
get_all_bins(Particle* p, int estimator, TallyFilterMatch& match)
const override
{
auto search = map_.find(std::abs(p->surface)-1);
if (search != map_.end()) {
match.bins.push_back(search->second + 1);
if (p->surface < 0) {
match.weights.push_back(-1);
} else {
match.weights.push_back(1);
}
}
}
virtual void
to_statepoint(hid_t filter_group) const override
{
TallyFilter::to_statepoint(filter_group);
std::vector<int32_t> surface_ids;
for (auto c : surfaces_) surface_ids.push_back(surfaces[c]->id_);
write_dataset(filter_group, "bins", surface_ids);
}
virtual std::string
text_label(int bin) const override
{
return "Surface " + std::to_string(surfaces[surfaces_[bin-1]]->id_);
}
std::vector<int32_t> surfaces_;
std::unordered_map<int32_t, int> map_;
};
} // namespace openmc
#endif // OPENMC_TALLY_FILTER_SURFACE_H

View file

@ -468,9 +468,6 @@ contains
surfaces(i) % ptr = surface_pointer(i - 1);
if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true.
! Add surface to dictionary
call surface_dict % set(surfaces(i) % id(), i)
end do
! Check to make sure a boundary condition was applied to at least one
@ -628,8 +625,6 @@ contains
do i = 1, n_surfaces
surfaces(i) % ptr = surface_pointer(i - 1);
! Add surface to dictionary
call surface_dict % set(surfaces(i) % id(), i)
end do
end subroutine allocate_surfaces

View file

@ -2,7 +2,6 @@ module surface_header
use, intrinsic :: ISO_C_BINDING
use dict_header, only: DictIntInt
use hdf5_interface
implicit none
@ -83,9 +82,6 @@ module surface_header
type(Surface), allocatable, target :: surfaces(:)
! Dictionary that maps user IDs to indices in 'surfaces'
type(DictIntInt) :: surface_dict
contains
pure function surface_id(this) result(id)
@ -128,7 +124,6 @@ contains
subroutine free_memory_surfaces()
if (allocated(surfaces)) deallocate(surfaces)
call surface_dict % clear()
call free_memory_surfaces_c()
end subroutine free_memory_surfaces

View file

@ -13,6 +13,7 @@
#include "openmc/tallies/tally_filter_meshsurface.h"
#include "openmc/tallies/tally_filter_mu.h"
#include "openmc/tallies/tally_filter_polar.h"
#include "openmc/tallies/tally_filter_surface.h"
#include "openmc/tallies/tally_filter_universe.h"
@ -105,6 +106,8 @@ extern "C" {
tally_filters.push_back(new MuFilter());
} else if (type_ == "polar") {
tally_filters.push_back(new PolarFilter());
} else if (type_ == "surface") {
tally_filters.push_back(new SurfaceFilter());
} else if (type_ == "universe") {
tally_filters.push_back(new UniverseFilter());
} else {

View file

@ -1,109 +1,16 @@
module tally_filter_surface
use, intrinsic :: ISO_C_BINDING
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error
use hdf5_interface
use surface_header
use particle_header, only: Particle
use string, only: to_str
use tally_filter_header
use xml_interface
implicit none
private
!===============================================================================
! SURFACEFILTER specifies which surface particles are crossing
!===============================================================================
type, public, extends(TallyFilter) :: SurfaceFilter
integer, allocatable :: surfaces(:)
type, extends(CppTallyFilter) :: SurfaceFilter
! True if this filter is used for surface currents
logical :: current = .false.
contains
procedure :: from_xml
procedure :: get_all_bins => get_all_bins_surface
procedure :: to_statepoint => to_statepoint_surface
procedure :: text_label => text_label_surface
procedure :: initialize => initialize_surface
end type SurfaceFilter
contains
subroutine from_xml(this, node)
class(SurfaceFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
integer :: n
n = node_word_count(node, "bins")
! Allocate and store bins
this % n_bins = n
allocate(this % surfaces(n))
call get_node_array(node, "bins", this % surfaces)
end subroutine from_xml
subroutine get_all_bins_surface(this, p, estimator, match)
class(SurfaceFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
integer :: i
do i = 1, this % n_bins
if (abs(p % surface) == this % surfaces(i)) then
call match % bins_push_back(i)
if (p % surface < 0) then
call match % weights_push_back(-ONE)
else
call match % weights_push_back(ONE)
end if
exit
end if
end do
end subroutine get_all_bins_surface
subroutine to_statepoint_surface(this, filter_group)
class(SurfaceFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
call write_dataset(filter_group, "type", "surface")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "bins", this % surfaces)
end subroutine to_statepoint_surface
subroutine initialize_surface(this)
class(SurfaceFilter), intent(inout) :: this
integer :: i, id
integer :: val
! Convert ids to indices.
do i = 1, this % n_bins
id = this % surfaces(i)
val = surface_dict % get(id)
if (val /= EMPTY) then
this % surfaces(i) = val
else
call fatal_error("Could not find surface " // trim(to_str(id)) &
&// " specified on tally filter.")
end if
end do
end subroutine initialize_surface
function text_label_surface(this, bin) result(label)
class(SurfaceFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
label = "Surface " // to_str(surfaces(this % surfaces(bin)) % id())
end function text_label_surface
end module tally_filter_surface