Merged with develop and tally triggers

This commit is contained in:
Will Boyd 2015-05-07 17:05:19 -04:00
commit f49749afb3
47 changed files with 2280 additions and 70 deletions

View file

@ -139,6 +139,8 @@ default. This element has the following attributes/sub-elements:
*Default*: 1.0
.. _eigenvalue:
``<eigenvalue>`` Element
------------------------
@ -170,6 +172,34 @@ should be performed. It has the following attributes/sub-elements:
*Default*: None
:keff_trigger:
This tag specifies a precision trigger on the combined :math:`k_{eff}`. The
trigger is a convergence criterion on the uncertainty of the estimated
eigenvalue. It has the following attributes/sub-elements:
:type:
The type of precision trigger. Accepted options are "variance", "std_dev",
and "rel_err".
:variance:
Variance of the batch mean :math:`\sigma^2`
:std_dev:
Standard deviation of the batch mean :math:`\sigma`
:rel_err:
Relative error of the batch mean :math:`\frac{\sigma}{\mu}`
*Default*: None
:threshold:
The precision trigger's convergence criterion for the
combined :math:`k_{eff}`.
*Default*: None
.. note:: See section on the :ref:`trigger` for more information.
``<energy_grid>`` Element
-------------------------
@ -214,8 +244,8 @@ problem. It has the following attributes/sub-elements:
``<fixed_source>`` Element
--------------------------
The ``<fixed_source>`` element indicates that a fixed source calculation should be
performed. It has the following attributes/sub-elements:
The ``<fixed_source>`` element indicates that a fixed source calculation should
be performed. It has the following attributes/sub-elements:
:batches:
The total number of batches. For fixed source calculations, each batch
@ -296,7 +326,8 @@ out the file and "false" will not.
*Default*: true
.. note:: The tally results will always be written to a binary/HDF5 state point file.
.. note:: The tally results will always be written to a binary/HDF5 state
point file.
``<output_path>`` Element
-------------------------
@ -545,8 +576,8 @@ attributes/sub-elements:
*Default*: None
:separate:
If this element is set to "true", a separate binary source point file will be
written. Otherwise, the source sites will be written in the state point
If this element is set to "true", a separate binary source point file will
be written. Otherwise, the source sites will be written in the state point
directly.
*Default*: false
@ -601,10 +632,56 @@ integers: the batch number, generation number, and particle number.
``<track>`` Element
-------------------
The ``<track>`` element specifies particles for which OpenMC will output binary files describing particle position at every step of its transport. This element should be followed by triplets of integers. Each triplet describes one particle. The integers in each triplet specify the batch number, generation number, and particle number, respectively.
The ``<track>`` element specifies particles for which OpenMC will output binary
files describing particle position at every step of its transport. This element
should be followed by triplets of integers. Each triplet describes one particle
. The integers in each triplet specify the batch number, generation number, and
particle number, respectively.
*Default*: None
.. _trigger:
``<trigger>`` Element
-------------------------
OpenMC includes tally precision triggers which allow the user to define
uncertainty thresholds on :math:`k_{eff}` in the ``<eigenvalue>`` subelement of
``settings.xml``, and/or tallies in ``tallies.xml``. When using triggers,
OpenMC will run until it completes as many batches as defined by ``<batches>``.
At this point, the uncertainties on all tallied values are computed and
compared with their corresponding trigger thresholds. If any triggers have not
been met, OpenMC will continue until either all trigger thresholds have been
satisfied or ``<max_batches>`` has been reached.
The ``<trigger>`` element provides an active "toggle switch" for tally
precision trigger(s), the maximum number of batches and the batch interval. It
has the following attributes/sub-elements:
:active:
This determines whether or not to use trigger(s). Trigger(s) are used when
this tag is set to "true".
:max_batches:
This describes the maximum number of batches allowed when using trigger(s).
.. note:: When max_batches is set, the number of ``batches`` shown in
``<eigenvalue>`` element represents minimum number of batches to
simulate when using the trigger(s).
:batch_interval:
This tag describes the number of batches in between convergence checks.
OpenMC will check if the trigger has been reached at each batch defined
by ``batch_interval`` after the minimum number of batches is reached.
.. note:: If this tag is not present, the ``batch_interval`` is predicted
dynamically by OpenMC for each convergence check. The predictive
model assumes no correlation between fission sources
distributions from batch-to-batch. This assumption is reasonable
for fixed source and small criticality calculations, but is very
optimistic for highly coupled full-core reactor problems.
``<uniform_fs>`` Element
------------------------
@ -1116,6 +1193,8 @@ post-collision energy, and an arbitrary structured mesh.
The three valid elements in the tallies.xml file are ``<tally>``, ``<mesh>``,
and ``<assume_separate>``.
.. _tally:
``<tally>`` Element
-------------------
@ -1139,11 +1218,13 @@ The ``<tally>`` element accepts the following sub-elements:
The ``filter`` element has the following attributes/sub-elements:
:type:
The type of the filter. Accepted options are "cell", "cellborn", "material",
"universe", "energy", "energyout", "mesh", and "distribcell".
The type of the filter. Accepted options are "cell", "cellborn",
"material", "universe", "energy", "energyout", "mesh", and
"distribcell".
:bins:
For each filter type, the corresponding ``bins`` entry is given as follows:
For each filter type, the corresponding ``bins`` entry is given as
follows:
:cell:
A list of cells in which the tally should be accumulated.
@ -1257,20 +1338,20 @@ The ``<tally>`` element accepts the following sub-elements:
:scatter-PN:
Tally all of the scattering moments from order 0 to N, where N is the
Legendre expansion order of the change in particle angle :math:`\left(\mu\right)`.
That is, ``scatter-P1`` is equivalent to requesting tallies of
``scatter-0`` and ``scatter-1``. Like for ``scatter-N``,
N must be between 0 and 10. As an example, tallying up to the
2\ :sup:`nd` \ scattering moment would be specified as
Legendre expansion order of the change in particle angle
:math:`\left(\mu\right)`. That is, ``scatter-P1`` is equivalent to
requesting tallies of ``scatter-0`` and ``scatter-1``. Like for
``scatter-N``, N must be between 0 and 10. As an example, tallying up
to the 2\ :sup:`nd` \ scattering moment would be specified as
``<scores> scatter-P2 </scores>``.
:scatter-YN:
``scatter-YN`` is similar to ``scatter-PN`` except an additional
expansion is performed for the incoming particle direction
:math:`\left(\Omega\right)` using the real spherical harmonics. This is useful
for performing angular flux moment weighting of the scattering moments.
Like ``scatter-PN``, ``scatter-YN`` will tally all of the moments from
order 0 to N; N again must be between 0 and 10.
:math:`\left(\Omega\right)` using the real spherical harmonics. This is
useful for performing angular flux moment weighting of the scattering
moments. Like ``scatter-PN``, ``scatter-YN`` will tally all of the
moments from order 0 to N; N again must be between 0 and 10.
:nu-scatter, nu-scatter-N, nu-scatter-PN, nu-scatter-YN:
These scores are similar in functionality to their ``scatter*``
@ -1300,6 +1381,40 @@ The ``<tally>`` element accepts the following sub-elements:
:events:
Number of scoring events
:trigger:
Precision trigger applied to all filter bins and nuclides for this tally.
It must specify the trigger's type, threshold and scores to which it will
be applied. It has the following attributes/sub-elements:
:type:
The type of the trigger. Accepted options are "variance", "std_dev",
and "rel_err".
:variance:
Variance of the batch mean :math:`\sigma^2`
:std_dev:
Standard deviation of the batch mean :math:`\sigma`
:rel_err:
Relative error of the batch mean :math:`\frac{\sigma}{\mu}`
*Default*: None
:threshold:
The precision trigger's convergence criterion for tallied values.
*Default*: None
:scores:
The score(s) in this tally to which the trigger should be applied.
.. note:: The ``scores`` in ``trigger`` must have been defined in
``scores`` in ``tally``. An optional "all" may be used to
select all scores in this tally.
*Default*: "all"
``<mesh>`` Element
------------------
@ -1338,8 +1453,8 @@ overhead. The effect of assuming all tallies are spatially separate is that once
one tally is scored to, the same event is assumed not to score to any other
tallies. This element should be followed by "true" or "false".
.. warning:: If used incorrectly, the assumption that all tallies are spatially
separate can lead to incorrect results.
.. warning:: If used incorrectly, the assumption that all tallies are
spatially separate can lead to incorrect results.
*Default*: false
@ -1355,8 +1470,10 @@ element of the plots.xml is simply ``<plots>`` and any number output plots can
be defined with ``<plot>`` sub-elements. Two plot types are currently
implemented in openMC:
* ``slice`` 2D pixel plot along one of the major axes. Produces a PPM image file.
* ``voxel`` 3D voxel data dump. Produces a binary file containing voxel xyz position and cell or material id.
* ``slice`` 2D pixel plot along one of the major axes. Produces a PPM image
file.
* ``voxel`` 3D voxel data dump. Produces a binary file containing voxel xyz
position and cell or material id.
``<plot>`` Element
@ -1540,9 +1657,9 @@ attributes or sub-elements. These are not used in "voxel" plots:
CMFD Specification -- cmfd.xml
------------------------------
Coarse mesh finite difference acceleration method has been implemented in OpenMC.
Currently, it allows users to accelerate fission source convergence during
inactive neutron batches. To run CMFD, the ``<run_cmfd>`` element in
Coarse mesh finite difference acceleration method has been implemented in
OpenMC. Currently, it allows users to accelerate fission source convergence
during inactive neutron batches. To run CMFD, the ``<run_cmfd>`` element in
``settings.xml`` should be set to "true".
``<begin>`` Element
@ -1557,8 +1674,8 @@ The ``<begin>`` element controls what batch CMFD calculations should begin.
The ``<display>`` element sets one additional CMFD output column. Options are:
* "balance" - prints the RMS [%] of the resdiual from the neutron balance equation
on CMFD tallies.
* "balance" - prints the RMS [%] of the resdiual from the neutron balance
equation on CMFD tallies.
* "dominance" - prints the estimated dominance ratio from the CMFD iterations.
**This will only work for power iteration eigensolver**.
* "entropy" - prints the *entropy* of the CMFD predicted fission source.
@ -1694,8 +1811,8 @@ not impact the calculation.
``<power_monitor>`` Element
---------------------------
The ``<power_monitor>`` element is used to view the convergence of power iteration.
This option can be turned on with "true" and turned off with "false".
The ``<power_monitor>`` element is used to view the convergence of power
iteration. This option can be turned on with "true" and turned off with "false".
*Default*: false

View file

@ -129,6 +129,9 @@ settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.keff_trigger = {'type' : 'std_dev', 'threshold' : 5E-4}
settings_file.trigger_active = True
settings_file.trigger_max_batches = 100
settings_file.export_to_xml()

View file

@ -130,6 +130,8 @@ settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
settings_file.set_source_space('box', [-1, -1, -1, 1, 1, 1])
settings_file.trigger_active = True
settings_file.trigger_max_batches = 100
settings_file.export_to_xml()
@ -164,10 +166,15 @@ mesh.width = [1, 1]
mesh_filter = openmc.Filter()
mesh_filter.mesh = mesh
# Instantiate tally Trigger
trigger = openmc.Trigger(trigger_type='rel_err', threshold=1E-2)
trigger.add_score('all')
# Instantiate the Tally
tally = openmc.Tally(tally_id=1)
tally.add_filter(mesh_filter)
tally.add_score('total')
tally.add_trigger(trigger)
# Instantiate a TalliesFile, register Tally/Mesh, and export to XML
tallies_file = openmc.TalliesFile()

View file

@ -317,7 +317,13 @@ module constants
OUT_FRONT = 4, &
IN_TOP = 5, &
OUT_TOP = 6
! 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 :: &

View file

@ -21,6 +21,7 @@ module eigenvalue
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
use trigger, only: check_triggers
use tracking, only: transport
implicit none
@ -53,7 +54,7 @@ contains
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_batches
BATCH_LOOP: do current_batch = 1, n_max_batches
call initialize_batch()
@ -95,6 +96,8 @@ contains
end do GENERATION_LOOP
call finalize_batch()
if (satisfy_triggers) exit BATCH_LOOP
end do BATCH_LOOP
@ -215,12 +218,22 @@ contains
! Display output
if (master) call print_batch_keff()
! Calculate combined estimate of k-effective
if (master) call calculate_combined_keff()
! Check_triggers
if (master) call check_triggers()
#ifdef MPI
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
MPI_COMM_WORLD, mpi_err)
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
call statepoint_batch % add(current_batch)
end if
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
! Calculate combined estimate of k-effective
if (master) call calculate_combined_keff()
! Create state point file
call write_state_point()
end if
@ -230,7 +243,7 @@ contains
call write_source_point()
end if
if (master .and. current_batch == n_batches) then
if (master .and. current_batch == n_max_batches) then
! Make sure combined estimate of k-effective is calculated at the last
! batch in case no state point is written
call calculate_combined_keff()
@ -491,7 +504,7 @@ contains
call MPI_WAITALL(n_request, request, MPI_STATUSES_IGNORE, mpi_err)
! Deallocate space for bank_position on the very last generation
if (current_batch == n_batches .and. current_gen == gen_per_batch) &
if (current_batch == n_max_batches .and. current_gen == gen_per_batch) &
deallocate(bank_position)
#else
source_bank = temp_sites(1:n_particles)
@ -500,7 +513,7 @@ contains
call time_bank_sendrecv % stop()
! Deallocate space for the temporary source bank on the last generation
if (current_batch == n_batches .and. current_gen == gen_per_batch) &
if (current_batch == n_max_batches .and. current_gen == gen_per_batch) &
deallocate(temp_sites)
end subroutine synchronize_bank

View file

@ -1,5 +1,9 @@
module fixed_source
#ifdef MPI
use mpi
#endif
use constants, only: ZERO, MAX_LINE_LEN
use global
use output, only: write_message, header
@ -9,6 +13,7 @@ module fixed_source
use state_point, only: write_state_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies
use trigger, only: check_triggers
use tracking, only: transport
implicit none
@ -36,7 +41,7 @@ contains
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_batches
BATCH_LOOP: do current_batch = 1, n_max_batches
! In a restart run, skip any batches that have already been simulated
if (restart_run .and. current_batch <= restart_batch) then
@ -79,6 +84,8 @@ contains
call finalize_batch()
if (satisfy_triggers) exit BATCH_LOOP
end do BATCH_LOOP
call time_active % stop()
@ -115,6 +122,17 @@ contains
call synchronize_tallies()
call time_tallies % stop()
! Check_triggers
if (master) call check_triggers()
#ifdef MPI
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
MPI_COMM_WORLD, mpi_err)
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
call statepoint_batch % add(current_batch)
end if
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
call write_state_point()

View file

@ -13,6 +13,7 @@ module global
use set_header, only: SetInt
use source_header, only: ExtSource
use tally_header, only: TallyObject, TallyMap, TallyResult
use trigger_header, only: KTrigger
use timer_header, only: Timer
#ifdef HDF5
@ -150,6 +151,16 @@ module global
integer :: current_gen = 0 ! current generation within a batch
integer :: overall_gen = 0 ! overall generation in the run
! ============================================================================
! TALLY PRECISION TRIGGER VARIABLES
integer :: n_max_batches ! max # of batches
integer :: n_batch_interval = 1 ! batch interval for triggers
logical :: pred_batches = .false. ! predict batches for triggers
logical :: trigger_on = .false. ! flag for turning triggers on/off
type(KTrigger) :: keff_trigger ! trigger for k-effective
logical :: satisfy_triggers = .false. ! whether triggers are satisfied
! External source
type(ExtSource), target :: external_source

View file

@ -74,6 +74,8 @@ contains
type(Node), pointer :: node_verb => null()
type(Node), pointer :: node_res_scat => null()
type(Node), pointer :: node_scatterer => null()
type(Node), pointer :: node_trigger => null()
type(Node), pointer :: node_keff_trigger => null()
type(NodeList), pointer :: node_scat_list => null()
! Display output message
@ -125,6 +127,43 @@ contains
path_output = trim(path_output) // "/"
end if
! Check for a trigger node and get trigger information
if (check_for_node(doc, "trigger")) then
call get_node_ptr(doc, "trigger", node_trigger)
! Check if trigger(s) are to be turned on
call get_node_value(node_trigger, "active", temp_str)
temp_str = trim(to_lower(temp_str))
if (temp_str == 'true' .or. temp_str == '1') then
trigger_on = .true.
elseif (temp_str == 'false' .or. temp_str == '0') then
trigger_on = .false.
else
call fatal_error("Unrecognized trigger active: " // temp_str)
end if
if (trigger_on) then
if (check_for_node(node_trigger, "max_batches") )then
call get_node_value(node_trigger, "max_batches", n_max_batches)
else
call fatal_error("The max_batches must be specified with triggers")
end if
! Get the batch interval to check triggers
if (.not. check_for_node(node_trigger, "batch_interval"))then
pred_batches = .true.
else
call get_node_value(node_trigger, "batch_interval", temp_int)
n_batch_interval = temp_int
if (n_batch_interval <= 0) then
call fatal_error("The batch interval must be greater than zero")
end if
end if
end if
end if
! Make sure that either eigenvalue or fixed source was specified
if (.not.check_for_node(doc, "eigenvalue") .and. &
.not.check_for_node(doc, "fixed_source")) then
@ -151,8 +190,13 @@ contains
! don't set it here
if (n_particles == 0) n_particles = temp_long
! Copy batch and generation information
! Get number of basic batches
call get_node_value(node_mode, "batches", n_batches)
if (.not. trigger_on) then
n_max_batches = n_batches
end if
! Get number of inactive batches
call get_node_value(node_mode, "inactive", n_inactive)
n_active = n_batches - n_inactive
if (check_for_node(node_mode, "generations_per_batch")) then
@ -160,9 +204,40 @@ contains
end if
! Allocate array for batch keff and entropy
allocate(k_generation(n_batches*gen_per_batch))
allocate(entropy(n_batches*gen_per_batch))
allocate(k_generation(n_max_batches*gen_per_batch))
allocate(entropy(n_max_batches*gen_per_batch))
entropy = ZERO
! Get the trigger information for keff
if (check_for_node(node_mode, "keff_trigger")) then
call get_node_ptr(node_mode, "keff_trigger", node_keff_trigger)
if (check_for_node(node_keff_trigger, "type")) then
call get_node_value(node_keff_trigger, "type", temp_str)
temp_str = trim(to_lower(temp_str))
select case (temp_str)
case ('std_dev')
keff_trigger % trigger_type = STANDARD_DEVIATION
case ('variance')
keff_trigger % trigger_type = VARIANCE
case ('rel_err')
keff_trigger % trigger_type = RELATIVE_ERROR
case default
call fatal_error("Unrecognized keff trigger type " // temp_str)
end select
else
call fatal_error("Specify keff trigger type in settings XML")
end if
if (check_for_node(node_keff_trigger, "threshold")) then
call get_node_value(node_keff_trigger, "threshold", &
keff_trigger % threshold)
else
call fatal_error("Specify keff trigger threshold in settings XML")
end if
end if
end if
! Fixed source calculation information
@ -187,6 +262,9 @@ contains
! Copy batch information
call get_node_value(node_mode, "batches", n_batches)
if (.not. trigger_on) then
n_max_batches = n_batches
end if
n_active = n_batches
n_inactive = 0
gen_per_batch = 1
@ -1995,9 +2073,13 @@ contains
integer :: n_filters ! number of filters
integer :: n_new ! number of new scores to add based on Yn/Pn tally
integer :: n_scores ! number of tot scores after adjusting for Yn/Pn tally
integer :: n_bins ! Total new bins for this score
integer :: n_order ! Moment order requested
integer :: n_order_pos ! Position of Scattering order in score name string
integer :: n_bins ! total new bins for this score
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
real(8) :: threshold ! trigger convergence threshold
integer :: n_order ! moment order requested
integer :: n_order_pos ! oosition of Scattering order in score name string
integer :: MT ! user-specified MT for score
integer :: iarray3(3) ! temporary integer array
integer :: imomstr ! Index of MOMENT_STRS & MOMENT_N_STRS
@ -2008,6 +2090,7 @@ contains
character(MAX_WORD_LEN) :: score_name
character(MAX_WORD_LEN) :: temp_str
character(MAX_WORD_LEN), allocatable :: sarray(:)
type(DictCharInt) :: trigger_scores
type(ElemKeyValueCI), pointer :: pair_list => null()
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
@ -2016,9 +2099,13 @@ contains
type(Node), pointer :: node_mesh => null()
type(Node), pointer :: node_tal => null()
type(Node), pointer :: node_filt => null()
type(Node), pointer :: node_trigger=>null()
type(NodeList), pointer :: node_mesh_list => null()
type(NodeList), pointer :: node_tal_list => null()
type(NodeList), pointer :: node_filt_list => null()
type(NodeList), pointer :: node_trigger_list => null()
type(ElemKeyValueCI), pointer :: scores
type(ElemKeyValueCI), pointer :: next
! Check if tallies.xml exists
filename = trim(path_input) // "tallies.xml"
@ -2774,9 +2861,8 @@ contains
! 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 scoring functions in the &
&same tally as surface currents. Separate other scoring &
&functions into a distinct tally.")
call fatal_error("Cannot tally other scores in the &
&same tally as surface currents")
end if
! Since the number of bins for the mesh filter was already set
@ -2846,6 +2932,10 @@ contains
end if
end select
! Append the score to the list of possible trigger scores
call trigger_scores % add_key(trim(score_name), l)
end do
t % n_score_bins = n_scores
t % n_user_score_bins = n_words
@ -2856,6 +2946,184 @@ contains
call fatal_error("No <scores> specified on tally " &
&// trim(to_str(t % id)) // ".")
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 = get_list_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
call get_list_item(node_trigger_list, user_trig_ind, node_trigger)
! Get scores for this trigger
if (check_for_node(node_trigger, "scores")) then
n_words = get_arraysize_string(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
scores => trigger_scores % keys()
do while (associated(scores))
next => scores % next
deallocate(scores)
scores => next
t % n_triggers = t % n_triggers + 1
end do
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
call get_list_item(node_trigger_list, user_trig_ind, node_trigger)
! 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 = get_arraysize_string(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
scores => trigger_scores % keys()
! Loop over all tally scores
do while (associated(scores))
score_name = trim(scores % key)
! Store the score name and index in the trigger
t % triggers(trig_ind) % score_name = trim(score_name)
t % triggers(trig_ind) % score_index = &
trigger_scores % get_key(trim(score_name))
! 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
! Move to next score
next => scores % next
deallocate(scores)
scores => next
! 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_key(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()
end if
! =======================================================================
! SET TALLY ESTIMATOR

View file

@ -1951,7 +1951,7 @@ contains
#endif
#ifdef HDF5
# ifdef MPI
#ifdef MPI
! Set size of total dataspace for all procs and rank
dims1(1) = n_particles

View file

@ -10,7 +10,12 @@ element settings {
(element particles { xsd:positiveInteger } |
attribute particles { xsd:positiveInteger }) &
(element generations_per_batch { xsd:positiveInteger } |
attribute generations_per_batch { xsd:positiveInteger })?
attribute generations_per_batch { xsd:positiveInteger })? &
(element keff_trigger {
(element type { xsd:string } | attribute type { xsd:string }) &
(element threshold { xsd:double} | attribute threshold { xsd:double })
}
)?
} |
element fixed_source {
(element batches { xsd:positiveInteger } |
@ -123,6 +128,12 @@ element settings {
element track { list { xsd:positiveInteger+ } }? &
element trigger {
(element active { xsd:boolean } | attribute active { xsd:boolean }) &
(element max_batches { xsd:positiveInteger } | attribute max_batches { xsd:positiveInteger }) &
(element batch_interval { xsd:positiveInteger } | attribute batch_interval { xsd:positiveInteger })?
}? &
element verbosity { xsd:positiveInteger }? &
element uniform_fs{

View file

@ -43,6 +43,28 @@
</attribute>
</choice>
</optional>
<optional>
<element name="keff_trigger">
<interleave>
<choice>
<element name="type">
<data type="string"/>
</element>
<attribute name="type">
<data type="string"/>
</attribute>
</choice>
<choice>
<element name="threshold">
<data type="double"/>
</element>
<attribute name="threshold">
<data type="double"/>
</attribute>
</choice>
</interleave>
</element>
</optional>
</interleave>
</element>
<element name="fixed_source">
@ -567,6 +589,38 @@
</list>
</element>
</optional>
<optional>
<element name="trigger">
<interleave>
<choice>
<element name="active">
<data type="boolean"/>
</element>
<attribute name="active">
<data type="boolean"/>
</attribute>
</choice>
<choice>
<element name="max_batches">
<data type="positiveInteger"/>
</element>
<attribute name="max_batches">
<data type="positiveInteger"/>
</attribute>
</choice>
<optional>
<choice>
<element name="batch_interval">
<data type="positiveInteger"/>
</element>
<attribute name="batch_interval">
<data type="positiveInteger"/>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
<optional>
<element name="verbosity">
<data type="positiveInteger"/>

View file

@ -34,7 +34,12 @@ element tallies {
}? &
element scores {
list { xsd:string { maxLength = "20" }+ }
}
} &
element trigger {
(element type { xsd:string } | attribute type { xsd:string }) &
(element threshold { xsd:double} | attribute threshold { xsd:double }) &
(element scores { list { xsd:string { maxLength = "20" }+ } } | attribute scores { list { xsd:string { maxLength = "20"}+ } } )?
}?
}* &
element assume_separate { xsd:boolean }?

View file

@ -204,6 +204,50 @@
</oneOrMore>
</list>
</element>
<optional>
<element name="trigger">
<interleave>
<choice>
<element name="type">
<data type="string"/>
</element>
<attribute name="type">
<data type="string"/>
</attribute>
</choice>
<choice>
<element name="threshold">
<data type="double"/>
</element>
<attribute name="threshold">
<data type="double"/>
</attribute>
</choice>
<optional>
<choice>
<element name="scores">
<list>
<oneOrMore>
<data type="string">
<param name="maxLength">20</param>
</data>
</oneOrMore>
</list>
</element>
<attribute name="scores">
<list>
<oneOrMore>
<data type="string">
<param name="maxLength">20</param>
</data>
</oneOrMore>
</list>
</attribute>
</choice>
</optional>
</interleave>
</element>
</optional>
</interleave>
</element>
</zeroOrMore>

View file

@ -52,7 +52,7 @@ contains
! Set filename for state point
filename = trim(path_output) // 'statepoint.' // &
& zero_padded(current_batch, count_digits(n_batches))
& zero_padded(current_batch, count_digits(n_max_batches))
! Append appropriate extension
#ifdef HDF5
@ -412,7 +412,7 @@ contains
! Set filename
filename = trim(path_output) // 'source.' // &
& zero_padded(current_batch, count_digits(n_batches))
& zero_padded(current_batch, count_digits(n_max_batches))
#ifdef HDF5
filename = trim(filename) // '.h5'
@ -434,7 +434,7 @@ contains
! Set filename for state point
filename = trim(path_output) // 'statepoint.' // &
& zero_padded(current_batch, count_digits(n_batches))
& zero_padded(current_batch, count_digits(n_max_batches))
#ifdef HDF5
filename = trim(filename) // '.h5'
#else
@ -530,7 +530,7 @@ contains
#endif
! Transfer values to value on master
if (current_batch == n_batches) then
if (current_batch == n_max_batches .or. satisfy_triggers) then
global_tallies(:) % sum = global_temp(1,:)
global_tallies(:) % sum_sq = global_temp(2,:)
end if
@ -602,7 +602,7 @@ contains
! At the end of the simulation, store the results back in the
! regular TallyResults array
if (current_batch == n_batches) then
if (current_batch == n_max_batches .or. satisfy_triggers) then
tally % results(:,:) % sum = tally_temp(1,:,:)
tally % results(:,:) % sum_sq = tally_temp(2,:,:)
end if

View file

@ -1,6 +1,7 @@
module tally_header
use constants, only: NONE, N_FILTER_TYPES
use constants, only: NONE, N_FILTER_TYPES
use trigger_header, only: TriggerObject
implicit none
@ -123,6 +124,10 @@ module tally_header
! Number of realizations of tally random variables
integer :: n_realizations = 0
! Tally precision triggers
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
! Type-Bound procedures
contains
@ -192,6 +197,11 @@ module tally_header
this % reset = .false.
this % n_realizations = 0
if (allocated(this % triggers)) &
deallocate (this % triggers)
this % n_triggers = 0
end subroutine tallyobject_clear

550
src/trigger.F90 Normal file
View file

@ -0,0 +1,550 @@
module trigger
#ifdef MPI
use mpi
#endif
use global
use string, only: to_str
use output, only: warning, write_message
use mesh, only: mesh_indices_to_bin
use trigger_header, only: TriggerObject
use tally, only: TallyObject
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()
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)))
! 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))
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)))
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)))
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 ! level in tally hierarchy
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 :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
real(8) :: uncertainty ! trigger uncertainty
real(8) :: std_dev = 0.0 ! trigger standard deviation
real(8) :: rel_err = 0.0 ! trigger relative error
real(8) :: ratio ! ratio of the uncertainty/trigger threshold
type(TallyObject), pointer :: t ! tally pointer
type(TriggerObject), pointer :: trigger ! tally trigger
! 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
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
t => tallies(i)
! 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
trigger => t % triggers(s)
! Initialize trigger uncertainties to zero
trigger % std_dev = 0.
trigger % rel_err = 0.
trigger % variance = 0.
! Surface current tally triggers require special treatment
if (t % type == TALLY_SURFACE_CURRENT) then
call compute_tally_current(t, trigger)
else
! Initialize bins, filter level
matching_bins(1:t % n_filters) = 0
j = 1
! Find filter index
FILTER_LOOP: do
find_bin: do
if (t % n_filters == 0) exit find_bin
matching_bins(j) = matching_bins(j) + 1
if (matching_bins(j) > t % filters(j) % n_bins) then
if (j == 1) exit FILTER_LOOP
matching_bins(j) = 0
j = j - 1
else
if (j == t % n_filters) exit find_bin
end if
end do find_bin
if (t % n_filters > 0) then
filter_index = sum((max(matching_bins(1:t%n_filters),1) - 1) * &
t % stride) + 1
else
filter_index = 1
end if
! Initialize score index
score_index = trigger % score_index
! Initialize score bin index
NUCLIDE_LOOP: do n = 1, t % n_nuclide_bins
select case(t % score_bins(trigger % score_index))
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
score_index = score_index - 1
do n_order = 0, t % moment_order(trigger % score_index)
score_index = score_index + 1
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
end do
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
SCORE_TOTAL_YN)
score_index = score_index - 1
do n_order = 0, t % moment_order(trigger % score_index)
do nm_order = -n_order, n_order
score_index = score_index + 1
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
end do
end do
case default
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
end select
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 (t % n_filters == 0) exit FILTER_LOOP
end do FILTER_LOOP
end if
end do TRIGGER_LOOP
end do TALLY_LOOP
end if
end subroutine check_tally_triggers
!===============================================================================
! COMPUTE_TALLY_CURRENT computes the current for a surface current tally with
! precision trigger(s).
!===============================================================================
subroutine compute_tally_current(t, trigger)
integer :: i ! mesh index for x
integer :: j ! mesh index for y
integer :: k ! mesh index for z
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 = 0.0 ! temporary relative error of result
real(8) :: std_dev = 0.0 ! temporary standard deviration of result
type(TallyObject), pointer :: t ! surface current tally
type(TriggerObject) :: trigger ! surface current tally trigger
type(StructuredMesh), pointer :: m ! surface current mesh
! Get pointer to mesh
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_surf = t % find_filter(FILTER_SURFACE)
m => meshes(t % filters(i_filter_mesh) % int_bins(1))
! initialize bins array
matching_bins(1:t % n_filters) = 1
! 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 = t % filters(i_filter_ein) % n_bins
else
print_ebin = .false.
n = 1
end if
do i = 1, m % dimension(1)
do j = 1, m % dimension(2)
do k = 1, m % dimension(3)
do l = 1, n
if (print_ebin) then
matching_bins(i_filter_ein) = l
end if
! Left Surface
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i-1, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
! Right Surface
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_surf) = OUT_RIGHT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j-1, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_surf) = OUT_FRONT
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k-1 /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_surf) = OUT_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_mesh) = &
mesh_indices_to_bin(m, (/ i, j, k /) + 1, .true.)
matching_bins(i_filter_surf) = IN_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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
matching_bins(i_filter_surf) = OUT_TOP
filter_index = &
sum((matching_bins(1:t % n_filters) - 1) * t % stride) + 1
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 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
integer :: n ! number of realizations
real(8) :: mean ! tally mean
type(TallyResult) :: tally_result ! pointer to TallyResult
type(TallyObject), pointer :: t ! tally pointer
n = t % n_realizations
tally_result = t % results(score_index, filter_index)
! Compute the tally mean and standard deviation
mean = tally_result % sum / n
std_dev = sqrt((tally_result % sum_sq / n - mean * mean) / (n - 1))
! Compute the relative error if the mean is non-zero
if (mean == 0.) then
rel_err = 0.
else
rel_err = std_dev / mean
end if
end subroutine get_trigger_uncertainty
end module trigger

29
src/trigger_header.F90 Normal file
View file

@ -0,0 +1,29 @@
module trigger_header
use constants, only: NONE, N_FILTER_TYPES
implicit none
!===============================================================================
! TRIGGEROBJECT stores the variance, relative error and standard deviation
! for some user-specified trigger.
!===============================================================================
type 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=0.0 ! temp variance container
real(8) :: std_dev =0.0 ! temp std. dev. container
real(8) :: rel_err =0.0 ! temp rel. err. container
end type TriggerObject
!===============================================================================
! KTRIGGER describes a user-specified precision trigger for k-effective
!===============================================================================
type KTrigger
integer :: trigger_type = 0
real(8) :: threshold = 0
end type KTrigger
end module trigger_header

View file

@ -17,6 +17,7 @@ class SettingsFile(object):
self._generations_per_batch = None
self._inactive = None
self._particles = None
self._keff_trigger = None
# Source subelement
self._source_subelement = None
@ -41,6 +42,12 @@ class SettingsFile(object):
self._entropy_lower_left = None
self._entropy_upper_right = None
# Trigger subelement
self._trigger_subelement = None
self._trigger_active = None
self._trigger_max_batches = None
self._trigger_batch_interval = None
self._output = None
self._output_path = None
@ -79,7 +86,7 @@ class SettingsFile(object):
self._dd_count_interactions = False
self._settings_file = ET.Element("settings")
self._eigenvalue_element = None
self._eigenvalue_subelement = None
self._source_element = None
@ -103,6 +110,11 @@ class SettingsFile(object):
return self._particles
@property
def keff_trigger(self):
return self._keff_trigger
@property
def source_file(self):
return self._source_file
@ -188,6 +200,21 @@ class SettingsFile(object):
return self._entropy_upper_right
@property
def trigger_active(self):
return self._trigger_active
@property
def trigger_max_batches(self):
return self._trigger_max_batches
@property
def trigger_batch_interval(self):
return self._batch_interval
@property
def output(self):
return self._output
@ -376,6 +403,37 @@ class SettingsFile(object):
self._particles = particles
@keff_trigger.setter
def keff_trigger(self, keff_trigger):
if not isinstance(keff_trigger, dict):
msg = 'Unable to set a trigger on keff from {0} which ' \
'is not a Python dictionary'.format(keff_trigger)
raise ValueError(msg)
elif not 'type' in keff_trigger:
msg = 'Unable to set a trigger on keff from {0} which ' \
'does not have a "type" key'.format(keff_trigger)
raise ValueError(msg)
elif not keff_trigger['type'] in ['variance', 'std_dev', 'rel_err']:
msg = 'Unable to set a trigger on keff with ' \
'type {0}'.format(keff_trigger['type'])
raise ValueError(msg)
elif not 'threshold' in keff_trigger:
msg = 'Unable to set a trigger on keff from {0} which ' \
'does not have a "threshold" key'.format(keff_trigger)
raise ValueError(msg)
elif not is_float(keff_trigger['threshold']):
msg = 'Unable to set a trigger on keff with ' \
'threshold {0}'.format(keff_trigger['threshold'])
raise ValueError(msg)
self._keff_trigger = keff_trigger
@source_file.setter
def source_file(self, source_file):
@ -836,6 +894,50 @@ class SettingsFile(object):
self._entropy_upper_right = upper_right
@trigger_active.setter
def trigger_active(self, trigger_active):
if not isinstance(trigger_active, bool):
msg = 'Unable to set trigger active to a ' \
'non-boolean value {0}'.format(trigger_active)
raise ValueError(msg)
self._trigger_active = trigger_active
@trigger_max_batches.setter
def trigger_max_batches(self, trigger_max_batches):
if not is_integer(trigger_max_batches):
msg = 'Unable to set trigger max batches to a non-integer ' \
'value {0}'.format(trigger_max_batches)
raise ValueError(msg)
elif trigger_max_batches <= 0:
msg = 'Unable to set trigger max batches to a non-positive ' \
'value {0}'.format(trigger_max_batches)
raise ValueError(msg)
self._trigger_max_batches = trigger_max_batches
@trigger_batch_interval.setter
def trigger_batch_interval(self, trigger_batch_interval):
if not is_integer(trigger_batch_interval):
msg = 'Unable to set trigger batch interval to a non-integer ' \
'value {0}'.format(trigger_batch_interval)
raise ValueError(msg)
elif trigger_batch_interval <= 0:
msg = 'Unable to set trigger batch interval to a non-positive ' \
'value {0}'.format(trigger_batch_interval)
raise ValueError(msg)
self._trigger_batch_interval = trigger_batch_interval
@no_reduce.setter
def no_reduce(self, no_reduce):
@ -1103,17 +1205,18 @@ class SettingsFile(object):
self.create_batches_subelement()
self.create_inactive_subelement()
self.create_generations_per_batch_subelement()
self.create_keff_trigger_subelement()
def create_batches_subelement(self):
if not self._batches is None:
if self._eigenvalue_element is None:
self._eigenvalue_element = ET.SubElement(self._settings_file,
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_element, "batches")
element = ET.SubElement(self._eigenvalue_subelement, "batches")
element.text = '{0}'.format(self._batches)
@ -1121,11 +1224,11 @@ class SettingsFile(object):
if not self._generations_per_batch is None:
if self._eigenvalue_element is None:
self._eigenvalue_element = ET.SubElement(self._settings_file,
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_element,
element = ET.SubElement(self._eigenvalue_subelement,
"generations_per_batch")
element.text = '{0}'.format(self._generations_per_batch)
@ -1134,11 +1237,11 @@ class SettingsFile(object):
if not self._inactive is None:
if self._eigenvalue_element is None:
self._eigenvalue_element = ET.SubElement(self._settings_file,
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_element, "inactive")
element = ET.SubElement(self._eigenvalue_subelement, "inactive")
element.text = '{0}'.format(self._inactive)
@ -1146,14 +1249,29 @@ class SettingsFile(object):
if not self._particles is None:
if self._eigenvalue_element is None:
self._eigenvalue_element = ET.SubElement(self._settings_file,
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_element, "particles")
element = ET.SubElement(self._eigenvalue_subelement, "particles")
element.text = '{0}'.format(self._particles)
def create_keff_trigger_subelement(self):
if not self._keff_trigger is None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement, "keff_trigger")
for key in self._keff_trigger:
subelement = ET.SubElement(element, key)
subelement.text = str(self._keff_trigger[key]).lower()
def create_source_subelement(self):
self.create_source_space_subelement()
@ -1380,6 +1498,49 @@ class SettingsFile(object):
self._entropy_upper_right[2])
def create_trigger_subelement(self):
self.create_trigger_active_subelement()
self.create_trigger_max_batches_subelement()
self.create_trigger_batch_interval_subelement()
def create_trigger_active_subelement(self):
if not self._trigger_active is None:
if self._trigger_subelement is None:
self._trigger_subelement = ET.SubElement(self._settings_file,
"trigger")
element = ET.SubElement(self._trigger_subelement, "active")
element.text = '{0}'.format(str(self._trigger_active).lower())
def create_trigger_max_batches_subelement(self):
if not self._trigger_max_batches is None:
if self._trigger_subelement is None:
self._trigger_subelement = ET.SubElement(self._settings_file,
"trigger")
element = ET.SubElement(self._trigger_subelement, "max_batches")
element.text = '{0}'.format(self._trigger_max_batches)
def create_trigger_batch_interval_subelement(self):
if not self._trigger_batch_interval is None:
if self._trigger_subelement is None:
self._trigger_subelement = ET.SubElement(self._settings_file,
"trigger")
element = ET.SubElement(self._trigger_subelement, "batch_interval")
element.text = '{0}'.format(self._trigger_batch_interval)
def create_no_reduce_subelement(self):
if not self._no_reduce is None:
@ -1497,6 +1658,7 @@ class SettingsFile(object):
self.create_survival_biasing_subelement()
self.create_cutoff_subelement()
self.create_entropy_subelement()
self.create_trigger_subelement()
self.create_no_reduce_subelement()
self.create_threads_subelement()
self.create_verbosity_subelement()

View file

@ -581,6 +581,116 @@ class Mesh(object):
return element
class Trigger(object):
def __init__(self, trigger_type, threshold):
# Initialize Mesh class attributes
self.trigger_type = trigger_type
self.threshold = threshold
self._scores = []
def __deepcopy__(self, memo):
existing = memo.get(id(self))
# If this is first time we have tried to copy this object, create a copy
if existing is None:
clone = type(self).__new__(type(self))
clone._trigger_type = self._trigger_type
clone._threshold = self._threshold
clone._scores = []
for score in self._scores:
clone.add_score(score)
memo[id(self)] = clone
return clone
# If this object has been copied before, return the first copy made
else:
return existing
@property
def trigger_type(self):
return self._trigger_type
@property
def threshold(self):
return self._threshold
@property
def scores(self):
return self._scores
@trigger_type.setter
def trigger_type(self, trigger_type):
if not trigger_type in ['variance', 'std_dev', 'rel_err']:
msg = 'Unable to create a tally trigger with ' \
'type {0}'.format(trigger_type)
raise ValueError(msg)
self._trigger_type = trigger_type
@threshold.setter
def threshold(self, threshold):
if not is_float(threshold):
msg = 'Unable to set a tally trigger threshold with ' \
'threshold {0}'.format(threshold)
raise ValueError(msg)
self._threshold = threshold
def add_score(self, score):
if not is_string(score):
msg = 'Unable to add score {0} to tally trigger since ' \
'it is not a string'.format(score)
raise ValueError(msg)
# If the score is already in the Tally, don't add it again
if score in self._scores:
return
else:
self._scores.append(score)
def __repr__(self):
string = 'Trigger\n'
string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self._trigger_type)
string += '{0: <16}{1}{2}\n'.format('\tThreshold', '=\t', self._threshold)
string += '{0: <16}{1}{2}\n'.format('\tScores', '=\t', self._scores)
return string
def get_trigger_xml(self, element):
subelement = ET.SubElement(element, "trigger")
subelement.set("type", self._trigger_type)
subelement.set("threshold", str(self._threshold))
# Scores
if len(self._scores) != 0:
scores = ''
for score in self._scores:
scores += '{0} '.format(score)
scores.rstrip(' ')
subelement.set("scores", scores)
class Tally(object):
@ -593,6 +703,7 @@ class Tally(object):
self._nuclides = []
self._scores = []
self._estimator = None
self._triggers = []
self._num_score_bins = 0
self._num_realizations = 0
@ -633,6 +744,10 @@ class Tally(object):
for score in self._scores:
clone.add_score(score)
clone._triggers = []
for trigger in self._triggers:
trigger.add_trigger(trigger)
memo[id(self)] = clone
return clone
@ -762,6 +877,11 @@ class Tally(object):
return self._estimator
@property
def triggers(self):
return self._triggers
@property
def num_realizations(self):
return self._num_realizations
@ -789,6 +909,7 @@ class Tally(object):
@estimator.setter
def estimator(self, estimator):
if not estimator in ['analog', 'tracklength']:
msg = 'Unable to set the estimator for Tally ID={0} to {1} since ' \
'it is not a valid estimator type'.format(self._id, estimator)
@ -797,6 +918,16 @@ class Tally(object):
self._estimator = estimator
def add_trigger(self, trigger):
if not isinstance(trigger, Trigger):
msg = 'Unable to add a tally trigger for for Tally ID={0} to ' \
'{1} since it is not a valid estimator type'.format(trigger)
raise ValueError(msg)
self._triggers.append(trigger)
@id.setter
def id(self, tally_id):
@ -969,6 +1100,7 @@ class Tally(object):
return string
def get_tally_xml(self):
element = ET.Element("tally")
@ -1027,6 +1159,10 @@ class Tally(object):
subelement = ET.SubElement(element, "estimator")
subelement.text = self._estimator
# Optional Triggers
for trigger in self._triggers:
trigger.get_trigger_xml(element)
return element

View file

@ -0,0 +1,10 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="sphere" coeffs="0. 0. 0. 1.4" />
<surface id="2" type="sphere" coeffs="0. 0. 0. 5.35" boundary="vacuum" />
<cell id="1" material="void" surfaces="-1" /> <!-- cavity -->
<cell id="2" material="1" surfaces=" 1 -2" /> <!-- Pu Core -->
</geometry>

View file

@ -0,0 +1,13 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Plutonium -->
<material id="1">
<density value="18.4757" units="g/cm3" />
<nuclide name="Pu-239" ao="1.0" />
</material>
</materials>

View file

@ -0,0 +1,43 @@
#!/usr/bin/env python
import sys
import numpy as np
sys.path.insert(0, '../../src/utils')
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = StatePoint(sys.argv[1])
else:
sp = StatePoint('statepoint.19.binary')
sp.read_results()
# extract tally results and convert to vector
sum1 = sp._tallies[1]._sum
sum_sq1 = sp._tallies[1]._sum_sq
results1 = np.concatenate([sum1, sum_sq1])
sum2 = sp._tallies[2]._sum
sum_sq2 = sp._tallies[2]._sum_sq
results2 = np.concatenate([sum2, sum_sq2])
results = np.concatenate([results1.flatten(), results2.flatten()])
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write out tally results
outstr += 'tallies:\n'
for item in results:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

View file

@ -0,0 +1,27 @@
k-combined:
9.851940E-01 4.283951E-03
tallies:
1.972289E+01
4.489781E+00
4.356175E+00
1.523311E+01
1.972289E+01
4.489781E+00
4.356175E+00
1.523311E+01
2.779778E+01
1.440247E+00
1.355797E+00
1.658397E+01
2.779778E+01
1.440247E+00
1.355797E+00
1.658397E+01
1.972289E+01
4.489781E+00
4.356175E+00
1.523311E+01
2.779778E+01
1.440247E+00
1.355797E+00
1.658397E+01

View file

@ -0,0 +1,26 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<keff_trigger>
<type>std_dev</type>
<threshold>0.00445</threshold>
</keff_trigger>
</eigenvalue>
<trigger>
<active>true</active>
<max_batches>30</max_batches>
<batch_interval>1</batch_interval>
</trigger>
<source>
<space>
<type>point</type>
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>

View file

@ -0,0 +1,14 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<nuclides>all</nuclides>
<scores>total absorption fission scatter</scores>
</tally>
<tally id="2">
<nuclides>Pu-239</nuclides>
<scores>total absorption fission scatter</scores>
</tally>
</tallies>

View file

@ -0,0 +1,59 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from optparse import OptionParser
parser = OptionParser()
parser.add_option('--mpi_exec', dest='mpi_exec', default='')
parser.add_option('--mpi_np', dest='mpi_np', default='3')
parser.add_option('--exe', dest='exe')
(opts, args) = parser.parse_args()
cwd = os.getcwd()
def test_run():
if opts.mpi_exec != '':
proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE)
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0, 'OpenMC did not exit successfully.'
def test_created_statepoint():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.19.*'))
assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.'
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\
'Statepoint file is not a binary or hdf5 file.'
def test_results():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.19.*'))
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare, 'Results do not agree.'
def teardown():
output = glob.glob(os.path.join(cwd, 'statepoint.19.*'))
output.append(os.path.join(cwd, 'results_test.dat'))
for f in output:
if os.path.exists(f):
os.remove(f)
if __name__ == '__main__':
# test for openmc executable
if opts.exe is None:
raise Exception('Must specify OpenMC executable from command line with --exe.')
# run tests
try:
test_run()
test_created_statepoint()
test_results()
finally:
teardown()

View file

@ -0,0 +1,10 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="sphere" coeffs="0. 0. 0. 1.4" />
<surface id="2" type="sphere" coeffs="0. 0. 0. 5.35" boundary="vacuum" />
<cell id="1" material="void" surfaces="-1" /> <!-- cavity -->
<cell id="2" material="1" surfaces=" 1 -2" /> <!-- Pu Core -->
</geometry>

View file

@ -0,0 +1,13 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Plutonium -->
<material id="1">
<density value="18.4757" units="g/cm3" />
<nuclide name="Pu-239" ao="1.0" />
</material>
</materials>

View file

@ -0,0 +1,42 @@
#!/usr/bin/env python
import sys
import numpy as np
sys.path.insert(0, '../../src/utils')
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = StatePoint(sys.argv[1])
else:
sp = StatePoint('statepoint.19.binary')
sp.read_results()
# extract tally results and convert to vector
sum1 = sp._tallies[1]._sum
sum_sq1 = sp._tallies[1]._sum_sq
results1 = np.concatenate([sum1, sum_sq1])
sum2 = sp._tallies[2]._sum
sum_sq2 = sp._tallies[2]._sum_sq
results2 = np.concatenate([sum2, sum_sq2])
results = np.concatenate([results1.flatten(), results2.flatten()])
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write out tally results
outstr += 'tallies:\n'
for item in results:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

View file

@ -0,0 +1,27 @@
k-combined:
9.851940E-01 4.283951E-03
tallies:
1.972289E+01
4.489781E+00
4.356175E+00
1.523311E+01
1.972289E+01
4.489781E+00
4.356175E+00
1.523311E+01
2.779778E+01
1.440247E+00
1.355797E+00
1.658397E+01
2.779778E+01
1.440247E+00
1.355797E+00
1.658397E+01
1.972289E+01
4.489781E+00
4.356175E+00
1.523311E+01
2.779778E+01
1.440247E+00
1.355797E+00
1.658397E+01

View file

@ -0,0 +1,25 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<keff_trigger>
<type>std_dev</type>
<threshold>0.00445</threshold>
</keff_trigger>
</eigenvalue>
<trigger>
<active>true</active>
<max_batches>30</max_batches>
</trigger>
<source>
<space>
<type>point</type>
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>

View file

@ -0,0 +1,14 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<nuclides>all</nuclides>
<scores>total absorption fission scatter</scores>
</tally>
<tally id="2">
<nuclides>Pu-239</nuclides>
<scores>total absorption fission scatter</scores>
</tally>
</tallies>

View file

@ -0,0 +1,59 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from optparse import OptionParser
parser = OptionParser()
parser.add_option('--mpi_exec', dest='mpi_exec', default='')
parser.add_option('--mpi_np', dest='mpi_np', default='3')
parser.add_option('--exe', dest='exe')
(opts, args) = parser.parse_args()
cwd = os.getcwd()
def test_run():
if opts.mpi_exec != '':
proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE)
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0, 'OpenMC did not exit successfully.'
def test_created_statepoint():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.19.*'))
assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.'
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\
'Statepoint file is not a binary or hdf5 file.'
def test_results():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.19.*'))
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare, 'Results do not agree.'
def teardown():
output = glob.glob(os.path.join(cwd, 'statepoint.19.*'))
output.append(os.path.join(cwd, 'results_test.dat'))
for f in output:
if os.path.exists(f):
os.remove(f)
if __name__ == '__main__':
# test for openmc executable
if opts.exe is None:
raise Exception('Must specify OpenMC executable from command line with --exe.')
# run tests
try:
test_run()
test_created_statepoint()
test_results()
finally:
teardown()

View file

@ -0,0 +1,10 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="sphere" coeffs="0. 0. 0. 1.4" />
<surface id="2" type="sphere" coeffs="0. 0. 0. 5.35" boundary="vacuum" />
<cell id="1" material="void" surfaces="-1" /> <!-- cavity -->
<cell id="2" material="1" surfaces=" 1 -2" /> <!-- Pu Core -->
</geometry>

View file

@ -0,0 +1,13 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Plutonium -->
<material id="1">
<density value="18.4757" units="g/cm3" />
<nuclide name="Pu-239" ao="1.0" />
</material>
</materials>

View file

@ -0,0 +1,42 @@
#!/usr/bin/env python
import sys
import numpy as np
sys.path.insert(0, '../../src/utils')
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = StatePoint(sys.argv[1])
else:
sp = StatePoint('statepoint.10.binary')
sp.read_results()
# extract tally results and convert to vector
sum1 = sp._tallies[1]._sum
sum_sq1 = sp._tallies[1]._sum_sq
results1 = np.concatenate([sum1, sum_sq1])
sum2 = sp._tallies[2]._sum
sum_sq2 = sp._tallies[2]._sum_sq
results2 = np.concatenate([sum2, sum_sq2])
results = np.concatenate([results1.flatten(), results2.flatten()])
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write out tally results
outstr += 'tallies:\n'
for item in results:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

View file

@ -0,0 +1,27 @@
k-combined:
9.809303E-01 7.264435E-03
tallies:
7.031241E+00
1.599216E+00
1.550710E+00
5.432025E+00
7.031241E+00
1.599216E+00
1.550710E+00
5.432025E+00
9.892232E+00
5.115987E-01
4.810241E-01
5.904811E+00
9.892232E+00
5.115987E-01
4.810241E-01
5.904811E+00
7.031241E+00
1.599216E+00
1.550710E+00
5.432025E+00
9.892232E+00
5.115987E-01
4.810241E-01
5.904811E+00

View file

@ -0,0 +1,26 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<keff_trigger>
<type>std_dev</type>
<threshold>0.009</threshold>
</keff_trigger>
</eigenvalue>
<trigger>
<active>false</active>
<max_batches>15</max_batches>
<batch_interval>1</batch_interval>
</trigger>
<source>
<space>
<type>point</type>
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>

View file

@ -0,0 +1,14 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<nuclides>all</nuclides>
<scores>total absorption fission scatter</scores>
</tally>
<tally id="2">
<nuclides>Pu-239</nuclides>
<scores>total absorption fission scatter</scores>
</tally>
</tallies>

View file

@ -0,0 +1,59 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from optparse import OptionParser
parser = OptionParser()
parser.add_option('--mpi_exec', dest='mpi_exec', default='')
parser.add_option('--mpi_np', dest='mpi_np', default='3')
parser.add_option('--exe', dest='exe')
(opts, args) = parser.parse_args()
cwd = os.getcwd()
def test_run():
if opts.mpi_exec != '':
proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE)
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0, 'OpenMC did not exit successfully.'
def test_created_statepoint():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.'
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\
'Statepoint file is not a binary or hdf5 file.'
def test_results():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare, 'Results do not agree.'
def teardown():
output = glob.glob(os.path.join(cwd, 'statepoint.10.*'))
output.append(os.path.join(cwd, 'results_test.dat'))
for f in output:
if os.path.exists(f):
os.remove(f)
if __name__ == '__main__':
# test for openmc executable
if opts.exe is None:
raise Exception('Must specify OpenMC executable from command line with --exe.')
# run tests
try:
test_run()
test_created_statepoint()
test_results()
finally:
teardown()

View file

@ -0,0 +1,10 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="sphere" coeffs="0. 0. 0. 1.4" />
<surface id="2" type="sphere" coeffs="0. 0. 0. 5.35" boundary="vacuum" />
<cell id="1" material="void" surfaces="-1" /> <!-- cavity -->
<cell id="2" material="1" surfaces=" 1 -2" /> <!-- Pu Core -->
</geometry>

View file

@ -0,0 +1,13 @@
<?xml version="1.0"?>
<materials>
<default_xs>71c</default_xs>
<!-- Plutonium -->
<material id="1">
<density value="18.4757" units="g/cm3" />
<nuclide name="Pu-239" ao="1.0" />
</material>
</materials>

View file

@ -0,0 +1,43 @@
#!/usr/bin/env python
import sys
import numpy as np
sys.path.insert(0, '../../src/utils')
from openmc.statepoint import StatePoint
# read in statepoint file
if len(sys.argv) > 1:
sp = StatePoint(sys.argv[1])
else:
sp = StatePoint('statepoint.15.binary')
sp.read_results()
# extract tally results and convert to vector
sum1 = sp._tallies[1]._sum
sum_sq1 = sp._tallies[1]._sum_sq
results1 = np.concatenate([sum1, sum_sq1])
sum2 = sp._tallies[2]._sum
sum_sq2 = sp._tallies[2]._sum_sq
results2 = np.concatenate([sum2, sum_sq2])
results = np.concatenate([results1.flatten(), results2.flatten()])
# set up output string
outstr = ''
# write out k-combined
outstr += 'k-combined:\n'
outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1])
# write out tally results
outstr += 'tallies:\n'
for item in results:
outstr += "{0:12.6E}\n".format(item)
# write results to file
with open('results_test.dat','w') as fh:
fh.write(outstr)

View file

@ -0,0 +1,27 @@
k-combined:
9.824135E-01 6.844128E-03
tallies:
1.408810E+01
3.199600E+00
3.102648E+00
1.088850E+01
1.408810E+01
3.199600E+00
3.102648E+00
1.088850E+01
1.985836E+01
1.024050E+00
9.629173E-01
1.186391E+01
1.985836E+01
1.024050E+00
9.629173E-01
1.186391E+01
1.408810E+01
3.199600E+00
3.102648E+00
1.088850E+01
1.985836E+01
1.024050E+00
9.629173E-01
1.186391E+01

View file

@ -0,0 +1,26 @@
<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<keff_trigger>
<type>std_dev</type>
<threshold>0.02</threshold>
</keff_trigger>
</eigenvalue>
<trigger>
<active>true</active>
<max_batches>15</max_batches>
<batch_interval>1</batch_interval>
</trigger>
<source>
<space>
<type>point</type>
<parameters>0.0 0.0 0.0</parameters>
</space>
</source>
</settings>

View file

@ -0,0 +1,15 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<nuclides>all</nuclides>
<scores>total absorption fission scatter</scores>
</tally>
<tally id="2">
<nuclides>Pu-239</nuclides>
<scores>total absorption fission scatter</scores>
<trigger type = 'rel_err' threshold = '0.0095' scores = 'all'/>
</tally>
</tallies>

View file

@ -0,0 +1,59 @@
#!/usr/bin/env python
import os
from subprocess import Popen, STDOUT, PIPE, call
import filecmp
import glob
from optparse import OptionParser
parser = OptionParser()
parser.add_option('--mpi_exec', dest='mpi_exec', default='')
parser.add_option('--mpi_np', dest='mpi_np', default='3')
parser.add_option('--exe', dest='exe')
(opts, args) = parser.parse_args()
cwd = os.getcwd()
def test_run():
if opts.mpi_exec != '':
proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd],
stderr=STDOUT, stdout=PIPE)
else:
proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE)
print(proc.communicate()[0])
returncode = proc.returncode
assert returncode == 0, 'OpenMC did not exit successfully.'
def test_created_statepoint():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.15.*'))
assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.'
assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\
'Statepoint file is not a binary or hdf5 file.'
def test_results():
statepoint = glob.glob(os.path.join(cwd, 'statepoint.15.*'))
call(['python', 'results.py', statepoint[0]])
compare = filecmp.cmp('results_test.dat', 'results_true.dat')
if not compare:
os.rename('results_test.dat', 'results_error.dat')
assert compare, 'Results do not agree.'
def teardown():
output = glob.glob(os.path.join(cwd, 'statepoint.15.*'))
output.append(os.path.join(cwd, 'results_test.dat'))
for f in output:
if os.path.exists(f):
os.remove(f)
if __name__ == '__main__':
# test for openmc executable
if opts.exe is None:
raise Exception('Must specify OpenMC executable from command line with --exe.')
# run tests
try:
test_run()
test_created_statepoint()
test_results()
finally:
teardown()