crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge branch 'develop' into cpp-volume

Browse source
This commit is contained in:
Paul Romano 2019-02-13 13:42:20 -06:00
commit 3dfe29b84f
51 changed files with 5262 additions and 6541 deletions
Download patch file Download diff file Expand all files Collapse all files

9
CMakeLists.txt
View file

@ -346,16 +346,10 @@ add_library(libopenmc SHARED
src/tallies/tally_derivative_header.F90
src/tallies/tally_filter.F90
src/tallies/tally_filter_header.F90
src/tallies/tally_filter_delayedgroup.F90
src/tallies/tally_filter_distribcell.F90
src/tallies/tally_filter_energy.F90
src/tallies/tally_filter_mesh.F90
src/tallies/tally_filter_meshsurface.F90
src/tallies/tally_filter_particle.F90
src/tallies/tally_filter_sph_harm.F90
src/tallies/tally_header.F90
src/tallies/trigger.F90
src/tallies/trigger_header.F90
src/bank.cpp
src/bremsstrahlung.cpp
src/dagmc.cpp
@ -409,6 +403,7 @@ add_library(libopenmc SHARED
src/string_utils.cpp
src/summary.cpp
src/surface.cpp
src/tallies/derivative.cpp
src/tallies/filter.cpp
src/tallies/filter_azimuthal.cpp
src/tallies/filter_cellborn.cpp
@ -431,6 +426,8 @@ add_library(libopenmc SHARED
src/tallies/filter_universe.cpp
src/tallies/filter_zernike.cpp
src/tallies/tally.cpp
src/tallies/tally_scoring.cpp
src/tallies/trigger.cpp
src/timer.cpp
src/thermal.cpp
src/volume_calc.cpp

2
include/openmc/capi.h
View file

@ -105,7 +105,7 @@ extern "C" {
int openmc_tally_get_active(int32_t index, bool* active);
int openmc_tally_get_estimator(int32_t index, int32_t* estimator);
int openmc_tally_get_id(int32_t index, int32_t* id);
int openmc_tally_get_filters(int32_t index, int32_t** indices, int* n);
int openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n);
int openmc_tally_get_n_realizations(int32_t index, int32_t* n);
int openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n);
int openmc_tally_get_scores(int32_t index, int** scores, int* n);

28
include/openmc/constants.h
View file

@ -368,24 +368,6 @@ constexpr int NO_BIN_FOUND {-1};
constexpr int FILTER_UNIVERSE {1};
constexpr int FILTER_MATERIAL {2};
constexpr int FILTER_CELL {3};
constexpr int FILTER_CELLBORN {4};
constexpr int FILTER_SURFACE {5};
constexpr int FILTER_MESH {6};
constexpr int FILTER_ENERGYIN {7};
constexpr int FILTER_ENERGYOUT {8};
constexpr int FILTER_DISTRIBCELL {9};
constexpr int FILTER_MU {10};
constexpr int FILTER_POLAR {11};
constexpr int FILTER_AZIMUTHAL {12};
constexpr int FILTER_DELAYEDGROUP {13};
constexpr int FILTER_ENERGYFUNCTION {14};
constexpr int FILTER_CELLFROM {15};
constexpr int FILTER_MESHSURFACE {16};
constexpr int FILTER_LEGENDRE {17};
constexpr int FILTER_SPH_HARMONICS {18};
constexpr int FILTER_SPTL_LEGENDRE {19};
constexpr int FILTER_ZERNIKE {20};
constexpr int FILTER_PARTICLE {21};
// Mesh types
constexpr int MESH_REGULAR {1};
@ -404,11 +386,6 @@ constexpr int IN_BOTTOM {10}; // z min
constexpr int OUT_TOP {11}; // z max
constexpr int IN_TOP {12}; // z max
// Tally trigger types and threshold
constexpr int VARIANCE {1};
constexpr int RELATIVE_ERROR {2};
constexpr int STANDARD_DEVIATION {3};
// Global tally parameters
constexpr int N_GLOBAL_TALLIES {4};
constexpr int K_COLLISION {0};
@ -416,11 +393,6 @@ constexpr int K_ABSORPTION {1};
constexpr int K_TRACKLENGTH {2};
constexpr int LEAKAGE {3};
// Differential tally independent variables
constexpr int DIFF_DENSITY {1};
constexpr int DIFF_NUCLIDE_DENSITY {2};
constexpr int DIFF_TEMPERATURE {3};
// Miscellaneous
constexpr int C_NONE {-1};
constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE

3
include/openmc/mgxs.h
View file

@ -134,7 +134,8 @@ class Mgxs {
//! @param dg delayed group index; use nullptr if irrelevant.
//! @return Requested cross section value.
double
get_xs(int xstype, int gin, int* gout, double* mu, int* dg);
get_xs(int xstype, int gin, const int* gout, const double* mu,
const int* dg);
//! \brief Samples the fission neutron energy and if prompt or delayed.
//!

23
include/openmc/mgxs_interface.h
View file

@ -51,20 +51,21 @@ extern "C" void
calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs);
extern "C" double
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
double
get_nuclide_xs(int index, int xstype, int gin, const int* gout,
const double* mu, const int* dg);
extern "C" double
get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
inline double
get_nuclide_xs(int index, int xstype, int gin)
{return get_nuclide_xs(index, xstype, gin, nullptr, nullptr, nullptr);}
extern "C" void
set_nuclide_angle_index_c(int index, const double uvw[3]);
double
get_macro_xs(int index, int xstype, int gin, const int* gout,
const double* mu, const int* dg);
extern "C" void
set_macro_angle_index_c(int index, const double uvw[3]);
extern "C" void
set_nuclide_temperature_index_c(int index, double sqrtkT);
inline double
get_macro_xs(int index, int xstype, int gin)
{return get_macro_xs(index, xstype, gin, nullptr, nullptr, nullptr);}
//==============================================================================
// General Mgxs methods

14
include/openmc/nuclide.h
View file

@ -26,10 +26,10 @@ namespace openmc {
constexpr double CACHE_INVALID {-1.0};
//===============================================================================
//==============================================================================
//! Cached microscopic cross sections for a particular nuclide at the current
//! energy
//===============================================================================
//==============================================================================
struct NuclideMicroXS {
// Microscopic cross sections in barns
@ -64,10 +64,10 @@ struct NuclideMicroXS {
//!< * temperature (eV))
};
//===============================================================================
//==============================================================================
// MATERIALMACROXS contains cached macroscopic cross sections for the material a
// particle is traveling through
//===============================================================================
//==============================================================================
struct MaterialMacroXS {
double total; //!< macroscopic total xs
@ -83,9 +83,9 @@ struct MaterialMacroXS {
double pair_production; //!< macroscopic pair production xs
};
//===============================================================================
//==============================================================================
// Data for a nuclide
//===============================================================================
//==============================================================================
class Nuclide {
public:
@ -176,6 +176,8 @@ private:
//! Checks for the right version of nuclear data within HDF5 files
void check_data_version(hid_t file_id);
extern "C" bool multipole_in_range(const Nuclide* nuc, double E);
//==============================================================================
// Global variables
//==============================================================================

64
include/openmc/tallies/derivative.h Normal file
View file

@ -0,0 +1,64 @@
#ifndef OPENMC_TALLIES_DERIVATIVE_H
#define OPENMC_TALLIES_DERIVATIVE_H
#include "openmc/particle.h"
#include <unordered_map>
#include <vector>
#include "pugixml.hpp"
//==============================================================================
//! Describes a first-order derivative that can be applied to tallies.
//==============================================================================
namespace openmc {
struct TallyDerivative {
int id; //!< User-defined identifier
int variable; //!< Independent variable (like temperature)
int diff_material; //!< Material this derivative is applied to
int diff_nuclide; //!< Nuclide this material is applied to
double flux_deriv; //!< Derivative of the current particle's weight
TallyDerivative() {}
explicit TallyDerivative(pugi::xml_node node);
};
//==============================================================================
// Non-method functions
//==============================================================================
//! Scale the given score by its logarithmic derivative
void
apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
double atom_density, int score_bin, double& score);
} // namespace openmc
//==============================================================================
// Global variables
//==============================================================================
// Explicit vector template specialization declaration of threadprivate variable
// outside of the openmc namespace for the picky Intel compiler.
extern template class std::vector<openmc::TallyDerivative>;
namespace openmc {
namespace model {
extern std::vector<TallyDerivative> tally_derivs;
#pragma omp threadprivate(tally_derivs)
extern std::unordered_map<int, int> tally_deriv_map;
} // namespace model
// Independent variables
//TODO: convert to enum
constexpr int DIFF_DENSITY {1};
constexpr int DIFF_NUCLIDE_DENSITY {2};
constexpr int DIFF_TEMPERATURE {3};
} // namespace openmc
#endif // OPENMC_TALLIES_DERIVATIVE_H

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

@ -22,6 +22,8 @@ class FilterMatch
public:
std::vector<int> bins_;
std::vector<double> weights_;
int i_bin_;
bool bins_present_ {false};
};
} // namespace openmc
@ -70,6 +72,8 @@ public:
virtual void initialize() {}
int32_t id_;
int n_bins_;
};
@ -93,10 +97,6 @@ extern std::vector<std::unique_ptr<Filter>> tally_filters;
//==============================================================================
extern "C" void free_memory_tally_c();
//==============================================================================
// Filter-related Fortran functions that will be called from C++
extern "C" int verify_filter(int32_t index);
extern "C" Filter* filter_from_f(int32_t index);

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

@ -2,17 +2,120 @@
#define OPENMC_TALLIES_TALLY_H
#include "openmc/constants.h"
#include "openmc/tallies/trigger.h"
#include "pugixml.hpp"
#include "xtensor/xtensor.hpp"
#include <memory> // for unique_ptr
#include <string>
#include <vector>
namespace openmc {
//==============================================================================
//! A user-specified flux-weighted (or current) measurement.
//==============================================================================
class Tally {
public:
Tally() {}
void init_from_xml(pugi::xml_node node);
void set_scores(pugi::xml_node node);
void set_scores(std::vector<std::string> scores);
void set_nuclides(pugi::xml_node node);
//----------------------------------------------------------------------------
// Methods for getting and setting filter/stride data.
const std::vector<int32_t>& filters() const {return filters_;}
int32_t filters(int i) const {return filters_[i];}
void set_filters(const int32_t filter_indices[], int n);
int32_t strides(int i) const {return strides_[i];}
int32_t n_filter_bins() const {return n_filter_bins_;}
//----------------------------------------------------------------------------
// Other methods.
void init_triggers(pugi::xml_node node, int i_tally);
//----------------------------------------------------------------------------
// Major public data members.
int id_; //!< User-defined identifier
std::string name_; //!< User-defined name
int type_ {TALLY_VOLUME}; //!< e.g. volume, surface current
//! Event type that contributes to this tally
int estimator_ {ESTIMATOR_TRACKLENGTH};
//! Whether this tally is currently being updated
bool active_ {false};
std::vector<int> scores_; //!< Filter integrands (e.g. flux, fission)
//! Index of each nuclide to be tallied. -1 indicates total material.
std::vector<int> nuclides_ {-1};
//! True if this tally has a bin for every nuclide in the problem
bool all_nuclides_ {false};
//----------------------------------------------------------------------------
// Miscellaneous public members.
// We need to have quick access to some filters. The following gives indices
// for various filters that could be in the tally or C_NONE if they are not
// present.
int energyout_filter_ {C_NONE};
int delayedgroup_filter_ {C_NONE};
bool depletion_rx_ {false}; //!< Has depletion reactions (e.g. (n,2n))
std::vector<Trigger> triggers_;
int deriv_ {C_NONE}; //!< Index of a TallyDerivative object for diff tallies.
private:
//----------------------------------------------------------------------------
// Private data.
std::vector<int32_t> filters_; //!< Filter indices in global filters array
//! Index strides assigned to each filter to support 1D indexing.
std::vector<int32_t> strides_;
int32_t n_filter_bins_ {0};
};
//==============================================================================
// Global variable declarations
//==============================================================================
extern "C" double total_weight;
namespace model {
extern std::vector<std::unique_ptr<Tally>> tallies;
extern std::vector<int> active_tallies;
extern std::vector<int> active_analog_tallies;
extern std::vector<int> active_tracklength_tallies;
extern std::vector<int> active_collision_tallies;
extern std::vector<int> active_meshsurf_tallies;
extern std::vector<int> active_surface_tallies;
} // namespace model
// Threadprivate variables
extern "C" double global_tally_absorption;
extern "C" double global_tally_collision;
@ -44,6 +147,8 @@ adaptor_type<3> tally_results(int idx);
extern "C" void reduce_tally_results();
#endif
extern "C" void free_memory_tally_c();
} // namespace openmc
#endif // OPENMC_TALLIES_TALLY_H

55
include/openmc/tallies/tally_scoring.h Normal file
View file

@ -0,0 +1,55 @@
#ifndef OPENMC_TALLIES_TALLY_SCORING_H
#define OPENMC_TALLIES_TALLY_SCORING_H
#include "openmc/particle.h"
#include "openmc/tallies/tally.h"
namespace openmc {
//==============================================================================
//! An iterator over all combinations of a tally's matching filter bins.
//
//! This iterator handles two distinct tasks. First, it maps the N-dimensional
//! space created by the indices of N filters onto a 1D sequence. In other
//! words, it provides a single number that uniquely identifies a combination of
//! bins for many filters. Second, it handles the task of finding each valid
//! combination of filter bins given that each filter can have 1 or 2 or many
//! bins that are valid for the current tally event.
//==============================================================================
class FilterBinIter
{
public:
//! Construct an iterator over bins that match a given particle's state.
FilterBinIter(const Tally& tally, const Particle* p);
//! Construct an iterator over all filter bin combinations.
//
//! \param end if true, the returned iterator indicates the end of a loop.
FilterBinIter(const Tally& tally, bool end);
bool operator==(const FilterBinIter& other) const
{return index_ == other.index_;}
bool operator!=(const FilterBinIter& other) const
{return !(*this == other);}
FilterBinIter& operator++();
int index_ {1};
double weight_ {1.};
private:
void compute_index_weight();
const Tally& tally_;
};
//==============================================================================
// Non-member functions
//==============================================================================
} // namespace openmc
#endif // OPENMC_TALLIES_TALLY_SCORING_H

50
include/openmc/tallies/trigger.h Normal file
View file

@ -0,0 +1,50 @@
#ifndef OPENMC_TALLIES_TRIGGER_H
#define OPENMC_TALLIES_TRIGGER_H
#include <string>
#include "pugixml.hpp"
namespace openmc {
//==============================================================================
// Type definitions
//==============================================================================
enum class TriggerMetric {
variance, relative_error, standard_deviation, not_active
};
//! Stops the simulation early if a desired tally uncertainty is reached.
struct Trigger {
TriggerMetric metric; //!< The type of uncertainty (e.g. std dev) measured
double threshold; //!< Uncertainty value below which trigger is satisfied
int score_index; //!< Index of the relevant score in the tally's arrays
};
//! Stops the simulation early if a desired k-effective uncertainty is reached.
struct KTrigger
{
TriggerMetric metric {TriggerMetric::not_active};
double threshold {0.};
};
//==============================================================================
// Global variable declarations
//==============================================================================
//TODO: consider a different namespace
namespace settings {
extern KTrigger keff_trigger;
}
//==============================================================================
// Non-memeber functions
//==============================================================================
void check_triggers();
} // namespace openmc
#endif // OPENMC_TALLIES_TRIGGER_H

6
openmc/_utils.py
View file

@ -7,7 +7,7 @@ from urllib.request import urlopen
_BLOCK_SIZE = 16384
def download(url, checksum=None):
def download(url, checksum=None, **kwargs):
"""Download file from a URL
Parameters
@ -16,6 +16,7 @@ def download(url, checksum=None):
URL from which to download
checksum : str or None
MD5 checksum to check against
Keyword arguments passed to :func:urllib.request.urlopen
Returns
-------
@ -23,8 +24,7 @@ def download(url, checksum=None):
Name of file written locally
"""
req = urlopen(url)
req = urlopen(url, **kwargs)
# Get file size from header
file_size = req.length

6
openmc/capi/tally.py
View file

@ -254,13 +254,13 @@ class Tally(_FortranObjectWithID):
filt_idx = POINTER(c_int32)()
n = c_int()
_dll.openmc_tally_get_filters(self._index, filt_idx, n)
return [_get_filter(filt_idx[i]) for i in range(n.value)]
return [_get_filter(filt_idx[i]+1) for i in range(n.value)]
@filters.setter
def filters(self, filters):
# Get filter indices as int32_t[]
n = len(filters)
indices = (c_int32*n)(*(f._index for f in filters))
indices = (c_int32*n)(*(f._index-1 for f in filters))
_dll.openmc_tally_set_filters(self._index, n, indices)
@ -278,7 +278,7 @@ class Tally(_FortranObjectWithID):
nucs = POINTER(c_int)()
n = c_int()
_dll.openmc_tally_get_nuclides(self._index, nucs, n)
return [Nuclide(nucs[i]).name if nucs[i] > 0 else 'total'
return [Nuclide(nucs[i]+1).name if nucs[i] >= 0 else 'total'
for i in range(n.value)]
@nuclides.setter

3
src/api.F90
View file

@ -65,10 +65,8 @@ contains
use settings
use simulation_header
use surface_header
use tally_derivative_header
use tally_filter_header
use tally_header
use trigger_header
interface
subroutine free_memory_source() bind(C)
@ -109,7 +107,6 @@ contains
call free_memory_mesh()
call free_memory_tally()
call free_memory_tally_filter()
call free_memory_tally_derivative()
call free_memory_bank()
call free_memory_cmfd()
#ifdef DAGMC

28
src/constants.F90
View file

@ -201,30 +201,10 @@ module constants
integer, parameter :: NO_BIN_FOUND = -1
! Tally filter and map types
integer, parameter :: N_FILTER_TYPES = 22
integer, parameter :: &
FILTER_UNIVERSE = 1, &
FILTER_MATERIAL = 2, &
FILTER_CELL = 3, &
FILTER_CELLBORN = 4, &
FILTER_SURFACE = 5, &
FILTER_MESH = 6, &
FILTER_ENERGYIN = 7, &
FILTER_ENERGYOUT = 8, &
FILTER_DISTRIBCELL = 9, &
FILTER_MU = 10, &
FILTER_POLAR = 11, &
FILTER_AZIMUTHAL = 12, &
FILTER_DELAYEDGROUP = 13, &
FILTER_ENERGYFUNCTION = 14, &
FILTER_CELLFROM = 15, &
FILTER_MESHSURFACE = 16, &
FILTER_LEGENDRE = 17, &
FILTER_SPH_HARMONICS = 18, &
FILTER_SPTL_LEGENDRE = 19, &
FILTER_ZERNIKE = 20, &
FILTER_ZERNIKE_RADIAL = 21, &
FILTER_PARTICLE = 22
FILTER_CELL = 3
! Tally surface current directions
integer, parameter :: &
@ -241,12 +221,6 @@ module constants
OUT_TOP = 11, & ! z max
IN_TOP = 12 ! z max
! Tally trigger types and threshold
integer, parameter :: &
VARIANCE = 1, &
RELATIVE_ERROR = 2, &
STANDARD_DEVIATION = 3
! Global tally parameters
integer, parameter :: N_GLOBAL_TALLIES = 4
integer, parameter :: &

699
src/input_xml.F90
View file

@ -29,7 +29,6 @@ module input_xml
use tally_derivative_header
use tally_filter_header
use tally_filter
use trigger_header
use xml_interface
implicit none
@ -414,53 +413,67 @@ contains
integer :: i ! loop over user-specified tallies
integer :: j ! loop over words
integer :: k ! another loop index
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 :: i_elem ! index of entry in dictionary
integer :: n ! size of arrays in mesh specification
integer :: n_words ! number of words read
integer :: n_filter ! number of filters
integer :: n_scores ! number of scores
integer :: n_user_trig ! number of user-specified tally triggers
integer :: trig_ind ! index of triggers array for each tally
integer :: user_trig_ind ! index of user-specified triggers for each tally
integer :: i_start, i_end
integer(C_INT) :: err
real(8) :: threshold ! trigger convergence threshold
integer :: MT ! user-specified MT for score
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) :: word
character(MAX_WORD_LEN) :: score_name
character(MAX_WORD_LEN) :: temp_str
character(MAX_WORD_LEN), allocatable :: sarray(:)
type(DictCharInt) :: trigger_scores
type(TallyFilterContainer), pointer :: f
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_tal
type(XMLNode) :: node_filt
type(XMLNode) :: node_trigger
type(XMLNode), allocatable :: node_tal_list(:)
type(XMLNode), allocatable :: node_filt_list(:)
type(XMLNode), allocatable :: node_trigger_list(:)
type(XMLNode), allocatable :: node_deriv_list(:)
type(DictEntryCI) :: elem
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
end interface
! Check if tallies.xml exists
filename = trim(path_input) // "tallies.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
! We need to allocate tally_derivs to avoid segfaults. Also needs to be
! done in parallel because tally derivs are threadprivate.
!$omp parallel
allocate(tally_derivs(0))
!$omp end parallel
! Since a tallies.xml file is optional, no error is issued here
return
end if
@ -501,29 +514,7 @@ contains
! ==========================================================================
! READ DATA FOR DERIVATIVES
! Get pointer list to XML <derivative> nodes and allocate global array.
! The array is threadprivate so it must be allocated in parallel.
call get_node_list(root, "derivative", node_deriv_list)
!$omp parallel
allocate(tally_derivs(size(node_deriv_list)))
!$omp end parallel
! Make sure this is not an MG run.
if (.not. run_CE .and. size(node_deriv_list) > 0) then
call fatal_error("Differential tallies not supported in multi-group mode")
end if
! Read derivative attributes.
do i = 1, size(node_deriv_list)
!$omp parallel
!$omp critical (ReadTallyDeriv)
call tally_derivs(i) % from_xml(node_deriv_list(i))
!$omp end critical (ReadTallyDeriv)
!$omp end parallel
! Update tally derivative dictionary
call tally_deriv_dict % set(tally_derivs(i) % id, i)
end do
call read_tally_derivatives(root % ptr)
! ==========================================================================
! READ FILTER DATA
@ -611,6 +602,8 @@ contains
! 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)
@ -656,11 +649,11 @@ contains
else
call fatal_error("Could not find filter " &
// trim(to_str(temp_filter(j))) // " specified on tally " &
// trim(to_str(t % id)))
// trim(to_str(t % id())))
end if
! Store the index of the filter
temp_filter(j) = i_filt
temp_filter(j) = i_filt - 1
end do
! Set the filters
@ -668,368 +661,53 @@ contains
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
if (check_for_node(node_tal, "nuclides")) then
! Allocate a temporary string array for nuclides and copy values over
allocate(sarray(node_word_count(node_tal, "nuclides")))
call get_node_array(node_tal, "nuclides", sarray)
if (trim(sarray(1)) == 'all') then
! Handle special case <nuclides>all</nuclides>
allocate(t % nuclide_bins(n_nuclides + 1))
! Set bins to 1, 2, 3, ..., n_nuclides, -1
t % nuclide_bins(1:n_nuclides) = &
(/ (j, j=1, n_nuclides) /)
t % nuclide_bins(n_nuclides + 1) = -1
! Set number of nuclide bins
t % n_nuclide_bins = n_nuclides + 1
! Set flag so we can treat this case specially
t % all_nuclides = .true.
else
! Any other case, e.g. <nuclides>U-235 Pu-239</nuclides>
n_words = node_word_count(node_tal, "nuclides")
allocate(t % nuclide_bins(n_words))
do j = 1, n_words
! Check if total material was specified
if (trim(sarray(j)) == 'total') then
t % nuclide_bins(j) = -1
cycle
end if
! If a specific nuclide was specified
word = sarray(j)
! Search through nuclides
k = nuclide_map_get(to_c_string(word))
if (k == -1) then
call fatal_error("Could not find the nuclide " &
// trim(word) // " specified in tally " &
// trim(to_str(t % id)) // " in any material.")
end if
! Set bin to index in nuclides array
t % nuclide_bins(j) = k
end do
! Set number of nuclide bins
t % n_nuclide_bins = n_words
end if
! Deallocate temporary string array
deallocate(sarray)
else
! No <nuclides> were specified -- create only one bin will be added
! for the total material.
allocate(t % nuclide_bins(1))
t % nuclide_bins(1) = -1
t % n_nuclide_bins = 1
end if
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")
allocate(sarray(n_words))
call get_node_array(node_tal, "scores", sarray)
! Append the score to the list of possible trigger scores
do j = 1, n_words
sarray(j) = to_lower(sarray(j))
score_name = trim(sarray(j))
if (trigger_on) call trigger_scores % set(trim(score_name), j)
end do
n_scores = n_words
! Allocate score storage accordingly
allocate(t % score_bins(n_scores))
! Check the validity of the scores and their filters
do j = 1, n_scores
score_name = sarray(j)
! Check if delayed group filter is used with any score besides
! delayed-nu-fission or decay-rate
if ((score_name /= 'delayed-nu-fission' .and. &
score_name /= 'decay-rate') .and. &
t % find_filter(FILTER_DELAYEDGROUP) > 0) then
call fatal_error("Cannot tally " // trim(score_name) // " with a &
&delayedgroup filter.")
end if
select case (trim(score_name))
case ('flux')
! Prohibit user from tallying flux for an individual nuclide
if (.not. (t % n_nuclide_bins == 1 .and. &
t % nuclide_bins(1) == -1)) then
call fatal_error("Cannot tally flux for an individual nuclide.")
end if
t % score_bins(j) = SCORE_FLUX
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
call fatal_error("Cannot tally flux with an outgoing energy &
&filter.")
end if
case ('total', '(n,total)')
t % score_bins(j) = SCORE_TOTAL
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
call fatal_error("Cannot tally total reaction rate with an &
&outgoing energy filter.")
end if
case ('scatter')
t % score_bins(j) = SCORE_SCATTER
if (t % find_filter(FILTER_ENERGYOUT) > 0 .or. &
t % find_filter(FILTER_LEGENDRE) > 0) then
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
case ('nu-scatter')
t % score_bins(j) = SCORE_NU_SCATTER
! Set tally estimator to analog for CE mode
! (MG mode has all data available without a collision being
! necessary)
if (run_CE) then
t % estimator = ESTIMATOR_ANALOG
else
if (t % find_filter(FILTER_ENERGYOUT) > 0 .or. &
t % find_filter(FILTER_LEGENDRE) > 0) then
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
end if
case ('n2n', '(n,2n)')
t % score_bins(j) = N_2N
t % depletion_rx = .true.
case ('n3n', '(n,3n)')
t % score_bins(j) = N_3N
t % depletion_rx = .true.
case ('n4n', '(n,4n)')
t % score_bins(j) = N_4N
t % depletion_rx = .true.
case ('absorption')
t % score_bins(j) = SCORE_ABSORPTION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
call fatal_error("Cannot tally absorption rate with an outgoing &
&energy filter.")
end if
case ('fission', '18')
t % score_bins(j) = SCORE_FISSION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
call fatal_error("Cannot tally fission rate with an outgoing &
&energy filter.")
end if
case ('nu-fission')
t % score_bins(j) = SCORE_NU_FISSION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
case ('decay-rate')
t % score_bins(j) = SCORE_DECAY_RATE
case ('delayed-nu-fission')
t % score_bins(j) = SCORE_DELAYED_NU_FISSION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
case ('prompt-nu-fission')
t % score_bins(j) = SCORE_PROMPT_NU_FISSION
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Set tally estimator to analog
t % estimator = ESTIMATOR_ANALOG
end if
case ('kappa-fission')
t % score_bins(j) = SCORE_KAPPA_FISSION
case ('inverse-velocity')
t % score_bins(j) = SCORE_INVERSE_VELOCITY
case ('fission-q-prompt')
t % score_bins(j) = SCORE_FISS_Q_PROMPT
case ('fission-q-recoverable')
t % score_bins(j) = SCORE_FISS_Q_RECOV
case ('current')
! Check which type of current is desired: mesh currents or
! surface currents
if (t % find_filter(FILTER_SURFACE) > 0 .or. &
&t % find_filter(FILTER_CELL) > 0 .or. &
&t % find_filter(FILTER_CELLFROM) > 0) then
! Check to make sure that mesh surface currents are not desired as well
if (t % find_filter(FILTER_MESHSURFACE) > 0) then
call fatal_error("Cannot tally mesh surface currents &
&in the same tally as normal surface currents")
end if
t % type = TALLY_SURFACE
t % score_bins(j) = SCORE_CURRENT
else if (t % find_filter(FILTER_MESHSURFACE) > 0) then
t % score_bins(j) = SCORE_CURRENT
t % type = TALLY_MESH_SURFACE
! Check to make sure that current is the only desired response
! for this tally
if (n_words > 1) then
call fatal_error("Cannot tally other scores in the &
&same tally as surface currents")
end if
else
call fatal_error("Cannot tally currents without surface &
&type filters")
end if
case ('events')
t % score_bins(j) = SCORE_EVENTS
case ('elastic', '(n,elastic)')
t % score_bins(j) = ELASTIC
case ('(n,2nd)')
t % score_bins(j) = N_2ND
case ('(n,na)')
t % score_bins(j) = N_2NA
case ('(n,n3a)')
t % score_bins(j) = N_N3A
case ('(n,2na)')
t % score_bins(j) = N_2NA
case ('(n,3na)')
t % score_bins(j) = N_3NA
case ('(n,np)')
t % score_bins(j) = N_NP
case ('(n,n2a)')
t % score_bins(j) = N_N2A
case ('(n,2n2a)')
t % score_bins(j) = N_2N2A
case ('(n,nd)')
t % score_bins(j) = N_ND
case ('(n,nt)')
t % score_bins(j) = N_NT
case ('(n,nHe-3)')
t % score_bins(j) = N_N3HE
case ('(n,nd2a)')
t % score_bins(j) = N_ND2A
case ('(n,nt2a)')
t % score_bins(j) = N_NT2A
case ('(n,3nf)')
t % score_bins(j) = N_3NF
case ('(n,2np)')
t % score_bins(j) = N_2NP
case ('(n,3np)')
t % score_bins(j) = N_3NP
case ('(n,n2p)')
t % score_bins(j) = N_N2P
case ('(n,npa)')
t % score_bins(j) = N_NPA
case ('(n,n1)')
t % score_bins(j) = N_N1
case ('(n,nc)')
t % score_bins(j) = N_NC
case ('(n,gamma)')
t % score_bins(j) = N_GAMMA
t % depletion_rx = .true.
case ('(n,p)')
t % score_bins(j) = N_P
t % depletion_rx = .true.
case ('(n,d)')
t % score_bins(j) = N_D
case ('(n,t)')
t % score_bins(j) = N_T
case ('(n,3He)')
t % score_bins(j) = N_3HE
case ('(n,a)')
t % score_bins(j) = N_A
t % depletion_rx = .true.
case ('(n,2a)')
t % score_bins(j) = N_2A
case ('(n,3a)')
t % score_bins(j) = N_3A
case ('(n,2p)')
t % score_bins(j) = N_2P
case ('(n,pa)')
t % score_bins(j) = N_PA
case ('(n,t2a)')
t % score_bins(j) = N_T2A
case ('(n,d2a)')
t % score_bins(j) = N_D2A
case ('(n,pd)')
t % score_bins(j) = N_PD
case ('(n,pt)')
t % score_bins(j) = N_PT
case ('(n,da)')
t % score_bins(j) = N_DA
case default
! First look for deprecated scores
if (starts_with(trim(score_name), 'scatter-') .or. &
starts_with(trim(score_name), 'nu-scatter-') .or. &
starts_with(trim(score_name), 'total-y') .or. &
starts_with(trim(score_name), 'flux-y')) then
call fatal_error(trim(score_name) // " is no longer available.")
end if
! Assume that user has specified an MT number
MT = int(str_to_int(score_name))
if (MT /= ERROR_INT) then
! Specified score was an integer
if (MT > 1) then
t % score_bins(j) = MT
else
call fatal_error("Invalid MT on <scores>: " // trim(score_name))
end if
else
! Specified score was not an integer
call fatal_error("Unknown scoring function: " // trim(score_name))
end if
end select
! Do a check at the end (instead of for every case) to make sure
! the tallies are compatible with MG mode where we have less detailed
! nuclear data
if (.not. run_CE .and. t % score_bins(j) > 0) then
call fatal_error("Cannot tally " // trim(score_name) // &
" reaction rate in multi-group mode")
end if
end do
t % n_score_bins = n_scores
! Deallocate temporary string array of scores
deallocate(sarray)
! Check that no duplicate scores exist
do j = 1, n_scores - 1
do k = j + 1, n_scores
if (t % score_bins(j) == t % score_bins(k)) then
call fatal_error("Duplicate score of type '" // trim(&
reaction_name(t % score_bins(j))) // "' found in tally " &
// trim(to_str(t % id)))
end if
end do
end do
! Check if tally is compatible with particle type
if (photon_transport) then
if (t % find_filter(FILTER_PARTICLE) == 0) then
do j = 1, n_scores
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 &
@ -1043,18 +721,18 @@ contains
end select
end do
else
select type(filt => filters(t % find_filter(FILTER_PARTICLE)) % obj)
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
t % estimator = ESTIMATOR_ANALOG
call t % set_estimator(ESTIMATOR_ANALOG)
end if
end do
end select
end if
else
if (t % find_filter(FILTER_PARTICLE) > 0) then
select type(filt => filters(t % find_filter(FILTER_PARTICLE)) % obj)
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
@ -1068,38 +746,39 @@ contains
end if
else
call fatal_error("No <scores> specified on tally " &
// trim(to_str(t % id)) // ".")
// trim(to_str(t % id())) // ".")
end if
! Check for a tally derivative.
if (check_for_node(node_tal, "derivative")) then
! Temporarily store the derivative id.
call get_node_value(node_tal, "derivative", t % deriv)
call get_node_value(node_tal, "derivative", deriv_id)
! Find the derivative with the given id, and store it's index.
do j = 1, size(tally_derivs)
if (tally_derivs(j) % id == t % deriv) then
t % deriv = j
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
t % estimator = ESTIMATOR_COLLISION
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 == size(tally_derivs)) then
if (j == n_tally_derivs()) then
call fatal_error("Could not find derivative " &
// trim(to_str(t % deriv)) // " specified on tally " &
// trim(to_str(t % id)))
// trim(to_str(deriv_id)) // " specified on tally " &
// trim(to_str(t % id())))
end if
end do
if (tally_derivs(t % deriv) % variable == DIFF_NUCLIDE_DENSITY &
.or. tally_derivs(t % deriv) % variable == DIFF_TEMPERATURE) then
if (any(t % nuclide_bins == -1)) then
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
call fatal_error("Error on tally " // trim(to_str(t % id)) &
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 &
@ -1109,176 +788,14 @@ contains
! (e.g. pertrubing moderator but only tallying fuel), but this
! case would be hard to check for by only reading inputs.
end if
end if
end do
end if
end if
! If settings.xml trigger is turned on, create tally triggers
if (trigger_on) then
! Get list of trigger nodes for this tally
call get_node_list(node_tal, "trigger", node_trigger_list)
! Initialize the number of triggers
n_user_trig = size(node_trigger_list)
! Count the number of triggers needed for all scores including "all"
t % n_triggers = 0
COUNT_TRIGGERS: do user_trig_ind = 1, n_user_trig
! Get pointer to trigger node
node_trigger = node_trigger_list(user_trig_ind)
! Get scores for this trigger
if (check_for_node(node_trigger, "scores")) then
n_words = node_word_count(node_trigger, "scores")
allocate(sarray(n_words))
call get_node_array(node_trigger, "scores", sarray)
else
n_words = 1
allocate(sarray(n_words))
sarray(1) = "all"
end if
! Count the number of scores for this trigger
do j = 1, n_words
score_name = trim(to_lower(sarray(j)))
if (score_name == "all") then
t % n_triggers = t % n_triggers + trigger_scores % size()
else
t % n_triggers = t % n_triggers + 1
end if
end do
deallocate(sarray)
end do COUNT_TRIGGERS
! Allocate array of triggers for this tally
if (t % n_triggers > 0) then
allocate(t % triggers(t % n_triggers))
end if
! Initialize overall trigger index for this tally to zero
trig_ind = 1
! Create triggers for all scores specified on each trigger
TRIGGER_LOOP: do user_trig_ind = 1, n_user_trig
! Get pointer to trigger node
node_trigger = node_trigger_list(user_trig_ind)
! Get the trigger type - "variance", "std_dev" or "rel_err"
if (check_for_node(node_trigger, "type")) then
call get_node_value(node_trigger, "type", temp_str)
temp_str = to_lower(temp_str)
else
call fatal_error("Must specify trigger type for tally " // &
trim(to_str(t % id)) // " in tally XML file.")
end if
! Get the convergence threshold for the trigger
if (check_for_node(node_trigger, "threshold")) then
call get_node_value(node_trigger, "threshold", threshold)
else
call fatal_error("Must specify trigger threshold for tally " // &
trim(to_str(t % id)) // " in tally XML file.")
end if
! Get list scores for this trigger
if (check_for_node(node_trigger, "scores")) then
n_words = node_word_count(node_trigger, "scores")
allocate(sarray(n_words))
call get_node_array(node_trigger, "scores", sarray)
else
n_words = 1
allocate(sarray(n_words))
sarray(1) = "all"
end if
! Create a trigger for each score
SCORE_LOOP: do j = 1, n_words
score_name = trim(to_lower(sarray(j)))
! Expand "all" to include TriggerObjects for each score in tally
if (score_name == "all") then
! Loop over all tally scores
i_elem = 0
do
! Move to next score
call trigger_scores % next_entry(elem, i_elem)
if (i_elem == 0) exit
score_name = trim(elem % key)
! Store the score name and index in the trigger
t % triggers(trig_ind) % score_name = score_name
t % triggers(trig_ind) % score_index = elem % value
! Set the trigger convergence threshold type
select case (temp_str)
case ('std_dev')
t % triggers(trig_ind) % type = STANDARD_DEVIATION
case ('variance')
t % triggers(trig_ind) % type = VARIANCE
case ('rel_err')
t % triggers(trig_ind) % type = RELATIVE_ERROR
case default
call fatal_error("Unknown trigger type " // &
trim(temp_str) // " in tally " // trim(to_str(t % id)))
end select
! Store the trigger convergence threshold
t % triggers(trig_ind) % threshold = threshold
! Increment the overall trigger index
trig_ind = trig_ind + 1
end do
! Scores other than the "all" placeholder
else
! Store the score name and index
t % triggers(trig_ind) % score_name = trim(score_name)
t % triggers(trig_ind) % score_index = &
trigger_scores % get(trim(score_name))
! Check if an invalid score was set for the trigger
if (t % triggers(trig_ind) % score_index == 0) then
call fatal_error("The trigger score " // trim(score_name) // &
" is not set for tally " // trim(to_str(t % id)))
end if
! Store the trigger convergence threshold
t % triggers(trig_ind) % threshold = threshold
! Set the trigger convergence threshold type
select case (temp_str)
case ('std_dev')
t % triggers(trig_ind) % type = STANDARD_DEVIATION
case ('variance')
t % triggers(trig_ind) % type = VARIANCE
case ('rel_err')
t % triggers(trig_ind) % type = RELATIVE_ERROR
case default
call fatal_error("Unknown trigger type " // trim(temp_str) // &
" in tally " // trim(to_str(t % id)))
end select
! Increment the overall trigger index
trig_ind = trig_ind + 1
end if
end do SCORE_LOOP
! Deallocate the list of tally scores used to create triggers
deallocate(sarray)
end do TRIGGER_LOOP
! Deallocate dictionary of scores/indices used to populate triggers
call trigger_scores % clear()
!TODO: off-by-one
call tally_init_triggers(t % ptr, i_start + i - 1 - 1, node_tal % ptr)
end if
! =======================================================================
@ -1290,33 +807,33 @@ contains
call get_node_value(node_tal, "estimator", temp_str)
select case(trim(temp_str))
case ('analog')
t % estimator = ESTIMATOR_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
if (t % estimator() == ESTIMATOR_ANALOG) then
call fatal_error("Cannot use track-length estimator for tally " &
// to_str(t % id))
// to_str(t % id()))
end if
! Set estimator to track-length estimator
t % estimator = ESTIMATOR_TRACKLENGTH
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
if (t % estimator() == ESTIMATOR_ANALOG) then
call fatal_error("Cannot use collision estimator for tally " &
// to_str(t % id))
// to_str(t % id()))
end if
! Set estimator to collision estimator
t % estimator = ESTIMATOR_COLLISION
call t % set_estimator(ESTIMATOR_COLLISION)
case default
call fatal_error("Invalid estimator '" // trim(temp_str) &
// "' on tally " // to_str(t % id))
// "' on tally " // to_str(t % id()))
end select
end if

7
src/material.cpp
View file

@ -636,7 +636,7 @@ void Material::init_nuclide_index()
int n = settings::run_CE ?
data::nuclides.size() : data::nuclides_MG.size();
mat_nuclide_index_.resize(n);
std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), -1);
std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), C_NONE);
for (int i = 0; i < nuclide_.size(); ++i) {
mat_nuclide_index_[nuclide_[i]] = i;
}
@ -1135,11 +1135,6 @@ extern "C" {
return model::materials[i_mat - 1]->density_gpcc_;
}
int material_nuclide_index(int32_t i_mat, int i_nuc)
{
return model::materials[i_mat - 1]->mat_nuclide_index_[i_nuc - 1] + 1;
}
void material_calculate_xs(const Particle* p)
{
model::materials[p->material - 1]->calculate_xs(*p);

8
src/material_header.F90
View file

@ -11,7 +11,6 @@ module material_header
public :: material_id
public :: material_nuclide
public :: material_nuclide_size
public :: material_nuclide_index
public :: material_atom_density
public :: material_density_gpcc
@ -40,13 +39,6 @@ module material_header
integer(C_INT) :: n
end function
function material_nuclide_index(i_mat, i_nuc) bind(C) result(idx)
import C_INT32_T, C_INT
integer(C_INT32_T), value :: i_mat
integer(C_INT), value :: i_nuc
integer(C_INT) :: idx
end function
function material_atom_density(i_mat, idx) bind(C) result(density)
import C_INT32_T, C_INT, C_DOUBLE
integer(C_INT32_T), value :: i_mat

3
src/mgxs.cpp
View file

@ -424,7 +424,8 @@ Mgxs::combine(const std::vector<Mgxs*>& micros, const std::vector<double>& scala
//==============================================================================
double
Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
const int* dg)
{
// This method assumes that the temperature and angle indices are set
#ifdef _OPENMP

47
src/mgxs_interface.F90
View file

@ -35,53 +35,6 @@ module mgxs_interface
real(C_DOUBLE) :: awr
end function get_awr_c
function get_nuclide_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: xstype
integer(C_INT), value, intent(in) :: gin
integer(C_INT), optional, intent(in) :: gout
real(C_DOUBLE), optional, intent(in) :: mu
integer(C_INT), optional, intent(in) :: dg
real(C_DOUBLE) :: val
end function get_nuclide_xs_c
function get_macro_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: xstype
integer(C_INT), value, intent(in) :: gin
integer(C_INT), optional, intent(in) :: gout
real(C_DOUBLE), optional, intent(in) :: mu
integer(C_INT), optional, intent(in) :: dg
real(C_DOUBLE) :: val
end function get_macro_xs_c
subroutine set_nuclide_angle_index_c(index, uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE), intent(in) :: uvw(1:3)
end subroutine set_nuclide_angle_index_c
subroutine set_macro_angle_index_c(index, uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE), intent(in) :: uvw(1:3)
end subroutine set_macro_angle_index_c
subroutine set_nuclide_temperature_index_c(index, sqrtkT) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
real(C_DOUBLE), value, intent(in) :: sqrtkT
end subroutine set_nuclide_temperature_index_c
end interface
! Number of energy groups

41
src/mgxs_interface.cpp
View file

@ -231,12 +231,13 @@ calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
//==============================================================================
double
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
get_nuclide_xs(int index, int xstype, int gin, const int* gout,
const double* mu, const int* dg)
{
int gout_c;
int* gout_c_p;
const int* gout_c_p;
int dg_c;
int* dg_c_p;
const int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
@ -255,12 +256,13 @@ get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
//==============================================================================
double
get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
get_macro_xs(int index, int xstype, int gin, const int* gout,
const double* mu, const int* dg)
{
int gout_c;
int* gout_c_p;
const int* gout_c_p;
int dg_c;
int* dg_c_p;
const int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
@ -276,33 +278,6 @@ get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
return data::macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
}
//==============================================================================
void
set_nuclide_angle_index_c(int index, const double uvw[3])
{
// Update the values
data::nuclides_MG[index - 1].set_angle_index(uvw);
}
//==============================================================================
void
set_macro_angle_index_c(int index, const double uvw[3])
{
// Update the values
data::macro_xs[index - 1].set_angle_index(uvw);
}
//==============================================================================
void
set_nuclide_temperature_index_c(int index, double sqrtkT)
{
// Update the values
data::nuclides_MG[index - 1].set_temperature_index(sqrtkT);
}
//==============================================================================
// General Mgxs methods
//==============================================================================

4
src/nuclide.cpp
View file

@ -274,7 +274,7 @@ void Nuclide::create_derived()
xs_.emplace_back(shape, 0.0);
}
reaction_index_.fill(-1);
reaction_index_.fill(C_NONE);
for (int i = 0; i < reactions_.size(); ++i) {
const auto& rx {reactions_[i]};
@ -997,7 +997,7 @@ extern "C" void multipole_deriv_eval(Nuclide* nuc, double E, double sqrtkT,
std::tie(*sig_s, *sig_a, *sig_f) = nuc->multipole_->evaluate_deriv(E, sqrtkT);
}
extern "C" bool multipole_in_range(Nuclide* nuc, double E)
extern "C" bool multipole_in_range(const Nuclide* nuc, double E)
{
return nuc->multipole_ && E >= nuc->multipole_->E_min_&&
E <= nuc->multipole_->E_max_;

250
src/output.F90
View file

@ -19,8 +19,6 @@ module output
use tally_header
use tally_derivative_header
use tally_filter
use tally_filter_mesh, only: MeshFilter
use tally_filter_header, only: TallyFilterMatch
implicit none
@ -33,6 +31,9 @@ module output
import Particle
type(Particle), intent(in) :: p
end subroutine
subroutine write_tallies() bind(C)
end subroutine
end interface
contains
@ -75,249 +76,4 @@ contains
end subroutine header
!===============================================================================
! WRITE_TALLIES creates an output file and writes out the mean values of all
! tallies and their standard deviations
!===============================================================================
subroutine write_tallies() bind(C)
integer :: i ! index in tallies array
integer :: j ! level in tally hierarchy
integer :: k ! loop index for scoring bins
integer :: n ! loop index for nuclides
integer :: h ! loop index for tally filters
integer :: indent ! number of spaces to preceed output
integer :: filter_index ! index in results array for filters
integer :: score_index ! scoring bin index
integer :: i_nuclide ! index in nuclides array
integer :: unit_tally ! tallies.out file unit
integer :: nr ! number of realizations
real(8) :: t_value ! t-values for confidence intervals
real(8) :: alpha ! significance level for CI
real(8) :: x(2) ! mean and standard deviation
character(MAX_FILE_LEN) :: filename ! name of output file
character(36) :: score_names(N_SCORE_TYPES) ! names of scoring function
character(36) :: score_name ! names of scoring function
! to be applied at write-time
integer, allocatable :: filter_bins(:)
character(MAX_WORD_LEN) :: temp_name
! Skip if there are no tallies
if (n_tallies == 0) return
allocate(filter_bins(n_filters))
! Initialize names for scores
score_names(abs(SCORE_FLUX)) = "Flux"
score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate"
score_names(abs(SCORE_SCATTER)) = "Scattering Rate"
score_names(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate"
score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate"
score_names(abs(SCORE_FISSION)) = "Fission Rate"
score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate"
score_names(abs(SCORE_KAPPA_FISSION)) = "Kappa-Fission Rate"
score_names(abs(SCORE_EVENTS)) = "Events"
score_names(abs(SCORE_DECAY_RATE)) = "Decay Rate"
score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate"
score_names(abs(SCORE_PROMPT_NU_FISSION)) = "Prompt-Nu-Fission Rate"
score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity"
score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power"
score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power"
score_names(abs(SCORE_CURRENT)) = "Current"
! Create filename for tally output
filename = trim(path_output) // "tallies.out"
! Open tally file for writing
open(FILE=filename, NEWUNIT=unit_tally, STATUS='replace', ACTION='write')
! Calculate t-value for confidence intervals
if (confidence_intervals) then
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1)
else
t_value = ONE
end if
TALLY_LOOP: do i = 1, n_tallies
associate (t => tallies(i) % obj)
nr = t % n_realizations
if (confidence_intervals) then
! Calculate t-value for confidence intervals
alpha = ONE - CONFIDENCE_LEVEL
t_value = t_percentile(ONE - alpha/TWO, nr - 1)
else
t_value = ONE
end if
! Write header block
if (t % name == "") then
call header("TALLY " // trim(to_str(t % id)), 1, unit=unit_tally)
else
call header("TALLY " // trim(to_str(t % id)) // ": " &
// trim(t % name), 1, unit=unit_tally)
endif
! Write derivative information.
if (t % deriv /= NONE) then
associate(deriv => tally_derivs(t % deriv))
select case (deriv % variable)
case (DIFF_DENSITY)
write(unit=unit_tally, fmt="(' Density derivative Material ',A)") &
to_str(deriv % diff_material)
case (DIFF_NUCLIDE_DENSITY)
write(unit=unit_tally, fmt="(' Nuclide density derivative &
&Material ',A,' Nuclide ',A)") &
trim(to_str(deriv % diff_material)), &
trim(nuclides(deriv % diff_nuclide) % name)
case (DIFF_TEMPERATURE)
write(unit=unit_tally, fmt="(' Temperature derivative Material ',&
&A)") to_str(deriv % diff_material)
case default
call fatal_error("Differential tally dependent variable for tally "&
// trim(to_str(t % id)) // " not defined in output.F90.")
end select
end associate
end if
! WARNING: Admittedly, the logic for moving for printing results is
! extremely confusing and took quite a bit of time to get correct. The
! logic is structured this way since it is not practical to have a do
! loop for each filter variable (given that only a few filters are likely
! to be used for a given tally.
! Initialize bins, filter level, and indentation
do h = 1, size(t % filter)
filter_bins(t % filter(h)) = 0
end do
j = 1
indent = 0
print_bin: do
find_bin: do
! Check for no filters
if (size(t % filter) == 0) exit find_bin
! Increment bin combination
filter_bins(t % filter(j)) = filter_bins(t % filter(j)) + 1
! =================================================================
! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER
if (filter_bins(t % filter(j)) > &
filters(t % filter(j)) % obj % n_bins) then
! If this is the first filter, then exit
if (j == 1) exit print_bin
filter_bins(t % filter(j)) = 0
j = j - 1
indent = indent - 2
! =================================================================
! VALID BIN -- WRITE FILTER INFORMATION OR EXIT TO WRITE RESULTS
else
! Check if this is last filter
if (j == size(t % filter)) exit find_bin
! Print current filter information
write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
trim(filters(t % filter(j)) % obj % &
text_label(filter_bins(t % filter(j))))
indent = indent + 2
j = j + 1
end if
end do find_bin
! Print filter information
if (size(t % filter) > 0) then
write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
trim(filters(t % filter(j)) % obj % &
text_label(filter_bins(t % filter(j))))
end if
! Determine scoring index for this bin combination -- note that unlike
! in the score_tally subroutine, we have to use max(bins,1) since all
! bins below the lowest filter level will be zeros
filter_index = 1
do h = 1, size(t % filter)
filter_index = filter_index &
+ (max(filter_bins(t % filter(h)) ,1) - 1) * t % stride(h)
end do
! Write results for this filter bin combination
score_index = 0
if (size(t % filter) > 0) indent = indent + 2
do n = 1, t % n_nuclide_bins
! Write label for nuclide
i_nuclide = t % nuclide_bins(n)
if (i_nuclide == -1) then
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
"Total Material"
else
if (run_CE) then
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(nuclides(i_nuclide) % name)
else
call get_name_c(i_nuclide, len(temp_name), temp_name)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), &
trim(temp_name)
end if
end if
indent = indent + 2
do k = 1, t % n_score_bins
score_index = score_index + 1
associate(r => t % results(RESULT_SUM:RESULT_SUM_SQ, :, :))
if (t % score_bins(k) > 0) then
score_name = reaction_name(t % score_bins(k))
else
score_name = score_names(abs(t % score_bins(k)))
end if
x(:) = mean_stdev(r(:, score_index, filter_index), nr)
write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(x(1)), trim(to_str(t_value * x(2)))
end associate
end do
indent = indent - 2
end do
indent = indent - 2
if (size(t % filter) == 0) exit print_bin
end do print_bin
end associate
end do TALLY_LOOP
close(UNIT=unit_tally)
end subroutine write_tallies
!===============================================================================
! MEAN_STDEV computes the sample mean and standard deviation of the mean of a
! single tally score
!===============================================================================
pure function mean_stdev(result_, n) result(x)
real(8), intent(in) :: result_(2) ! sum and sum-of-squares
integer, intent(in) :: n ! number of realizations
real(8) :: x(2) ! mean and standard deviation
x(1) = result_(1) / n
if (n > 1) then
x(2) = sqrt((result_(2) / n - x(1)*x(1))/(n - 1))
else
x(2) = ZERO
end if
end function mean_stdev
end module output

165
src/output.cpp
View file

@ -6,7 +6,9 @@
#include <iomanip> // for setw, setprecision, put_time
#include <ios> // for fixed, scientific, left
#include <iostream>
#include <fstream>
#include <sstream>
#include <unordered_map>
#include <utility> // for pair
#include <omp.h>
@ -21,12 +23,18 @@
#include "openmc/lattice.h"
#include "openmc/math_functions.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/plot.h"
#include "openmc/reaction.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/surface.h"
#include "openmc/timer.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/timer.h"
namespace openmc {
@ -87,8 +95,8 @@ void title()
//==============================================================================
void
header(const char* msg, int level) {
std::string
header(const char* msg) {
// Determine how many times to repeat the '=' character.
int n_prefix = (63 - strlen(msg)) / 2;
int n_suffix = n_prefix;
@ -105,10 +113,18 @@ header(const char* msg, int level) {
out << "> " << upper << " <";
for (int i = 0; i < n_suffix; i++) out << '=';
return out.str();
}
std::string header(const std::string& msg) {return header(msg.c_str());}
void
header(const char* msg, int level) {
auto out = header(msg);
// Print header based on verbosity level.
if (settings::verbosity >= level) {
std::cout << '\n' << out.str() << "\n\n";
}
if (settings::verbosity >= level)
std::cout << '\n' << out << "\n\n";
}
//==============================================================================
@ -502,7 +518,8 @@ void print_runtime()
//==============================================================================
std::pair<double, double> mean_stdev(const double* x, int n)
std::pair<double, double>
mean_stdev(const double* x, int n)
{
double mean = x[RESULT_SUM] / n;
double stdev = n > 1 ? std::sqrt((x[RESULT_SUM_SQ]/n
@ -580,4 +597,138 @@ void print_results()
std::cout.flags(f);
}
//==============================================================================
const std::unordered_map<int, const char*> score_names = {
{SCORE_FLUX, "Flux"},
{SCORE_TOTAL, "Total Reaction Rate"},
{SCORE_SCATTER, "Scattering Rate"},
{SCORE_NU_SCATTER, "Scattering Production Rate"},
{SCORE_ABSORPTION, "Absorption Rate"},
{SCORE_FISSION, "Fission Rate"},
{SCORE_NU_FISSION, "Nu-Fission Rate"},
{SCORE_KAPPA_FISSION, "Kappa-Fission Rate"},
{SCORE_EVENTS, "Events"},
{SCORE_DECAY_RATE, "Decay Rate"},
{SCORE_DELAYED_NU_FISSION, "Delayed-Nu-Fission Rate"},
{SCORE_PROMPT_NU_FISSION, "Prompt-Nu-Fission Rate"},
{SCORE_INVERSE_VELOCITY, "Flux-Weighted Inverse Velocity"},
{SCORE_FISS_Q_PROMPT, "Prompt fission power"},
{SCORE_FISS_Q_RECOV, "Recoverable fission power"},
{SCORE_CURRENT, "Current"},
};
//! Create an ASCII output file showing all tally results.
extern "C" void
write_tallies()
{
if (model::tallies.empty()) return;
// Open the tallies.out file.
std::ofstream tallies_out;
tallies_out.open("tallies.out", std::ios::out | std::ios::trunc);
tallies_out << std::setprecision(6);
// Loop over each tally.
for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
const auto& tally {*model::tallies[i_tally]};
auto results = tally_results(i_tally);
// TODO: get this directly from the tally object when it's been translated
int32_t n_realizations;
auto err = openmc_tally_get_n_realizations(i_tally+1, &n_realizations);
// Calculate t-value for confidence intervals
double t_value = 1;
if (settings::confidence_intervals) {
auto alpha = 1 - CONFIDENCE_LEVEL;
t_value = t_percentile(1 - alpha*0.5, n_realizations - 1);
}
// Write header block.
std::string tally_header("TALLY " + std::to_string(tally.id_));
if (!tally.name_.empty()) tally_header += ": " + tally.name_;
tallies_out << header(tally_header) << "\n\n";
// Write derivative information.
if (tally.deriv_ != C_NONE) {
const auto& deriv {model::tally_derivs[tally.deriv_]};
switch (deriv.variable) {
case DIFF_DENSITY:
tallies_out << " Density derivative Material "
<< std::to_string(deriv.diff_material) << "\n";
break;
case DIFF_NUCLIDE_DENSITY:
tallies_out << " Nuclide density derivative Material "
<< std::to_string(deriv.diff_material) << " Nuclide "
// TODO: off-by-one
<< data::nuclides[deriv.diff_nuclide-1]->name_ << "\n";
break;
case DIFF_TEMPERATURE:
tallies_out << " Temperature derivative Material "
<< std::to_string(deriv.diff_material) << "\n";
break;
default:
fatal_error("Differential tally dependent variable for tally "
+ std::to_string(tally.id_) + " not defined in output.cpp");
}
}
// Loop over all filter bin combinations.
auto filter_iter = FilterBinIter(tally, false);
auto end = FilterBinIter(tally, true);
for (; filter_iter != end; ++filter_iter) {
auto filter_index = filter_iter.index_;
// Print info about this combination of filter bins. The stride check
// prevents redundant output.
int indent = 0;
for (auto i = 0; i < tally.filters().size(); ++i) {
if ((filter_index-1) % tally.strides(i) == 0) {
auto i_filt = tally.filters(i);
const auto& filt {*model::tally_filters[i_filt]};
auto& match {simulation::filter_matches[i_filt]};
tallies_out << std::string(indent+1, ' ')
// TODO: off-by-one
<< filt.text_label(match.i_bin_+1) << "\n";
}
indent += 2;
}
// Loop over all nuclide and score combinations.
int score_index = 0;
for (auto i_nuclide : tally.nuclides_) {
// Write label for this nuclide bin.
if (i_nuclide == -1) {
tallies_out << std::string(indent+1, ' ') << "Total Material\n";
} else {
if (settings::run_CE) {
tallies_out << std::string(indent+1, ' ')
<< data::nuclides[i_nuclide]->name_ << "\n";
} else {
tallies_out << std::string(indent+1, ' ')
<< data::nuclides_MG[i_nuclide].name << "\n";
}
}
// Write the score, mean, and uncertainty.
indent += 2;
for (auto score : tally.scores_) {
std::string score_name = score > 0 ? reaction_name(score)
: score_names.at(score);
double mean, stdev;
//TODO: off-by-one
std::tie(mean, stdev) = mean_stdev(
&results(filter_index-1, score_index, 0), n_realizations);
tallies_out << std::string(indent+1, ' ') << std::left
<< std::setw(36) << score_name << " " << mean << " +/- "
<< t_value * stdev << "\n";
score_index += 1;
}
indent -= 2;
}
}
}
}
} // namespace openmc

2
src/settings.F90
View file

@ -51,8 +51,6 @@ module settings
integer(C_INT32_T), bind(C) :: gen_per_batch ! # of generations per batch
integer(C_INT), bind(C) :: n_max_batches ! max # of batches
integer(C_INT), bind(C, name='trigger_batch_interval') :: n_batch_interval ! batch interval for triggers
logical(C_BOOL), bind(C, name='trigger_predict') :: pred_batches ! predict batches for triggers
logical(C_BOOL), bind(C) :: trigger_on ! flag for turning triggers on/off
logical(C_BOOL), bind(C) :: entropy_on

14
src/settings.cpp
View file

@ -25,6 +25,7 @@
#include "openmc/simulation.h"
#include "openmc/source.h"
#include "openmc/string_utils.h"
#include "openmc/tallies/trigger.h"
#include "openmc/volume_calc.h"
#include "openmc/xml_interface.h"
@ -106,13 +107,6 @@ int verbosity {7};
double weight_cutoff {0.25};
double weight_survive {1.0};
// TODO: Move to separate file
struct KTrigger {
int type;
double threshold;
};
extern "C" KTrigger keff_trigger;
} // namespace settings
//==============================================================================
@ -161,11 +155,11 @@ void get_run_parameters(pugi::xml_node node_base)
if (check_for_node(node_keff_trigger, "type")) {
auto temp = get_node_value(node_keff_trigger, "type", true, true);
if (temp == "std_dev") {
keff_trigger.type = STANDARD_DEVIATION;
keff_trigger.metric = TriggerMetric::standard_deviation;
} else if (temp == "variance") {
keff_trigger.type = VARIANCE;
keff_trigger.metric = TriggerMetric::variance;
} else if (temp == "rel_err") {
keff_trigger.type = RELATIVE_ERROR;
keff_trigger.metric = TriggerMetric::relative_error;
} else {
fatal_error("Unrecognized keff trigger type " + temp);
}

13
src/simulation.F90
View file

@ -4,7 +4,6 @@ module simulation
use nuclide_header, only: micro_xs, n_nuclides
use photon_header, only: micro_photon_xs, n_elements
use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer
implicit none
private
@ -17,18 +16,10 @@ contains
subroutine simulation_init_f() bind(C)
integer :: i
!$omp parallel
! Allocate array for microscopic cross section cache
allocate(micro_xs(n_nuclides))
allocate(micro_photon_xs(n_elements))
! Allocate array for matching filter bins
allocate(filter_matches(n_filters))
do i = 1, n_filters
filter_matches(i) % ptr = filter_match_pointer(i - 1)
end do
!$omp end parallel
end subroutine
@ -39,12 +30,10 @@ contains
!===============================================================================
subroutine simulation_finalize_f() bind(C)
! Free up simulation-specific memory
!$omp parallel
deallocate(micro_xs, micro_photon_xs, filter_matches)
deallocate(micro_xs, micro_photon_xs)
!$omp end parallel
end subroutine
end module simulation

14
src/simulation.cpp
View file

@ -17,6 +17,7 @@
#include "openmc/timer.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/trigger.h"
#include <omp.h>
@ -28,7 +29,6 @@ namespace openmc {
// data/functions from Fortran side
extern "C" void accumulate_tallies();
extern "C" void allocate_tally_results();
extern "C" void check_triggers();
extern "C" void init_tally_routines();
extern "C" void setup_active_tallies();
extern "C" void simulation_init_f();
@ -138,11 +138,6 @@ int openmc_simulation_finalize()
simulation::time_active.stop();
simulation::time_finalize.start();
#pragma omp parallel
{
simulation::filter_matches.clear();
}
// Deallocate Fortran variables, set tallies to inactive
for (auto& mat : model::materials) {
mat->mat_nuclide_index_.clear();
@ -159,6 +154,11 @@ int openmc_simulation_finalize()
// Write tally results to tallies.out
if (settings::output_tallies && mpi::master) write_tallies();
#pragma omp parallel
{
simulation::filter_matches.clear();
}
// Deactivate all tallies
for (int i = 1; i <= n_tallies; ++i) {
openmc_tally_set_active(i, false);
@ -540,7 +540,7 @@ void calculate_work()
#ifdef OPENMC_MPI
void broadcast_results() {
// Broadcast tally results so that each process has access to results
for (int i = 1; i <= n_tallies; ++i) {
for (int i = 0; i < n_tallies; ++i) {
// Create a new datatype that consists of all values for a given filter
// bin and then use that to broadcast. This is done to minimize the
// chance of the 'count' argument of MPI_BCAST exceeding 2**31

5
src/simulation_header.F90
View file

@ -23,11 +23,6 @@ module simulation_header
integer(C_INT), bind(C) :: total_gen ! total number of generations simulated
logical(C_BOOL), bind(C) :: need_depletion_rx ! need to calculate depletion reaction rx?
! ============================================================================
! TALLY PRECISION TRIGGER VARIABLES
logical(C_BOOL), bind(C) :: satisfy_triggers ! whether triggers are satisfied
integer(C_INT64_T), bind(C) :: work ! number of particles per processor
! ============================================================================

175
src/state_point.F90
View file

@ -13,19 +13,10 @@ module state_point
use, intrinsic :: ISO_C_BINDING
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message
use hdf5_interface
use message_passing
use mgxs_interface
use nuclide_header, only: nuclides
use settings
use simulation_header
use string, only: to_str
use tally_header
use tally_filter_header
use tally_derivative_header, only: tally_derivs
implicit none
@ -38,174 +29,18 @@ contains
subroutine statepoint_write_f(file_id) bind(C)
integer(HID_T), value :: file_id
integer :: i, j
integer :: i_xs
integer, allocatable :: id_array(:)
integer(HID_T) :: tallies_group, tally_group, &
filters_group, filter_group, derivs_group, &
deriv_group
character(MAX_WORD_LEN), allocatable :: str_array(:)
character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name
integer :: i
integer(HID_T) :: tallies_group, tally_group
! Open tallies group
tallies_group = open_group(file_id, "tallies")
! Write information for derivatives.
if (size(tally_derivs) > 0) then
derivs_group = create_group(tallies_group, "derivatives")
do i = 1, size(tally_derivs)
associate(deriv => tally_derivs(i))
deriv_group = create_group(derivs_group, "derivative " &
// trim(to_str(deriv % id)))
select case (deriv % variable)
case (DIFF_DENSITY)
call write_dataset(deriv_group, "independent variable", "density")
call write_dataset(deriv_group, "material", deriv % diff_material)
case (DIFF_NUCLIDE_DENSITY)
call write_dataset(deriv_group, "independent variable", &
"nuclide_density")
call write_dataset(deriv_group, "material", deriv % diff_material)
call write_dataset(deriv_group, "nuclide", &
nuclides(deriv % diff_nuclide) % name)
case (DIFF_TEMPERATURE)
call write_dataset(deriv_group, "independent variable", &
"temperature")
call write_dataset(deriv_group, "material", deriv % diff_material)
case default
call fatal_error("Independent variable for derivative " &
// trim(to_str(deriv % id)) // " not defined in &
&state_point.F90.")
end select
call close_group(deriv_group)
end associate
end do
call close_group(derivs_group)
end if
! Write number of filters
filters_group = create_group(tallies_group, "filters")
call write_attribute(filters_group, "n_filters", n_filters)
if (n_filters > 0) then
! Write IDs of filters
allocate(id_array(n_filters))
do i = 1, n_filters
id_array(i) = filters(i) % obj % id
end do
call write_attribute(filters_group, "ids", id_array)
deallocate(id_array)
! Write filter information
FILTER_LOOP: do i = 1, n_filters
filter_group = create_group(filters_group, "filter " // &
trim(to_str(filters(i) % obj % id)))
call filters(i) % obj % to_statepoint(filter_group)
call close_group(filter_group)
end do FILTER_LOOP
end if
call close_group(filters_group)
! Write number of tallies
call write_attribute(tallies_group, "n_tallies", n_tallies)
if (n_tallies > 0) then
! Write array of tally IDs
allocate(id_array(n_tallies))
do i = 1, n_tallies
id_array(i) = tallies(i) % obj % id
end do
call write_attribute(tallies_group, "ids", id_array)
deallocate(id_array)
! Write all tally information except results
TALLY_METADATA: do i = 1, n_tallies
! Get pointer to tally
associate (tally => tallies(i) % obj)
tally_group = create_group(tallies_group, "tally " // &
trim(to_str(tally % id)))
! Write the name for this tally
call write_dataset(tally_group, "name", tally % name)
select case(tally % estimator)
case (ESTIMATOR_ANALOG)
call write_dataset(tally_group, "estimator", "analog")
case (ESTIMATOR_TRACKLENGTH)
call write_dataset(tally_group, "estimator", "tracklength")
case (ESTIMATOR_COLLISION)
call write_dataset(tally_group, "estimator", "collision")
end select
call write_dataset(tally_group, "n_realizations", &
tally % n_realizations)
call write_dataset(tally_group, "n_filters", size(tally % filter))
if (size(tally % filter) > 0) then
! Write IDs of filters
allocate(id_array(size(tally % filter)))
do j = 1, size(tally % filter)
id_array(j) = filters(tally % filter(j)) % obj % id
end do
call write_dataset(tally_group, "filters", id_array)
deallocate(id_array)
end if
! Set up nuclide bin array and then write
allocate(str_array(tally % n_nuclide_bins))
NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins
if (tally % nuclide_bins(j) > 0) then
if (run_CE) then
i_xs = index(nuclides(tally % nuclide_bins(j)) % name, '.')
if (i_xs > 0) then
str_array(j) = nuclides(tally % nuclide_bins(j)) % name(1 : i_xs-1)
else
str_array(j) = nuclides(tally % nuclide_bins(j)) % name
end if
else
call get_name_c(tally % nuclide_bins(j), len(temp_name), &
temp_name)
i_xs = index(temp_name, '.')
if (i_xs > 0) then
str_array(j) = trim(temp_name(1 : i_xs-1))
else
str_array(j) = trim(temp_name)
end if
end if
else
str_array(j) = 'total'
end if
end do NUCLIDE_LOOP
call write_dataset(tally_group, "nuclides", str_array)
deallocate(str_array)
! Write derivative information.
if (tally % deriv /= NONE) then
call write_dataset(tally_group, "derivative", &
tally_derivs(tally % deriv) % id)
end if
! Write scores.
call write_dataset(tally_group, "n_score_bins", tally % n_score_bins)
allocate(str_array(size(tally % score_bins)))
do j = 1, size(tally % score_bins)
str_array(j) = reaction_name(tally % score_bins(j))
end do
call write_dataset(tally_group, "score_bins", str_array)
deallocate(str_array)
call close_group(tally_group)
end associate
end do TALLY_METADATA
end if
if (reduce_tallies) then
! Write global tallies
call write_dataset(file_id, "global_tallies", global_tallies)
! Write tallies
if (active_tallies % size() > 0) then
if (active_tallies_size() > 0) then
! Indicate that tallies are on
call write_attribute(file_id, "tallies_present", 1)
@ -214,7 +49,7 @@ contains
associate (tally => tallies(i) % obj)
! Write sum and sum_sq for each bin
tally_group = open_group(tallies_group, "tally " &
// to_str(tally % id))
// to_str(tally % id()))
call tally % write_results_hdf5(tally_group)
call close_group(tally_group)
end associate
@ -254,7 +89,7 @@ contains
associate (t => tallies(i) % obj)
! Read sum, sum_sq, and N for each bin
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(t % id)))
trim(to_str(t % id())))
call t % read_results_hdf5(tally_group)
call read_dataset(t % n_realizations, tally_group, &
"n_realizations")

139
src/state_point.cpp
View file

@ -17,9 +17,12 @@
#include "openmc/hdf5_interface.h"
#include "openmc/mesh.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/output.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/timer.h"
@ -101,13 +104,135 @@ openmc_statepoint_write(const char* filename, bool* write_source)
write_attribute(file_id, "source_present", write_source_);
// Write out information for eigenvalue run
if (settings::run_mode == RUN_MODE_EIGENVALUE) write_eigenvalue_hdf5(file_id);
if (settings::run_mode == RUN_MODE_EIGENVALUE)
write_eigenvalue_hdf5(file_id);
// Write tally-related data
hid_t tallies_group = create_group(file_id, "tallies");
// Write meshes
meshes_to_hdf5(tallies_group);
statepoint_write_f(file_id);
// Write information for derivatives
if (!model::tally_derivs.empty()) {
hid_t derivs_group = create_group(tallies_group, "derivatives");
for (const auto& deriv : model::tally_derivs) {
hid_t deriv_group = create_group(derivs_group,
"derivative " + std::to_string(deriv.id));
write_dataset(deriv_group, "material", deriv.diff_material);
if (deriv.variable == DIFF_DENSITY) {
write_dataset(deriv_group, "independent variable", "density");
} else if (deriv.variable == DIFF_NUCLIDE_DENSITY) {
write_dataset(deriv_group, "independent variable", "nuclide_density");
//TODO: off-by-one
write_dataset(deriv_group, "nuclide",
data::nuclides[deriv.diff_nuclide-1]->name_);
} else if (deriv.variable == DIFF_TEMPERATURE) {
write_dataset(deriv_group, "independent variable", "temperature");
} else {
fatal_error("Independent variable for derivative "
+ std::to_string(deriv.id) + " not defined in state_point.cpp");
}
close_group(deriv_group);
}
close_group(derivs_group);
}
// Write information for filters
hid_t filters_group = create_group(tallies_group, "filters");
write_attribute(filters_group, "n_filters", model::tally_filters.size());
if (!model::tally_filters.empty()) {
// Write filter IDs
std::vector<int32_t> filter_ids;
filter_ids.reserve(model::tally_filters.size());
for (const auto& filt : model::tally_filters)
filter_ids.push_back(filt->id_);
write_attribute(filters_group, "ids", filter_ids);
// Write info for each filter
for (const auto& filt : model::tally_filters) {
hid_t filter_group = create_group(filters_group,
"filter " + std::to_string(filt->id_));
filt->to_statepoint(filter_group);
close_group(filter_group);
}
}
close_group(filters_group);
// Write information for tallies
write_attribute(tallies_group, "n_tallies", model::tallies.size());
if (!model::tallies.empty()) {
// Write tally IDs
std::vector<int32_t> tally_ids;
tally_ids.reserve(model::tallies.size());
for (const auto& tally : model::tallies)
tally_ids.push_back(tally->id_);
write_attribute(tallies_group, "ids", tally_ids);
// Write all tally information except results
//TODO: use these two lines when openmc_get_n_realizations isn't needed
//for (const auto& tally_ptr : model::tallies) {
// const auto& tally {*tally_ptr};
for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
const auto& tally {*model::tallies[i_tally]};
hid_t tally_group = create_group(tallies_group,
"tally " + std::to_string(tally.id_));
write_dataset(tally_group, "name", tally.name_);
if (tally.estimator_ == ESTIMATOR_ANALOG) {
write_dataset(tally_group, "estimator", "analog");
} else if (tally.estimator_ == ESTIMATOR_TRACKLENGTH) {
write_dataset(tally_group, "estimator", "tracklength");
} else if (tally.estimator_ == ESTIMATOR_COLLISION) {
write_dataset(tally_group, "estimator", "collision");
}
// TODO: get this directly from the tally object when it's moved to C++
int32_t n_realizations;
auto err = openmc_tally_get_n_realizations(i_tally+1, &n_realizations);
write_dataset(tally_group, "n_realizations", n_realizations);
// Write the ID of each filter attached to this tally
write_dataset(tally_group, "n_filters", tally.filters().size());
if (!tally.filters().empty()) {
std::vector<int32_t> filter_ids;
filter_ids.reserve(tally.filters().size());
for (auto i_filt : tally.filters())
filter_ids.push_back(model::tally_filters[i_filt]->id_);
write_dataset(tally_group, "filters", filter_ids);
}
// Write the nuclides this tally scores
std::vector<std::string> nuclides;
for (auto i_nuclide : tally.nuclides_) {
if (i_nuclide == -1) {
nuclides.push_back("total");
} else {
if (settings::run_CE) {
nuclides.push_back(data::nuclides[i_nuclide]->name_);
} else {
nuclides.push_back(data::nuclides_MG[i_nuclide].name);
}
}
}
write_dataset(tally_group, "nuclides", nuclides);
if (tally.deriv_ != C_NONE) write_dataset(tally_group, "derivative",
model::tally_derivs[tally.deriv_].id);
// Write the tally score bins
std::vector<std::string> scores;
for (auto sc : tally.scores_) scores.push_back(reaction_name(sc));
write_dataset(tally_group, "n_score_bins", scores.size());
write_dataset(tally_group, "score_bins", scores);
close_group(tally_group);
}
}
close_group(tallies_group);
statepoint_write_f(file_id);
}
// Check for the no-tally-reduction method
@ -529,10 +654,11 @@ void write_tally_results_nr(hid_t file_id)
write_dataset(file_id, "global_tallies", gt);
}
for (int i = 1; i <= n_tallies; ++i) {
for (int i = 0; i < n_tallies; ++i) {
// Skip any tallies that are not active
bool active;
openmc_tally_get_active(i, &active);
// TODO: off-by-one
openmc_tally_get_active(i+1, &active);
if (!active) continue;
if (mpi::master && !object_exists(file_id, "tallies_present")) {
@ -550,7 +676,8 @@ void write_tally_results_nr(hid_t file_id)
if (mpi::master) {
// Open group for tally
int id;
openmc_tally_get_id(i, &id);
// TODO: off-by-one
openmc_tally_get_id(i+1, &id);
std::string groupname {"tally " + std::to_string(id)};
hid_t tally_group = open_group(tallies_group, groupname.c_str());

720
src/tallies/derivative.cpp Normal file
View file

@ -0,0 +1,720 @@
#include "openmc/tallies/derivative.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/nuclide.h"
#include "openmc/settings.h"
#include "openmc/tallies/tally.h"
#include "openmc/xml_interface.h"
#include <sstream>
template class std::vector<openmc::TallyDerivative>;
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
namespace model {
std::vector<TallyDerivative> tally_derivs;
std::unordered_map<int, int> tally_deriv_map;
}
//==============================================================================
// TallyDerivative implementation
//==============================================================================
TallyDerivative::TallyDerivative(pugi::xml_node node)
{
if (check_for_node(node, "id")) {
id = std::stoi(get_node_value(node, "id"));
} else {
fatal_error("Must specify an ID for <derivative> elements in the tally "
"XML file");
}
if (id <= 0)
fatal_error("<derivative> IDs must be an integer greater than zero");
std::string variable_str = get_node_value(node, "variable");
if (variable_str == "density") {
variable = DIFF_DENSITY;
} else if (variable_str == "nuclide_density") {
variable = DIFF_NUCLIDE_DENSITY;
std::string nuclide_name = get_node_value(node, "nuclide");
bool found = false;
for (auto i = 0; i < data::nuclides.size(); ++i) {
if (data::nuclides[i]->name_ == nuclide_name) {
found = true;
//TODO: off-by-one
diff_nuclide = i + 1;
}
}
if (!found) {
std::stringstream out;
out << "Could not find the nuclide \"" << nuclide_name
<< "\" specified in derivative " << id << " in any material.";
fatal_error(out);
}
} else if (variable_str == "temperature") {
variable = DIFF_TEMPERATURE;
} else {
std::stringstream out;
out << "Unrecognized variable \"" << variable_str
<< "\" on derivative " << id;
fatal_error(out);
}
diff_material = std::stoi(get_node_value(node, "material"));
}
//==============================================================================
// Non-method functions
//==============================================================================
extern "C" 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"))
model::tally_derivs.emplace_back(deriv_node);
}
// Fill the derivative map.
for (auto i = 0; i < model::tally_derivs.size(); ++i) {
auto id = model::tally_derivs[i].id;
auto search = model::tally_deriv_map.find(id);
if (search == model::tally_deriv_map.end()) {
model::tally_deriv_map[id] = i;
} else {
fatal_error("Two or more derivatives use the same unique ID: "
+ std::to_string(id));
}
}
// Make sure derivatives were not requested for an MG run.
if (!settings::run_CE && !model::tally_derivs.empty())
fatal_error("Differential tallies not supported in multi-group mode");
}
void
apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
double atom_density, int score_bin, double& score)
{
const Tally& tally {*model::tallies[i_tally]};
if (score == 0.0) return;
// If our score was previously c then the new score is
// c * (1/f * d_f/d_p + 1/c * d_c/d_p)
// where (1/f * d_f/d_p) is the (logarithmic) flux derivative and p is the
// perturbated variable.
const auto& deriv {model::tally_derivs[tally.deriv_]};
auto flux_deriv = deriv.flux_deriv;
// Handle special cases where we know that d_c/d_p must be zero.
if (score_bin == SCORE_FLUX) {
score *= flux_deriv;
return;
} else if (p->material == MATERIAL_VOID) {
score *= flux_deriv;
return;
}
//TODO: off-by-one
const Material& material {*model::materials[p->material-1]};
if (material.id_ != deriv.diff_material) {
score *= flux_deriv;
return;
}
switch (deriv.variable) {
//============================================================================
// Density derivative:
// c = Sigma_MT
// c = sigma_MT * N
// c = sigma_MT * rho * const
// d_c / d_rho = sigma_MT * const
// (1 / c) * (d_c / d_rho) = 1 / rho
case DIFF_DENSITY:
switch (tally.estimator_) {
case ESTIMATOR_ANALOG:
case ESTIMATOR_COLLISION:
switch (score_bin) {
case SCORE_TOTAL:
case SCORE_SCATTER:
case SCORE_ABSORPTION:
case SCORE_FISSION:
case SCORE_NU_FISSION:
score *= flux_deriv + 1. / material.density_gpcc_;
break;
default:
fatal_error("Tally derivative not defined for a score on tally "
+ std::to_string(tally.id_));
}
break;
default:
fatal_error("Differential tallies are only implemented for analog and "
"collision estimators.");
}
break;
//============================================================================
// Nuclide density derivative:
// If we are scoring a reaction rate for a single nuclide then
// c = Sigma_MT_i
// c = sigma_MT_i * N_i
// d_c / d_N_i = sigma_MT_i
// (1 / c) * (d_c / d_N_i) = 1 / N_i
// If the score is for the total material (i_nuclide = -1)
// c = Sum_i(Sigma_MT_i)
// d_c / d_N_i = sigma_MT_i
// (1 / c) * (d_c / d_N) = sigma_MT_i / Sigma_MT
// where i is the perturbed nuclide.
case DIFF_NUCLIDE_DENSITY:
//TODO: off-by-one throughout on diff_nuclide
switch (tally.estimator_) {
case ESTIMATOR_ANALOG:
if (p->event_nuclide != deriv.diff_nuclide) {
score *= flux_deriv;
return;
}
switch (score_bin) {
case SCORE_TOTAL:
case SCORE_SCATTER:
case SCORE_ABSORPTION:
case SCORE_FISSION:
case SCORE_NU_FISSION:
{
// Find the index of the perturbed nuclide.
int i;
for (i = 0; i < material.nuclide_.size(); ++i)
if (material.nuclide_[i] == deriv.diff_nuclide - 1) break;
score *= flux_deriv + 1. / material.atom_density_(i);
}
break;
default:
fatal_error("Tally derivative not defined for a score on tally "
+ std::to_string(tally.id_));
}
break;
case ESTIMATOR_COLLISION:
switch (score_bin) {
case SCORE_TOTAL:
if (i_nuclide == -1 && simulation::material_xs.total > 0.0) {
score *= flux_deriv
+ simulation::micro_xs[deriv.diff_nuclide-1].total
/ simulation::material_xs.total;
} else if (i_nuclide == deriv.diff_nuclide-1
&& simulation::micro_xs[i_nuclide].total) {
score *= flux_deriv + 1. / atom_density;
} else {
score *= flux_deriv;
}
break;
case SCORE_SCATTER:
if (i_nuclide == -1 && (simulation::material_xs.total
- simulation::material_xs.absorption) > 0.0) {
score *= flux_deriv
+ (simulation::micro_xs[deriv.diff_nuclide-1].total
- simulation::micro_xs[deriv.diff_nuclide-1].absorption)
/ (simulation::material_xs.total
- simulation::material_xs.absorption);
} else if (i_nuclide == deriv.diff_nuclide-1) {
score *= flux_deriv + 1. / atom_density;
} else {
score *= flux_deriv;
}
break;
case SCORE_ABSORPTION:
if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) {
score *= flux_deriv
+ simulation::micro_xs[deriv.diff_nuclide-1].absorption
/ simulation::material_xs.absorption;
} else if (i_nuclide == deriv.diff_nuclide-1
&& simulation::micro_xs[i_nuclide].absorption) {
score *= flux_deriv + 1. / atom_density;
} else {
score *= flux_deriv;
}
break;
case SCORE_FISSION:
if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) {
score *= flux_deriv
+ simulation::micro_xs[deriv.diff_nuclide-1].fission
/ simulation::material_xs.fission;
} else if (i_nuclide == deriv.diff_nuclide-1
&& simulation::micro_xs[i_nuclide].fission) {
score *= flux_deriv + 1. / atom_density;
} else {
score *= flux_deriv;
}
break;
case SCORE_NU_FISSION:
if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) {
score *= flux_deriv
+ simulation::micro_xs[deriv.diff_nuclide-1].nu_fission
/ simulation::material_xs.nu_fission;
} else if (i_nuclide == deriv.diff_nuclide-1
&& simulation::micro_xs[i_nuclide].nu_fission) {
score *= flux_deriv + 1. / atom_density;
} else {
score *= flux_deriv;
}
break;
default:
fatal_error("Tally derivative not defined for a score on tally "
+ std::to_string(tally.id_));
}
break;
default:
fatal_error("Differential tallies are only implemented for analog and "
"collision estimators.");
}
break;
//============================================================================
// Temperature derivative:
// If we are scoring a reaction rate for a single nuclide then
// c = Sigma_MT_i
// c = sigma_MT_i * N_i
// d_c / d_T = (d_sigma_Mt_i / d_T) * N_i
// (1 / c) * (d_c / d_T) = (d_sigma_MT_i / d_T) / sigma_MT_i
// If the score is for the total material (i_nuclide = -1)
// (1 / c) * (d_c / d_T) = Sum_i((d_sigma_MT_i / d_T) * N_i) / Sigma_MT_i
// where i is the perturbed nuclide. The d_sigma_MT_i / d_T term is
// computed by multipole_deriv_eval. It only works for the resolved
// resonance range and requires multipole data.
case DIFF_TEMPERATURE:
switch (tally.estimator_) {
case ESTIMATOR_ANALOG:
{
// Find the index of the event nuclide.
int i;
for (i = 0; i < material.nuclide_.size(); ++i)
if (material.nuclide_[i] == p->event_nuclide-1) break;
const auto& nuc {*data::nuclides[p->event_nuclide-1]};
if (!multipole_in_range(&nuc, p->last_E)) {
score *= flux_deriv;
break;
}
switch (score_bin) {
case SCORE_TOTAL:
if (simulation::micro_xs[p->event_nuclide-1].total) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i)
/ simulation::material_xs.total;
} else {
score *= flux_deriv;
}
break;
case SCORE_SCATTER:
if (simulation::micro_xs[p->event_nuclide-1].total
- simulation::micro_xs[p->event_nuclide-1].absorption) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv + dsig_s * material.atom_density_(i)
/ (simulation::material_xs.total
- simulation::material_xs.absorption);
} else {
score *= flux_deriv;
}
break;
case SCORE_ABSORPTION:
if (simulation::micro_xs[p->event_nuclide-1].absorption) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv + dsig_a * material.atom_density_(i)
/ simulation::material_xs.absorption;
} else {
score *= flux_deriv;
}
break;
case SCORE_FISSION:
if (simulation::micro_xs[p->event_nuclide-1].fission) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv + dsig_f * material.atom_density_(i)
/ simulation::material_xs.fission;
} else {
score *= flux_deriv;
}
break;
case SCORE_NU_FISSION:
if (simulation::micro_xs[p->event_nuclide-1].fission) {
double nu = simulation::micro_xs[p->event_nuclide-1].nu_fission
/ simulation::micro_xs[p->event_nuclide-1].fission;
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv + nu * dsig_f * material.atom_density_(i)
/ simulation::material_xs.nu_fission;
} else {
score *= flux_deriv;
}
break;
default:
fatal_error("Tally derivative not defined for a score on tally "
+ std::to_string(tally.id_));
}
}
break;
case ESTIMATOR_COLLISION:
if (i_nuclide != -1) {
const auto& nuc {data::nuclides[i_nuclide]};
if (!multipole_in_range(nuc.get(), p->last_E)) {
score *= flux_deriv;
return;
}
}
switch (score_bin) {
case SCORE_TOTAL:
if (i_nuclide == -1 && simulation::material_xs.total > 0.0) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
&& simulation::micro_xs[i_nuc].total) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / simulation::material_xs.total;
} else if (simulation::micro_xs[i_nuclide].total) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv
+ (dsig_s + dsig_a) / simulation::micro_xs[i_nuclide].total;
} else {
score *= flux_deriv;
}
break;
case SCORE_SCATTER:
if (i_nuclide == -1 && (simulation::material_xs.total
- simulation::material_xs.absorption)) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
&& (simulation::micro_xs[i_nuc].total
- simulation::micro_xs[i_nuc].absorption)) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
cum_dsig += dsig_s * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / (simulation::material_xs.total
- simulation::material_xs.absorption);
} else if (simulation::micro_xs[i_nuclide].total
- simulation::micro_xs[i_nuclide].absorption) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total
- simulation::micro_xs[i_nuclide].absorption);
} else {
score *= flux_deriv;
}
break;
case SCORE_ABSORPTION:
if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
&& simulation::micro_xs[i_nuc].absorption) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
cum_dsig += dsig_a * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / simulation::material_xs.absorption;
} else if (simulation::micro_xs[i_nuclide].absorption) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv
+ dsig_a / simulation::micro_xs[i_nuclide].absorption;
} else {
score *= flux_deriv;
}
break;
case SCORE_FISSION:
if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
&& simulation::micro_xs[i_nuc].fission) {
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
cum_dsig += dsig_f * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / simulation::material_xs.fission;
} else if (simulation::micro_xs[i_nuclide].fission) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv
+ dsig_f / simulation::micro_xs[i_nuclide].fission;
} else {
score *= flux_deriv;
}
break;
case SCORE_NU_FISSION:
if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) {
double cum_dsig = 0;
for (auto i = 0; i < material.nuclide_.size(); ++i) {
auto i_nuc = material.nuclide_[i];
const auto& nuc {*data::nuclides[i_nuc]};
if (multipole_in_range(&nuc, p->last_E)
&& simulation::micro_xs[i_nuc].fission) {
double nu = simulation::micro_xs[i_nuc].nu_fission
/ simulation::micro_xs[i_nuc].fission;
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
cum_dsig += nu * dsig_f * material.atom_density_(i);
}
}
score *= flux_deriv + cum_dsig / simulation::material_xs.nu_fission;
} else if (simulation::micro_xs[i_nuclide].fission) {
const auto& nuc {*data::nuclides[i_nuclide]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
score *= flux_deriv
+ dsig_f / simulation::micro_xs[i_nuclide].fission;
} else {
score *= flux_deriv;
}
break;
default:
break;
}
break;
default:
fatal_error("Differential tallies are only implemented for analog and "
"collision estimators.");
}
break;
}
}
//! Adjust diff tally flux derivatives for a particle tracking event.
extern "C" void
score_track_derivative(const Particle* p, double distance)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material == MATERIAL_VOID) return;
//TODO: off-by-one
const Material& material {*model::materials[p->material-1]};
for (auto& deriv : model::tally_derivs) {
if (deriv.diff_material != material.id_) continue;
switch (deriv.variable) {
case DIFF_DENSITY:
// phi is proportional to e^(-Sigma_tot * dist)
// (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
// (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
deriv.flux_deriv -= distance * simulation::material_xs.total
/ material.density_gpcc_;
break;
case DIFF_NUCLIDE_DENSITY:
// phi is proportional to e^(-Sigma_tot * dist)
// (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
// (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
//TODO: off-by-one
deriv.flux_deriv -= distance
* simulation::micro_xs[deriv.diff_nuclide-1].total;
break;
case DIFF_TEMPERATURE:
for (auto i = 0; i < material.nuclide_.size(); ++i) {
const auto& nuc {*data::nuclides[material.nuclide_[i]]};
if (multipole_in_range(&nuc, p->last_E)) {
// phi is proportional to e^(-Sigma_tot * dist)
// (1 / phi) * (d_phi / d_T) = - (d_Sigma_tot / d_T) * dist
// (1 / phi) * (d_phi / d_T) = - N (d_sigma_tot / d_T) * dist
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->E, p->sqrtkT);
deriv.flux_deriv -= distance * (dsig_s + dsig_a)
* material.atom_density_(i);
}
}
break;
}
}
}
//! Adjust diff tally flux derivatives for a particle scattering event.
//
//! Note that this subroutine will be called after absorption events in
//! addition to scattering events, but any flux derivatives scored after an
//! absorption will never be tallied. The paricle will be killed before any
//! further tallies are scored.
extern "C" void
score_collision_derivative(const Particle* p)
{
// A void material cannot be perturbed so it will not affect flux derivatives.
if (p->material == MATERIAL_VOID) return;
//TODO: off-by-one
const Material& material {*model::materials[p->material-1]};
for (auto& deriv : model::tally_derivs) {
if (deriv.diff_material != material.id_) continue;
switch (deriv.variable) {
case DIFF_DENSITY:
// phi is proportional to Sigma_s
// (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s
// (1 / phi) * (d_phi / d_rho) = 1 / rho
deriv.flux_deriv += 1. / material.density_gpcc_;
break;
case DIFF_NUCLIDE_DENSITY:
//TODO: off-by-one throughout on diff_nuclide
if (p->event_nuclide != deriv.diff_nuclide) continue;
// Find the index in this material for the diff_nuclide.
int i;
for (i = 0; i < material.nuclide_.size(); ++i)
if (material.nuclide_[i] == deriv.diff_nuclide - 1) break;
// Make sure we found the nuclide.
if (material.nuclide_[i] != deriv.diff_nuclide - 1) {
std::stringstream err_msg;
err_msg << "Could not find nuclide "
<< data::nuclides[deriv.diff_nuclide-1]->name_ << " in material "
<< material.id_ << " for tally derivative " << deriv.id;
fatal_error(err_msg);
}
// phi is proportional to Sigma_s
// (1 / phi) * (d_phi / d_N) = (d_Sigma_s / d_N) / Sigma_s
// (1 / phi) * (d_phi / d_N) = sigma_s / Sigma_s
// (1 / phi) * (d_phi / d_N) = 1 / N
deriv.flux_deriv += 1. / material.atom_density_(i);
break;
case DIFF_TEMPERATURE:
// Loop over the material's nuclides until we find the event nuclide.
for (auto i_nuc : material.nuclide_) {
const auto& nuc {*data::nuclides[i_nuc]};
//TODO: off-by-one
if (i_nuc == p->event_nuclide - 1
&& multipole_in_range(&nuc, p->last_E)) {
// phi is proportional to Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
// (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
const auto& micro_xs {simulation::micro_xs[i_nuc]};
double dsig_s, dsig_a, dsig_f;
std::tie(dsig_s, dsig_a, dsig_f)
= nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
deriv.flux_deriv += dsig_s / (micro_xs.total - micro_xs.absorption);
// Note that this is an approximation! The real scattering cross
// section is
// Sigma_s(E'->E, uvw'->uvw) = Sigma_s(E') * P(E'->E, uvw'->uvw).
// We are assuming that d_P(E'->E, uvw'->uvw) / d_T = 0 and only
// computing d_S(E') / d_T. Using this approximation in the vicinity
// of low-energy resonances causes errors (~2-5% for PWR pincell
// eigenvalue derivatives).
}
}
break;
}
}
}
//! Set the flux derivatives on differential tallies to zero.
extern "C" void
zero_flux_derivs()
{
for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.;
}
//==============================================================================
// Fortran interop
//==============================================================================
extern "C" int n_tally_derivs() {return model::tally_derivs.size();}
extern "C" TallyDerivative*
tally_deriv_c(int i)
{
return &model::tally_derivs[i];
}
}// namespace openmc

51
src/tallies/filter.cpp
View file

@ -43,21 +43,6 @@ namespace model {
std::vector<std::unique_ptr<Filter>> tally_filters;
} // namespace model
//==============================================================================
// Non-member functions
//==============================================================================
void
free_memory_tally_c()
{
#pragma omp parallel
{
simulation::filter_matches.clear();
}
model::tally_filters.clear();
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
@ -65,38 +50,6 @@ free_memory_tally_c()
extern "C" {
// filter_match_point moved to simulation.cpp
void
filter_match_bins_push_back(FilterMatch* match, int val)
{match->bins_.push_back(val);}
void
filter_match_weights_push_back(FilterMatch* match, double val)
{match->weights_.push_back(val);}
void
filter_match_bins_clear(FilterMatch* match)
{match->bins_.clear();}
void
filter_match_weights_clear(FilterMatch* match)
{match->weights_.clear();}
int
filter_match_bins_size(FilterMatch* match)
{return match->bins_.size();}
int
filter_match_bins_data(FilterMatch* match, int indx)
{return match->bins_[indx-1];}
double
filter_match_weights_data(FilterMatch* match, int indx)
{return match->weights_[indx-1];}
void
filter_match_bins_set_data(FilterMatch* match, int indx, int val)
{match->bins_[indx-1] = val;}
Filter*
allocate_filter(const char* type)
{
@ -151,6 +104,10 @@ extern "C" {
return model::tally_filters.back().get();
}
int32_t filter_get_id(Filter* filt) {return filt->id_;}
void filter_set_id(Filter* filt, int32_t id) {filt->id_ = id;}
void filter_from_xml(Filter* filt, pugi::xml_node* node)
{filt->from_xml(*node);}

7
src/tallies/filter_delayedgroup.cpp
View file

@ -48,11 +48,4 @@ DelayedGroupFilter::text_label(int bin) const
return "Delayed Group " + std::to_string(groups_[bin-1]);
}
//==============================================================================
// Fortran interoperability
//==============================================================================
extern "C" int delayedgroup_filter_groups(DelayedGroupFilter* filt, int i)
{return filt->groups_[i-1];}
} // namespace openmc

18
src/tallies/filter_energy.cpp
View file

@ -112,24 +112,6 @@ EnergyoutFilter::text_label(int bin) const
return out.str();
}
//==============================================================================
// Fortran interoperability
//==============================================================================
extern "C" bool energy_filter_matches_transport_groups(EnergyFilter* filt)
{return filt->matches_transport_groups_;}
extern "C" int
energy_filter_search(EnergyFilter* filt, double val)
{
if (val < filt->bins_.front() || val > filt->bins_.back()) {
return -1;
} else {
//TODO: off-by-one
return lower_bound_index(filt->bins_.begin(), filt->bins_.end(), val) + 1;
}
}
//==============================================================================
// C-API functions
//==============================================================================

3733
src/tallies/tally.F90
View file

File diff suppressed because it is too large Load diff

820
src/tallies/tally.cpp
View file

@ -2,7 +2,24 @@
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/reaction_product.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/filter_cell.h"
#include "openmc/tallies/filter_cellfrom.h"
#include "openmc/tallies/filter_delayedgroup.h"
#include "openmc/tallies/filter_energy.h"
#include "openmc/tallies/filter_legendre.h"
#include "openmc/tallies/filter_mesh.h"
#include "openmc/tallies/filter_meshsurface.h"
#include "openmc/tallies/filter_surface.h"
#include "openmc/xml_interface.h"
#include "xtensor/xadapt.hpp"
#include "xtensor/xbuilder.hpp" // for empty_like
@ -10,6 +27,8 @@
#include <array>
#include <cstddef>
#include <sstream>
#include <string>
namespace openmc {
@ -17,11 +36,473 @@ namespace openmc {
// Global variable definitions
//==============================================================================
namespace model {
std::vector<std::unique_ptr<Tally>> tallies;
std::vector<int> active_tallies;
std::vector<int> active_analog_tallies;
std::vector<int> active_tracklength_tallies;
std::vector<int> active_collision_tallies;
std::vector<int> active_meshsurf_tallies;
std::vector<int> active_surface_tallies;
}
double global_tally_absorption;
double global_tally_collision;
double global_tally_tracklength;
double global_tally_leakage;
int
score_str_to_int(std::string score_str)
{
if (score_str == "flux")
return SCORE_FLUX;
if (score_str == "total" || score_str == "(n,total)")
return SCORE_TOTAL;
if (score_str == "scatter")
return SCORE_SCATTER;
if (score_str == "nu-scatter")
return SCORE_NU_SCATTER;
if (score_str == "absorption")
return SCORE_ABSORPTION;
if (score_str == "fission" || score_str == "18")
return SCORE_FISSION;
if (score_str == "nu-fission")
return SCORE_NU_FISSION;
if (score_str == "decay-rate")
return SCORE_DECAY_RATE;
if (score_str == "delayed-nu-fission")
return SCORE_DELAYED_NU_FISSION;
if (score_str == "prompt-nu-fission")
return SCORE_PROMPT_NU_FISSION;
if (score_str == "kappa-fission")
return SCORE_KAPPA_FISSION;
if (score_str == "inverse-velocity")
return SCORE_INVERSE_VELOCITY;
if (score_str == "fission-q-prompt")
return SCORE_FISS_Q_PROMPT;
if (score_str == "fission-q-recoverable")
return SCORE_FISS_Q_RECOV;
if (score_str == "current")
return SCORE_CURRENT;
if (score_str == "events")
return SCORE_EVENTS;
if (score_str == "elastic" || score_str == "(n,elastic)")
return ELASTIC;
if (score_str == "n2n" || score_str == "(n,2n)")
return N_2N;
if (score_str == "n3n" || score_str == "(n,3n)")
return N_3N;
if (score_str == "n4n" || score_str == "(n,4n)")
return N_4N;
if (score_str == "(n,2nd)")
return N_2ND;
if (score_str == "(n,na)")
return N_2NA;
if (score_str == "(n,n3a)")
return N_N3A;
if (score_str == "(n,2na)")
return N_2NA;
if (score_str == "(n,3na)")
return N_3NA;
if (score_str == "(n,np)")
return N_NP;
if (score_str == "(n,n2a)")
return N_N2A;
if (score_str == "(n,2n2a)")
return N_2N2A;
if (score_str == "(n,nd)")
return N_ND;
if (score_str == "(n,nt)")
return N_NT;
if (score_str == "(n,nHe-3)")
return N_N3HE;
if (score_str == "(n,nd2a)")
return N_ND2A;
if (score_str == "(n,nt2a)")
return N_NT2A;
if (score_str == "(n,3nf)")
return N_3NF;
if (score_str == "(n,2np)")
return N_2NP;
if (score_str == "(n,3np)")
return N_3NP;
if (score_str == "(n,n2p)")
return N_N2P;
if (score_str == "(n,npa)")
return N_NPA;
if (score_str == "(n,n1)")
return N_N1;
if (score_str == "(n,nc)")
return N_NC;
if (score_str == "(n,gamma)")
return N_GAMMA;
if (score_str == "(n,p)")
return N_P;
if (score_str == "(n,d)")
return N_D;
if (score_str == "(n,t)")
return N_T;
if (score_str == "(n,3He)")
return N_3HE;
if (score_str == "(n,a)")
return N_A;
if (score_str == "(n,2a)")
return N_2A;
if (score_str == "(n,3a)")
return N_3A;
if (score_str == "(n,2p)")
return N_2P;
if (score_str == "(n,pa)")
return N_PA;
if (score_str == "(n,t2a)")
return N_T2A;
if (score_str == "(n,d2a)")
return N_D2A;
if (score_str == "(n,pd)")
return N_PD;
if (score_str == "(n,pt)")
return N_PT;
if (score_str == "(n,da)")
return N_DA;
// So far we have not identified this score string. Check to see if it is a
// deprecated score.
if (score_str.rfind("scatter-", 0) == 0
|| score_str.rfind("nu-scatter-", 0) == 0
|| score_str.rfind("total-y", 0) == 0
|| score_str.rfind("flux-y", 0) == 0)
fatal_error(score_str + " is no longer an available score");
// Assume the given string is a reaction MT number. Make sure it's a natural
// number then return.
int MT;
try {
MT = std::stoi(score_str);
} catch (const std::invalid_argument& ex) {
throw std::invalid_argument("Invalid tally score \"" + score_str + "\"");
}
if (MT < 1)
throw std::invalid_argument("Invalid tally score \"" + score_str + "\"");
return MT;
}
//==============================================================================
// Tally object implementation
//==============================================================================
void
Tally::init_from_xml(pugi::xml_node node)
{
if (check_for_node(node, "name")) name_ = get_node_value(node, "name");
}
void
Tally::set_filters(const int32_t filter_indices[], int n)
{
// Clear old data.
filters_.clear();
strides_.clear();
// Copy in the given filter indices.
filters_.assign(filter_indices, filter_indices + n);
for (int i = 0; i < n; ++i) {
auto i_filt = filters_[i];
if (i_filt < 0 || i_filt >= model::tally_filters.size())
throw std::out_of_range("Index in tally filter array out of bounds.");
// Keep track of indices for special filters.
const auto* filt = model::tally_filters[i_filt].get();
if (dynamic_cast<const EnergyoutFilter*>(filt)) {
energyout_filter_ = i;
} else if (dynamic_cast<const DelayedGroupFilter*>(filt)) {
delayedgroup_filter_ = i;
}
}
// Set the strides. Filters are traversed in reverse so that the last filter
// has the shortest stride in memory and the first filter has the longest
// stride.
strides_.resize(n, 0);
int stride = 1;
for (int i = n-1; i >= 0; --i) {
strides_[i] = stride;
stride *= model::tally_filters[filters_[i]]->n_bins_;
}
n_filter_bins_ = stride;
}
void
Tally::set_scores(pugi::xml_node node)
{
if (!check_for_node(node, "scores"))
fatal_error("No scores specified on tally " + std::to_string(id_));
auto scores = get_node_array<std::string>(node, "scores");
set_scores(scores);
}
void
Tally::set_scores(std::vector<std::string> scores)
{
// Reset state and prepare for the new scores.
scores_.clear();
depletion_rx_ = false;
scores_.reserve(scores.size());
// Check for the presence of certain restrictive filters.
bool energyout_present = energyout_filter_ != C_NONE;
bool legendre_present = false;
bool cell_present = false;
bool cellfrom_present = false;
bool surface_present = false;
bool meshsurface_present = false;
for (auto i_filt : filters_) {
const auto* filt {model::tally_filters[i_filt].get()};
if (dynamic_cast<const LegendreFilter*>(filt)) {
legendre_present = true;
} else if (dynamic_cast<const CellFromFilter*>(filt)) {
cellfrom_present = true;
} else if (dynamic_cast<const CellFilter*>(filt)) {
cell_present = true;
} else if (dynamic_cast<const SurfaceFilter*>(filt)) {
surface_present = true;
} else if (dynamic_cast<const MeshSurfaceFilter*>(filt)) {
meshsurface_present = true;
}
}
// Iterate over the given scores.
for (auto score_str : scores) {
// Make sure a delayed group filter wasn't used with an incompatible score.
bool has_delayedgroup = delayedgroup_filter_ != C_NONE;
if (delayedgroup_filter_ != C_NONE) {
if (score_str != "delayed-nu-fission" && score_str != "decay-rate")
fatal_error("Cannot tally " + score_str + "with a delayedgroup filter");
}
auto score = score_str_to_int(score_str);
switch (score) {
case SCORE_FLUX:
if (!nuclides_.empty())
if (!(nuclides_.size() == 1 && nuclides_[0] == -1))
fatal_error("Cannot tally flux for an individual nuclide.");
if (energyout_present)
fatal_error("Cannot tally flux with an outgoing energy filter.");
break;
case SCORE_TOTAL:
case SCORE_ABSORPTION:
case SCORE_FISSION:
if (energyout_present)
fatal_error("Cannot tally " + score_str + " reaction rate with an "
"outgoing energy filter");
break;
case SCORE_SCATTER:
if (legendre_present)
estimator_ = ESTIMATOR_ANALOG;
case SCORE_NU_FISSION:
case SCORE_DELAYED_NU_FISSION:
case SCORE_PROMPT_NU_FISSION:
if (energyout_present)
estimator_ = ESTIMATOR_ANALOG;
break;
case SCORE_NU_SCATTER:
if (settings::run_CE) {
estimator_ = ESTIMATOR_ANALOG;
} else {
if (energyout_present || legendre_present)
estimator_ = ESTIMATOR_ANALOG;
}
break;
case N_2N:
case N_3N:
case N_4N:
case N_GAMMA:
case N_P:
case N_A:
depletion_rx_ = true;
break;
case SCORE_CURRENT:
// Check which type of current is desired: mesh or surface currents.
if (surface_present || cell_present || cellfrom_present) {
if (meshsurface_present)
fatal_error("Cannot tally mesh surface currents in the same tally as "
"normal surface currents");
type_ = TALLY_SURFACE;
} else if (meshsurface_present) {
type_ = TALLY_MESH_SURFACE;
} else {
fatal_error("Cannot tally currents without surface type filters");
}
break;
}
scores_.push_back(score);
}
// Make sure that no duplicate scores exist.
for (auto it1 = scores_.begin(); it1 != scores_.end(); ++it1) {
for (auto it2 = it1 + 1; it2 != scores_.end(); ++it2) {
if (*it1 == *it2)
fatal_error("Duplicate score of type \"" + reaction_name(*it1)
+ "\" found in tally " + std::to_string(id_));
}
}
// Make sure all scores are compatible with multigroup mode.
if (!settings::run_CE) {
for (auto sc : scores_)
if (sc > 0)
fatal_error("Cannot tally " + reaction_name(sc) + " reaction rate "
"in multi-group mode");
}
// Make sure current scores are not mixed in with volumetric scores.
if (type_ == TALLY_SURFACE || type_ == TALLY_MESH_SURFACE) {
if (scores_.size() != 1)
fatal_error("Cannot tally other scores in the same tally as surface "
"currents");
}
}
void
Tally::set_nuclides(pugi::xml_node node)
{
nuclides_.clear();
// By default, we tally just the total material rates.
if (!check_for_node(node, "nuclides")) {
nuclides_.push_back(-1);
return;
}
if (get_node_value(node, "nuclides") == "all") {
// This tally should bin every nuclide in the problem. It should also bin
// the total material rates. To achieve this, set the nuclides_ vector to
// 0, 1, 2, ..., -1.
nuclides_.reserve(data::nuclides.size() + 1);
for (auto i = 0; i < data::nuclides.size(); ++i)
nuclides_.push_back(i);
nuclides_.push_back(-1);
all_nuclides_ = true;
} else {
// The user provided specifics nuclides. Parse it as an array with either
// "total" or a nuclide name like "U-235" in each position.
auto words = get_node_array<std::string>(node, "nuclides");
for (auto word : words) {
if (word == "total") {
nuclides_.push_back(-1);
} else {
auto search = data::nuclide_map.find(word);
if (search == data::nuclide_map.end())
fatal_error("Could not find the nuclide " + word
+ " specified in tally " + std::to_string(id_)
+ " in any material");
nuclides_.push_back(search->second);
}
}
}
}
void
Tally::init_triggers(pugi::xml_node node, int i_tally)
{
for (auto trigger_node: node.children("trigger")) {
// Read the trigger type.
TriggerMetric metric;
if (check_for_node(trigger_node, "type")) {
auto type_str = get_node_value(trigger_node, "type");
if (type_str == "std_dev") {
metric = TriggerMetric::standard_deviation;
} else if (type_str == "variance") {
metric = TriggerMetric::variance;
} else if (type_str == "rel_err") {
metric = TriggerMetric::relative_error;
} else {
std::stringstream msg;
msg << "Unknown trigger type \"" << type_str << "\" in tally " << id_;
fatal_error(msg);
}
} else {
std::stringstream msg;
msg << "Must specify trigger type for tally " << id_
<< " in tally XML file";
fatal_error(msg);
}
// Read the trigger threshold.
double threshold;
if (check_for_node(trigger_node, "threshold")) {
threshold = std::stod(get_node_value(trigger_node, "threshold"));
} else {
std::stringstream msg;
msg << "Must specify trigger threshold for tally " << id_
<< " in tally XML file";
fatal_error(msg);
}
// Read the trigger scores.
std::vector<std::string> trigger_scores;
if (check_for_node(trigger_node, "scores")) {
trigger_scores = get_node_array<std::string>(trigger_node, "scores");
} else {
trigger_scores.push_back("all");
}
// 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_.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;
}
if (i_score == tally.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";
fatal_error(msg);
}
triggers_.push_back({metric, threshold, i_score});
}
}
}
}
//==============================================================================
// Non-member functions
//==============================================================================
@ -44,7 +525,8 @@ adaptor_type<3> tally_results(int idx)
// Get pointer to tally results
double* results;
std::array<std::size_t, 3> shape;
openmc_tally_results(idx, &results, shape.data());
// TODO: off-by-one
openmc_tally_results(idx+1, &results, shape.data());
// Adapt array into xtensor with no ownership
std::size_t size {shape[0] * shape[1] * shape[2]};
@ -54,10 +536,11 @@ adaptor_type<3> tally_results(int idx)
#ifdef OPENMC_MPI
void reduce_tally_results()
{
for (int i = 1; i <= n_tallies; ++i) {
for (int i = 0; i < n_tallies; ++i) {
// Skip any tallies that are not active
bool active;
openmc_tally_get_active(i, &active);
// TODO: off-by-one
openmc_tally_get_active(i+1, &active);
if (!active) continue;
// Get view of accumulated tally values
@ -108,4 +591,335 @@ void reduce_tally_results()
}
#endif
extern "C" void
setup_active_tallies_c()
{
model::active_tallies.clear();
model::active_analog_tallies.clear();
model::active_tracklength_tallies.clear();
model::active_collision_tallies.clear();
model::active_meshsurf_tallies.clear();
model::active_surface_tallies.clear();
for (auto i = 0; i < model::tallies.size(); ++i) {
const auto& tally {*model::tallies[i]};
if (tally.active_) {
model::active_tallies.push_back(i);
switch (tally.type_) {
case TALLY_VOLUME:
switch (tally.estimator_) {
case ESTIMATOR_ANALOG:
model::active_analog_tallies.push_back(i);
break;
case ESTIMATOR_TRACKLENGTH:
model::active_tracklength_tallies.push_back(i);
break;
case ESTIMATOR_COLLISION:
model::active_collision_tallies.push_back(i);
}
break;
case TALLY_MESH_SURFACE:
model::active_meshsurf_tallies.push_back(i);
break;
case TALLY_SURFACE:
model::active_surface_tallies.push_back(i);
}
}
}
}
extern "C" void
free_memory_tally_c()
{
#pragma omp parallel
{
model::tally_derivs.clear();
}
model::tally_filters.clear();
model::tallies.clear();
model::active_tallies.clear();
model::active_analog_tallies.clear();
model::active_tracklength_tallies.clear();
model::active_collision_tallies.clear();
model::active_meshsurf_tallies.clear();
model::active_surface_tallies.clear();
}
//==============================================================================
// C-API functions
//==============================================================================
extern "C" int
openmc_tally_get_type(int32_t index, int32_t* type)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
//TODO: off-by-one
*type = model::tallies[index-1]->type_;
return 0;
}
extern "C" int
openmc_tally_set_type(int32_t index, const char* type)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
if (strcmp(type, "volume") == 0) {
model::tallies[index-1]->type_ = TALLY_VOLUME;
} else if (strcmp(type, "mesh-surface") == 0) {
model::tallies[index-1]->type_ = TALLY_MESH_SURFACE;
} else if (strcmp(type, "surface") == 0) {
model::tallies[index-1]->type_ = TALLY_SURFACE;
} else {
std::stringstream errmsg;
errmsg << "Unknown tally type: " << type;
set_errmsg(errmsg);
return OPENMC_E_INVALID_ARGUMENT;
}
return 0;
}
extern "C" int
openmc_tally_get_active(int32_t index, bool* active)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
//TODO: off-by-one
*active = model::tallies[index-1]->active_;
return 0;
}
extern "C" int
openmc_tally_set_active(int32_t index, bool active)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
//TODO: off-by-one
model::tallies[index-1]->active_ = active;
return 0;
}
extern "C" int
openmc_tally_get_scores(int32_t index, int** scores, int* n)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
//TODO: off-by-one
*scores = model::tallies[index-1]->scores_.data();
*n = model::tallies[index-1]->scores_.size();
return 0;
}
extern "C" int
openmc_tally_set_scores(int32_t index, int n, const char** scores)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
std::vector<std::string> scores_str(scores, scores+n);
try {
//TODO: off-by-one
model::tallies[index-1]->set_scores(scores_str);
} catch (const std::invalid_argument& ex) {
set_errmsg(ex.what());
return OPENMC_E_INVALID_ARGUMENT;
}
return 0;
}
extern "C" int
openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n)
{
// Make sure the index fits in the array bounds.
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
//TODO: off-by-one
*n = model::tallies[index-1]->nuclides_.size();
*nuclides = model::tallies[index-1]->nuclides_.data();
return 0;
}
extern "C" int
openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides)
{
// Make sure the index fits in the array bounds.
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
std::vector<std::string> words(nuclides, nuclides+n);
std::vector<int> nucs;
for (auto word : words){
if (word == "total") {
nucs.push_back(-1);
} else {
auto search = data::nuclide_map.find(word);
if (search == data::nuclide_map.end()) {
set_errmsg("Nuclide \"" + word + "\" has not been loaded yet");
return OPENMC_E_DATA;
}
nucs.push_back(search->second);
}
}
//TODO: off-by-one
model::tallies[index-1]->nuclides_ = nucs;
return 0;
}
extern "C" int
openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n)
{
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
//TODO: off-by-one
*indices = model::tallies[index-1]->filters().data();
*n = model::tallies[index-1]->filters().size();
return 0;
}
extern "C" int
openmc_tally_set_filters(int32_t index, int n, const int32_t* indices)
{
// Make sure the index fits in the array bounds.
if (index < 1 || index > model::tallies.size()) {
set_errmsg("Index in tallies array is out of bounds.");
return OPENMC_E_OUT_OF_BOUNDS;
}
// Set the filters.
try {
//TODO: off-by-one
model::tallies[index-1]->set_filters(indices, n);
} catch (const std::out_of_range& ex) {
set_errmsg(ex.what());
return OPENMC_E_OUT_OF_BOUNDS;
}
return 0;
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
Tally* tally_pointer(int indx) {return model::tallies[indx].get();}
void
extend_tallies_c(int n)
{
for (int i = 0; i < n; ++i)
model::tallies.push_back(std::make_unique<Tally>());
}
int active_tallies_data(int i)
{return model::active_tallies[i-1];}
int active_tallies_size()
{return model::active_tallies.size();}
int active_analog_tallies_size()
{return model::active_analog_tallies.size();}
int active_tracklength_tallies_size()
{return model::active_tracklength_tallies.size();}
int active_collision_tallies_size()
{return model::active_collision_tallies.size();}
int active_meshsurf_tallies_size()
{return model::active_meshsurf_tallies.size();}
int active_surface_tallies_size()
{return model::active_surface_tallies.size();}
void tally_init_from_xml(Tally* tally, pugi::xml_node* node)
{tally->init_from_xml(*node);}
int tally_get_id_c(Tally* tally) {return tally->id_;}
void tally_set_id_c(Tally* tally, int id) {tally->id_ = id;}
int tally_get_type_c(Tally* tally) {return tally->type_;}
void tally_set_type_c(Tally* tally, int type) {tally->type_ = type;}
int tally_get_estimator_c(Tally* tally) {return tally->estimator_;}
void tally_set_estimator_c(Tally* tally, int e) {tally->estimator_ = e;}
bool tally_get_depletion_rx_c(Tally* tally) {return tally->depletion_rx_;}
int tally_get_n_scores_c(Tally* tally) {return tally->scores_.size();}
int tally_get_score_c(Tally* tally, int i) {return tally->scores_[i];}
void tally_set_filters_c(Tally* tally, int n, int32_t filter_indices[])
{tally->set_filters(filter_indices, n);}
int tally_get_n_filters_c(Tally* tally) {return tally->filters().size();}
int32_t tally_get_filter_c(Tally* tally, int i) {return tally->filters(i);}
int32_t tally_get_n_filter_bins_c(Tally* tally)
{return tally->n_filter_bins();}
int tally_get_n_nuclide_bins_c(Tally* tally)
{return tally->nuclides_.size();}
int tally_get_nuclide_bins_c(Tally* tally, int i)
{return tally->nuclides_[i-1];}
int tally_get_energyout_filter_c(Tally* tally)
{return tally->energyout_filter_;}
void tally_set_scores(Tally* tally, pugi::xml_node* node)
{tally->set_scores(*node);}
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);}
int tally_get_deriv_c(Tally* tally) {return tally->deriv_;}
int tally_set_deriv_c(Tally* tally, int deriv) {tally->deriv_ = deriv;}
}
} // namespace openmc

106
src/tallies/tally_derivative_header.F90
View file

@ -1,103 +1,33 @@
module tally_derivative_header
use constants
use dict_header, only: DictIntInt
use error, only: fatal_error
use nuclide_header, only: nuclide_map_get
use string, only: to_str, to_lower, to_c_string
use xml_interface
use, intrinsic :: ISO_C_BINDING
implicit none
private
public :: free_memory_tally_derivative
!===============================================================================
! TALLYDERIVATIVE describes a first-order derivative that can be applied to
! tallies.
!===============================================================================
type, public :: TallyDerivative
integer :: id
integer :: variable
integer :: diff_material
integer :: diff_nuclide
real(8) :: flux_deriv
contains
procedure :: from_xml
type, bind(C), public :: TallyDerivative
integer(C_INT) :: id
integer(C_INT) :: variable
integer(C_INT) :: diff_material
integer(C_INT) :: diff_nuclide
real(C_DOUBLE) :: flux_deriv
end type TallyDerivative
type(TallyDerivative), public, allocatable :: tally_derivs(:)
!$omp threadprivate(tally_derivs)
interface
function n_tally_derivs() result(n) bind(C)
import C_INT
integer(C_INT) :: n
end function
! Dictionary that maps user IDs to indices in 'tally_derivs'
type(DictIntInt), public :: tally_deriv_dict
contains
subroutine from_xml(this, node)
class(TallyDerivative), intent(inout) :: this
type(XMLNode), intent(in) :: node
character(MAX_WORD_LEN) :: temp_str
character(MAX_WORD_LEN) :: word
integer :: val
! Copy the derivative id.
if (check_for_node(node, "id")) then
call get_node_value(node, "id", this % id)
else
call fatal_error("Must specify an ID for <derivative> elements in the&
& tally XML file")
end if
! Make sure the id is > 0.
if (this % id <= 0) then
call fatal_error("<derivative> IDs must be an integer greater than &
&zero")
end if
! Make sure this id has not already been used.
if (tally_deriv_dict % has(this % id)) then
call fatal_error("Two or more <derivative>'s use the same unique &
&ID: " // trim(to_str(this % id)))
end if
! Read the independent variable name.
call get_node_value(node, "variable", temp_str)
temp_str = to_lower(temp_str)
select case(temp_str)
case("density")
this % variable = DIFF_DENSITY
case("nuclide_density")
this % variable = DIFF_NUCLIDE_DENSITY
call get_node_value(node, "nuclide", word)
val = nuclide_map_get(to_c_string(word))
if (val == -1) then
call fatal_error("Could not find the nuclide " &
// trim(word) // " specified in derivative " &
// trim(to_str(this % id)) // " in any material.")
end if
this % diff_nuclide = val
case("temperature")
this % variable = DIFF_TEMPERATURE
end select
call get_node_value(node, "material", this % diff_material)
end subroutine from_xml
!===============================================================================
! FREE_MEMORY_TALLY_DERIVATIVE deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_tally_derivative()
!$omp parallel
if (allocated(tally_derivs)) deallocate(tally_derivs)
!$omp end parallel
end subroutine free_memory_tally_derivative
function tally_deriv_c(i) result(deriv) bind(C)
import C_INT, TallyDerivative
integer(C_INT), value :: i
type(TallyDerivative), pointer :: deriv
end function
end interface
end module tally_derivative_header

4
src/tallies/tally_filter.F90
View file

@ -7,11 +7,7 @@ module tally_filter
use tally_filter_header
! Inherit other filters
use tally_filter_delayedgroup
use tally_filter_distribcell
use tally_filter_energy
use tally_filter_mesh
use tally_filter_meshsurface
use tally_filter_particle
use tally_filter_sph_harm

38
src/tallies/tally_filter_delayedgroup.F90
View file

@ -1,38 +0,0 @@
module tally_filter_delayedgroup
use, intrinsic :: ISO_C_BINDING
use tally_filter_header
implicit none
private
!===============================================================================
! DELAYEDGROUPFILTER bins outgoing fission neutrons in their delayed groups.
! The get_all_bins functionality is not actually used. The bins are manually
! iterated over in the scoring subroutines.
!===============================================================================
type, public, extends(TallyFilter) :: DelayedGroupFilter
contains
procedure :: groups
end type DelayedGroupFilter
contains
function groups(this, i) result(group)
class(DelayedGroupFilter), intent(in) :: this
integer, intent(in) :: i
integer :: group
interface
function delayedgroup_filter_groups(filt, i) result(group) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: filt
integer(C_INT), value :: i
integer(C_INT) :: group
end function
end interface
group = delayedgroup_filter_groups(this % ptr, i)
end function groups
end module tally_filter_delayedgroup

91
src/tallies/tally_filter_energy.F90
View file

@ -1,91 +0,0 @@
module tally_filter_energy
use, intrinsic :: ISO_C_BINDING
use tally_filter_header
use xml_interface
implicit none
private
public :: openmc_energy_filter_get_bins
public :: openmc_energy_filter_set_bins
interface
function openmc_energy_filter_get_bins(index, energies, n) result(err) bind(C)
import C_INT32_T, C_PTR, C_INT
integer(C_INT32_T), value :: index
type(C_PTR), intent(out) :: energies
integer(C_INT32_T), intent(out) :: n
integer(C_INT) :: err
end function
function openmc_energy_filter_set_bins(index, n, energies) result(err) bind(C)
import C_INT32_T, C_DOUBLE, C_INT
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT32_T), value, intent(in) :: n
real(C_DOUBLE), intent(in) :: energies(n)
integer(C_INT) :: err
end function
end interface
!===============================================================================
! ENERGYFILTER bins the incident neutron energy.
!===============================================================================
type, public, extends(TallyFilter) :: EnergyFilter
! 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 :: search
end type EnergyFilter
!===============================================================================
! ENERGYOUTFILTER bins the outgoing neutron energy. Only scattering events use
! the get_all_bins functionality. Nu-fission tallies manually iterate over the
! filter bins.
!===============================================================================
type, public, extends(EnergyFilter) :: EnergyoutFilter
end type EnergyoutFilter
contains
!===============================================================================
! EnergyFilter methods
!===============================================================================
subroutine from_xml_energy(this, node)
class(EnergyFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
interface
function energy_filter_matches_transport_groups(filt) result(matches) &
bind(C)
import C_PTR, C_BOOL
type(C_PTR), value :: filt
logical(C_BOOL) :: matches
end function
end interface
call this % from_xml_cpp(node)
this % n_bins = this % n_bins_cpp()
this % matches_transport_groups = &
energy_filter_matches_transport_groups(this % ptr)
end subroutine from_xml_energy
function search(this, val) result(bin)
class(EnergyFilter), intent(in) :: this
real(8), intent(in) :: val
integer :: bin
interface
function energy_filter_search(filt, val) result(bin) bind(C)
import C_PTR, C_DOUBLE, C_INT
type(C_PTR), value :: filt
real(C_DOUBLE), value :: val
integer(C_INT) :: bin
end function
end interface
bin = energy_filter_search(this % ptr, val)
end function search
end module tally_filter_energy

228
src/tallies/tally_filter_header.F90
View file

@ -19,39 +19,6 @@ module tally_filter_header
public :: openmc_filter_set_id
public :: openmc_get_filter_index
public :: openmc_get_filter_next_id
public :: filter_match_pointer
interface
function filter_match_pointer(indx) bind(C) result(ptr)
import C_PTR, C_INT
integer(C_INT), intent(in), value :: indx
type(C_PTR) :: ptr
end function filter_match_pointer
end interface
!===============================================================================
! TALLYFILTERMATCH stores every valid bin and weight for a filter
!===============================================================================
type, public :: TallyFilterMatch
type(C_PTR) :: ptr
! Index of the bin and weight being used in the current filter combination
integer :: i_bin
! Indicates whether all valid bins for this filter have been found
logical :: bins_present = .false.
contains
procedure :: bins_push_back
procedure :: weights_push_back
procedure :: bins_clear
procedure :: weights_clear
procedure :: bins_size
procedure :: bins_data
procedure :: weights_data
procedure :: bins_set_data
end type TallyFilterMatch
!===============================================================================
! TALLYFILTER describes a filter that limits what events score to a tally. For
@ -60,14 +27,12 @@ module tally_filter_header
!===============================================================================
type, public, abstract :: TallyFilter
integer :: id
integer :: n_bins = 0
type(C_PTR) :: ptr
contains
procedure :: id => filter_get_id_f90
procedure :: from_xml
procedure :: get_all_bins
procedure :: to_statepoint
procedure :: text_label
procedure :: initialize
procedure :: n_bins_cpp
procedure :: from_xml_cpp
@ -90,6 +55,15 @@ module tally_filter_header
type, public, extends(CellFilter) :: CellFromFilter
end type
type, public, extends(TallyFilter) :: DelayedGroupFilter
end type
type, public, extends(TallyFilter) :: EnergyFilter
end type
type, public, extends(EnergyFilter) :: EnergyoutFilter
end type
type, public, extends(TallyFilter) :: EnergyFunctionFilter
end type
@ -99,6 +73,12 @@ module tally_filter_header
type, public, extends(TallyFilter) :: MaterialFilter
end type
type, public, extends(TallyFilter) :: MeshFilter
end type
type, public, extends(TallyFilter) :: MeshSurfaceFilter
end type
type, public, extends(TallyFilter) :: MuFilter
end type
@ -134,8 +114,6 @@ module tally_filter_header
integer(C_INT32_T), public, bind(C) :: n_filters = 0 ! # of filters
type(TallyFilterContainer), public, allocatable, target :: filters(:)
type(TallyFilterMatch), public, allocatable :: filter_matches(:)
!$omp threadprivate(filter_matches)
! Dictionary that maps user IDs to indices in 'filters'
type(DictIntInt), public :: filter_dict
@ -146,120 +124,23 @@ module tally_filter_header
contains
!===============================================================================
! TallyFilterMatch implementation
!===============================================================================
subroutine bins_push_back(this, val)
class(TallyFilterMatch), intent(inout) :: this
integer, intent(in) :: val
interface
subroutine filter_match_bins_push_back(ptr, val) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), intent(in), value :: val
end subroutine
end interface
call filter_match_bins_push_back(this % ptr, val)
end subroutine bins_push_back
subroutine weights_push_back(this, val)
class(TallyFilterMatch), intent(inout) :: this
real(8), intent(in) :: val
interface
subroutine filter_match_weights_push_back(ptr, val) bind(C)
import C_PTR, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE), intent(in), value :: val
end subroutine
end interface
call filter_match_weights_push_back(this % ptr, val)
end subroutine weights_push_back
subroutine bins_clear(this)
class(TallyFilterMatch), intent(inout) :: this
interface
subroutine filter_match_bins_clear(ptr) bind(C)
import C_PTR
type(C_PTR), value :: ptr
end subroutine
end interface
call filter_match_bins_clear(this % ptr)
end subroutine bins_clear
subroutine weights_clear(this)
class(TallyFilterMatch), intent(inout) :: this
interface
subroutine filter_match_weights_clear(ptr) bind(C)
import C_PTR
type(C_PTR), value :: ptr
end subroutine
end interface
call filter_match_weights_clear(this % ptr)
end subroutine weights_clear
function bins_size(this) result(len)
class(TallyFilterMatch), intent(inout) :: this
integer :: len
interface
function filter_match_bins_size(ptr) bind(C) result(len)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT) :: len
end function
end interface
len = filter_match_bins_size(this % ptr)
end function bins_size
function bins_data(this, indx) result(val)
class(TallyFilterMatch), intent(inout) :: this
integer, intent(in) :: indx
integer :: val
interface
function filter_match_bins_data(ptr, indx) bind(C) result(val)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), intent(in), value :: indx
integer(C_INT) :: val
end function
end interface
val = filter_match_bins_data(this % ptr, indx)
end function bins_data
function weights_data(this, indx) result(val)
class(TallyFilterMatch), intent(inout) :: this
integer, intent(in) :: indx
real(8) :: val
interface
function filter_match_weights_data(ptr, indx) bind(C) result(val)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), intent(in), value :: indx
real(C_DOUBLE) :: val
end function
end interface
val = filter_match_weights_data(this % ptr, indx)
end function weights_data
subroutine bins_set_data(this, indx, val)
class(TallyFilterMatch), intent(inout) :: this
integer, intent(in) :: indx
integer, intent(in) :: val
interface
subroutine filter_match_bins_set_data(ptr, indx, val) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), value, intent(in) :: indx
integer(C_INT), value, intent(in) :: val
end subroutine
end interface
call filter_match_bins_set_data(this % ptr, indx, val)
end subroutine bins_set_data
!===============================================================================
! TallyFilter implementation
!===============================================================================
function filter_get_id_f90(this) result(id)
class(TallyFilter) :: this
integer(C_INT32_T) :: id
interface
function filter_get_id(filt) result(id) bind(C)
import C_PTR, C_INT32_T
type(C_PTR), value :: filt
integer(C_INT32_T) :: id
end function
end interface
id = filter_get_id(this % ptr)
end function
subroutine from_xml(this, node)
class(TallyFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
@ -267,23 +148,6 @@ contains
this % n_bins = this % n_bins_cpp()
end subroutine from_xml
subroutine get_all_bins(this, p, estimator, match)
class(TallyFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
interface
subroutine filter_get_all_bins(filt, p, estimator, match) bind(C)
import C_PTR, Particle, C_INT
type(C_PTR), value :: filt
type(Particle), intent(in) :: p
integer(C_INT), intent(in), value :: estimator
type(C_PTR), value :: match
end subroutine filter_get_all_bins
end interface
call filter_get_all_bins(this % ptr, p, estimator, match % ptr)
end subroutine get_all_bins
subroutine to_statepoint(this, filter_group)
class(TallyFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
@ -297,28 +161,6 @@ contains
call filter_to_statepoint(this % ptr, filter_group)
end subroutine to_statepoint
function text_label(this, bin) result(label)
class(TallyFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
character(kind=C_CHAR) :: label_(MAX_LINE_LEN+1)
integer :: i
interface
subroutine filter_text_label(filt, bin, label) bind(C)
import C_PTR, C_INT, C_CHAR
type(C_PTR), value :: filt
integer(C_INT), value :: bin
character(kind=C_CHAR) :: label(*)
end subroutine filter_text_label
end interface
call filter_text_label(this % ptr, bin, label_)
label = " "
do i = 1, MAX_LINE_LEN
if (label_(i) == C_NULL_CHAR) exit
label(i:i) = label_(i)
end do
end function text_label
subroutine initialize(this)
class(TallyFilter), intent(inout) :: this
call this % initialize_cpp()
@ -456,7 +298,7 @@ contains
integer(C_INT) :: err
if (index >= 1 .and. index <= n_filters) then
id = filters(index) % obj % id
id = filters(index) % obj % id()
err = 0
else
err = E_OUT_OF_BOUNDS
@ -471,9 +313,17 @@ contains
integer(C_INT32_T), value, intent(in) :: id
integer(C_INT) :: err
interface
subroutine filter_set_id(filt, id) bind(C)
import C_PTR, C_INT32_T
type(C_PTR), value :: filt
integer(C_INT32_T), value :: id
end subroutine
end interface
if (index >= 1 .and. index <= n_filters) then
if (allocated(filters(index) % obj)) then
filters(index) % obj % id = id
call filter_set_id(filters(index) % obj % ptr, id)
call filter_dict % set(id, index)
if (id > largest_filter_id) largest_filter_id = id

38
src/tallies/tally_filter_mesh.F90
View file

@ -1,38 +0,0 @@
module tally_filter_mesh
use, intrinsic :: ISO_C_BINDING
use tally_filter_header
implicit none
interface
function openmc_mesh_filter_set_mesh(index, index_mesh) result(err) bind(C)
import C_INT32_T, C_INT
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT32_T), value, intent(in) :: index_mesh
integer(C_INT) :: err
end function
end interface
type, public, extends(TallyFilter) :: MeshFilter
contains
procedure :: mesh => get_mesh
end type MeshFilter
contains
function get_mesh(this) result(mesh)
class(MeshFilter) :: this
integer :: mesh
interface
function mesh_filter_get_mesh(filt) result(index_mesh) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: filt
integer(C_INT) :: index_mesh
end function mesh_filter_get_mesh
end interface
mesh = mesh_filter_get_mesh(this % ptr)
end function get_mesh
end module tally_filter_mesh

22
src/tallies/tally_filter_meshsurface.F90
View file

@ -1,22 +0,0 @@
module tally_filter_meshsurface
use, intrinsic :: ISO_C_BINDING
use tally_filter_header
implicit none
interface
function openmc_meshsurface_filter_set_mesh(index, index_mesh) result(err) &
bind(C)
import C_INT32_T, C_INT
integer(C_INT32_T), value, intent(in) :: index
integer(C_INT32_T), value, intent(in) :: index_mesh
integer(C_INT) :: err
end function
end interface
type, extends(TallyFilter) :: MeshSurfaceFilter
end type MeshSurfaceFilter
end module tally_filter_meshsurface

916
src/tallies/tally_header.F90
View file

File diff suppressed because it is too large Load diff

2308
src/tallies/tally_scoring.cpp Normal file
View file

File diff suppressed because it is too large Load diff

443
src/tallies/trigger.F90
View file

@ -1,443 +0,0 @@
module trigger
use, intrinsic :: ISO_C_BINDING
#ifdef OPENMC_MPI
use message_passing
#endif
use constants
use eigenvalue, only: openmc_get_keff
use error, only: warning, write_message
use string, only: to_str
use mesh_header, only: RegularMesh, meshes
use message_passing, only: master
use settings
use simulation_header
use trigger_header
use tally, only: TallyObject
use tally_filter_mesh, only: MeshFilter
use tally_filter_header, only: filter_matches, filters
use tally_header, only: tallies, n_tallies
implicit none
contains
!===============================================================================
! CHECK_TRIGGERS checks any user-specified precision triggers' for convergence
! and predicts the number of remainining batches to convergence.
!===============================================================================
subroutine check_triggers() bind(C)
implicit none
! Variables to reflect distance to trigger convergence criteria
real(8) :: max_ratio ! max uncertainty/thresh ratio
integer :: tally_id ! id for tally with max ratio
character(len=52) :: name ! "eigenvalue" or tally score
integer :: n_pred_batches ! predicted # batches to satisfy all triggers
! Checks if current_batch is one for which the triggers must be checked
if (current_batch < n_batches .or. (.not. trigger_on)) return
if (mod((current_batch - n_batches), n_batch_interval) /= 0 .and. &
current_batch /= n_max_batches) return
! Check the trigger and output the result
call check_tally_triggers(max_ratio, tally_id, name)
! When trigger threshold is reached, write information
if (satisfy_triggers) then
call write_message("Triggers satisfied for batch " // &
trim(to_str(current_batch)), 7)
! When trigger is not reached write convergence info for user
elseif (name == "eigenvalue") then
call write_message("Triggers unsatisfied, max unc./thresh. is " // &
trim(to_str(max_ratio)) // " for " // trim(name), 7)
else
call write_message("Triggers unsatisfied, max unc./thresh. is " // &
trim(to_str(max_ratio)) // " for " // trim(name) // &
" in tally " // trim(to_str(tally_id)), 7)
end if
! If batch_interval is not set, estimate batches till triggers are satisfied
if (pred_batches .and. .not. satisfy_triggers) then
! Estimate the number of remaining batches to convergence
! The prediction uses the fact that tally variances are proportional
! to 1/N where N is the number of the batches/particles
n_batch_interval = int((current_batch-n_inactive) * &
(max_ratio ** 2)) + n_inactive-n_batches + 1
n_pred_batches = n_batch_interval + n_batches
! Write the predicted number of batches for the user
if (n_pred_batches > n_max_batches) then
call warning("The estimated number of batches is " // &
trim(to_str(n_pred_batches)) // &
" -- greater than max batches. ")
else
call write_message("The estimated number of batches is " // &
trim(to_str(n_pred_batches)), 7)
end if
end if
end subroutine check_triggers
!===============================================================================
! CHECK_TALLY_TRIGGERS checks whether uncertainties are below the threshold,
! and finds the maximum uncertainty/threshold ratio for all triggers
!===============================================================================
subroutine check_tally_triggers(max_ratio, tally_id, name)
! Variables to reflect distance to trigger convergence criteria
real(8), intent(inout) :: max_ratio ! max uncertainty/thresh ratio
integer, intent(inout) :: tally_id ! id for tally with max ratio
character(len=52), intent(inout) :: name ! "eigenvalue" or tally score
integer :: i ! index in tallies array
integer :: j ! index in tally filters
integer :: n ! loop index for nuclides
integer :: s ! loop index for triggers
integer :: filter_index ! index in results array for filters
integer :: score_index ! scoring bin index
integer(C_INT) :: err
real(8) :: uncertainty ! trigger uncertainty
real(8) :: std_dev = ZERO ! trigger standard deviation
real(8) :: rel_err = ZERO ! trigger relative error
real(8) :: ratio ! ratio of the uncertainty/trigger threshold
real(C_DOUBLE) :: k_combined(2)
! Initialize tally trigger maximum uncertainty ratio to zero
max_ratio = 0
if (master) then
! By default, assume all triggers are satisfied
satisfy_triggers = .true.
! Check eigenvalue trigger
if (run_mode == MODE_EIGENVALUE) then
if (keff_trigger % trigger_type /= 0) then
err = openmc_get_keff(k_combined)
select case (keff_trigger % trigger_type)
case(VARIANCE)
uncertainty = k_combined(2) ** 2
case(STANDARD_DEVIATION)
uncertainty = k_combined(2)
case default
uncertainty = k_combined(2) / k_combined(1)
end select
! If uncertainty is above threshold, store uncertainty ratio
if (uncertainty > keff_trigger % threshold) then
satisfy_triggers = .false.
if (keff_trigger % trigger_type == VARIANCE) then
ratio = sqrt(uncertainty / keff_trigger % threshold)
else
ratio = uncertainty / keff_trigger % threshold
end if
if (max_ratio < ratio) then
max_ratio = ratio
name = "eigenvalue"
end if
end if
end if
end if
! Compute uncertainties for all tallies, scores with triggers
TALLY_LOOP: do i = 1, n_tallies
associate (t => tallies(i) % obj)
! Cycle through if only one batch has been simumlate
if (t % n_realizations == 1) then
cycle TALLY_LOOP
end if
TRIGGER_LOOP: do s = 1, t % n_triggers
associate (trigger => t % triggers(s))
! Initialize trigger uncertainties to zero
trigger % std_dev = ZERO
trigger % rel_err = ZERO
trigger % variance = ZERO
! Mesh current tally triggers require special treatment
if (t % type == TALLY_MESH_SURFACE) then
call compute_tally_current(t, trigger)
else
! Initialize bins, filter level
do j = 1, size(t % filter)
call filter_matches(t % filter(j)) % bins_clear()
call filter_matches(t % filter(j)) % bins_push_back(0)
end do
FILTER_LOOP: do filter_index = 1, t % n_filter_bins
! Initialize score index
score_index = trigger % score_index
! Initialize score bin index
NUCLIDE_LOOP: do n = 1, t % n_nuclide_bins
call get_trigger_uncertainty(std_dev, rel_err, &
score_index, filter_index, t)
if (trigger % variance < variance) then
trigger % variance = std_dev ** 2
end if
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
end if
if (trigger % rel_err < rel_err) then
trigger % rel_err = rel_err
end if
select case (t % triggers(s) % type)
case(VARIANCE)
uncertainty = trigger % variance
case(STANDARD_DEVIATION)
uncertainty = trigger % std_dev
case default
uncertainty = trigger % rel_err
end select
if (uncertainty > t % triggers(s) % threshold) then
satisfy_triggers = .false.
if (t % triggers(s) % type == VARIANCE) then
ratio = sqrt(uncertainty / t % triggers(s) % threshold)
else
ratio = uncertainty / t % triggers(s) % threshold
end if
if (max_ratio < ratio) then
max_ratio = ratio
name = t % triggers(s) % score_name
tally_id = t % id
end if
end if
end do NUCLIDE_LOOP
if (size(t % filter) == 0) exit FILTER_LOOP
end do FILTER_LOOP
end if
end associate
end do TRIGGER_LOOP
end associate
end do TALLY_LOOP
end if
end subroutine check_tally_triggers
!===============================================================================
! COMPUTE_TALLY_CURRENT computes the current for a mesh current tally with
! precision trigger(s).
!===============================================================================
subroutine compute_tally_current(t, trigger)
type(TallyObject), intent(in) :: t ! mesh current tally
type(TriggerObject), intent(inout) :: trigger ! mesh current tally trigger
integer :: i ! mesh index
integer :: j ! loop index for tally filters
integer :: ijk(3) ! indices of mesh cells
integer :: n_dim ! number of mesh dimensions
integer :: n_cells ! number of mesh cells
integer :: l ! index for energy
integer :: i_filter_mesh ! index for mesh filter
integer :: i_filter_ein ! index for incoming energy filter
integer :: i_filter_surf ! index for surface filter
integer :: n ! number of incoming energy bins
integer :: filter_index ! index in results array for filters
logical :: print_ebin ! should incoming energy bin be displayed?
real(8) :: rel_err = ZERO ! temporary relative error of result
real(8) :: std_dev = ZERO ! temporary standard deviration of result
type(RegularMesh) :: m ! surface current mesh
! Get pointer to mesh
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
select type(filt => filters(i_filter_mesh) % obj)
type is (MeshFilter)
m = meshes(filt % mesh())
end select
! initialize bins array
do j = 1, size(t % filter)
call filter_matches(t % filter(j)) % bins_clear()
call filter_matches(t % filter(j)) % bins_push_back(1)
end do
! determine how many energyin bins there are
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
if (i_filter_ein > 0) then
print_ebin = .true.
n = filters(t % filter(i_filter_ein)) % obj % n_bins
i_filter_ein = t % filter(i_filter_ein)
else
print_ebin = .false.
n = 1
end if
! Get the dimensions and number of cells in the mesh
n_dim = m % n_dimension()
n_cells = 1
do j = 1, n_dim
n_cells = n_cells * m % dimension(j)
end do
! Loop over all the mesh cells
do i = 1, n_cells
! Get the indices for this cell
call m % get_indices_from_bin(i, ijk)
call filter_matches(i_filter_mesh) % bins_set_data(1, i)
do l = 1, n
if (print_ebin) then
call filter_matches(i_filter_ein) % bins_set_data(1, l)
end if
! Left Surface
call filter_matches(i_filter_surf) % bins_set_data(1, OUT_LEFT)
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins_data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
end if
if (trigger % rel_err < rel_err) then
trigger % rel_err = rel_err
end if
trigger % variance = std_dev**2
! Right Surface
call filter_matches(i_filter_surf) % bins_set_data(1, OUT_RIGHT)
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins_data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
end if
if (trigger % rel_err < rel_err) then
trigger % rel_err = rel_err
end if
trigger % variance = trigger % std_dev**2
! Back Surface
call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BACK)
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins_data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
end if
if (trigger % rel_err < rel_err) then
trigger % rel_err = rel_err
end if
trigger % variance = trigger % std_dev**2
! Front Surface
call filter_matches(i_filter_surf) % bins_set_data(1, OUT_FRONT)
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins_data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
end if
if (trigger % rel_err < rel_err) then
trigger % rel_err = rel_err
end if
trigger % variance = trigger % std_dev**2
! Bottom Surface
call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BOTTOM)
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins_data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
end if
if (trigger % rel_err < rel_err) then
trigger % rel_err = rel_err
end if
trigger % variance = trigger % std_dev**2
! Top Surface
call filter_matches(i_filter_surf) % bins_set_data(1, OUT_TOP)
filter_index = 1
do j = 1, size(t % filter)
filter_index = filter_index + (filter_matches(t % filter(j)) % &
bins_data(1) - 1) * t % stride(j)
end do
call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
if (trigger % std_dev < std_dev) then
trigger % std_dev = std_dev
end if
if (trigger % rel_err < rel_err) then
trigger % rel_err = rel_err
end if
trigger % variance = trigger % std_dev**2
end do
end do
end subroutine compute_tally_current
!===============================================================================
! GET_TRIGGER_UNCERTAINTY computes the standard deviation and relative error
! for a single tally bin for CHECK_TALLY_TRIGGERS.
!===============================================================================
subroutine get_trigger_uncertainty(std_dev, rel_err, score_index, &
filter_index, t)
real(8), intent(inout) :: std_dev ! tally standard deviation
real(8), intent(inout) :: rel_err ! tally relative error
integer, intent(in) :: score_index ! tally results score index
integer, intent(in) :: filter_index ! tally results filter index
type(TallyObject), intent(in) :: t ! tally
integer :: n ! number of realizations
real(8) :: mean ! tally mean
real(8) :: tally_sum, tally_sum_sq ! results for a single tally bin
n = t % n_realizations
tally_sum = t % results(RESULT_SUM, score_index, filter_index)
tally_sum_sq = t % results(RESULT_SUM_SQ, score_index, filter_index)
! Compute the tally mean and standard deviation
mean = tally_sum / n
std_dev = sqrt((tally_sum_sq / n - mean * mean) / (n - 1))
! Compute the relative error if the mean is non-zero
if (mean == ZERO) then
rel_err = ZERO
else
rel_err = std_dev / mean
end if
end subroutine get_trigger_uncertainty
end module trigger

201
src/tallies/trigger.cpp Normal file
View file

@ -0,0 +1,201 @@
#include "openmc/tallies/trigger.h"
#include <cmath>
#include <sstream>
#include <utility> // for std::pair
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/reaction.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/tally.h"
namespace openmc {
//==============================================================================
// Global variable definitions
//==============================================================================
namespace settings {
KTrigger keff_trigger;
}
//==============================================================================
// Non-member functions
//==============================================================================
std::pair<double, double>
get_tally_uncertainty(int i_tally, int score_index, int filter_index)
{
int n;
//TODO: off-by-one
int err = openmc_tally_get_n_realizations(i_tally+1, &n);
auto results = tally_results(i_tally);
auto sum = results(filter_index, score_index, RESULT_SUM);
auto sum_sq = results(filter_index, score_index, RESULT_SUM_SQ);
auto mean = sum / n;
double std_dev = std::sqrt((sum_sq/n - mean*mean) / (n-1));
double rel_err = (mean != 0.) ? std_dev / std::abs(mean) : 0.;
return {std_dev, rel_err};
}
//! Find the limiting limiting tally trigger.
//
//! param[out] ratio The uncertainty/threshold ratio for the most limiting
//! tally trigger
//! param[out] tally_id The ID number of the most limiting tally
//! param[out] score The most limiting tally score bin
void
check_tally_triggers(double& ratio, int& tally_id, int& score)
{
ratio = 0.;
for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) {
const Tally& t {*model::tallies[i_tally]};
// Ignore tallies with less than two realizations.
int n_reals;
//TODO: off-by-one
int err = openmc_tally_get_n_realizations(i_tally+1, &n_reals);
if (n_reals < 2) continue;
for (const auto& trigger : t.triggers_) {
const auto& results = tally_results(i_tally);
for (auto filter_index = 0; filter_index < results.shape()[0];
++filter_index) {
for (auto score_index = 0; score_index < results.shape()[1];
++score_index) {
// Compute the tally uncertainty metrics.
auto uncert_pair = get_tally_uncertainty(i_tally, score_index,
filter_index);
double std_dev = uncert_pair.first;
double rel_err = uncert_pair.second;
// Pick out the relevant uncertainty metric for this trigger.
double uncertainty;
switch (trigger.metric) {
case TriggerMetric::variance:
uncertainty = std_dev * std_dev;
break;
case TriggerMetric::standard_deviation:
uncertainty = std_dev;
break;
case TriggerMetric::relative_error:
uncertainty = rel_err;
break;
}
// Compute the uncertainty / threshold ratio.
double this_ratio = uncertainty / trigger.threshold;
if (trigger.metric == TriggerMetric::variance) {
this_ratio = std::sqrt(ratio);
}
// If this is the most uncertain value, set the output variables.
if (this_ratio > ratio) {
ratio = this_ratio;
score = t.scores_[trigger.score_index];
tally_id = t.id_;
}
}
}
}
}
}
//! Compute the uncertainty/threshold ratio for the eigenvalue trigger.
double
check_keff_trigger()
{
if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.;
if (settings::keff_trigger.metric == TriggerMetric::not_active) return 0.;
double k_combined[2];
int err = openmc_get_keff(k_combined);
double uncertainty = 0.;
switch (settings::keff_trigger.metric) {
case TriggerMetric::variance:
uncertainty = k_combined[1] * k_combined[1];
break;
case TriggerMetric::standard_deviation:
uncertainty = k_combined[1];
break;
case TriggerMetric::relative_error:
uncertainty = k_combined[1] / k_combined[0];
}
double ratio = uncertainty / settings::keff_trigger.threshold;
if (settings::keff_trigger.metric == TriggerMetric::variance)
ratio = std::sqrt(ratio);
return ratio;
}
//! See if tally and eigenvalue uncertainties are under trigger thresholds.
void
check_triggers()
{
// Make some aliases.
const auto current_batch {simulation::current_batch};
const auto n_batches {settings::n_batches};
const auto interval {settings::trigger_batch_interval};
// See if the current batch is one for which the triggers must be checked.
if (!settings::trigger_on) return;
if (current_batch < n_batches) return;
if (((current_batch - n_batches) % interval) != 0) return;
// Check the eigenvalue and tally triggers.
double keff_ratio = check_keff_trigger();
double tally_ratio;
int tally_id, score;
check_tally_triggers(tally_ratio, tally_id, score);
// If all the triggers are satisfied, alert the user and return.
if (std::max(keff_ratio, tally_ratio) <= 1.) {
simulation::satisfy_triggers = true;
write_message("Triggers satisfied for batch "
+ std::to_string(current_batch), 7);
return;
}
// At least one trigger is unsatisfied. Let the user know which one.
simulation::satisfy_triggers = false;
std::stringstream msg;
msg << "Triggers unsatisfied, max unc./thresh. is ";
if (keff_ratio >= tally_ratio) {
msg << keff_ratio << " for eigenvalue";
} else {
msg << tally_ratio << " for " << reaction_name(score) << " in tally "
<< tally_id;
}
write_message(msg, 7);
// Estimate batches til triggers are satisfied.
if (settings::trigger_predict) {
// This calculation assumes tally variance is proportional to 1/N where N is
// the number of batches.
auto max_ratio = std::max(keff_ratio, tally_ratio);
auto n_active = current_batch - settings::n_inactive;
auto n_pred_batches = static_cast<int>(n_active * max_ratio * max_ratio)
+ settings::n_inactive + 1;
std::stringstream msg;
msg << "The estimated number of batches is " << n_pred_batches;
if (n_pred_batches > settings::n_max_batches) {
msg << " --- greater than max batches";
warning(msg);
} else {
write_message(msg, 7);
}
}
}
} // namespace openmc

34
src/tallies/trigger_header.F90
View file

@ -1,34 +0,0 @@
module trigger_header
use, intrinsic :: ISO_C_BINDING
use constants, only: NONE, N_FILTER_TYPES, ZERO
implicit none
private
!===============================================================================
! TRIGGEROBJECT stores the variance, relative error and standard deviation
! for some user-specified trigger.
!===============================================================================
type, public :: TriggerObject
integer :: type ! "variance", "std_dev" or "rel_err"
real(8) :: threshold ! a convergence threshold
character(len=52) :: score_name ! the name of the score
integer :: score_index ! the index of the score
real(8) :: variance = ZERO ! temp variance container
real(8) :: std_dev = ZERO ! temp std. dev. container
real(8) :: rel_err = ZERO ! temp rel. err. container
end type TriggerObject
!===============================================================================
! KTRIGGER describes a user-specified precision trigger for k-effective
!===============================================================================
type, public, bind(C) :: KTrigger
integer(C_INT) :: trigger_type = 0
real(C_DOUBLE) :: threshold = ZERO
end type KTrigger
type(KTrigger), public, bind(C) :: keff_trigger ! trigger for k-effective
end module trigger_header

34
src/tracking.F90
View file

@ -24,6 +24,7 @@ module tracking
use tally_header
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_tally, &
score_meshsurface_tally, &
score_track_derivative, zero_flux_derivs, &
score_collision_derivative
use track_output, only: initialize_particle_track, write_particle_track, &
@ -86,7 +87,7 @@ contains
endif
! Every particle starts with no accumulated flux derivative.
if (active_tallies % size() > 0) call zero_flux_derivs()
if (active_tallies_size() > 0) call zero_flux_derivs()
EVENT_LOOP: do
! Set the random number stream
@ -172,7 +173,7 @@ contains
end do
! Score track-length tallies
if (active_tracklength_tallies % size() > 0) then
if (active_tracklength_tallies_size() > 0) then
call score_tracklength_tally(p, distance)
end if
@ -183,7 +184,7 @@ contains
end if
! Score flux derivative accumulators for differential tallies.
if (active_tallies % size() > 0) call score_track_derivative(p, distance)
if (active_tallies_size() > 0) call score_track_derivative(p, distance)
if (d_collision > d_boundary) then
! ====================================================================
@ -210,8 +211,7 @@ contains
p % event = EVENT_SURFACE
end if
! Score cell to cell partial currents
if(active_surface_tallies % size() > 0) &
call score_surface_tally(p, active_surface_tallies)
if (active_surface_tallies_size() > 0) call score_surface_tally(p)
else
! ====================================================================
! PARTICLE HAS COLLISION
@ -226,8 +226,7 @@ contains
! since the direction of the particle will change and we need to use the
! pre-collision direction to figure out what mesh surfaces were crossed
if (active_meshsurf_tallies % size() > 0) &
call score_surface_tally(p, active_meshsurf_tallies)
if (active_meshsurf_tallies_size() > 0) call score_meshsurface_tally(p)
! Clear surface component
p % surface = ERROR_INT
@ -241,8 +240,8 @@ contains
! Score collision estimator tallies -- this is done after a collision
! has occurred rather than before because we need information on the
! outgoing energy for any tallies with an outgoing energy filter
if (active_collision_tallies % size() > 0) call score_collision_tally(p)
if (active_analog_tallies % size() > 0) call score_analog_tally(p)
if (active_collision_tallies_size() > 0) call score_collision_tally(p)
if (active_analog_tallies_size() > 0) call score_analog_tally(p)
! Reset banked weight during collision
p % n_bank = 0
@ -273,7 +272,7 @@ contains
end do
! Score flux derivative accumulators for differential tallies.
if (active_tallies % size() > 0) call score_collision_derivative(p)
if (active_tallies_size() > 0) call score_collision_derivative(p)
end if
! If particle has too many events, display warning and kill it
@ -346,12 +345,12 @@ contains
! forward slightly so that if the mesh boundary is on the surface, it is
! still processed
if (active_meshsurf_tallies % size() > 0) then
if (active_meshsurf_tallies_size() > 0) then
! TODO: Find a better solution to score surface currents than
! physically moving the particle forward slightly
p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw
call score_surface_tally(p, active_meshsurf_tallies)
call score_meshsurface_tally(p)
end if
! Score to global leakage tally
@ -382,14 +381,13 @@ contains
! the particle slightly back in case the surface crossing is coincident
! with a mesh boundary
if(active_surface_tallies % size() > 0) &
call score_surface_tally(p, active_surface_tallies)
if (active_surface_tallies_size() > 0) call score_surface_tally(p)
if (active_meshsurf_tallies % size() > 0) then
if (active_meshsurf_tallies_size() > 0) then
xyz = p % coord(1) % xyz
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
call score_surface_tally(p, active_meshsurf_tallies)
call score_meshsurface_tally(p)
p % coord(1) % xyz = xyz
end if
@ -443,10 +441,10 @@ contains
! Score surface currents since reflection causes the direction of the
! particle to change -- artificially move the particle slightly back in
! case the surface crossing is coincident with a mesh boundary
if (active_meshsurf_tallies % size() > 0) then
if (active_meshsurf_tallies_size() > 0) then
xyz = p % coord(1) % xyz
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
call score_surface_tally(p, active_meshsurf_tallies)
call score_meshsurface_tally(p)
p % coord(1) % xyz = xyz
end if

1
tests/unit_tests/test_capi.py
View file

@ -156,6 +156,7 @@ def test_tally_mapping(capi_init):
def test_tally(capi_init):
t = openmc.capi.tallies[1]
assert t.type == 'volume'
assert len(t.filters) == 2
assert isinstance(t.filters[0], openmc.capi.MaterialFilter)
assert isinstance(t.filters[1], openmc.capi.EnergyFilter)