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

Move read_tallies_xml to C++

Browse source
This commit is contained in:
Paul Romano 2019-02-20 23:34:36 -06:00
parent 627ee79f0c
commit d0b6d9fbd6
14 changed files with 413 additions and 562 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -310,7 +310,6 @@ add_library(libopenmc SHARED
src/geometry.F90
src/hdf5_interface.F90
src/initialize.F90
src/input_xml.F90
src/material_header.F90
src/message_passing.F90
src/particle_header.F90

2
include/openmc/mesh.h
View file

@ -124,7 +124,7 @@ private:
//! Read meshes from either settings/tallies
//! \param[in] root XML node
extern "C" void read_meshes(pugi::xml_node* root);
void read_meshes(pugi::xml_node root);
//! Write mesh data to an HDF5 group
//! \param[in] group HDF5 group

2
include/openmc/settings.h
View file

@ -23,7 +23,7 @@ namespace openmc {
namespace settings {
// Boolean flags
extern "C" bool assume_separate; //!< assume tallies are spatially separate?
extern bool assume_separate; //!< assume tallies are spatially separate?
extern bool check_overlaps; //!< check overlaps in geometry?
extern bool confidence_intervals; //!< use confidence intervals for results?
extern bool create_fission_neutrons; //!< create fission neutrons (fixed source)?

3
include/openmc/tallies/derivative.h
View file

@ -29,6 +29,9 @@ struct TallyDerivative {
// Non-method functions
//==============================================================================
//! Read tally derivatives from a tallies.xml file
void read_tally_derivatives(pugi::xml_node node);
//! Scale the given score by its logarithmic derivative
void

14
include/openmc/tallies/filter.h
View file

@ -4,6 +4,7 @@
#include <cstdint>
#include <memory>
#include <string>
#include <unordered_map>
#include <vector>
#include "openmc/hdf5_interface.h"
@ -89,13 +90,16 @@ extern std::vector<FilterMatch> filter_matches;
} // namespace simulation
namespace model {
extern "C" int32_t n_filters;
extern std::vector<std::unique_ptr<Filter>> tally_filters;
} // namespace model
extern "C" int32_t n_filters;
extern std::vector<std::unique_ptr<Filter>> tally_filters;
extern std::unordered_map<int, int> filter_map;
}
//==============================================================================
// Non-member functions
//==============================================================================
Filter* allocate_filter(const std::string& type);
// Filter-related Fortran functions that will be called from C++
extern "C" int verify_filter(int32_t index);

8
include/openmc/tallies/tally.h
View file

@ -9,6 +9,7 @@
#include "xtensor/xtensor.hpp"
#include <memory> // for unique_ptr
#include <unordered_map>
#include <string>
#include <vector>
@ -46,7 +47,7 @@ public:
//----------------------------------------------------------------------------
// Other methods.
void init_triggers(pugi::xml_node node, int i_tally);
void init_triggers(pugi::xml_node node);
void init_results();
@ -123,6 +124,8 @@ namespace model {
extern std::vector<int> active_collision_tallies;
extern std::vector<int> active_meshsurf_tallies;
extern std::vector<int> active_surface_tallies;
extern std::unordered_map<int, int> tally_map;
}
namespace simulation {
@ -149,6 +152,9 @@ extern double global_tally_leakage;
// Non-member functions
//==============================================================================
//! Read tally specification from tallies.xml
void read_tallies_xml();
//! \brief Accumulate the sum of the contributions from each history within the
//! batch to a new random variable
void accumulate_tallies();

2
src/initialize.cpp
View file

@ -28,12 +28,12 @@
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/summary.h"
#include "openmc/tallies/tally.h"
#include "openmc/thermal.h"
#include "openmc/timer.h"
// data/functions from Fortran side
extern "C" void read_command_line();
extern "C" void read_tallies_xml();
// Paths to various files
extern "C" {

471
src/input_xml.F90
View file

@ -1,471 +0,0 @@
module input_xml
use, intrinsic :: ISO_C_BINDING
use constants
use error, only: fatal_error, warning, write_message, openmc_err_msg
use material_header
use message_passing
use settings
use stl_vector, only: VectorInt, VectorReal, VectorChar
use string, only: to_lower, to_str, str_to_int, &
starts_with, ends_with, to_c_string
use tally
use tally_header, only: openmc_extend_tallies
use tally_derivative_header
use tally_filter_header
use tally_filter
use xml_interface
implicit none
save
contains
!===============================================================================
! READ_TALLIES_XML reads data from a tallies.xml file and parses it, checking
! for errors and placing properly-formatted data in the right data structures
!===============================================================================
subroutine read_tallies_xml() bind(C)
integer :: i ! loop over user-specified tallies
integer :: j ! loop over words
integer :: l ! loop over bins
integer :: filter_id ! user-specified identifier for filter
integer :: tally_id ! user-specified identifier for filter
integer :: deriv_id
integer :: i_filt ! index in filters array
integer :: n ! size of arrays in mesh specification
integer :: n_words ! number of words read
integer :: n_filter ! number of filters
integer :: i_start, i_end
integer(C_INT) :: err
logical :: file_exists ! does tallies.xml file exist?
integer, allocatable :: temp_filter(:) ! temporary filter indices
logical :: has_energyout
integer :: particle_filter_index
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: temp_str
type(TallyFilterContainer), pointer :: f
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_tal
type(XMLNode) :: node_filt
type(XMLNode), allocatable :: node_tal_list(:)
type(XMLNode), allocatable :: node_filt_list(:)
type(TallyDerivative), pointer :: deriv
interface
subroutine tally_init_from_xml(tally_ptr, xml_node) bind(C)
import C_PTR
type(C_PTR), value :: tally_ptr
type(C_PTR) :: xml_node
end subroutine
subroutine tally_set_scores(tally_ptr, xml_node) bind(C)
import C_PTR
type(C_PTR), value :: tally_ptr
type(C_PTR) :: xml_node
end subroutine
subroutine tally_set_nuclides(tally_ptr, xml_node) bind(C)
import C_PTR
type(C_PTR), value :: tally_ptr
type(C_PTR) :: xml_node
end subroutine
subroutine tally_init_triggers(tally_ptr, i_tally, xml_node) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: tally_ptr
integer(C_INT), value :: i_tally
type(C_PTR) :: xml_node
end subroutine
subroutine read_tally_derivatives(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine
subroutine read_meshes(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine
end interface
! Check if tallies.xml exists
filename = trim(path_input) // "tallies.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
! Since a tallies.xml file is optional, no error is issued here
return
end if
! Display output message
call write_message("Reading tallies XML file...", 5)
! Parse tallies.xml file
call doc % load_file(filename)
root = doc % document_element()
! ==========================================================================
! DETERMINE SIZE OF ARRAYS AND ALLOCATE
! Get pointer list to XML <filter>
call get_node_list(root, "filter", node_filt_list)
! Get pointer list to XML <tally>
call get_node_list(root, "tally", node_tal_list)
! Check for <assume_separate> setting
if (check_for_node(root, "assume_separate")) then
call get_node_value(root, "assume_separate", assume_separate)
end if
! ==========================================================================
! READ MESH DATA
! Check for user meshes and allocate
call read_meshes(root % ptr)
! We only need the mesh info for plotting
if (run_mode == MODE_PLOTTING) then
call doc % clear()
return
end if
! ==========================================================================
! READ DATA FOR DERIVATIVES
call read_tally_derivatives(root % ptr)
! ==========================================================================
! READ FILTER DATA
! Check for user filters and allocate
n = size(node_filt_list)
if (n > 0) then
err = openmc_extend_filters(n, i_start, i_end)
end if
READ_FILTERS: do i = 1, n
f => filters(i_start + i - 1)
! Get pointer to filter xml node
node_filt = node_filt_list(i)
! Copy filter id
if (check_for_node(node_filt, "id")) then
call get_node_value(node_filt, "id", filter_id)
else
call fatal_error("Must specify id for filter in tally XML file.")
end if
! Check to make sure 'id' hasn't been used
if (filter_dict % has(filter_id)) then
call fatal_error("Two or more filters use the same unique ID: " &
// to_str(filter_id))
end if
! Convert filter type to lower case
temp_str = ''
if (check_for_node(node_filt, "type")) &
call get_node_value(node_filt, "type", temp_str)
temp_str = to_lower(temp_str)
! Make sure bins have been set
select case(temp_str)
case ("energy", "energyout", "mu", "polar", "azimuthal")
if (.not. check_for_node(node_filt, "bins")) then
call fatal_error("Bins not set in filter " // trim(to_str(filter_id)))
end if
case ("mesh", "meshsurface", "universe", "material", "cell", "distribcell", &
"cellborn", "cellfrom", "surface", "delayedgroup")
if (.not. check_for_node(node_filt, "bins")) then
call fatal_error("Bins not set in filter " // trim(to_str(filter_id)))
end if
end select
! Allocate according to the filter type
err = openmc_filter_set_type(i_start + i - 1, to_c_string(temp_str))
if (err /= 0) call fatal_error(to_f_string(openmc_err_msg))
! Read filter data from XML
call f % obj % from_xml(node_filt)
! Set filter id
err = openmc_filter_set_id(i_start + i - 1, filter_id)
! Initialize filter
call f % obj % initialize()
end do READ_FILTERS
! ==========================================================================
! READ TALLY DATA
! Check for user tallies
n = size(node_tal_list)
if (n == 0) then
if (master) call warning("No tallies present in tallies.xml file!")
end if
! Allocate user tallies
if (n > 0 .and. run_mode /= MODE_PLOTTING) then
err = openmc_extend_tallies(n, i_start, i_end)
end if
READ_TALLIES: do i = 1, n
! Allocate tally
err = openmc_tally_allocate(i_start + i - 1, &
C_CHAR_'generic' // C_NULL_CHAR)
! Get pointer to tally
associate (t => tallies(i_start + i - 1) % obj)
! Get pointer to tally xml node
node_tal = node_tal_list(i)
call tally_init_from_xml(t % ptr, node_tal % ptr)
! Copy and set tally id
if (check_for_node(node_tal, "id")) then
call get_node_value(node_tal, "id", tally_id)
err = openmc_tally_set_id(i_start + i - 1, tally_id)
if (err /= 0) call fatal_error(to_f_string(openmc_err_msg))
else
call fatal_error("Must specify id for tally in tally XML file.")
end if
! Copy tally name
if (check_for_node(node_tal, "name")) &
call get_node_value(node_tal, "name", t % name)
! =======================================================================
! READ DATA FOR FILTERS
! Check if user is using old XML format and throw an error if so
if (check_for_node(node_tal, "filter")) then
call fatal_error("Tally filters must be specified independently of &
&tallies in a <filter> element. The <tally> element itself should &
&have a list of filters that apply, e.g., <filters>1 2</filters> &
&where 1 and 2 are the IDs of filters specified outside of &
&<tally>.")
end if
! Determine number of filters
if (check_for_node(node_tal, "filters")) then
n_filter = node_word_count(node_tal, "filters")
else
n_filter = 0
end if
! Allocate and store filter user ids
allocate(temp_filter(n_filter))
if (n_filter > 0) then
call get_node_array(node_tal, "filters", temp_filter)
do j = 1, n_filter
! Get pointer to filter
if (filter_dict % has(temp_filter(j))) then
i_filt = filter_dict % get(temp_filter(j))
f => filters(i_filt)
else
call fatal_error("Could not find filter " &
// trim(to_str(temp_filter(j))) // " specified on tally " &
// trim(to_str(t % id())))
end if
! Store the index of the filter
temp_filter(j) = i_filt - 1
end do
! Set the filters
err = openmc_tally_set_filters(i_start + i - 1, n_filter, temp_filter)
end if
deallocate(temp_filter)
! Check for the presence of certain filter types
has_energyout = (t % energyout_filter() > 0)
particle_filter_index = 0
do j = 1, t % n_filters()
select type (filt => filters(t % filter(j) + 1) % obj)
type is (ParticleFilter)
particle_filter_index = j
end select
end do
! Change the tally estimator if a filter demands it
do j = 1, t % n_filters()
select type (filt => filters(t % filter(j) + 1) % obj)
type is (EnergyoutFilter)
call t % set_estimator(ESTIMATOR_ANALOG)
type is (LegendreFilter)
call t % set_estimator(ESTIMATOR_ANALOG)
type is (SphericalHarmonicsFilter)
if (filt % cosine() == COSINE_SCATTER) then
call t % set_estimator(ESTIMATOR_ANALOG)
end if
type is (SpatialLegendreFilter)
call t % set_estimator(ESTIMATOR_COLLISION)
type is (ZernikeFilter)
call t % set_estimator(ESTIMATOR_COLLISION)
type is (ZernikeRadialFilter)
call t % set_estimator(ESTIMATOR_COLLISION)
end select
end do
! =======================================================================
! READ DATA FOR NUCLIDES
call tally_set_nuclides(t % ptr, node_tal % ptr)
! =======================================================================
! READ DATA FOR SCORES
call tally_set_scores(t % ptr, node_tal % ptr)
if (check_for_node(node_tal, "scores")) then
n_words = node_word_count(node_tal, "scores")
! Check if tally is compatible with particle type
if (photon_transport) then
if (particle_filter_index == 0) then
do j = 1, t % n_score_bins()
select case (t % score_bins(j))
case (SCORE_INVERSE_VELOCITY)
call fatal_error("Particle filter must be used with photon &
&transport on and inverse velocity score")
case (SCORE_FLUX, SCORE_TOTAL, SCORE_SCATTER, SCORE_NU_SCATTER, &
SCORE_ABSORPTION, SCORE_FISSION, SCORE_NU_FISSION, &
SCORE_CURRENT, SCORE_EVENTS, SCORE_DELAYED_NU_FISSION, &
SCORE_PROMPT_NU_FISSION, SCORE_DECAY_RATE)
call warning("Particle filter is not used with photon transport&
& on and " // trim(to_str(t % score_bins(j))) // " score")
end select
end do
else
select type(filt => filters(particle_filter_index) % obj)
type is (ParticleFilter)
do l = 1, filt % n_bins
if (filt % particles(l) == ELECTRON .or. filt % particles(l) == POSITRON) then
call t % set_estimator(ESTIMATOR_ANALOG)
end if
end do
end select
end if
else
if (particle_filter_index > 0) then
select type(filt => filters(particle_filter_index) % obj)
type is (ParticleFilter)
do l = 1, filt % n_bins
if (filt % particles(l) /= NEUTRON) then
call warning("Particle filter other than NEUTRON used with &
&photon transport turned off. All tallies for particle &
&type " // trim(to_str(filt % particles(l))) // " will have no scores")
end if
end do
end select
end if
end if
else
call fatal_error("No <scores> specified on tally " &
// trim(to_str(t % id())) // ".")
end if
! Check for a tally derivative.
if (check_for_node(node_tal, "derivative")) then
call get_node_value(node_tal, "derivative", deriv_id)
! Find the derivative with the given id, and store it's index.
do j = 0, n_tally_derivs() - 1
deriv => tally_deriv_c(j)
if (deriv % id == deriv_id) then
call t % set_deriv(j)
! Only analog or collision estimators are supported for differential
! tallies.
if (t % estimator() == ESTIMATOR_TRACKLENGTH) then
call t % set_estimator(ESTIMATOR_COLLISION)
end if
! We found the derivative we were looking for; exit the do loop.
exit
end if
if (j == n_tally_derivs()) then
call fatal_error("Could not find derivative " &
// trim(to_str(deriv_id)) // " specified on tally " &
// trim(to_str(t % id())))
end if
end do
deriv => tally_deriv_c(t % deriv())
if (deriv % variable == DIFF_NUCLIDE_DENSITY &
.or. deriv % variable == DIFF_TEMPERATURE) then
do j = 1, t % n_nuclide_bins()
if (has_energyout .and. t % nuclide_bins(j) == -1) then
call fatal_error("Error on tally " // trim(to_str(t % id())) &
// ": Cannot use a 'nuclide_density' or 'temperature' &
&derivative on a tally with an outgoing energy filter and &
&'total' nuclide rate. Instead, tally each nuclide in the &
&material individually.")
! Note that diff tallies with these characteristics would work
! correctly if no tally events occur in the perturbed material
! (e.g. pertrubing moderator but only tallying fuel), but this
! case would be hard to check for by only reading inputs.
end if
end do
end if
end if
! If settings.xml trigger is turned on, create tally triggers
if (trigger_on) then
!TODO: off-by-one
call tally_init_triggers(t % ptr, i_start + i - 1 - 1, node_tal % ptr)
end if
! =======================================================================
! SET TALLY ESTIMATOR
! Check if user specified estimator
if (check_for_node(node_tal, "estimator")) then
temp_str = ''
call get_node_value(node_tal, "estimator", temp_str)
select case(trim(temp_str))
case ('analog')
call t % set_estimator(ESTIMATOR_ANALOG)
case ('tracklength', 'track-length', 'pathlength', 'path-length')
! If the estimator was set to an analog estimator, this means the
! tally needs post-collision information
if (t % estimator() == ESTIMATOR_ANALOG) then
call fatal_error("Cannot use track-length estimator for tally " &
// to_str(t % id()))
end if
! Set estimator to track-length estimator
call t % set_estimator(ESTIMATOR_TRACKLENGTH)
case ('collision')
! If the estimator was set to an analog estimator, this means the
! tally needs post-collision information
if (t % estimator() == ESTIMATOR_ANALOG) then
call fatal_error("Cannot use collision estimator for tally " &
// to_str(t % id()))
end if
! Set estimator to collision estimator
call t % set_estimator(ESTIMATOR_COLLISION)
case default
call fatal_error("Invalid estimator '" // trim(temp_str) &
// "' on tally " // to_str(t % id()))
end select
end if
end associate
end do READ_TALLIES
! Close XML document
call doc % clear()
end subroutine read_tallies_xml
end module input_xml

4
src/mesh.cpp
View file

@ -863,9 +863,9 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
// Non-member functions
//==============================================================================
void read_meshes(pugi::xml_node* root)
void read_meshes(pugi::xml_node root)
{
for (auto node : root->children("mesh")) {
for (auto node : root.children("mesh")) {
// Read mesh and add to vector
model::meshes.emplace_back(new RegularMesh{node});

2
src/settings.cpp
View file

@ -490,7 +490,7 @@ void read_settings_xml()
}
// Read meshes
read_meshes(&root);
read_meshes(root);
// Shannon Entropy mesh
if (check_for_node(root, "entropy_mesh")) {

6
src/tallies/derivative.cpp
View file

@ -79,14 +79,14 @@ TallyDerivative::TallyDerivative(pugi::xml_node node)
// Non-method functions
//==============================================================================
extern "C" void
read_tally_derivatives(pugi::xml_node* node)
void
read_tally_derivatives(pugi::xml_node node)
{
// Populate the derivatives array. This must be done in parallel because
// the derivatives are threadprivate.
#pragma omp parallel
{
for (auto deriv_node : node->children("derivative"))
for (auto deriv_node : node.children("derivative"))
model::tally_derivs.emplace_back(deriv_node);
}

116
src/tallies/filter.cpp
View file

@ -36,12 +36,67 @@ namespace openmc {
//==============================================================================
namespace simulation {
std::vector<FilterMatch> filter_matches;
} // namespace simulation
std::vector<FilterMatch> filter_matches;
}
namespace model {
std::vector<std::unique_ptr<Filter>> tally_filters;
} // namespace model
std::vector<std::unique_ptr<Filter>> tally_filters;
std::unordered_map<int, int> filter_map;
}
Filter*
allocate_filter(const std::string& type)
{
if (type == "azimuthal") {
model::tally_filters.push_back(std::make_unique<AzimuthalFilter>());
} else if (type == "cell") {
model::tally_filters.push_back(std::make_unique<CellFilter>());
} else if (type == "cellborn") {
model::tally_filters.push_back(std::make_unique<CellbornFilter>());
} else if (type == "cellfrom") {
model::tally_filters.push_back(std::make_unique<CellFromFilter>());
} else if (type == "distribcell") {
model::tally_filters.push_back(std::make_unique<DistribcellFilter>());
} else if (type == "delayedgroup") {
model::tally_filters.push_back(std::make_unique<DelayedGroupFilter>());
} else if (type == "energyfunction") {
model::tally_filters.push_back(std::make_unique<EnergyFunctionFilter>());
} else if (type == "energy") {
model::tally_filters.push_back(std::make_unique<EnergyFilter>());
} else if (type == "energyout") {
model::tally_filters.push_back(std::make_unique<EnergyoutFilter>());
} else if (type == "legendre") {
model::tally_filters.push_back(std::make_unique<LegendreFilter>());
} else if (type == "material") {
model::tally_filters.push_back(std::make_unique<MaterialFilter>());
} else if (type == "mesh") {
model::tally_filters.push_back(std::make_unique<MeshFilter>());
} else if (type == "meshsurface") {
model::tally_filters.push_back(std::make_unique<MeshSurfaceFilter>());
} else if (type == "mu") {
model::tally_filters.push_back(std::make_unique<MuFilter>());
} else if (type == "particle") {
model::tally_filters.push_back(std::make_unique<ParticleFilter>());
} else if (type == "polar") {
model::tally_filters.push_back(std::make_unique<PolarFilter>());
} else if (type == "surface") {
model::tally_filters.push_back(std::make_unique<SurfaceFilter>());
} else if (type == "spatiallegendre") {
model::tally_filters.push_back(std::make_unique<SpatialLegendreFilter>());
} else if (type == "sphericalharmonics") {
model::tally_filters.push_back(std::make_unique<SphericalHarmonicsFilter>());
} else if (type == "universe") {
model::tally_filters.push_back(std::make_unique<UniverseFilter>());
} else if (type == "zernike") {
model::tally_filters.push_back(std::make_unique<ZernikeFilter>());
} else if (type == "zernikeradial") {
model::tally_filters.push_back(std::make_unique<ZernikeRadialFilter>());
} else {
throw std::runtime_error{"Unknown filter type: " + type};
}
return model::tally_filters.back().get();
}
//==============================================================================
// Fortran compatibility functions
@ -50,59 +105,6 @@ std::vector<std::unique_ptr<Filter>> tally_filters;
extern "C" {
// filter_match_point moved to simulation.cpp
Filter*
allocate_filter(const char* type)
{
std::string type_ {type};
if (type_ == "azimuthal") {
model::tally_filters.push_back(std::make_unique<AzimuthalFilter>());
} else if (type_ == "cell") {
model::tally_filters.push_back(std::make_unique<CellFilter>());
} else if (type_ == "cellborn") {
model::tally_filters.push_back(std::make_unique<CellbornFilter>());
} else if (type_ == "cellfrom") {
model::tally_filters.push_back(std::make_unique<CellFromFilter>());
} else if (type_ == "distribcell") {
model::tally_filters.push_back(std::make_unique<DistribcellFilter>());
} else if (type_ == "delayedgroup") {
model::tally_filters.push_back(std::make_unique<DelayedGroupFilter>());
} else if (type_ == "energyfunction") {
model::tally_filters.push_back(std::make_unique<EnergyFunctionFilter>());
} else if (type_ == "energy") {
model::tally_filters.push_back(std::make_unique<EnergyFilter>());
} else if (type_ == "energyout") {
model::tally_filters.push_back(std::make_unique<EnergyoutFilter>());
} else if (type_ == "legendre") {
model::tally_filters.push_back(std::make_unique<LegendreFilter>());
} else if (type_ == "material") {
model::tally_filters.push_back(std::make_unique<MaterialFilter>());
} else if (type_ == "mesh") {
model::tally_filters.push_back(std::make_unique<MeshFilter>());
} else if (type_ == "meshsurface") {
model::tally_filters.push_back(std::make_unique<MeshSurfaceFilter>());
} else if (type_ == "mu") {
model::tally_filters.push_back(std::make_unique<MuFilter>());
} else if (type_ == "particle") {
model::tally_filters.push_back(std::make_unique<ParticleFilter>());
} else if (type_ == "polar") {
model::tally_filters.push_back(std::make_unique<PolarFilter>());
} else if (type_ == "surface") {
model::tally_filters.push_back(std::make_unique<SurfaceFilter>());
} else if (type_ == "spatiallegendre") {
model::tally_filters.push_back(std::make_unique<SpatialLegendreFilter>());
} else if (type_ == "sphericalharmonics") {
model::tally_filters.push_back(std::make_unique<SphericalHarmonicsFilter>());
} else if (type_ == "universe") {
model::tally_filters.push_back(std::make_unique<UniverseFilter>());
} else if (type_ == "zernike") {
model::tally_filters.push_back(std::make_unique<ZernikeFilter>());
} else if (type_ == "zernikeradial") {
model::tally_filters.push_back(std::make_unique<ZernikeRadialFilter>());
} else {
return nullptr;
}
return model::tally_filters.back().get();
}
int32_t filter_get_id(Filter* filt) {return filt->id_;}

334
src/tallies/tally.cpp
View file

@ -3,9 +3,12 @@
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/message_passing.h"
#include "openmc/mesh.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/particle.h"
#include "openmc/reaction_product.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
@ -19,6 +22,8 @@
#include "openmc/tallies/filter_legendre.h"
#include "openmc/tallies/filter_mesh.h"
#include "openmc/tallies/filter_meshsurface.h"
#include "openmc/tallies/filter_particle.h"
#include "openmc/tallies/filter_sph_harm.h"
#include "openmc/tallies/filter_surface.h"
#include "openmc/xml_interface.h"
@ -45,6 +50,7 @@ namespace model {
std::vector<int> active_collision_tallies;
std::vector<int> active_meshsurf_tallies;
std::vector<int> active_surface_tallies;
std::unordered_map<int, int> tally_map;
}
namespace simulation {
@ -439,7 +445,7 @@ Tally::set_nuclides(pugi::xml_node node)
}
void
Tally::init_triggers(pugi::xml_node node, int i_tally)
Tally::init_triggers(pugi::xml_node node)
{
for (auto trigger_node: node.children("trigger")) {
// Read the trigger type.
@ -484,19 +490,18 @@ Tally::init_triggers(pugi::xml_node node, int i_tally)
}
// Parse the trigger scores and populate the triggers_ vector.
const auto& tally {*model::tallies[i_tally]};
for (auto score_str : trigger_scores) {
if (score_str == "all") {
triggers_.reserve(triggers_.size() + tally.scores_.size());
for (auto i_score = 0; i_score < tally.scores_.size(); ++i_score) {
triggers_.reserve(triggers_.size() + this->scores_.size());
for (auto i_score = 0; i_score < this->scores_.size(); ++i_score) {
triggers_.push_back({metric, threshold, i_score});
}
} else {
int i_score = 0;
for (; i_score < tally.scores_.size(); ++i_score) {
if (reaction_name(tally.scores_[i_score]) == score_str) break;
for (; i_score < this->scores_.size(); ++i_score) {
if (reaction_name(this->scores_[i_score]) == score_str) break;
}
if (i_score == tally.scores_.size()) {
if (i_score == this->scores_.size()) {
std::stringstream msg;
msg << "Could not find the score \"" << score_str << "\" in tally "
<< id_ << " but it was listed in a trigger on that tally";
@ -547,6 +552,317 @@ void Tally::accumulate()
// Non-member functions
//==============================================================================
void read_tallies_xml()
{
// Check if tallies.xml exists. If not, just return since it is optional
std::string filename = settings::path_input + "tallies.xml";
if (!file_exists(filename)) return;
write_message("Reading tallies XML file...", 5);
// Parse tallies.xml file
pugi::xml_document doc;
doc.load_file(filename.c_str());
pugi::xml_node root = doc.document_element();
// Check for <assume_separate> setting
if (check_for_node(root, "assume_separate")) {
settings::assume_separate = get_node_value_bool(root, "assume_separate");
}
// Check for user meshes and allocate
read_meshes(root);
// We only need the mesh info for plotting
if (settings::run_mode == RUN_MODE_PLOTTING) return;
// Read data for tally derivatives
read_tally_derivatives(root);
// ==========================================================================
// READ FILTER DATA
// Check for user filters and allocate
for (auto node_filt : root.children("filter")) {
// Copy filter id
if (!check_for_node(node_filt, "id")) {
fatal_error("Must specify id for filter in tally XML file.");
}
int filter_id = std::stoi(get_node_value(node_filt, "id"));
// Check to make sure 'id' hasn't been used
if (model::filter_map.find(filter_id) != model::filter_map.end()) {
fatal_error("Two or more filters use the same unique ID: "
+ std::to_string(filter_id));
}
// Convert filter type to lower case
std::string s;
if (check_for_node(node_filt, "type")) {
s = get_node_value(node_filt, "type", true);
}
// Make sure bins have been set
if (s == "energy" || s == "energyout" || s == "mu" || s == "polar"
|| s == "azimuthal" || s == "mesh" || s == "meshsurface" || s == "universe"
|| s == "material" || s == "cell" || s == "distribcell" || s == "cellborn"
|| s == "cellfrom" || s == "surface" || s == "delayedgroup") {
if (!check_for_node(node_filt, "bins")) {
fatal_error("Bins not set in filter " + std::to_string(filter_id));
}
}
// Allocate according to the filter type
Filter* f = allocate_filter(s);
// Read filter data from XML
f->from_xml(node_filt);
// Set filter id
f->id_ = filter_id;
model::filter_map[filter_id] = model::tally_filters.size() - 1;
// Initialize filter
f->initialize();
}
// ==========================================================================
// READ TALLY DATA
// Check for user tallies
int n = 0;
for (auto node : root.children("tally")) ++n;
if (n == 0 && mpi::master) {
warning("No tallies present in tallies.xml file.");
}
for (auto node_tal : root.children("tally")) {
model::tallies.push_back(std::make_unique<Tally>());
auto& t {model::tallies.back()};
t->init_from_xml(node_tal);
// Copy and set tally id
if (!check_for_node(node_tal, "id")) {
fatal_error("Must specify id for tally in tally XML file.");
}
t->id_ = std::stoi(get_node_value(node_tal, "id"));
model::tally_map[t->id_] = model::tallies.size() - 1;
// Copy tally name
if (check_for_node(node_tal, "name")) {
t->name_ = get_node_value(node_tal, "name");
}
// =======================================================================
// READ DATA FOR FILTERS
// Check if user is using old XML format and throw an error if so
if (check_for_node(node_tal, "filter")) {
fatal_error("Tally filters must be specified independently of "
"tallies in a <filter> element. The <tally> element itself should "
"have a list of filters that apply, e.g., <filters>1 2</filters> "
"where 1 and 2 are the IDs of filters specified outside of "
"<tally>.");
}
// Determine number of filters
std::vector<int> filters;
if (check_for_node(node_tal, "filters")) {
filters = get_node_array<int>(node_tal, "filters");
}
// Allocate and store filter user ids
if (!filters.empty()) {
std::vector<int> filter_indices;
for (int filter_id : filters) {
// Determine if filter ID is valid
auto it = model::filter_map.find(filter_id);
if (it == model::filter_map.end()) {
fatal_error("Could not find filter " + std::to_string(filter_id)
+ " specified on tally " + std::to_string(t->id_));
}
// Store the index of the filter
filter_indices.push_back(it->second);
}
// Set the filters
t->set_filters(filter_indices.data(), filter_indices.size());
}
// Check for the presence of certain filter types
bool has_energyout = t->energyout_filter_ >= 0;
int particle_filter_index = C_NONE;
for (int j = 0; j < t->filters().size(); ++j) {
int i_filter = t->filters(j);
const auto& f = model::tally_filters[i_filter].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
if (pf) particle_filter_index = j;
// Change the tally estimator if a filter demands it
std::string filt_type = f->type();
if (filt_type == "energyout" || filt_type == "legendre") {
t->estimator_ = ESTIMATOR_ANALOG;
} else if (filt_type == "sphericalharmonics") {
auto sf = dynamic_cast<SphericalHarmonicsFilter*>(f);
if (sf->cosine_ == SphericalHarmonicsCosine::scatter) {
t->estimator_ = ESTIMATOR_ANALOG;
}
} else if (filt_type == "spatiallegendre" || filt_type == "zernike"
|| filt_type == "zernikeradial") {
t->estimator_ = ESTIMATOR_COLLISION;
}
}
// =======================================================================
// READ DATA FOR NUCLIDES
t->set_nuclides(node_tal);
// =======================================================================
// READ DATA FOR SCORES
t->set_scores(node_tal);
if (!check_for_node(node_tal, "scores")) {
fatal_error("No scores specified on tally " + std::to_string(t->id_)
+ ".");
}
// Check if tally is compatible with particle type
if (settings::photon_transport) {
if (particle_filter_index == C_NONE) {
for (int score : t->scores_) {
switch (score) {
case SCORE_INVERSE_VELOCITY:
fatal_error("Particle filter must be used with photon "
"transport on and inverse velocity score");
break;
case SCORE_FLUX:
case SCORE_TOTAL:
case SCORE_SCATTER:
case SCORE_NU_SCATTER:
case SCORE_ABSORPTION:
case SCORE_FISSION:
case SCORE_NU_FISSION:
case SCORE_CURRENT:
case SCORE_EVENTS:
case SCORE_DELAYED_NU_FISSION:
case SCORE_PROMPT_NU_FISSION:
case SCORE_DECAY_RATE:
warning("Particle filter is not used with photon transport"
" on and " + std::to_string(score) + " score.");
break;
}
}
} else {
const auto& f = model::tally_filters[particle_filter_index].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
for (int p : pf->particles_) {
if (p == static_cast<int>(ParticleType::electron) ||
p == static_cast<int>(ParticleType::positron)) {
t->estimator_ = ESTIMATOR_ANALOG;
}
}
}
} else {
if (particle_filter_index >= 0) {
const auto& f = model::tally_filters[particle_filter_index].get();
auto pf = dynamic_cast<ParticleFilter*>(f);
for (int p : pf->particles_) {
if (p != static_cast<int>(ParticleType::neutron)) {
warning("Particle filter other than NEUTRON used with photon "
"transport turned off. All tallies for particle type " +
std::to_string(p) + " will have no scores");
}
}
}
}
// Check for a tally derivative.
if (check_for_node(node_tal, "derivative")) {
int deriv_id = std::stoi(get_node_value(node_tal, "derivative"));
// Find the derivative with the given id, and store it's index.
auto it = model::tally_deriv_map.find(deriv_id);
if (it == model::tally_deriv_map.end()) {
fatal_error("Could not find derivative " + std::to_string(deriv_id)
+ " specified on tally " + std::to_string(t->id_));
}
t->deriv_ = it->second;
// Only analog or collision estimators are supported for differential
// tallies.
if (t->estimator_ == ESTIMATOR_TRACKLENGTH) {
t->estimator_ = ESTIMATOR_COLLISION;
}
const auto& deriv = model::tally_derivs[t->deriv_];
if (deriv.variable == DIFF_NUCLIDE_DENSITY
|| deriv.variable == DIFF_TEMPERATURE) {
for (int i_nuc : t->nuclides_) {
if (has_energyout && i_nuc == -1) {
fatal_error("Error on tally " + std::to_string(t->id_)
+ ": Cannot use a 'nuclide_density' or 'temperature' "
"derivative on a tally with an outgoing energy filter and "
"'total' nuclide rate. Instead, tally each nuclide in the "
"material individually.");
// Note that diff tallies with these characteristics would work
// correctly if no tally events occur in the perturbed material
// (e.g. pertrubing moderator but only tallying fuel), but this
// case would be hard to check for by only reading inputs.
}
}
}
}
// If settings.xml trigger is turned on, create tally triggers
if (settings::trigger_on) {
t->init_triggers(node_tal);
}
// =======================================================================
// SET TALLY ESTIMATOR
// Check if user specified estimator
if (check_for_node(node_tal, "estimator")) {
std::string est = get_node_value(node_tal, "estimator");
if (est == "analog") {
t->estimator_ = ESTIMATOR_ANALOG;
} else if (est == "tracklength" || est == "track-length"
|| est == "pathlength" || est == "path-length") {
// If the estimator was set to an analog estimator, this means the
// tally needs post-collision information
if (t->estimator_ == ESTIMATOR_ANALOG) {
fatal_error("Cannot use track-length estimator for tally "
+ std::to_string(t->id_));
}
// Set estimator to track-length estimator
t->estimator_ = ESTIMATOR_TRACKLENGTH;
} else if (est == "collision") {
// If the estimator was set to an analog estimator, this means the
// tally needs post-collision information
if (t->estimator_ == ESTIMATOR_ANALOG) {
fatal_error("Cannot use collision estimator for tally " +
std::to_string(t->id_));
}
// Set estimator to collision estimator
t->estimator_ = ESTIMATOR_COLLISION;
} else {
fatal_error("Invalid estimator '" + est + "' on tally " +
std::to_string(t->id_));
}
}
}
}
#ifdef OPENMC_MPI
void reduce_tally_results()
{
@ -1041,8 +1357,8 @@ extern "C" {
void tally_set_nuclides(Tally* tally, pugi::xml_node* node)
{tally->set_nuclides(*node);}
void tally_init_triggers(Tally* tally, int i_tally, pugi::xml_node* node)
{tally->init_triggers(*node, i_tally);}
// void tally_init_triggers(Tally* tally, int i_tally, pugi::xml_node* node)
// {tally->init_triggers(*node, i_tally);}
int tally_get_deriv_c(Tally* tally) {return tally->deriv_;}

10
src/tallies/tally_filter.F90
View file

@ -97,14 +97,6 @@ contains
character(:), allocatable :: type_
interface
function allocate_filter(type) result(ptr) bind(C)
import C_CHAR, C_PTR
character(kind=C_CHAR), intent(in) :: type(*)
type(C_PTR) :: ptr
end function allocate_filter
end interface
! Convert C string to Fortran string
type_ = to_f_string(type)
@ -164,7 +156,7 @@ contains
call set_errmsg("Unknown filter type: " // trim(type_))
end select
filters(index) % obj % ptr = allocate_filter(type)
!filters(index) % obj % ptr = allocate_filter(type)
if (.not. c_associated(filters(index) % obj % ptr)) then
err = E_UNASSIGNED
call set_errmsg("Could not allocate C++ tally filter")