diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index c407019f2..7bf772e2b 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -139,6 +139,8 @@ default. This element has the following attributes/sub-elements: *Default*: 1.0 +.. _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. + ```` Element ------------------------- @@ -214,8 +244,8 @@ problem. It has the following attributes/sub-elements: ```` Element -------------------------- -The ```` element indicates that a fixed source calculation should be -performed. It has the following attributes/sub-elements: +The ```` 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. ```` 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. ```` Element ------------------- -The ```` 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 ```` 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: + +```` Element +------------------------- + +OpenMC includes tally precision triggers which allow the user to define +uncertainty thresholds on :math:`k_{eff}` in the ```` 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 ````. +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 ```` has been reached. + +The ```` 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 + ```` 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. + + ```` Element ------------------------ @@ -1116,6 +1193,8 @@ post-collision energy, and an arbitrary structured mesh. The three valid elements in the tallies.xml file are ````, ````, and ````. +.. _tally: + ```` Element ------------------- @@ -1139,11 +1218,13 @@ The ```` 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 ```` 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 `` scatter-P2 ``. :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 ```` 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" + ```` 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 ```` and any number output plots can be defined with ```` 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. ```` 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 ```` 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 ```` element in ``settings.xml`` should be set to "true". ```` Element @@ -1557,8 +1674,8 @@ The ```` element controls what batch CMFD calculations should begin. The ```` 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. ```` Element --------------------------- -The ```` element is used to view the convergence of power iteration. -This option can be turned on with "true" and turned off with "false". +The ```` 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 diff --git a/examples/python/lattice/hexagonal/build-xml.py b/examples/python/lattice/hexagonal/build-xml.py index c837e0194..5cf8eed06 100644 --- a/examples/python/lattice/hexagonal/build-xml.py +++ b/examples/python/lattice/hexagonal/build-xml.py @@ -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() diff --git a/examples/python/lattice/simple/build-xml.py b/examples/python/lattice/simple/build-xml.py index 98ee1fa66..675c7e08b 100644 --- a/examples/python/lattice/simple/build-xml.py +++ b/examples/python/lattice/simple/build-xml.py @@ -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() diff --git a/src/constants.F90 b/src/constants.F90 index e1c2fe3dc..0c9d15f6a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 :: & diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 860e2271f..3b7f149a8 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -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 diff --git a/src/fixed_source.F90 b/src/fixed_source.F90 index 4e7d9642f..78ae3d5a7 100644 --- a/src/fixed_source.F90 +++ b/src/fixed_source.F90 @@ -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() diff --git a/src/global.F90 b/src/global.F90 index 0f2cf0ed9..0ea4a1956 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8f5676647..32adfdc02 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 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 diff --git a/src/output_interface.F90 b/src/output_interface.F90 index d8134aa58..816ff2575 100644 --- a/src/output_interface.F90 +++ b/src/output_interface.F90 @@ -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 diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 092bf9983..a95b1f297 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -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{ diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index 10770f80c..8dd1a6feb 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -43,6 +43,28 @@ + + + + + + + + + + + + + + + + + + + + + + @@ -567,6 +589,38 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 6911dab2f..9cdd7a281 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -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 }? diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index a3668925a..e893b2fcf 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -204,6 +204,50 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + 20 + + + + + + + + + 20 + + + + + + + + + diff --git a/src/state_point.F90 b/src/state_point.F90 index 52492fc6c..b209d1707 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -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 diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 5c7efe84e..43ec868c5 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -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 diff --git a/src/trigger.F90 b/src/trigger.F90 new file mode 100644 index 000000000..d5efe96d5 --- /dev/null +++ b/src/trigger.F90 @@ -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 diff --git a/src/trigger_header.F90 b/src/trigger_header.F90 new file mode 100644 index 000000000..8289ca615 --- /dev/null +++ b/src/trigger_header.F90 @@ -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 diff --git a/src/utils/openmc/settings.py b/src/utils/openmc/settings.py index ce054487b..632565946 100644 --- a/src/utils/openmc/settings.py +++ b/src/utils/openmc/settings.py @@ -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() diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py index 1dc8e536e..14ca91f6f 100644 --- a/src/utils/openmc/tallies.py +++ b/src/utils/openmc/tallies.py @@ -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 diff --git a/tests/test_trigger_batch_interval/geometry.xml b/tests/test_trigger_batch_interval/geometry.xml new file mode 100644 index 000000000..65954a785 --- /dev/null +++ b/tests/test_trigger_batch_interval/geometry.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/tests/test_trigger_batch_interval/materials.xml b/tests/test_trigger_batch_interval/materials.xml new file mode 100644 index 000000000..2761be30c --- /dev/null +++ b/tests/test_trigger_batch_interval/materials.xml @@ -0,0 +1,13 @@ + + + + 71c + + + + + + + + + diff --git a/tests/test_trigger_batch_interval/results.py b/tests/test_trigger_batch_interval/results.py new file mode 100644 index 000000000..922fd4f74 --- /dev/null +++ b/tests/test_trigger_batch_interval/results.py @@ -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) diff --git a/tests/test_trigger_batch_interval/results_true.dat b/tests/test_trigger_batch_interval/results_true.dat new file mode 100644 index 000000000..95bbef6f6 --- /dev/null +++ b/tests/test_trigger_batch_interval/results_true.dat @@ -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 diff --git a/tests/test_trigger_batch_interval/settings.xml b/tests/test_trigger_batch_interval/settings.xml new file mode 100644 index 000000000..d7afd9231 --- /dev/null +++ b/tests/test_trigger_batch_interval/settings.xml @@ -0,0 +1,26 @@ + + + + 10 + 5 + 1000 + + std_dev + 0.00445 + + + + + true + 30 + 1 + + + + + point + 0.0 0.0 0.0 + + + + diff --git a/tests/test_trigger_batch_interval/tallies.xml b/tests/test_trigger_batch_interval/tallies.xml new file mode 100644 index 000000000..cf20668c8 --- /dev/null +++ b/tests/test_trigger_batch_interval/tallies.xml @@ -0,0 +1,14 @@ + + + + + all + total absorption fission scatter + + + + Pu-239 + total absorption fission scatter + + + diff --git a/tests/test_trigger_batch_interval/test_trigger_batch_interval.py b/tests/test_trigger_batch_interval/test_trigger_batch_interval.py new file mode 100644 index 000000000..b2c3a5ae7 --- /dev/null +++ b/tests/test_trigger_batch_interval/test_trigger_batch_interval.py @@ -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() diff --git a/tests/test_trigger_no_batch_interval/geometry.xml b/tests/test_trigger_no_batch_interval/geometry.xml new file mode 100644 index 000000000..65954a785 --- /dev/null +++ b/tests/test_trigger_no_batch_interval/geometry.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/tests/test_trigger_no_batch_interval/materials.xml b/tests/test_trigger_no_batch_interval/materials.xml new file mode 100644 index 000000000..2761be30c --- /dev/null +++ b/tests/test_trigger_no_batch_interval/materials.xml @@ -0,0 +1,13 @@ + + + + 71c + + + + + + + + + diff --git a/tests/test_trigger_no_batch_interval/results.py b/tests/test_trigger_no_batch_interval/results.py new file mode 100644 index 000000000..3665b48a9 --- /dev/null +++ b/tests/test_trigger_no_batch_interval/results.py @@ -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) diff --git a/tests/test_trigger_no_batch_interval/results_true.dat b/tests/test_trigger_no_batch_interval/results_true.dat new file mode 100644 index 000000000..95bbef6f6 --- /dev/null +++ b/tests/test_trigger_no_batch_interval/results_true.dat @@ -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 diff --git a/tests/test_trigger_no_batch_interval/settings.xml b/tests/test_trigger_no_batch_interval/settings.xml new file mode 100644 index 000000000..5ceff1082 --- /dev/null +++ b/tests/test_trigger_no_batch_interval/settings.xml @@ -0,0 +1,25 @@ + + + + 10 + 5 + 1000 + + std_dev + 0.00445 + + + + + true + 30 + + + + + point + 0.0 0.0 0.0 + + + + diff --git a/tests/test_trigger_no_batch_interval/tallies.xml b/tests/test_trigger_no_batch_interval/tallies.xml new file mode 100644 index 000000000..cf20668c8 --- /dev/null +++ b/tests/test_trigger_no_batch_interval/tallies.xml @@ -0,0 +1,14 @@ + + + + + all + total absorption fission scatter + + + + Pu-239 + total absorption fission scatter + + + diff --git a/tests/test_trigger_no_batch_interval/test_trigger_no_batch_interval.py b/tests/test_trigger_no_batch_interval/test_trigger_no_batch_interval.py new file mode 100644 index 000000000..b2c3a5ae7 --- /dev/null +++ b/tests/test_trigger_no_batch_interval/test_trigger_no_batch_interval.py @@ -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() diff --git a/tests/test_trigger_no_status/geometry.xml b/tests/test_trigger_no_status/geometry.xml new file mode 100644 index 000000000..65954a785 --- /dev/null +++ b/tests/test_trigger_no_status/geometry.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/tests/test_trigger_no_status/materials.xml b/tests/test_trigger_no_status/materials.xml new file mode 100644 index 000000000..2761be30c --- /dev/null +++ b/tests/test_trigger_no_status/materials.xml @@ -0,0 +1,13 @@ + + + + 71c + + + + + + + + + diff --git a/tests/test_trigger_no_status/results.py b/tests/test_trigger_no_status/results.py new file mode 100644 index 000000000..7754ec7e5 --- /dev/null +++ b/tests/test_trigger_no_status/results.py @@ -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) diff --git a/tests/test_trigger_no_status/results_true.dat b/tests/test_trigger_no_status/results_true.dat new file mode 100644 index 000000000..4ab52de53 --- /dev/null +++ b/tests/test_trigger_no_status/results_true.dat @@ -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 diff --git a/tests/test_trigger_no_status/settings.xml b/tests/test_trigger_no_status/settings.xml new file mode 100644 index 000000000..3d215a432 --- /dev/null +++ b/tests/test_trigger_no_status/settings.xml @@ -0,0 +1,26 @@ + + + + 10 + 5 + 1000 + + std_dev + 0.009 + + + + false + 15 + 1 + + + + + + point + 0.0 0.0 0.0 + + + + diff --git a/tests/test_trigger_no_status/tallies.xml b/tests/test_trigger_no_status/tallies.xml new file mode 100644 index 000000000..cf20668c8 --- /dev/null +++ b/tests/test_trigger_no_status/tallies.xml @@ -0,0 +1,14 @@ + + + + + all + total absorption fission scatter + + + + Pu-239 + total absorption fission scatter + + + diff --git a/tests/test_trigger_no_status/test_trigger_no_status.py b/tests/test_trigger_no_status/test_trigger_no_status.py new file mode 100644 index 000000000..6fdbf8745 --- /dev/null +++ b/tests/test_trigger_no_status/test_trigger_no_status.py @@ -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() diff --git a/tests/test_trigger_tallies/geometry.xml b/tests/test_trigger_tallies/geometry.xml new file mode 100644 index 000000000..65954a785 --- /dev/null +++ b/tests/test_trigger_tallies/geometry.xml @@ -0,0 +1,10 @@ + + + + + + + + + + diff --git a/tests/test_trigger_tallies/materials.xml b/tests/test_trigger_tallies/materials.xml new file mode 100644 index 000000000..2761be30c --- /dev/null +++ b/tests/test_trigger_tallies/materials.xml @@ -0,0 +1,13 @@ + + + + 71c + + + + + + + + + diff --git a/tests/test_trigger_tallies/results.py b/tests/test_trigger_tallies/results.py new file mode 100644 index 000000000..88aa4146f --- /dev/null +++ b/tests/test_trigger_tallies/results.py @@ -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) + diff --git a/tests/test_trigger_tallies/results_true.dat b/tests/test_trigger_tallies/results_true.dat new file mode 100644 index 000000000..86c3cfa03 --- /dev/null +++ b/tests/test_trigger_tallies/results_true.dat @@ -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 diff --git a/tests/test_trigger_tallies/settings.xml b/tests/test_trigger_tallies/settings.xml new file mode 100644 index 000000000..c6d986ce8 --- /dev/null +++ b/tests/test_trigger_tallies/settings.xml @@ -0,0 +1,26 @@ + + + + 10 + 5 + 1000 + + std_dev + 0.02 + + + + + true + 15 + 1 + + + + + point + 0.0 0.0 0.0 + + + + diff --git a/tests/test_trigger_tallies/tallies.xml b/tests/test_trigger_tallies/tallies.xml new file mode 100644 index 000000000..2dbc836ee --- /dev/null +++ b/tests/test_trigger_tallies/tallies.xml @@ -0,0 +1,15 @@ + + + + + all + total absorption fission scatter + + + + Pu-239 + total absorption fission scatter + + + + diff --git a/tests/test_trigger_tallies/test_trigger_tallies.py b/tests/test_trigger_tallies/test_trigger_tallies.py new file mode 100644 index 000000000..cb973cb4d --- /dev/null +++ b/tests/test_trigger_tallies/test_trigger_tallies.py @@ -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()