From 87c06cd334ccc61f96429151d479e6bc9ddefbab Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 28 Oct 2017 22:04:57 -0400 Subject: [PATCH 01/26] Update list of publications. Add note about citing OpenMC. --- docs/source/index.rst | 26 ++++++----- docs/source/publications.rst | 86 +++++++++++++++++++++++++----------- readme.rst | 12 +++++ 3 files changed, 88 insertions(+), 36 deletions(-) diff --git a/docs/source/index.rst b/docs/source/index.rst index 00b09db9e7..4071b9794c 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -10,21 +10,25 @@ interaction data is based on a native HDF5 format that can be generated from ACE files used by the MCNP and Serpent Monte Carlo codes. OpenMC was originally developed by members of the `Computational Reactor Physics -Group`_ at the `Massachusetts Institute of Technology`_ starting -in 2011. Various universities, laboratories, and other organizations now -contribute to the development of OpenMC. For more information on OpenMC, feel -free to send a message to the User's Group `mailing list`_. +Group `_ at the `Massachusetts Institute of Technology +`_ starting in 2011. Various universities, laboratories, and +other organizations now contribute to the development of OpenMC. For more +information on OpenMC, feel free to send a message to the User's Group `mailing +list `_. -.. _Computational Reactor Physics Group: http://crpg.mit.edu -.. _Massachusetts Institute of Technology: http://web.mit.edu -.. _mailing list: https://groups.google.com/forum/?fromgroups=#!forum/openmc-users -.. _Read the Docs: http://openmc.readthedocs.io/en/latest/ +.. admonition:: Recommended publication for citing + :class: tip + + Paul K. Romano, Nicholas E. Horelik, Bryan R. Herman, Adam G. Nelson, Benoit + Forget, and Kord Smith, "`OpenMC: A State-of-the-Art Monte Carlo Code for + Research and Development `_," + *Ann. Nucl. Energy*, **82**, 90--97 (2015). .. only:: html - -------- - Contents - -------- + -------- + Contents + -------- .. toctree:: :maxdepth: 1 diff --git a/docs/source/publications.rst b/docs/source/publications.rst index f50ec8a68a..f47ac704c1 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -10,7 +10,7 @@ Overviews - Paul K. Romano, Nicholas E. Horelik, Bryan R. Herman, Adam G. Nelson, Benoit Forget, and Kord Smith, "`OpenMC: A State-of-the-Art Monte Carlo Code for - Research and Development `_," + Research and Development `_," *Ann. Nucl. Energy*, **82**, 90--97 (2015). - Paul K. Romano, Bryan R. Herman, Nicholas E. Horelik, Benoit Forget, Kord @@ -19,16 +19,19 @@ Overviews Nuclear Science and Engineering*, Sun Valley, Idaho, May 5--9 (2013). - Paul K. Romano and Benoit Forget, "`The OpenMC Monte Carlo Particle Transport - Code `_," + Code `_," *Ann. Nucl. Energy*, **51**, 274--281 (2013). ------------ Benchmarking ------------ +- Travis J. Labossiere-Hickman and Benoit Forget, "Selected VERA Core Physics + Benchmarks in OpenMC," *Trans. Am. Nucl. Soc.*, **117**, 1520-1523 (2017). + - Khurrum S. Chaudri and Sikander M. Mirza, "`Burnup dependent Monte Carlo neutron physics calculations of IAEA MTR benchmark - `_," *Prog. Nucl. Energy*, + `_," *Prog. Nucl. Energy*, **81**, 43-52 (2015). - Daniel J. Kelly, Brian N. Aviles, Paul K. Romano, Bryan R. Herman, @@ -54,10 +57,15 @@ Benchmarking Coupling and Multi-physics -------------------------- +- Tianliang Hu, Liangzhu Cao, Hongchun Wu, Xianan Du, and Mingtao He, "`Coupled + neutrons and thermal-hydraulics simulation of molten salt reactors based on + OpenMC/TANSY `_," + *Ann. Nucl. Energy*, **109**, 260-276 (2017). + - Matthew Ellis, Derek Gaston, Benoit Forget, and Kord Smith, "`Preliminary Coupling of the Monte Carlo Code OpenMC and the Multiphysics Object-Oriented Simulation Environment for Analyzing Doppler Feedback in Monte Carlo - Simulations `_," *Nucl. Sci. Eng.*, + Simulations `_," *Nucl. Sci. Eng.*, **185**, 184-193 (2017). - Matthew Ellis, Benoit Forget, Kord Smith, and Derek Gaston, "Continuous @@ -79,7 +87,7 @@ Coupling and Multi-physics - Bryan R. Herman, Benoit Forget, and Kord Smith, "`Progress toward Monte Carlo-thermal hydraulic coupling using low-order nonlinear diffusion - acceleration methods `_," + acceleration methods `_," *Ann. Nucl. Energy*, **84**, 63-72 (2015). - Bryan R. Herman, Benoit Forget, and Kord Smith, "Utilizing CMFD in OpenMC to @@ -106,6 +114,15 @@ Geometry and Visualization Miscellaneous ------------- +- Adam G. Nelson, Samuel Shaner, William Boyd, and Paul K. Romano, + "Incorporation of a Multigroup Transport Capability in the OpenMC Monte Carlo + Particle Transport Code," *Trans. Am. Nucl. Soc.*, **117**, 679-681 (2017). + +- Youqi Zheng, Yunlong Xiao, and Hongchun Wu, "`Application of the virtual + density theory in fast reactor analysis based on the neutron transport + calculation `_," + *Nucl. Eng. Des.*, **320**, 200-206 (2017). + - Amanda L. Lund, Paul K. Romano, and Andrew R. Siegel, "Accelerating Source Convergence in Monte Carlo Criticality Calculations Using a Particle Ramp-Up Technique," *Proc. Int. Conf. Mathematics & Computational Methods Applied to @@ -126,7 +143,7 @@ Miscellaneous - Yunzhao Li, Qingming He, Liangzhi Cao, Hongchun Wu, and Tiejun Zu, "`Resonance Elastic Scattering and Interference Effects Treatments in Subgroup Method - `_," *Nucl. Eng. Tech.*, **48**, + `_," *Nucl. Eng. Tech.*, **48**, 339-350 (2016). - William Boyd, Sterling Harper, and Paul K. Romano, "Equipping OpenMC for the @@ -135,7 +152,7 @@ Miscellaneous - Michal Kostal, Vojtech Rypar, Jan Milcak, Vlastimil Juricek, Evzen Losa, Benoit Forget, and Sterling Harper, "`Study of graphite reactivity worth on well-defined cores assembled on LR-0 reactor - `_," *Ann. Nucl. Energy*, + `_," *Ann. Nucl. Energy*, **87**, 601-611 (2016). - Qicang Shen, William Boyd, Benoit Forget, and Kord Smith, "Tally precision @@ -154,9 +171,19 @@ Miscellaneous Multi-group Cross Section Generation ------------------------------------ +- Gang Yang, Tongkyu Park, and Won Sik Yang, "Effects of Fuel Salt Velocity + Field on Neutronics Performances in Molten Salt Reactors with Open Flow + Channels," *Trans. Am. Nucl. Soc.*, **117**, 1339-1342 (2017). + +- William Boyd, Nathan Gibson, Benoit Forget, and Kord Smith, "`An analysis of + condensation errors in multi-group cross section generation for fine-mesh + neutron transport calculations + `_," *Ann. Nucl. Energy*, + **112**, 267-276 (2018). + - Hong Shuang, Yang Yongwei, Zhang Lu, and Gao Yucui, "`Fabrication and validation of multigroup cross section library based on the OpenMC code - `_," + `_," *Nucl. Techniques* **40** (4), 040504 (2017). (in Mandarin) - Nicholas E. Stauff, Changho Lee, Paul K. Romano, and Taek K. Kim, @@ -207,7 +234,7 @@ Doppler Broadening - Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "`Windowed multipole for cross section Doppler broadening - `_," *J. Comput. Phys.*, **307**, + `_," *J. Comput. Phys.*, **307**, 715-727 (2016). - Jonathan A. Walsh, Benoit Forget, Kord S. Smith, and Forrest B. Brown, @@ -217,12 +244,12 @@ Doppler Broadening - Colin Josey, Benoit Forget, and Kord Smith, "`Windowed multipole sensitivity to target accuracy of the optimization procedure - `_," + `_," *J. Nucl. Sci. Technol.*, **52**, 987-992 (2015). - Paul K. Romano and Timothy H. Trumbull, "`Comparison of algorithms for Doppler broadening pointwise tabulated cross sections - `_," *Ann. Nucl. Energy*, + `_," *Ann. Nucl. Energy*, **75**, 358--364 (2015). - Tuomas Viitanen, Jaakko Leppanen, and Benoit Forget, "Target motion sampling @@ -231,13 +258,17 @@ Doppler Broadening - Benoit Forget, Sheng Xu, and Kord Smith, "`Direct Doppler broadening in Monte Carlo simulations using the multipole representation - `_," *Ann. Nucl. Energy*, + `_," *Ann. Nucl. Energy*, **64**, 78--85 (2014). ------------ Nuclear Data ------------ +- Jonathan A. Walsh, "Comparison of Unresolved Resonance Region Cross Section + Formalisms in Transport Simulations," *Trans. Am. Nucl. Soc.*, **117**, + 749-752 (2017). + - Jonathan A. Walsh, Benoit Forget, Kord S. Smith, and Forrest B. Brown, "`Uncertainty in Fast Reactor-Relevant Critical Benchmark Simulations Due to Unresolved Resonance Structure @@ -255,13 +286,13 @@ Nuclear Data - Jonathan A. Walsh, Benoit Froget, Kord S. Smith, and Forrest B. Brown, "`Neutron Cross Section Processing Methods for Improved Integral Benchmarking of Unresolved Resonance Region Evaluations - `_," *Eur. Phys. J. Web Conf.* + `_," *Eur. Phys. J. Web Conf.* **111**, 06001 (2016). - Jonathan A. Walsh, Paul K. Romano, Benoit Forget, and Kord S. Smith, "`Optimizations of the energy grid search algorithm in continuous-energy Monte Carlo particle transport codes - `_", *Comput. Phys. Commun.*, + `_", *Comput. Phys. Commun.*, **196**, 134-142 (2015). - Jonathan A. Walsh, Benoit Forget, Kord S. Smith, Brian C. Kiedrowski, and @@ -280,7 +311,7 @@ Nuclear Data - Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "`Accelerated sampling of the free gas resonance elastic scattering kernel - `_," *Ann. Nucl. Energy*, + `_," *Ann. Nucl. Energy*, **69**, 116--124 (2014). ----------- @@ -325,45 +356,45 @@ Parallelism - Nicholas Horelik, Andrew Siegel, Benoit Forget, and Kord Smith, "`Monte Carlo domain decomposition for robust nuclear reactor analysis - `_," *Parallel Comput.*, + `_," *Parallel Comput.*, **40**, 646--660 (2014). - Andrew Siegel, Kord Smith, Kyle Felker, Paul Romano, Benoit Forget, and Peter Beckman, "`Improved cache performance in Monte Carlo transport calculations - using energy banding `_," + using energy banding `_," *Comput. Phys. Commun.*, **185** (4), 1195--1199 (2014). - Paul K. Romano, Benoit Forget, Kord Smith, and Andrew Siegel, "`On the use of tally servers in Monte Carlo simulations of light-water reactors - `_," *Proc. Joint International + `_," *Proc. Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo*, Paris, France, Oct. 27--31 (2013). - Kyle G. Felker, Andrew R. Siegel, Kord S. Smith, Paul K. Romano, and Benoit Forget, "`The energy band memory server algorithm for parallel Monte Carlo - calculations `_," *Proc. Joint + calculations `_," *Proc. Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo*, Paris, France, Oct. 27--31 (2013). - John R. Tramm and Andrew R. Siegel, "`Memory Bottlenecks and Memory Contention in Multi-Core Monte Carlo Transport Codes - `_," *Proc. Joint International + `_," *Proc. Joint International Conference on Supercomputing in Nuclear Applications and Monte Carlo*, Paris, France, Oct. 27--31 (2013). - Andrew R. Siegel, Kord Smith, Paul K. Romano, Benoit Forget, and Kyle Felker, "`Multi-core performance studies of a Monte Carlo neutron transport code - `_," *Int. J. High + `_," *Int. J. High Perform. Comput. Appl.*, **28** (1), 87--96 (2014). - Paul K. Romano, Andrew R. Siegel, Benoit Forget, and Kord Smith, "`Data decomposition of Monte Carlo particle transport simulations via tally servers - `_," *J. Comput. Phys.*, **252**, + `_," *J. Comput. Phys.*, **252**, 20--36 (2013). - Andrew R. Siegel, Kord Smith, Paul K. Romano, Benoit Forget, and Kyle Felker, "`The effect of load imbalances on the performance of Monte Carlo codes in LWR - analysis `_," *J. Comput. Phys.*, + analysis `_," *J. Comput. Phys.*, **235**, 901--911 (2013). @@ -372,13 +403,18 @@ Parallelism 519--522 (2012). - Paul K. Romano and Benoit Forget, "`Parallel Fission Bank Algorithms in Monte - Carlo Criticality Calculations `_," + Carlo Criticality Calculations `_," *Nucl. Sci. Eng.*, **170**, 125--135 (2012). --------- Depletion --------- +- Colin Josey, Benoit Forget, and Kord Smith, "`High order methods for the + integration of the Bateman equations and other problems of the form of y' = + F(y,t)y `_," *J. Comput. Phys.*, + **350**, 296-313 (2017). + - Matthew S. Ellis, Colin Josey, Benoit Forget, and Kord Smith, "`Spatially Continuous Depletion Algorithm for Monte Carlo Simulations `_," *Trans. Am. Nucl. Soc.*, **115**, @@ -386,7 +422,7 @@ Depletion - Anas Gul, K. S. Chaudri, R. Khan, and M. Azeen, "`Development and verification of LOOP: A Linkage of ORIGEN2.2 and OpenMC - `_," *Ann. Nucl. Energy*, + `_," *Ann. Nucl. Energy*, **99**, 321--327 (2017). - Kai Huang, Hongchun Wu, Yunzhao Li, and Liangzhi Cao, "Generalized depletion diff --git a/readme.rst b/readme.rst index 75e4a1438f..b7366fe3e1 100644 --- a/readme.rst +++ b/readme.rst @@ -20,6 +20,18 @@ Installation Detailed `installation instructions`_ can be found in the User's Guide. +------ +Citing +------ + +If you use OpenMC in your research, please consider giving proper attribution by +citing the following publication: + +- Paul K. Romano, Nicholas E. Horelik, Bryan R. Herman, Adam G. Nelson, Benoit + Forget, and Kord Smith, "`OpenMC: A State-of-the-Art Monte Carlo Code for + Research and Development `_," + *Ann. Nucl. Energy*, **82**, 90--97 (2015). + --------------- Troubleshooting --------------- From 3af20e8faf4fc65ff5abae78c2c66ad2ef4bffa8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 29 Oct 2017 05:47:48 -0400 Subject: [PATCH 02/26] Fix openmc.capi.Tally.scores getter --- openmc/capi/tally.py | 31 ++++++++++++++++++++++++++++++- src/api.F90 | 1 + src/constants.F90 | 3 ++- src/tallies/tally_header.F90 | 26 ++++++++++++++++++++++++++ 4 files changed, 59 insertions(+), 2 deletions(-) diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 799cb42ee8..7037484aea 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -4,6 +4,7 @@ from weakref import WeakValueDictionary from numpy.ctypeslib import as_array +from openmc.data.reaction import REACTION_NAME from . import _dll, Nuclide from .core import _FortranObjectWithID from .error import _error_handler, AllocationError, InvalidIDError @@ -30,6 +31,10 @@ _dll.openmc_tally_get_nuclides.argtypes = [ c_int32, POINTER(POINTER(c_int)), POINTER(c_int)] _dll.openmc_tally_get_nuclides.restype = c_int _dll.openmc_tally_get_nuclides.errcheck = _error_handler +_dll.openmc_tally_get_scores.argtypes = [ + c_int32, POINTER(POINTER(c_int)), POINTER(c_int)] +_dll.openmc_tally_get_scores.restype = c_int +_dll.openmc_tally_get_scores.errcheck = _error_handler _dll.openmc_tally_results.argtypes = [ c_int32, POINTER(POINTER(c_double)), POINTER(c_int*3)] _dll.openmc_tally_results.restype = c_int @@ -51,6 +56,15 @@ _dll.openmc_tally_set_type.restype = c_int _dll.openmc_tally_set_type.errcheck = _error_handler +_SCORES = { + -1: 'flux', -2: 'total', -3: 'scatter', -4: 'nu-scatter', + -9: 'absorption', -10: 'fission', -11: 'nu-fission', -12: 'kappa-fission', + -13: 'current', -18: 'events', -19: 'delayed-nu-fission', + -20: 'prompt-nu-fission', -21: 'inverse-velocity', -22: 'fission-q-prompt', + -23: 'fission-q-recoverable', -24: 'decay-rate' +} + + class Tally(_FortranObjectWithID): """Tally stored internally. @@ -161,7 +175,22 @@ class Tally(_FortranObjectWithID): @property def scores(self): - pass + scores_as_int = POINTER(c_int)() + n = c_int() + try: + _dll.openmc_tally_get_scores(self._index, scores_as_int, n) + except AllocationError: + return [] + else: + scores = [] + for i in range(n.value): + if scores_as_int[i] in _SCORES: + scores.append(_SCORES[scores_as_int[i]]) + elif scores_as_int[i] in REACTION_NAME: + scores.append(REACTION_NAME[scores_as_int[i]]) + else: + scores.append(str(scores_as_int[i])) + return scores @scores.setter def scores(self, scores): diff --git a/src/api.F90 b/src/api.F90 index f5b27e7181..933faaa09e 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -74,6 +74,7 @@ module openmc_api public :: openmc_tally_get_id public :: openmc_tally_get_filters public :: openmc_tally_get_nuclides + public :: openmc_tally_get_scores public :: openmc_tally_results public :: openmc_tally_set_filters public :: openmc_tally_set_id diff --git a/src/constants.F90 b/src/constants.F90 index 73aca25d29..f6a7e8da92 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -304,7 +304,8 @@ module constants EVENT_SCATTER = 1, & EVENT_ABSORB = 2 - ! Tally score type + ! Tally score type -- if you change these, make sure you also update the + ! _SCORES dictionary in openmc/capi/tally.py integer, parameter :: N_SCORE_TYPES = 24 integer, parameter :: & SCORE_FLUX = -1, & ! flux diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index df9d935e94..e137a13274 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -25,6 +25,7 @@ module tally_header public :: openmc_tally_get_id public :: openmc_tally_get_filters public :: openmc_tally_get_nuclides + public :: openmc_tally_get_scores public :: openmc_tally_results public :: openmc_tally_set_filters public :: openmc_tally_set_id @@ -535,6 +536,31 @@ contains end function openmc_tally_get_nuclides + function openmc_tally_get_scores(index, scores, n) result(err) bind(C) + ! Return the list of nuclides assigned to a tally + integer(C_INT32_T), value :: index + type(C_PTR), intent(out) :: scores + integer(C_INT), intent(out) :: n + integer(C_INT) :: err + + if (index >= 1 .and. index <= size(tallies)) then + associate (t => tallies(index) % obj) + if (allocated(t % score_bins)) then + scores = C_LOC(t % score_bins(1)) + n = size(t % score_bins) + err = 0 + else + err = E_ALLOCATE + call set_errmsg("Tally scores have not been allocated yet.") + end if + end associate + else + err = E_OUT_OF_BOUNDS + call set_errmsg('Index in tallies array is out of bounds.') + end if + end function openmc_tally_get_scores + + function openmc_tally_results(index, ptr, shape_) result(err) bind(C) ! Returns a pointer to a tally results array along with its shape. This ! allows a user to obtain in-memory tally results from Python directly. From 9f66edb7fb01ea50cf11e2056c45a77e488c1257 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 3 Nov 2017 07:25:23 -0500 Subject: [PATCH 03/26] Move batch logic into openmc_next_batch --- src/simulation.F90 | 129 +++++++++++++++++++++++---------------------- 1 file changed, 65 insertions(+), 64 deletions(-) diff --git a/src/simulation.F90 b/src/simulation.F90 index a240f55b06..5fd6f60f70 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -43,6 +43,8 @@ module simulation implicit none private public :: openmc_run + public :: openmc_simulation_init + public :: openmc_simulation_finalize contains @@ -54,75 +56,71 @@ contains subroutine openmc_run() bind(C) - type(Particle) :: p - integer(8) :: i_work + call openmc_simulation_init() - call initialize_simulation() - - ! Turn on inactive timer - call time_inactive % start() - - ! ========================================================================== - ! LOOP OVER BATCHES BATCH_LOOP: do current_batch = 1, n_max_batches - - call initialize_batch() - - ! Handle restart runs - if (restart_run .and. current_batch <= restart_batch) then - call replay_batch_history() - cycle BATCH_LOOP - end if - - ! ======================================================================= - ! LOOP OVER GENERATIONS - GENERATION_LOOP: do current_gen = 1, gen_per_batch - - call initialize_generation() - - ! Start timer for transport - call time_transport % start() - - ! ==================================================================== - ! LOOP OVER PARTICLES -!$omp parallel do schedule(static) firstprivate(p) copyin(tally_derivs) - PARTICLE_LOOP: do i_work = 1, work - current_work = i_work - - ! grab source particle from bank - call initialize_history(p, current_work) - - ! transport particle - call transport(p) - - end do PARTICLE_LOOP -!$omp end parallel do - - ! Accumulate time for transport - call time_transport % stop() - - call finalize_generation() - - end do GENERATION_LOOP - - call finalize_batch() - + call openmc_next_batch() if (satisfy_triggers) exit BATCH_LOOP - end do BATCH_LOOP call time_active % stop() - ! ========================================================================== - ! END OF RUN WRAPUP - - call finalize_simulation() - - ! Clear particle - call p % clear() + call openmc_simulation_finalize() end subroutine openmc_run +!=============================================================================== +! OPENMC_NEXT_BATCH +!=============================================================================== + + subroutine openmc_next_batch() bind(C) + + type(Particle) :: p + integer(8) :: i_work + + call initialize_batch() + + ! Handle restart runs + if (restart_run .and. current_batch <= restart_batch) then + call replay_batch_history() + return + end if + + ! ======================================================================= + ! LOOP OVER GENERATIONS + GENERATION_LOOP: do current_gen = 1, gen_per_batch + + call initialize_generation() + + ! Start timer for transport + call time_transport % start() + + ! ==================================================================== + ! LOOP OVER PARTICLES +!$omp parallel do schedule(static) firstprivate(p) copyin(tally_derivs) + PARTICLE_LOOP: do i_work = 1, work + current_work = i_work + + ! grab source particle from bank + call initialize_history(p, current_work) + + ! transport particle + call transport(p) + + end do PARTICLE_LOOP +!$omp end parallel do + + ! Accumulate time for transport + call time_transport % stop() + + call finalize_generation() + + end do GENERATION_LOOP + + call finalize_batch() + + end subroutine openmc_next_batch + !=============================================================================== ! INITIALIZE_HISTORY !=============================================================================== @@ -184,7 +182,10 @@ contains ! Reset total starting particle weight used for normalizing tallies total_weight = ZERO - if (current_batch == n_inactive + 1) then + if (n_inactive > 0 .and. current_batch == 1) then + ! Turn on inactive timer + call time_inactive % start() + elseif (current_batch == n_inactive + 1) then ! Switch from inactive batch timer to active batch timer call time_inactive % stop() call time_active % start() @@ -385,7 +386,7 @@ contains ! INITIALIZE_SIMULATION !=============================================================================== - subroutine initialize_simulation() + subroutine openmc_simulation_init() bind(C) ! Set up tally procedure pointers call init_tally_routines() @@ -426,14 +427,14 @@ contains end if end if - end subroutine initialize_simulation + end subroutine openmc_simulation_init !=============================================================================== ! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays ! execution time and results !=============================================================================== - subroutine finalize_simulation() + subroutine openmc_simulation_finalize() bind(C) #ifdef MPI integer :: i ! loop index for tallies @@ -496,7 +497,7 @@ contains if (check_overlaps) call print_overlap_check() end if - end subroutine finalize_simulation + end subroutine openmc_simulation_finalize !=============================================================================== ! CALCULATE_WORK determines how many particles each processor should simulate From ac565724dcb06c82ef847f20c68ff3010f8bd1a9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 3 Nov 2017 09:46:48 -0500 Subject: [PATCH 04/26] Have openmc_next_batch return a status --- src/simulation.F90 | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/simulation.F90 b/src/simulation.F90 index 5fd6f60f70..e45e431268 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -58,10 +58,9 @@ contains call openmc_simulation_init() - BATCH_LOOP: do current_batch = 1, n_max_batches - call openmc_next_batch() - if (satisfy_triggers) exit BATCH_LOOP - end do BATCH_LOOP + do + if (openmc_next_batch() < 0) exit + end do call time_active % stop() @@ -73,7 +72,8 @@ contains ! OPENMC_NEXT_BATCH !=============================================================================== - subroutine openmc_next_batch() bind(C) + function openmc_next_batch() result(retval) bind(C) + integer(C_INT) :: retval type(Particle) :: p integer(8) :: i_work @@ -119,7 +119,16 @@ contains call finalize_batch() - end subroutine openmc_next_batch + ! Check simulation ending criteria + if (current_batch == n_max_batches) then + retval = -1 + elseif (satisfy_triggers) then + retval = -2 + else + retval = 0 + end if + + end function openmc_next_batch !=============================================================================== ! INITIALIZE_HISTORY @@ -174,6 +183,9 @@ contains integer :: i + ! Increment current batch + current_batch = current_batch + 1 + if (run_mode == MODE_FIXEDSOURCE) then call write_message("Simulating batch " // trim(to_str(current_batch)) & // "...", 6) @@ -427,6 +439,9 @@ contains end if end if + ! Reset current batch + current_batch = 0 + end subroutine openmc_simulation_init !=============================================================================== From f14957940dad27da3258ce0b11cdfa3e2ba666c2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 3 Nov 2017 21:45:49 -0500 Subject: [PATCH 05/26] Add openmc.capi.global_tallies() --- openmc/capi/tally.py | 33 +++++++++++++++++++++++++++++---- src/tallies/tally_header.F90 | 15 ++++++++++++++- 2 files changed, 43 insertions(+), 5 deletions(-) diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 7037484aea..14f9d1c1ce 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -2,6 +2,7 @@ from collections import Mapping from ctypes import c_int, c_int32, c_double, c_char_p, POINTER from weakref import WeakValueDictionary +import numpy as np from numpy.ctypeslib import as_array from openmc.data.reaction import REACTION_NAME @@ -11,15 +12,18 @@ from .error import _error_handler, AllocationError, InvalidIDError from .filter import _get_filter -__all__ = ['Tally', 'tallies'] +__all__ = ['Tally', 'tallies', 'global_tallies'] # Tally functions -_dll.openmc_get_tally_index.argtypes = [c_int32, POINTER(c_int32)] -_dll.openmc_get_tally_index.restype = c_int -_dll.openmc_get_tally_index.errcheck = _error_handler _dll.openmc_extend_tallies.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)] _dll.openmc_extend_tallies.restype = c_int _dll.openmc_extend_tallies.errcheck = _error_handler +_dll.openmc_get_tally_index.argtypes = [c_int32, POINTER(c_int32)] +_dll.openmc_get_tally_index.restype = c_int +_dll.openmc_get_tally_index.errcheck = _error_handler +_dll.openmc_global_tallies.argtypes = [POINTER(POINTER(c_double))] +_dll.openmc_global_tallies.restype = c_int +_dll.openmc_global_tallies.errcheck = _error_handler _dll.openmc_tally_get_id.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_tally_get_id.restype = c_int _dll.openmc_tally_get_id.errcheck = _error_handler @@ -65,6 +69,27 @@ _SCORES = { } +def global_tallies(): + ptr = POINTER(c_double)() + _dll.openmc_global_tallies(ptr) + array = as_array(ptr, (4, 3)) + + # Get sum, sum-of-squares, and number of realizations + sum_ = array[:, 1] + sum_sq = array[:, 2] + n = c_int32.in_dll(_dll, 'n_realizations').value + + # Determine mean + mean = sum_ / n + + # Determine standard deviation + nonzero = np.abs(mean) > 0 + stdev = np.zeros_like(mean) + stdev[nonzero] = np.sqrt((sum_sq[nonzero]/n - mean[nonzero]**2)/(n - 1)) + + return list(zip(mean, stdev)) + + class Tally(_FortranObjectWithID): """Tally stored internally. diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index e137a13274..5aa16c2d06 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -146,7 +146,7 @@ module tally_header type(VectorInt), public :: active_surface_tallies ! Normalization for statistics - integer, public :: n_realizations = 0 ! # of independent realizations + integer(C_INT32_T), public, bind(C) :: n_realizations = 0 ! # of independent realizations real(8), public :: total_weight ! total starting particle weight in realization contains @@ -470,6 +470,19 @@ contains end function openmc_get_tally_index + function openmc_global_tallies(ptr) result(err) bind(C) + type(C_PTR), intent(out) :: ptr + integer(C_INT) :: err + + if (.not. allocated(global_tallies)) then + err = E_ALLOCATE + else + err = 0 + ptr = C_LOC(global_tallies) + end if + end function openmc_global_tallies + + function openmc_tally_get_id(index, id) result(err) bind(C) ! Return the ID of a tally integer(C_INT32_T), value :: index From 9ddc7316e494a8a81490a543bb5ba9222f03c1f4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Nov 2017 16:26:28 -0600 Subject: [PATCH 06/26] Add simulation_init/finalize, iter_batches, source_bank in openmc.capi --- openmc/capi/core.py | 73 +++++++++++++++++++++++++++++++++++- src/bank_header.F90 | 22 +++++++++++ src/simulation.F90 | 4 +- src/state_point.F90 | 6 +-- src/tallies/tally_header.F90 | 2 + 5 files changed, 101 insertions(+), 6 deletions(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 0d8ff5c436..530ed50635 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -1,11 +1,23 @@ from contextlib import contextmanager -from ctypes import CDLL, c_int, c_int32, c_double, POINTER +from ctypes import CDLL, c_int, c_int32, c_int64, c_double, POINTER, Structure from warnings import warn +import numpy as np +from numpy.ctypeslib import as_array + from . import _dll from .error import _error_handler +class _Bank(Structure): + _fields_ = [('wgt', c_double), + ('xyz', c_double*3), + ('uvw', c_double*3), + ('E', c_double), + ('delayed_group', c_int)] + + + _dll.openmc_calculate_volumes.restype = None _dll.openmc_finalize.restype = None _dll.openmc_find.argtypes = [POINTER(c_double*3), c_int, POINTER(c_int32), @@ -18,9 +30,16 @@ _dll.openmc_init.restype = None _dll.openmc_get_keff.argtypes = [POINTER(c_double*2)] _dll.openmc_get_keff.restype = c_int _dll.openmc_get_keff.errcheck = _error_handler +_dll.openmc_next_batch.restype = c_int _dll.openmc_plot_geometry.restype = None _dll.openmc_run.restype = None _dll.openmc_reset.restype = None +_dll.openmc_source_bank.argtypes = [POINTER(POINTER(_Bank)), POINTER(c_int64)] +_dll.openmc_source_bank.restype = c_int +_dll.openmc_source_bank.errcheck = _error_handler +_dll.openmc_simulation_init.restype = None +_dll.openmc_simulation_finalize.restype = None +_dll.openmc_statepoint_write.restype = None def calculate_volumes(): @@ -102,6 +121,20 @@ def init(intracomm=None): _dll.openmc_init(None) +def iter_batches(): + """Iterator over batches.""" + while True: + # Run next batch + retval = _dll.openmc_next_batch() + + # Provide opportunity for user to perform action between batches + yield + + # End the iteration + if retval < 0: + break + + def keff(): """Return the calculated k-eigenvalue and its standard deviation. @@ -115,6 +148,10 @@ def keff(): _dll.openmc_get_keff(k) return tuple(k) +def next_batch(): + """Run next batch.""" + return _dll.openmc_next_batch() + def plot_geometry(): """Plot geometry""" @@ -131,6 +168,40 @@ def run(): _dll.openmc_run() +def simulation_init(): + """Initialize simulation""" + _dll.openmc_simulation_init() + + +def simulation_finalize(): + """Finalize simulation""" + _dll.openmc_simulation_finalize() + + +def source_bank(): + """Return source bank as NumPy array + + Returns + ------- + numpy.ndarray + Source sites + + """ + # Get pointer to source bank + ptr = POINTER(_Bank)() + n = c_int64() + _dll.openmc_source_bank(ptr, n) + + # Convert to numpy array with appropriate datatype + bank_dtype = np.dtype(_Bank) + return as_array(ptr, (n.value,)).view(bank_dtype) + + +def statepoint_write(): + """Write a statepoint.""" + _dll.openmc_statepoint_write() + + @contextmanager def run_in_memory(intracomm=None): """Provides context manager for calling OpenMC shared library functions. diff --git a/src/bank_header.F90 b/src/bank_header.F90 index 961824a406..10dea5b8eb 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -2,6 +2,8 @@ module bank_header use, intrinsic :: ISO_C_BINDING + use error, only: E_ALLOCATE, set_errmsg + implicit none !=============================================================================== @@ -48,4 +50,24 @@ contains end subroutine free_memory_bank +!=============================================================================== +! C API FUNCTIONS +!=============================================================================== + + function openmc_source_bank(ptr, n) result(err) bind(C) + ! Return a pointer to the source bank + type(C_PTR), intent(out) :: ptr + integer(C_INT64_T), intent(out) :: n + integer(C_INT) :: err + + if (.not. allocated(source_bank)) then + err = E_ALLOCATE + call set_errmsg("Source bank has not been allocated.") + else + err = 0 + ptr = C_LOC(source_bank) + n = size(source_bank) + end if + end function openmc_source_bank + end module bank_header diff --git a/src/simulation.F90 b/src/simulation.F90 index e45e431268..94a285688a 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -29,7 +29,7 @@ module simulation use settings use simulation_header use source, only: initialize_source, sample_external_source - use state_point, only: write_state_point, write_source_point, load_state_point + use state_point, only: openmc_statepoint_write, write_source_point, load_state_point use string, only: to_str use tally, only: accumulate_tallies, setup_active_tallies, & init_tally_routines @@ -349,7 +349,7 @@ contains ! Write out state point if it's been specified for this batch if (statepoint_batch % contains(current_batch)) then - call write_state_point() + call openmc_statepoint_write() end if ! Write out source point if it's been specified for this batch diff --git a/src/state_point.F90 b/src/state_point.F90 index 5dde708aa9..265cc49c02 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -40,10 +40,10 @@ module state_point contains !=============================================================================== -! WRITE_STATE_POINT +! OPENMC_STATEPOINT_WRITE writes an HDF5 statepoint file to disk !=============================================================================== - subroutine write_state_point() + subroutine openmc_statepoint_write() bind(C) integer :: i, j, k integer :: i_xs @@ -433,7 +433,7 @@ contains call file_close(file_id) end if - end subroutine write_state_point + end subroutine openmc_statepoint_write !=============================================================================== ! WRITE_SOURCE_POINT diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 5aa16c2d06..2812eef6ec 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -22,6 +22,7 @@ module tally_header public :: free_memory_tally public :: openmc_extend_tallies public :: openmc_get_tally_index + public :: openmc_global_tallies public :: openmc_tally_get_id public :: openmc_tally_get_filters public :: openmc_tally_get_nuclides @@ -476,6 +477,7 @@ contains if (.not. allocated(global_tallies)) then err = E_ALLOCATE + call set_errmsg("Global tallies have not been allocated yet.") else err = 0 ptr = C_LOC(global_tallies) From 7d00c0b5eb41cdd8ef7d7b665539a2459cb7a4b9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Nov 2017 23:30:02 -0600 Subject: [PATCH 07/26] Account for possibility of running beyond n_batches --- src/simulation.F90 | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/simulation.F90 b/src/simulation.F90 index 94a285688a..ad2841f54d 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -57,13 +57,8 @@ contains subroutine openmc_run() bind(C) call openmc_simulation_init() - - do - if (openmc_next_batch() < 0) exit + do while (openmc_next_batch() == 0) end do - - call time_active % stop() - call openmc_simulation_finalize() end subroutine openmc_run @@ -459,12 +454,15 @@ contains real(8) :: tempr(3) ! temporary array for communication #endif + ! Stop active batch timer + call time_active % stop() + !$omp parallel deallocate(micro_xs, filter_matches) !$omp end parallel ! Increment total number of generations - total_gen = total_gen + n_batches*gen_per_batch + total_gen = total_gen + current_batch*gen_per_batch ! Start finalization timer call time_finalize % start() From be28e1322031e04969f26f363d021ba4411de9d8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 Nov 2017 11:17:42 -0600 Subject: [PATCH 08/26] Allow k_generation and entropy to grow beyond n_max_batches --- src/eigenvalue.F90 | 28 ++++++++++++++-------------- src/input_xml.F90 | 7 +++---- src/output.F90 | 8 ++++---- src/simulation.F90 | 1 + src/simulation_header.F90 | 12 ++++++++---- src/state_point.F90 | 15 ++++++++++----- 6 files changed, 40 insertions(+), 31 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 5f99d1580a..27ae4063ea 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -301,8 +301,8 @@ contains subroutine shannon_entropy() - integer :: ent_idx ! entropy index integer :: i ! index for mesh elements + real(8) :: entropy_gen ! entropy at this generation logical :: sites_outside ! were there sites outside entropy box? associate (m => meshes(index_entropy_mesh)) @@ -320,14 +320,16 @@ contains ! Normalize to total weight of bank sites entropy_p = entropy_p / sum(entropy_p) - ent_idx = current_gen + gen_per_batch*(current_batch - 1) - entropy(ent_idx) = ZERO + entropy_gen = ZERO do i = 1, size(entropy_p, 2) if (entropy_p(1,i) > ZERO) then - entropy(ent_idx) = entropy(ent_idx) - & + entropy_gen = entropy_gen - & entropy_p(1,i) * log(entropy_p(1,i))/log(TWO) end if end do + + ! Add value to vector + call entropy % push_back(entropy_gen) end if end associate end subroutine shannon_entropy @@ -340,7 +342,7 @@ contains subroutine calculate_generation_keff() - integer :: i ! overall generation + real(8) :: keff_reduced #ifdef MPI integer :: mpi_err ! MPI error code #endif @@ -348,20 +350,18 @@ contains ! Get keff for this generation by subtracting off the starting value keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) - keff_generation - ! Determine overall generation - i = overall_generation() - #ifdef MPI ! Combine values across all processors - call MPI_ALLREDUCE(keff_generation, k_generation(i), 1, MPI_REAL8, & + call MPI_ALLREDUCE(keff_generation, keff_reduced, 1, MPI_REAL8, & MPI_SUM, mpi_intracomm, mpi_err) #else - k_generation(i) = keff_generation + keff_reduced = keff_generation #endif ! Normalize single batch estimate of k ! TODO: This should be normalized by total_weight, not by n_particles - k_generation(i) = k_generation(i) / n_particles + keff_reduced = keff_reduced / n_particles + call k_generation % push_back(keff_reduced) end subroutine calculate_generation_keff @@ -385,12 +385,12 @@ contains if (n <= 0) then ! For inactive generations, use current generation k as estimate for next ! generation - keff = k_generation(i) + keff = k_generation % data(i) else ! Sample mean of keff - k_sum(1) = k_sum(1) + k_generation(i) - k_sum(2) = k_sum(2) + k_generation(i)**2 + k_sum(1) = k_sum(1) + k_generation % data(i) + k_sum(2) = k_sum(2) + k_generation % data(i)**2 ! Determine mean keff = k_sum(1) / n diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 71ebe3ae6b..69b519ea5a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -861,10 +861,9 @@ contains call get_node_value(node_base, "generations_per_batch", gen_per_batch) end if - ! Allocate array for batch keff and entropy - allocate(k_generation(n_max_batches*gen_per_batch)) - allocate(entropy(n_max_batches*gen_per_batch)) - entropy = ZERO + ! Preallocate space for keff and entropy by generation + call k_generation % reserve(n_max_batches*gen_per_batch) + call entropy % initialize(n_max_batches*gen_per_batch) ! Get the trigger information for keff if (check_for_node(node_base, "keff_trigger")) then diff --git a/src/output.F90 b/src/output.F90 index f820ce4b9d..3d02489f78 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -384,11 +384,11 @@ contains ! write out information about batch and generation write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & trim(to_str(current_batch)) // "/" // trim(to_str(current_gen)) - write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation(i) + write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') k_generation % data(i) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - entropy(i) + entropy % data(i) if (n > 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & @@ -418,11 +418,11 @@ contains write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & trim(to_str(current_batch)) // "/" // trim(to_str(gen_per_batch)) write(UNIT=OUTPUT_UNIT, FMT='(3X,F8.5)', ADVANCE='NO') & - k_generation(i) + k_generation % data(i) ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - entropy(i) + entropy % data(i) ! write out accumulated k-effective if after first active batch if (n > 1) then diff --git a/src/simulation.F90 b/src/simulation.F90 index ad2841f54d..e3eba0d512 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -78,6 +78,7 @@ contains ! Handle restart runs if (restart_run .and. current_batch <= restart_batch) then call replay_batch_history() + retval = 0 return end if diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index cfb75b4492..9719cc72e3 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -3,6 +3,7 @@ module simulation_header use bank_header use constants use settings, only: gen_per_batch + use stl_vector, only: VectorReal implicit none @@ -31,7 +32,7 @@ module simulation_header integer(8) :: current_work ! index in source bank of current history simulated ! Temporary k-effective values - real(8), allocatable :: k_generation(:) ! single-generation estimates of k + type(VectorReal) :: k_generation ! single-generation estimates of k real(8) :: keff = ONE ! average k over active batches real(8) :: keff_std ! standard deviation of average k real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption @@ -39,7 +40,7 @@ module simulation_header real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength ! Shannon entropy - real(8), allocatable :: entropy(:) ! shannon entropy at each generation + type(VectorReal) :: entropy ! shannon entropy at each generation real(8), allocatable :: entropy_p(:,:) ! % of source sites in each cell ! Uniform fission source weighting @@ -85,11 +86,14 @@ contains subroutine free_memory_simulation() if (allocated(overlap_check_cnt)) deallocate(overlap_check_cnt) - if (allocated(k_generation)) deallocate(k_generation) - if (allocated(entropy)) deallocate(entropy) if (allocated(entropy_p)) deallocate(entropy_p) if (allocated(source_frac)) deallocate(source_frac) if (allocated(work_index)) deallocate(work_index) + + call k_generation % clear() + call k_generation % shrink_to_fit() + call entropy % clear() + call entropy % shrink_to_fit() end subroutine free_memory_simulation end module simulation_header diff --git a/src/state_point.F90 b/src/state_point.F90 index 265cc49c02..e93d6e4729 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -121,8 +121,11 @@ contains if (run_mode == MODE_EIGENVALUE) then call write_dataset(file_id, "n_inactive", n_inactive) call write_dataset(file_id, "generations_per_batch", gen_per_batch) - call write_dataset(file_id, "k_generation", k_generation) - call write_dataset(file_id, "entropy", entropy) + k = k_generation % size() + call write_dataset(file_id, "k_generation", k_generation % data(1:k)) + if (entropy_on) then + call write_dataset(file_id, "entropy", entropy % data(1:k)) + end if call write_dataset(file_id, "k_col_abs", k_col_abs) call write_dataset(file_id, "k_col_tra", k_col_tra) call write_dataset(file_id, "k_abs_tra", k_abs_tra) @@ -707,10 +710,12 @@ contains if (run_mode == MODE_EIGENVALUE) then call read_dataset(int_array(1), file_id, "n_inactive") call read_dataset(gen_per_batch, file_id, "generations_per_batch") - call read_dataset(k_generation(1:restart_batch*gen_per_batch), & + call read_dataset(k_generation % data(1:restart_batch*gen_per_batch), & file_id, "k_generation") - call read_dataset(entropy(1:restart_batch*gen_per_batch), & - file_id, "entropy") + if (entropy_on) then + call read_dataset(entropy % data(1:restart_batch*gen_per_batch), & + file_id, "entropy") + end if call read_dataset(k_col_abs, file_id, "k_col_abs") call read_dataset(k_col_tra, file_id, "k_col_tra") call read_dataset(k_abs_tra, file_id, "k_abs_tra") From 671db024db057c17350666429890a805515effdc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 Nov 2017 14:01:26 -0600 Subject: [PATCH 09/26] Fix bugs with statepoint loading. Add properties to capi.Tally. --- docs/source/conf.py | 5 +- docs/source/pythonapi/capi.rst | 6 ++ openmc/capi/core.py | 22 ++++++- openmc/capi/tally.py | 117 +++++++++++++++++++++++++-------- src/input_xml.F90 | 2 +- src/state_point.F90 | 12 ++-- src/tallies/tally_header.F90 | 16 +++++ 7 files changed, 143 insertions(+), 37 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 95930fedba..5f993b1a3b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -26,8 +26,9 @@ except ImportError: MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial', 'numpy.ctypeslib', 'scipy', 'scipy.sparse', 'scipy.interpolate', - 'scipy.integrate', 'scipy.optimize', 'scipy.special', 'h5py', - 'pandas', 'uncertainties', 'openmoc', 'openmc.data.reconstruct'] + 'scipy.integrate', 'scipy.optimize', 'scipy.special', + 'scipy.stats','h5py', 'pandas', 'uncertainties', 'openmoc', + 'openmc.data.reconstruct'] sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES) import numpy as np diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 457f053204..44a755bbd1 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -18,12 +18,18 @@ Functions openmc.capi.find_material openmc.capi.hard_reset openmc.capi.init + openmc.capi.iter_batches openmc.capi.keff openmc.capi.load_nuclide + openmc.capi.next_batch openmc.capi.plot_geometry openmc.capi.reset openmc.capi.run openmc.capi.run_in_memory + openmc.capi.simulation_init + openmc.capi.simulation_finalize + openmc.capi.source_bank + openmc.capi.statepoint_write Classes ------- diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 530ed50635..77c9ca7a46 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -122,7 +122,27 @@ def init(intracomm=None): def iter_batches(): - """Iterator over batches.""" + """Iterator over batches. + + This function returns a generator-iterator that allows Python code to be run + between batches in an OpenMC simulation. It should be used in conjunction + with :func:`openmc.capi.simulation_init` and + :func:`openmc.capi.simulation_finalize`. For example: + + .. code-block:: Python + + with openmc.capi.run_in_memory(): + openmc.capi.simulation_init() + for _ in openmc.capi.iter_batches(): + # Look at convergence of tallies, for example + ... + openmc.capi.simulation_finalize() + + See Also + -------- + openmc.capi.next_batch + + """ while True: # Run next batch retval = _dll.openmc_next_batch() diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 14f9d1c1ce..6725ee5c92 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -4,6 +4,7 @@ from weakref import WeakValueDictionary import numpy as np from numpy.ctypeslib import as_array +import scipy.stats from openmc.data.reaction import REACTION_NAME from . import _dll, Nuclide @@ -31,6 +32,9 @@ _dll.openmc_tally_get_filters.argtypes = [ c_int32, POINTER(POINTER(c_int32)), POINTER(c_int)] _dll.openmc_tally_get_filters.restype = c_int _dll.openmc_tally_get_filters.errcheck = _error_handler +_dll.openmc_tally_get_n_realizations.argtypes = [c_int32, POINTER(c_int32)] +_dll.openmc_tally_get_n_realizations.restype = c_int +_dll.openmc_tally_get_n_realizations.errcheck = _error_handler _dll.openmc_tally_get_nuclides.argtypes = [ c_int32, POINTER(POINTER(c_int)), POINTER(c_int)] _dll.openmc_tally_get_nuclides.restype = c_int @@ -70,6 +74,14 @@ _SCORES = { def global_tallies(): + """Mean and standard deviation of the mean for each global tally. + + Returns + ------- + list of tuple + For each global tally, a tuple of (mean, standard deviation) + + """ ptr = POINTER(c_double)() _dll.openmc_global_tallies(ptr) array = as_array(ptr, (4, 3)) @@ -113,10 +125,16 @@ class Tally(_FortranObjectWithID): ID of the tally filters : list List of tally filters + mean : numpy.ndarray + An array containing the sample mean for each bin nuclides : list of str List of nuclides to score results for + num_realizations : int + Number of realizations results : numpy.ndarray Array of tally results + std_dev : numpy.ndarray + An array containing the sample standard deviation for each bin """ __instances = WeakValueDictionary() @@ -169,21 +187,6 @@ class Tally(_FortranObjectWithID): _dll.openmc_tally_get_filters(self._index, filt_idx, n) return [_get_filter(filt_idx[i]) for i in range(n.value)] - @property - def nuclides(self): - nucs = POINTER(c_int)() - n = c_int() - _dll.openmc_tally_get_nuclides(self._index, nucs, n) - return [Nuclide(nucs[i]).name if nucs[i] > 0 else 'total' - for i in range(n.value)] - - @property - def results(self): - data = POINTER(c_double)() - shape = (c_int*3)() - _dll.openmc_tally_results(self._index, data, shape) - return as_array(data, tuple(shape[::-1])) - @filters.setter def filters(self, filters): # Get filter indices as int32_t[] @@ -192,12 +195,42 @@ class Tally(_FortranObjectWithID): _dll.openmc_tally_set_filters(self._index, n, indices) + @property + def mean(self): + n = self.num_realizations + sum_ = self.results[:, :, 1] + if n > 0: + return sum_ / n + else: + return sum_.copy() + + @property + def nuclides(self): + nucs = POINTER(c_int)() + n = c_int() + _dll.openmc_tally_get_nuclides(self._index, nucs, n) + return [Nuclide(nucs[i]).name if nucs[i] > 0 else 'total' + for i in range(n.value)] + @nuclides.setter def nuclides(self, nuclides): nucs = (c_char_p * len(nuclides))() nucs[:] = [x.encode() for x in nuclides] _dll.openmc_tally_set_nuclides(self._index, len(nuclides), nucs) + @property + def num_realizations(self): + n = c_int32() + _dll.openmc_tally_get_n_realizations(self._index, n) + return n.value + + @property + def results(self): + data = POINTER(c_double)() + shape = (c_int*3)() + _dll.openmc_tally_results(self._index, data, shape) + return as_array(data, tuple(shape[::-1])) + @property def scores(self): scores_as_int = POINTER(c_int)() @@ -223,21 +256,47 @@ class Tally(_FortranObjectWithID): scores_[:] = [x.encode() for x in scores] _dll.openmc_tally_set_scores(self._index, len(scores), scores_) - @classmethod - def new(cls, tally_id=None): - # Determine ID to assign - if tally_id is None: - try: - tally_id = max(tallies) + 1 - except ValueError: - tally_id = 1 + @property + def std_dev(self): + results = self.results + std_dev = np.empty(results.shape[:2]) + std_dev.fill(np.nan) - index = c_int32() - _dll.openmc_extend_tallies(1, index, None) - _dll.openmc_tally_set_type(index, b'generic') - tally = cls(index.value) - tally.id = tally_id - return tally + n = self.num_realizations + if n > 1: + # Get sum and sum-of-squares from results + sum_ = results[:, :, 1] + sum_sq = results[:, :, 2] + + # Determine non-zero entries + mean = sum_ / n + nonzero = np.abs(mean) > 0 + + # Calculate sample standard deviation of the mean + std_dev[nonzero] = np.sqrt( + (sum_sq[nonzero]/n - mean[nonzero]**2)/(n - 1)) + + return std_dev + + def ci_width(self, alpha=0.05): + """Confidence interval half-width based on a Student t distribution + + Parameters + ---------- + alpha : float + Significance level (one minus the confidence level!) + + Returns + ------- + float + Half-width of a two-sided (1 - :math:`alpha`) confidence interval + + """ + half_width = self.std_dev.copy() + n = self.num_realizations + if n > 1: + half_width *= scipy.stats.t.ppf(1 - alpha/2, n - 1) + return half_width class _TallyMapping(Mapping): diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 69b519ea5a..766dbfa9f2 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -863,7 +863,7 @@ contains ! Preallocate space for keff and entropy by generation call k_generation % reserve(n_max_batches*gen_per_batch) - call entropy % initialize(n_max_batches*gen_per_batch) + call entropy % reserve(n_max_batches*gen_per_batch) ! Get the trigger information for keff if (check_for_node(node_base, "keff_trigger")) then diff --git a/src/state_point.F90 b/src/state_point.F90 index e93d6e4729..0383150ac4 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -632,6 +632,7 @@ contains subroutine load_state_point() integer :: i + integer :: n integer :: int_array(3) integer, allocatable :: array(:) integer(HID_T) :: file_id @@ -710,11 +711,14 @@ contains if (run_mode == MODE_EIGENVALUE) then call read_dataset(int_array(1), file_id, "n_inactive") call read_dataset(gen_per_batch, file_id, "generations_per_batch") - call read_dataset(k_generation % data(1:restart_batch*gen_per_batch), & - file_id, "k_generation") + + n = restart_batch*gen_per_batch + call k_generation % resize(n) + call read_dataset(k_generation % data(1:n), file_id, "k_generation") + if (entropy_on) then - call read_dataset(entropy % data(1:restart_batch*gen_per_batch), & - file_id, "entropy") + call entropy % resize(n) + call read_dataset(entropy % data(1:n), file_id, "entropy") end if call read_dataset(k_col_abs, file_id, "k_col_abs") call read_dataset(k_col_tra, file_id, "k_col_tra") diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 2812eef6ec..1778474e04 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -526,6 +526,22 @@ contains end function openmc_tally_get_filters + function openmc_tally_get_n_realizations(index, n) result(err) bind(C) + ! Return the number of realizations for a tally + integer(C_INT32_T), value :: index + integer(C_INT32_T), intent(out) :: n + integer(C_INT) :: err + + if (index >= 1 .and. index <= size(tallies)) then + n = tallies(index) % obj % n_realizations + err = 0 + else + err = E_OUT_OF_BOUNDS + call set_errmsg('Index in tallies array is out of bounds.') + end if + end function openmc_tally_get_n_realizations + + function openmc_tally_get_nuclides(index, nuclides, n) result(err) bind(C) ! Return the list of nuclides assigned to a tally integer(C_INT32_T), value :: index From 198df07311cc5747e8eb5915c8dcb209b5fdeb63 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 Nov 2017 15:53:57 -0600 Subject: [PATCH 10/26] Add subdivide convenience function --- docs/source/pythonapi/index.rst | 8 ++++---- docs/source/pythonapi/model.rst | 10 ++++++++++ openmc/capi/core.py | 2 +- openmc/model/__init__.py | 1 + openmc/model/funcs.py | 25 +++++++++++++++++++++++++ 5 files changed, 41 insertions(+), 5 deletions(-) create mode 100644 openmc/model/funcs.py diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index bb99dd470f..862acb2384 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -40,13 +40,13 @@ Modules ------- .. toctree:: - :maxdepth: 2 + :maxdepth: 1 base - stats - mgxs model + examples + mgxs + stats data capi - examples openmoc diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 6e589b7eb1..55a5d17b12 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -2,6 +2,16 @@ :mod:`openmc.model` -- Model Building ------------------------------------- +Convenience Functions +--------------------- + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myfunction.rst + + openmc.model.subdivide + TRISO Fuel Modeling ------------------- diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 77c9ca7a46..1f7225f642 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -17,7 +17,6 @@ class _Bank(Structure): ('delayed_group', c_int)] - _dll.openmc_calculate_volumes.restype = None _dll.openmc_finalize.restype = None _dll.openmc_find.argtypes = [POINTER(c_double*3), c_int, POINTER(c_int32), @@ -168,6 +167,7 @@ def keff(): _dll.openmc_get_keff(k) return tuple(k) + def next_batch(): """Run next batch.""" return _dll.openmc_next_batch() diff --git a/openmc/model/__init__.py b/openmc/model/__init__.py index 557effcefa..9fa999dd4e 100644 --- a/openmc/model/__init__.py +++ b/openmc/model/__init__.py @@ -1,2 +1,3 @@ from .triso import * from .model import * +from .funcs import * diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py new file mode 100644 index 0000000000..d9755b0e8f --- /dev/null +++ b/openmc/model/funcs.py @@ -0,0 +1,25 @@ +def subdivide(surfaces): + """Create regions separated by a series of surfaces. + + This function allows regions to be constructed from a set of a surfaces that + are "in order". For example, if you had four instances of + :class:`openmc.ZPlane` at z=-10, z=-5, z=5, and z=10, this function would + return a list of regions corresponding to z < -10, -10 < z < -5, -5 < z < 5, + 5 < z < 10, and 10 < z. That is, for n surfaces, n+1 regions are returned. + + Parameters + ---------- + surfaces : sequence of openmc.Surface + Surfaces separating regions + + Returns + ------- + list of openmc.Region + Regions formed by the given surfaces + + """ + regions = [-surfaces[0]] + for s0, s1 in zip(surfaces[:-1], surfaces[1:]): + regions.append(+s0 & -s1) + regions.append(+surfaces[-1]) + return regions From 6e9b8a5ce52ae15d147bc503c908e947963ae0bf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 Nov 2017 16:06:58 -0600 Subject: [PATCH 11/26] Move get_*_prism functions to openmc.model --- docs/source/pythonapi/base.rst | 12 -- docs/source/pythonapi/model.rst | 7 + openmc/__init__.py | 3 + openmc/model/funcs.py | 254 ++++++++++++++++++++++++++++++++ openmc/surface.py | 253 +------------------------------ 5 files changed, 266 insertions(+), 263 deletions(-) diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index e6cf1e1ccc..f1633f3f9d 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -91,18 +91,6 @@ Many of the above classes are derived from several abstract classes: openmc.Region openmc.Lattice -Two helper function are also available to create rectangular and hexagonal -prisms defined by the intersection of four and six surface half-spaces, -respectively. - -.. autosummary:: - :toctree: generated - :nosignatures: - :template: myfunction.rst - - openmc.get_hexagonal_prism - openmc.get_rectangular_prism - .. _pythonapi_tallies: Constructing Tallies diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 55a5d17b12..4ba247468a 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -5,11 +5,18 @@ Convenience Functions --------------------- +Several helper functions are available here. Ther first two create rectangular +and hexagonal prisms defined by the intersection of four and six surface +half-spaces, respectively. The last function takes a sequence of surfaces and +returns the regions that separate them. + .. autosummary:: :toctree: generated :nosignatures: :template: myfunction.rst + openmc.model.get_hexagonal_prism + openmc.model.get_rectangular_prism openmc.model.subdivide TRISO Fuel Modeling diff --git a/openmc/__init__.py b/openmc/__init__.py index d692ebbae2..8fb9bcf37c 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -28,4 +28,7 @@ from openmc.mixin import * from openmc.plotter import * from openmc.search import * +# Import a few convencience functions that used to be here +from openmc.model import get_rectangular_prism, get_hexagonal_prism + __version__ = '0.9.0' diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index d9755b0e8f..589765d2c8 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -1,3 +1,257 @@ +from __future__ import division +from collections import Iterable, OrderedDict +from math import sqrt +from numbers import Real + +from openmc import XPlane, YPlane, Plane, ZCylinder +from openmc.checkvalue import check_type, check_value + + +def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), + boundary_type='transmission', corner_radius=0.): + """Get an infinite rectangular prism from four planar surfaces. + + Parameters + ---------- + width: float + Prism width in units of cm. The width is aligned with the y, x, + or x axes for prisms parallel to the x, y, or z axis, respectively. + height: float + Prism height in units of cm. The height is aligned with the z, z, + or y axes for prisms parallel to the x, y, or z axis, respectively. + axis : {'x', 'y', 'z'} + Axis with which the infinite length of the prism should be aligned. + Defaults to 'z'. + origin: Iterable of two floats + Origin of the prism. The two floats correspond to (y,z), (x,z) or + (x,y) for prisms parallel to the x, y or z axis, respectively. + Defaults to (0., 0.). + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + Boundary condition that defines the behavior for particles hitting the + surfaces comprising the rectangular prism (default is 'transmission'). + corner_radius: float + Prism corner radius in units of cm. Defaults to 0. + + Returns + ------- + openmc.Region + The inside of a rectangular prism + + """ + + check_type('width', width, Real) + check_type('height', height, Real) + check_type('corner_radius', corner_radius, Real) + check_value('axis', axis, ['x', 'y', 'z']) + check_type('origin', origin, Iterable, Real) + + # Define function to create a plane on given axis + def plane(axis, name, value): + cls = globals()['{}Plane'.format(axis.upper())] + return cls(name='{} {}'.format(name, axis), + boundary_type=boundary_type, + **{axis + '0': value}) + + if axis == 'x': + x1, x2 = 'y', 'z' + elif axis == 'y': + x1, x2 = 'x', 'z' + else: + x1, x2 = 'x', 'y' + + # Get cylinder class corresponding to given axis + cyl = globals()['{}Cylinder'.format(axis.upper())] + + # Create rectangular region + min_x1 = plane(x1, 'minimum', -width/2 + origin[0]) + max_x1 = plane(x1, 'maximum', width/2 + origin[0]) + min_x2 = plane(x2, 'minimum', -height/2 + origin[1]) + max_x2 = plane(x2, 'maximum', height/2 + origin[1]) + prism = +min_x1 & -max_x1 & +min_x2 & -max_x2 + + # Handle rounded corners if given + if corner_radius > 0.: + args = {'R': corner_radius, 'boundary_type': boundary_type} + + args[x1 + '0'] = origin[0] - width/2 + corner_radius + args[x2 + '0'] = origin[1] - height/2 + corner_radius + x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args) + + args[x1 + '0'] = origin[0] - width/2 + corner_radius + args[x2 + '0'] = origin[1] - height/2 + corner_radius + x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args) + + args[x1 + '0'] = origin[0] - width/2 + corner_radius + args[x2 + '0'] = origin[1] + height/2 - corner_radius + x1_min_x2_max = cyl(name='{} min {} max'.format(x1, x2), **args) + + args[x1 + '0'] = origin[0] + width/2 - corner_radius + args[x2 + '0'] = origin[1] - height/2 + corner_radius + x1_max_x2_min = cyl(name='{} max {} min'.format(x1, x2), **args) + + args[x1 + '0'] = origin[0] + width/2 - corner_radius + args[x2 + '0'] = origin[1] + height/2 - corner_radius + x1_max_x2_max = cyl(name='{} max {} max'.format(x1, x2), **args) + + x1_min = plane(x1, 'min', -width/2 + origin[0] + corner_radius) + x1_max = plane(x1, 'max', width/2 + origin[0] - corner_radius) + x2_min = plane(x2, 'min', -height/2 + origin[1] + corner_radius) + x2_max = plane(x2, 'max', height/2 + origin[1] - corner_radius) + + corners = (+x1_min_x2_min & -x1_min & -x2_min) | \ + (+x1_min_x2_max & -x1_min & +x2_max) | \ + (+x1_max_x2_min & +x1_max & -x2_min) | \ + (+x1_max_x2_max & +x1_max & +x2_max) + + prism = prism & ~corners + + return prism + + +def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), + boundary_type='transmission', corner_radius=0.): + """Create a hexagon region from six surface planes. + + Parameters + ---------- + edge_length : float + Length of a side of the hexagon in cm + orientation : {'x', 'y'} + An 'x' orientation means that two sides of the hexagon are parallel to + the x-axis and a 'y' orientation means that two sides of the hexagon are + parallel to the y-axis. + origin: Iterable of two floats + Origin of the prism. Defaults to (0., 0.). + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} + Boundary condition that defines the behavior for particles hitting the + surfaces comprising the hexagonal prism (default is 'transmission'). + corner_radius: float + Prism corner radius in units of cm. Defaults to 0. + + Returns + ------- + openmc.Region + The inside of a hexagonal prism + + """ + + l = edge_length + x, y = origin + + if orientation == 'y': + right = XPlane(x0=x + sqrt(3.)/2*l, boundary_type=boundary_type) + left = XPlane(x0=x - sqrt(3.)/2*l, boundary_type=boundary_type) + c = sqrt(3.)/3. + + # y = -x/sqrt(3) + a + upper_right = Plane(A=c, B=1., D=l+x*c+y, boundary_type=boundary_type) + + # y = x/sqrt(3) + a + upper_left = Plane(A=-c, B=1., D=l-x*c+y, boundary_type=boundary_type) + + # y = x/sqrt(3) - a + lower_right = Plane(A=-c, B=1., D=-l-x*c+y, boundary_type=boundary_type) + + # y = -x/sqrt(3) - a + lower_left = Plane(A=c, B=1., D=-l+x*c+y, boundary_type=boundary_type) + + prism = -right & +left & -upper_right & -upper_left & \ + +lower_right & +lower_left + + if boundary_type == 'periodic': + right.periodic_surface = left + upper_right.periodic_surface = lower_left + lower_right.periodic_surface = upper_left + + elif orientation == 'x': + top = YPlane(y0=y + sqrt(3.)/2*l, boundary_type=boundary_type) + bottom = YPlane(y0=y - sqrt(3.)/2*l, boundary_type=boundary_type) + c = sqrt(3.) + + # y = -sqrt(3)*(x - a) + upper_right = Plane(A=c, B=1., D=c*l+x*c+y, boundary_type=boundary_type) + + # y = sqrt(3)*(x + a) + lower_right = Plane(A=-c, B=1., D=-c*l-x*c+y, + boundary_type=boundary_type) + + # y = -sqrt(3)*(x + a) + lower_left = Plane(A=c, B=1., D=-c*l+x*c+y, boundary_type=boundary_type) + + # y = sqrt(3)*(x + a) + upper_left = Plane(A=-c, B=1., D=c*l-x*c+y, boundary_type=boundary_type) + + prism = -top & +bottom & -upper_right & +lower_right & \ + +lower_left & -upper_left + + if boundary_type == 'periodic': + top.periodic_surface = bottom + upper_right.periodic_surface = lower_left + lower_right.periodic_surface = upper_left + + # Handle rounded corners if given + if corner_radius > 0.: + if boundary_type == 'periodic': + raise ValueError('Periodic boundary conditions not permitted when ' + 'rounded corners are used.') + + c = sqrt(3.)/2 + t = l - corner_radius/c + + # Cylinder with corner radius and boundary type pre-applied + cyl1 = partial(ZCylinder, R=corner_radius, boundary_type=boundary_type) + cyl2 = partial(ZCylinder, R=corner_radius/(2*c), + boundary_type=boundary_type) + + if orientation == 'x': + x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t) + x_min_y_max_in = cyl1(name='x min y max in', x0=x+t/2, y0=y-c*t) + x_max_y_min_in = cyl1(name='x max y min in', x0=x-t/2, y0=y+c*t) + x_max_y_max_in = cyl1(name='x max y max in', x0=x+t/2, y0=y+c*t) + x_min_in = cyl1(name='x min in', x0=x-t, y0=y) + x_max_in = cyl1(name='x max in', x0=x+t, y0=y) + + x_min_y_min_out = cyl2(name='x min y min out', x0=x-l/2, y0=y-c*l) + x_min_y_max_out = cyl2(name='x min y max out', x0=x+l/2, y0=y-c*l) + x_max_y_min_out = cyl2(name='x max y min out', x0=x-l/2, y0=y+c*l) + x_max_y_max_out = cyl2(name='x max y max out', x0=x+l/2, y0=y+c*l) + x_min_out = cyl2(name='x min out', x0=x-l, y0=y) + x_max_out = cyl2(name='x max out', x0=x+l, y0=y) + + corners = (+x_min_y_min_in & -x_min_y_min_out | + +x_min_y_max_in & -x_min_y_max_out | + +x_max_y_min_in & -x_max_y_min_out | + +x_max_y_max_in & -x_max_y_max_out | + +x_min_in & -x_min_out | + +x_max_in & -x_max_out) + + elif orientation == 'y': + x_min_y_min_in = cyl1(name='x min y min in', x0=x-c*t, y0=y-t/2) + x_min_y_max_in = cyl1(name='x min y max in', x0=x-c*t, y0=y+t/2) + x_max_y_min_in = cyl1(name='x max y min in', x0=x+c*t, y0=y-t/2) + x_max_y_max_in = cyl1(name='x max y max in', x0=x+c*t, y0=y+t/2) + y_min_in = cyl1(name='y min in', x0=x, y0=y-t) + y_max_in = cyl1(name='y max in', x0=x, y0=y+t) + + x_min_y_min_out = cyl2(name='x min y min out', x0=x-c*l, y0=y-l/2) + x_min_y_max_out = cyl2(name='x min y max out', x0=x-c*l, y0=y+l/2) + x_max_y_min_out = cyl2(name='x max y min out', x0=x+c*l, y0=y-l/2) + x_max_y_max_out = cyl2(name='x max y max out', x0=x+c*l, y0=y+l/2) + y_min_out = cyl2(name='y min out', x0=x, y0=y-l) + y_max_out = cyl2(name='y max out', x0=x, y0=y+l) + + corners = (+x_min_y_min_in & -x_min_y_min_out | + +x_min_y_max_in & -x_min_y_max_out | + +x_max_y_min_in & -x_max_y_min_out | + +x_max_y_max_in & -x_max_y_max_out | + +y_min_in & -y_min_out | + +y_max_in & -y_max_out) + + prism = prism & ~corners + + return prism + + def subdivide(surfaces): """Create regions separated by a series of surfaces. diff --git a/openmc/surface.py b/openmc/surface.py index b415ecb93a..961ac4d459 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,23 +1,19 @@ from __future__ import division from abc import ABCMeta -from collections import Iterable, OrderedDict +from collections import OrderedDict from copy import deepcopy from functools import partial from numbers import Real, Integral from xml.etree import ElementTree as ET -from math import sqrt from six import add_metaclass, string_types import numpy as np -from openmc.checkvalue import check_type, check_value, check_greater_than +from openmc.checkvalue import check_type, check_value from openmc.region import Region, Intersection, Union from openmc.mixin import IDManagerMixin -# A static variable for auto-generated Surface IDs -AUTO_SURFACE_ID = 10000 - _BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic'] @@ -1944,248 +1940,3 @@ class Halfspace(Region): clone = deepcopy(self) clone.surface = self.surface.clone(memo) return clone - - -def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), - boundary_type='transmission', corner_radius=0.): - """Get an infinite rectangular prism from four planar surfaces. - - Parameters - ---------- - width: float - Prism width in units of cm. The width is aligned with the y, x, - or x axes for prisms parallel to the x, y, or z axis, respectively. - height: float - Prism height in units of cm. The height is aligned with the z, z, - or y axes for prisms parallel to the x, y, or z axis, respectively. - axis : {'x', 'y', 'z'} - Axis with which the infinite length of the prism should be aligned. - Defaults to 'z'. - origin: Iterable of two floats - Origin of the prism. The two floats correspond to (y,z), (x,z) or - (x,y) for prisms parallel to the x, y or z axis, respectively. - Defaults to (0., 0.). - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} - Boundary condition that defines the behavior for particles hitting the - surfaces comprising the rectangular prism (default is 'transmission'). - corner_radius: float - Prism corner radius in units of cm. Defaults to 0. - - Returns - ------- - openmc.Region - The inside of a rectangular prism - - """ - - check_type('width', width, Real) - check_type('height', height, Real) - check_type('corner_radius', corner_radius, Real) - check_value('axis', axis, ['x', 'y', 'z']) - check_type('origin', origin, Iterable, Real) - - # Define function to create a plane on given axis - def plane(axis, name, value): - cls = globals()['{}Plane'.format(axis.upper())] - return cls(name='{} {}'.format(name, axis), - boundary_type=boundary_type, - **{axis + '0': value}) - - if axis == 'x': - x1, x2 = 'y', 'z' - elif axis == 'y': - x1, x2 = 'x', 'z' - else: - x1, x2 = 'x', 'y' - - # Get cylinder class corresponding to given axis - cyl = globals()['{}Cylinder'.format(axis.upper())] - - # Create rectangular region - min_x1 = plane(x1, 'minimum', -width/2 + origin[0]) - max_x1 = plane(x1, 'maximum', width/2 + origin[0]) - min_x2 = plane(x2, 'minimum', -height/2 + origin[1]) - max_x2 = plane(x2, 'maximum', height/2 + origin[1]) - prism = +min_x1 & -max_x1 & +min_x2 & -max_x2 - - # Handle rounded corners if given - if corner_radius > 0.: - args = {'R': corner_radius, 'boundary_type': boundary_type} - - args[x1 + '0'] = origin[0] - width/2 + corner_radius - args[x2 + '0'] = origin[1] - height/2 + corner_radius - x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args) - - args[x1 + '0'] = origin[0] - width/2 + corner_radius - args[x2 + '0'] = origin[1] - height/2 + corner_radius - x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args) - - args[x1 + '0'] = origin[0] - width/2 + corner_radius - args[x2 + '0'] = origin[1] + height/2 - corner_radius - x1_min_x2_max = cyl(name='{} min {} max'.format(x1, x2), **args) - - args[x1 + '0'] = origin[0] + width/2 - corner_radius - args[x2 + '0'] = origin[1] - height/2 + corner_radius - x1_max_x2_min = cyl(name='{} max {} min'.format(x1, x2), **args) - - args[x1 + '0'] = origin[0] + width/2 - corner_radius - args[x2 + '0'] = origin[1] + height/2 - corner_radius - x1_max_x2_max = cyl(name='{} max {} max'.format(x1, x2), **args) - - x1_min = plane(x1, 'min', -width/2 + origin[0] + corner_radius) - x1_max = plane(x1, 'max', width/2 + origin[0] - corner_radius) - x2_min = plane(x2, 'min', -height/2 + origin[1] + corner_radius) - x2_max = plane(x2, 'max', height/2 + origin[1] - corner_radius) - - corners = (+x1_min_x2_min & -x1_min & -x2_min) | \ - (+x1_min_x2_max & -x1_min & +x2_max) | \ - (+x1_max_x2_min & +x1_max & -x2_min) | \ - (+x1_max_x2_max & +x1_max & +x2_max) - - prism = prism & ~corners - - return prism - - -def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), - boundary_type='transmission', corner_radius=0.): - """Create a hexagon region from six surface planes. - - Parameters - ---------- - edge_length : float - Length of a side of the hexagon in cm - orientation : {'x', 'y'} - An 'x' orientation means that two sides of the hexagon are parallel to - the x-axis and a 'y' orientation means that two sides of the hexagon are - parallel to the y-axis. - origin: Iterable of two floats - Origin of the prism. Defaults to (0., 0.). - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'} - Boundary condition that defines the behavior for particles hitting the - surfaces comprising the hexagonal prism (default is 'transmission'). - corner_radius: float - Prism corner radius in units of cm. Defaults to 0. - - Returns - ------- - openmc.Region - The inside of a hexagonal prism - - """ - - l = edge_length - x, y = origin - - if orientation == 'y': - right = XPlane(x0=x + sqrt(3.)/2*l, boundary_type=boundary_type) - left = XPlane(x0=x - sqrt(3.)/2*l, boundary_type=boundary_type) - c = sqrt(3.)/3. - - # y = -x/sqrt(3) + a - upper_right = Plane(A=c, B=1., D=l+x*c+y, boundary_type=boundary_type) - - # y = x/sqrt(3) + a - upper_left = Plane(A=-c, B=1., D=l-x*c+y, boundary_type=boundary_type) - - # y = x/sqrt(3) - a - lower_right = Plane(A=-c, B=1., D=-l-x*c+y, boundary_type=boundary_type) - - # y = -x/sqrt(3) - a - lower_left = Plane(A=c, B=1., D=-l+x*c+y, boundary_type=boundary_type) - - prism = -right & +left & -upper_right & -upper_left & \ - +lower_right & +lower_left - - if boundary_type == 'periodic': - right.periodic_surface = left - upper_right.periodic_surface = lower_left - lower_right.periodic_surface = upper_left - - elif orientation == 'x': - top = YPlane(y0=y + sqrt(3.)/2*l, boundary_type=boundary_type) - bottom = YPlane(y0=y - sqrt(3.)/2*l, boundary_type=boundary_type) - c = sqrt(3.) - - # y = -sqrt(3)*(x - a) - upper_right = Plane(A=c, B=1., D=c*l+x*c+y, boundary_type=boundary_type) - - # y = sqrt(3)*(x + a) - lower_right = Plane(A=-c, B=1., D=-c*l-x*c+y, - boundary_type=boundary_type) - - # y = -sqrt(3)*(x + a) - lower_left = Plane(A=c, B=1., D=-c*l+x*c+y, boundary_type=boundary_type) - - # y = sqrt(3)*(x + a) - upper_left = Plane(A=-c, B=1., D=c*l-x*c+y, boundary_type=boundary_type) - - prism = -top & +bottom & -upper_right & +lower_right & \ - +lower_left & -upper_left - - if boundary_type == 'periodic': - top.periodic_surface = bottom - upper_right.periodic_surface = lower_left - lower_right.periodic_surface = upper_left - - # Handle rounded corners if given - if corner_radius > 0.: - if boundary_type == 'periodic': - raise ValueError('Periodic boundary conditions not permitted when ' - 'rounded corners are used.') - - c = sqrt(3.)/2 - t = l - corner_radius/c - - # Cylinder with corner radius and boundary type pre-applied - cyl1 = partial(ZCylinder, R=corner_radius, boundary_type=boundary_type) - cyl2 = partial(ZCylinder, R=corner_radius/(2*c), - boundary_type=boundary_type) - - if orientation == 'x': - x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t) - x_min_y_max_in = cyl1(name='x min y max in', x0=x+t/2, y0=y-c*t) - x_max_y_min_in = cyl1(name='x max y min in', x0=x-t/2, y0=y+c*t) - x_max_y_max_in = cyl1(name='x max y max in', x0=x+t/2, y0=y+c*t) - x_min_in = cyl1(name='x min in', x0=x-t, y0=y) - x_max_in = cyl1(name='x max in', x0=x+t, y0=y) - - x_min_y_min_out = cyl2(name='x min y min out', x0=x-l/2, y0=y-c*l) - x_min_y_max_out = cyl2(name='x min y max out', x0=x+l/2, y0=y-c*l) - x_max_y_min_out = cyl2(name='x max y min out', x0=x-l/2, y0=y+c*l) - x_max_y_max_out = cyl2(name='x max y max out', x0=x+l/2, y0=y+c*l) - x_min_out = cyl2(name='x min out', x0=x-l, y0=y) - x_max_out = cyl2(name='x max out', x0=x+l, y0=y) - - corners = (+x_min_y_min_in & -x_min_y_min_out | - +x_min_y_max_in & -x_min_y_max_out | - +x_max_y_min_in & -x_max_y_min_out | - +x_max_y_max_in & -x_max_y_max_out | - +x_min_in & -x_min_out | - +x_max_in & -x_max_out) - - elif orientation == 'y': - x_min_y_min_in = cyl1(name='x min y min in', x0=x-c*t, y0=y-t/2) - x_min_y_max_in = cyl1(name='x min y max in', x0=x-c*t, y0=y+t/2) - x_max_y_min_in = cyl1(name='x max y min in', x0=x+c*t, y0=y-t/2) - x_max_y_max_in = cyl1(name='x max y max in', x0=x+c*t, y0=y+t/2) - y_min_in = cyl1(name='y min in', x0=x, y0=y-t) - y_max_in = cyl1(name='y max in', x0=x, y0=y+t) - - x_min_y_min_out = cyl2(name='x min y min out', x0=x-c*l, y0=y-l/2) - x_min_y_max_out = cyl2(name='x min y max out', x0=x-c*l, y0=y+l/2) - x_max_y_min_out = cyl2(name='x max y min out', x0=x+c*l, y0=y-l/2) - x_max_y_max_out = cyl2(name='x max y max out', x0=x+c*l, y0=y+l/2) - y_min_out = cyl2(name='y min out', x0=x, y0=y-l) - y_max_out = cyl2(name='y max out', x0=x, y0=y+l) - - corners = (+x_min_y_min_in & -x_min_y_min_out | - +x_min_y_max_in & -x_min_y_max_out | - +x_max_y_min_in & -x_max_y_min_out | - +x_max_y_max_in & -x_max_y_max_out | - +y_min_in & -y_min_out | - +y_max_in & -y_max_out) - - prism = prism & ~corners - - return prism From 7b5f2bb821b7346ada0ac84c9aed8936074f3009 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 10 Nov 2017 10:35:47 -0600 Subject: [PATCH 12/26] Simple optimizations to Universe._determine_paths --- openmc/universe.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index 26977ac0a2..cb6574f98f 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -553,14 +553,16 @@ class Universe(IDManagerMixin): for cell in self.cells.values(): cell_path = '{}->c{}'.format(univ_path, cell.id) + fill = cell._fill + fill_type = cell.fill_type # If universe-filled, recursively count cells in filling universe - if cell.fill_type == 'universe': - cell.fill._determine_paths(cell_path + '->', instances_only) + if fill_type == 'universe': + fill._determine_paths(cell_path + '->', instances_only) # If lattice-filled, recursively call for all universes in lattice - elif cell.fill_type == 'lattice': - latt = cell.fill + elif fill_type == 'lattice': + latt = fill # Count instances in each universe in the lattice for index in latt._natural_indices: @@ -570,10 +572,10 @@ class Universe(IDManagerMixin): univ._determine_paths(latt_path, instances_only) else: - if cell.fill_type == 'material': - mat = cell.fill - elif cell.fill_type == 'distribmat': - mat = cell.fill[cell._num_instances] + if fill_type == 'material': + mat = fill + elif fill_type == 'distribmat': + mat = fill[cell._num_instances] else: mat = None From 402b09ee3d81054bc97d2fff0d45d9ebc27fe358 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 17 Nov 2017 11:22:48 -0600 Subject: [PATCH 13/26] Add Python bindings for global num_realizations. Fix active rate. --- docs/source/pythonapi/capi.rst | 1 + openmc/capi/tally.py | 9 +++++++-- src/output.F90 | 2 +- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 44a755bbd1..3c6300bb60 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -22,6 +22,7 @@ Functions openmc.capi.keff openmc.capi.load_nuclide openmc.capi.next_batch + openmc.capi.num_realizations openmc.capi.plot_geometry openmc.capi.reset openmc.capi.run diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 6725ee5c92..4a76de41d5 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -13,7 +13,7 @@ from .error import _error_handler, AllocationError, InvalidIDError from .filter import _get_filter -__all__ = ['Tally', 'tallies', 'global_tallies'] +__all__ = ['Tally', 'tallies', 'global_tallies', 'num_realizations'] # Tally functions _dll.openmc_extend_tallies.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)] @@ -89,7 +89,7 @@ def global_tallies(): # Get sum, sum-of-squares, and number of realizations sum_ = array[:, 1] sum_sq = array[:, 2] - n = c_int32.in_dll(_dll, 'n_realizations').value + n = num_realizations() # Determine mean mean = sum_ / n @@ -102,6 +102,11 @@ def global_tallies(): return list(zip(mean, stdev)) +def num_realizations(): + """Number of realizations of global tallies.""" + return c_int32.in_dll(_dll, 'n_realizations').value + + class Tally(_FortranObjectWithID): """Tally stored internally. diff --git a/src/output.F90 b/src/output.F90 index 3d02489f78..c200775840 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -570,7 +570,7 @@ contains write(ou,100) "Total time elapsed", time_total % elapsed ! Calculate particle rate in active/inactive batches - n_active = n_batches - n_inactive + n_active = current_batch - n_inactive if (restart_run) then if (restart_batch < n_inactive) then speed_inactive = real(n_particles * (n_inactive - restart_batch) * & From ee8b7edd9c0a1a70dbcecd4d0511e3baf9ad01aa Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Nov 2017 14:17:45 -0600 Subject: [PATCH 14/26] Display message when writing summary file --- src/input_xml.F90 | 10 +++++----- src/summary.F90 | 5 ++++- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 766dbfa9f2..c76d66df7f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1920,10 +1920,7 @@ contains type(XMLDocument) :: doc type(XMLNode) :: root - ! Display output message - call write_message("Reading materials XML file...", 5) - - ! Check is materials.xml exists + ! Check if materials.xml exists filename = trim(path_input) // "materials.xml" inquire(FILE=filename, EXIST=file_exists) if (.not. file_exists) then @@ -2057,7 +2054,10 @@ contains type(XMLNode), allocatable :: node_macro_list(:) type(XMLNode), allocatable :: node_sab_list(:) - ! Check is materials.xml exists + ! Display output message + call write_message("Reading materials XML file...", 5) + + ! Check if materials.xml exists filename = trim(path_input) // "materials.xml" inquire(FILE=filename, EXIST=file_exists) if (.not. file_exists) then diff --git a/src/summary.F90 b/src/summary.F90 index e271fbaba6..248b48e865 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -11,7 +11,7 @@ module summary use message_passing use mgxs_header, only: nuclides_MG use nuclide_header - use output, only: time_stamp + use output, only: time_stamp, write_message use settings, only: run_CE use surface_header use string, only: to_str @@ -33,6 +33,9 @@ contains integer(HID_T) :: file_id + ! Display output message + call write_message("Writing summary.h5 file...", 5) + ! Create a new file using default properties. file_id = file_create("summary.h5") From 726d1a268f9419db35f3aaf21b2cafa976a98c4a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Nov 2017 23:04:12 -0600 Subject: [PATCH 15/26] Sort materials list before creating elements --- openmc/geometry.py | 2 +- openmc/material.py | 9 ++++----- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/openmc/geometry.py b/openmc/geometry.py index 0b02210b7e..88c9f18348 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -87,7 +87,7 @@ class Geometry(object): # Write the XML Tree to the geometry.xml file tree = ET.ElementTree(root_element) - tree.write(path, xml_declaration=True, encoding='utf-8', method="xml") + tree.write(path, xml_declaration=True, encoding='utf-8') def find(self, point): """Find cells/universes/lattices which contain a given point diff --git a/openmc/material.py b/openmc/material.py index 03ef7fbc42..3800d89c67 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -10,7 +10,7 @@ import numpy as np import openmc import openmc.data import openmc.checkvalue as cv -from openmc.clean_xml import sort_xml_elements, clean_xml_indentation +from openmc.clean_xml import clean_xml_indentation from .mixin import IDManagerMixin @@ -1129,7 +1129,7 @@ class Materials(cv.CheckedList): material.make_isotropic_in_lab() def _create_material_subelements(self, root_element): - for material in self: + for material in sorted(self, key=lambda x: x.id): root_element.append(material.to_xml_element(self.cross_sections)) def _create_cross_sections_subelement(self, root_element): @@ -1153,14 +1153,13 @@ class Materials(cv.CheckedList): """ root_element = ET.Element("materials") - self._create_material_subelements(root_element) self._create_cross_sections_subelement(root_element) self._create_multipole_library_subelement(root_element) + self._create_material_subelements(root_element) # Clean the indentation in the file to be user-readable - sort_xml_elements(root_element) clean_xml_indentation(root_element) # Write the XML Tree to the materials.xml file tree = ET.ElementTree(root_element) - tree.write(path, xml_declaration=True, encoding='utf-8', method="xml") + tree.write(path, xml_declaration=True, encoding='utf-8') From 38550f3bf87f6400ea577fcfab51545911e21e08 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Nov 2017 06:34:15 -0600 Subject: [PATCH 16/26] Use sorted() to sort XML elements in geometry.xml --- openmc/clean_xml.py | 67 --------------------------------------------- openmc/geometry.py | 7 +++-- 2 files changed, 5 insertions(+), 69 deletions(-) diff --git a/openmc/clean_xml.py b/openmc/clean_xml.py index 2497d0d9f2..d1002070b3 100644 --- a/openmc/clean_xml.py +++ b/openmc/clean_xml.py @@ -1,70 +1,3 @@ -def sort_xml_elements(tree): - - # Retrieve all children of the root XML node in the tree - elements = list(tree) - - # Initialize empty lists for the sorted and comment elements - sorted_elements = [] - - # Initialize an empty set of tags (e.g., Surface, Cell, and Lattice) - tags = set() - - # Find the unique tags in the tree - for element in elements: - tags.add(element.tag) - - # Initialize an empty list for the comment elements - comment_elements = [] - - # Find the comment elements and record their ordering within the - # tree using a precedence with respect to the subsequent nodes - for index, element in enumerate(elements): - next_element = None - - if 'Comment' in str(element.tag): - - if index < len(elements)-1: - next_element = elements[index+1] - - comment_elements.append((element, next_element)) - - # Now iterate over all tags and order the elements within each tag - for tag in sorted(list(tags)): - - # Retrieve all of the elements for this tag - try: - tagged_elements = tree.findall(tag) - except: - continue - - # Initialize an empty list of tuples to sort (id, element) - tagged_data = [] - - # Retrieve the IDs for each of the elements - for element in tagged_elements: - key = element.get('id') - - # If this element has an "ID" tag, append it to the list to sort - if key is not None: - tagged_data.append((int(key), element)) - - # Sort the elements according to the IDs for this tag - tagged_data.sort() - sorted_elements.extend(list(item[-1] for item in tagged_data)) - - # Add the comment elements while preserving the original precedence - for element, next_element in comment_elements: - index = sorted_elements.index(next_element) - sorted_elements.insert(index, element) - - # Remove all of the sorted elements from the tree - for element in sorted_elements: - tree.remove(element) - - # Add the sorted elements back to the tree in the proper order - tree.extend(sorted_elements) - - def clean_xml_indentation(element, level=0, spaces_per_level=2): """ copy and paste from http://effbot.org/zone/elementent-lib.htm#prettyprint diff --git a/openmc/geometry.py b/openmc/geometry.py index 88c9f18348..0738ae615a 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -5,7 +5,7 @@ from xml.etree import ElementTree as ET from six import string_types import openmc -from openmc.clean_xml import sort_xml_elements, clean_xml_indentation +from openmc.clean_xml import clean_xml_indentation from openmc.checkvalue import check_type @@ -81,8 +81,11 @@ class Geometry(object): root_element = ET.Element("geometry") self.root_universe.create_xml_subelement(root_element) + # Sort the elements in the file + root_element[:] = sorted(root_element, key=lambda x: ( + x.tag, int(x.get('id')))) + # Clean the indentation in the file to be user-readable - sort_xml_elements(root_element) clean_xml_indentation(root_element) # Write the XML Tree to the geometry.xml file From a9b831cc912d0858032b424a2b6e9a7cde8d976d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Nov 2017 07:19:18 -0600 Subject: [PATCH 17/26] Add one publication to list --- docs/source/publications.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/source/publications.rst b/docs/source/publications.rst index f47ac704c1..7578fee6b1 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -171,6 +171,12 @@ Miscellaneous Multi-group Cross Section Generation ------------------------------------ +- Zhaoyuan Liu, Kord Smith, Benoit Forget, and Javier Ortensi, "`Cumulative + migration method for computing rigorous diffusion coefficients and transport + cross sections from Monte Carlo + `_," *Ann. Nucl. Energy*, + **112**, 507-516 (2018). + - Gang Yang, Tongkyu Park, and Won Sik Yang, "Effects of Fuel Salt Velocity Field on Neutronics Performances in Molten Salt Reactors with Open Flow Channels," *Trans. Am. Nucl. Soc.*, **117**, 1339-1342 (2017). From ad0283d0166e5f9779dcfc1c65d8e4f2a7236e52 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Dec 2017 10:58:52 -0600 Subject: [PATCH 18/26] Account for 'ftn' as an MPI wrapper since it is common on Cray systems --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index dfd9221cec..791e31c44a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,7 +46,7 @@ add_definitions(-DMAX_COORD=${maxcoord}) #=============================================================================== set(MPI_ENABLED FALSE) -if($ENV{FC} MATCHES "mpi[^/]*$") +if($ENV{FC} MATCHES "(mpi[^/]*|ftn)$") message("-- Detected MPI wrapper: $ENV{FC}") add_definitions(-DMPI) set(MPI_ENABLED TRUE) From 4fbd9d75c22e4d6f22f0a6614f379a74a3fde169 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Dec 2017 11:21:30 -0600 Subject: [PATCH 19/26] Allow openmc.capi.statepoint_write to pass filename --- openmc/capi/core.py | 20 ++++++++++++++++---- src/state_point.F90 | 25 ++++++++++++++++--------- 2 files changed, 32 insertions(+), 13 deletions(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 1f7225f642..1e3b8a5102 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -1,5 +1,6 @@ from contextlib import contextmanager -from ctypes import CDLL, c_int, c_int32, c_int64, c_double, POINTER, Structure +from ctypes import (CDLL, c_int, c_int32, c_int64, c_double, c_char_p, + POINTER, Structure) from warnings import warn import numpy as np @@ -38,6 +39,7 @@ _dll.openmc_source_bank.restype = c_int _dll.openmc_source_bank.errcheck = _error_handler _dll.openmc_simulation_init.restype = None _dll.openmc_simulation_finalize.restype = None +_dll.openmc_statepoint_write.argtypes = [POINTER(c_char_p)] _dll.openmc_statepoint_write.restype = None @@ -217,9 +219,19 @@ def source_bank(): return as_array(ptr, (n.value,)).view(bank_dtype) -def statepoint_write(): - """Write a statepoint.""" - _dll.openmc_statepoint_write() +def statepoint_write(filename=None): + """Write a statepoint file. + + Parameters + ---------- + filename : str or None + Path to the statepoint to write. If None is passed, a default name that + contains the current batch will be written. + + """ + if filename is not None: + filename = c_char_p(filename.encode()) + _dll.openmc_statepoint_write(filename) @contextmanager diff --git a/src/state_point.F90 b/src/state_point.F90 index 0383150ac4..6d647b72e6 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -29,7 +29,7 @@ module state_point use random_lcg, only: seed use settings use simulation_header - use string, only: to_str, count_digits, zero_padded + use string, only: to_str, count_digits, zero_padded, to_f_string use tally_header use tally_filter_header use tally_derivative_header, only: tally_derivs @@ -43,7 +43,8 @@ contains ! OPENMC_STATEPOINT_WRITE writes an HDF5 statepoint file to disk !=============================================================================== - subroutine openmc_statepoint_write() bind(C) + subroutine openmc_statepoint_write(filename) bind(C) + type(C_PTR), intent(in), optional :: filename integer :: i, j, k integer :: i_xs @@ -57,19 +58,25 @@ contains integer(C_INT) :: err real(C_DOUBLE) :: k_combined(2) character(MAX_WORD_LEN), allocatable :: str_array(:) - character(MAX_FILE_LEN) :: filename + character(C_CHAR), pointer :: string(:) + character(len=:, kind=C_CHAR), allocatable :: filename_ - ! Set filename for state point - filename = trim(path_output) // 'statepoint.' // & - & zero_padded(current_batch, count_digits(n_max_batches)) - filename = trim(filename) // '.h5' + if (present(filename)) then + call c_f_pointer(filename, string, [MAX_FILE_LEN]) + filename_ = to_f_string(string) + else + ! Set filename for state point + filename_ = trim(path_output) // 'statepoint.' // & + & zero_padded(current_batch, count_digits(n_max_batches)) + filename_ = trim(filename_) // '.h5' + end if ! Write message - call write_message("Creating state point " // trim(filename) // "...", 5) + call write_message("Creating state point " // trim(filename_) // "...", 5) if (master) then ! Create statepoint file - file_id = file_create(filename) + file_id = file_create(filename_) ! Write file type call write_attribute(file_id, "filetype", "statepoint") From 5116916d7a21458e9e297912c0c0069a36eee1e6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 4 Dec 2017 17:57:27 -0600 Subject: [PATCH 20/26] Fix capi.Cell.__new__ --- openmc/capi/cell.py | 43 +++++++++++++++++++++++++++++++++-------- src/api.F90 | 1 + src/geometry_header.F90 | 32 ++++++++++++++++++++++++++++++ 3 files changed, 68 insertions(+), 8 deletions(-) diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index fe33825894..e14df3f4af 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -7,12 +7,15 @@ from numpy.ctypeslib import as_array from . import _dll from .core import _FortranObjectWithID -from .error import _error_handler +from .error import _error_handler, AllocationError, InvalidIDError from .material import Material __all__ = ['Cell', 'cells'] # Cell functions +_dll.openmc_extend_cells.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)] +_dll.openmc_extend_cells.restype = c_int +_dll.openmc_extend_cells.errcheck = _error_handler _dll.openmc_cell_get_id.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_cell_get_id.restype = c_int _dll.openmc_cell_get_id.errcheck = _error_handler @@ -51,11 +54,35 @@ class Cell(_FortranObjectWithID): """ __instances = WeakValueDictionary() - def __new__(cls, *args): - if args not in cls.__instances: - instance = super(Cell, self).__new__(cls) - cls.__instances[args] = instance - return cls.__instances[args] + def __new__(cls, uid=None, new=True, index=None): + mapping = cells + if index is None: + if new: + # Determine ID to assign + if uid is None: + try: + uid = max(mapping) + 1 + except ValueError: + uid = 1 + else: + if uid in mapping: + raise AllocationError('A cell with ID={} has already ' + 'been allocated.'.format(uid)) + + index = c_int32() + _dll.openmc_extend_cells(1, index, None) + index = index.value + else: + index = mapping[uid]._index + + if index not in cls.__instances: + instance = super(Cell, cls).__new__(cls) + instance._index = index + if uid is not None: + instance.id = uid + cls.__instances[index] = instance + + return cls.__instances[index] @property def id(self): @@ -113,11 +140,11 @@ class _CellMapping(Mapping): except (AllocationError, InvalidIDError) as e: # __contains__ expects a KeyError to work correctly raise KeyError(str(e)) - return Cell(index.value) + return Cell(index=index.value) def __iter__(self): for i in range(len(self)): - yield Cell(i + 1).id + yield Cell(index=i + 1).id def __len__(self): return c_int32.in_dll(_dll, 'n_cells').value diff --git a/src/api.F90 b/src/api.F90 index 933faaa09e..d0231c9a8f 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -40,6 +40,7 @@ module openmc_api public :: openmc_energy_filter_get_bins public :: openmc_energy_filter_set_bins public :: openmc_extend_filters + public :: openmc_extend_cells public :: openmc_extend_materials public :: openmc_extend_tallies public :: openmc_filter_get_id diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index dd2cad7f15..63b1bbe954 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -427,6 +427,38 @@ contains ! C API FUNCTIONS !=============================================================================== + function openmc_extend_cells(n, index_start, index_end) result(err) bind(C) + ! Extend the cells array by n elements + integer(C_INT32_T), value, intent(in) :: n + integer(C_INT32_T), optional, intent(out) :: index_start + integer(C_INT32_T), optional, intent(out) :: index_end + integer(C_INT) :: err + + type(Cell), allocatable :: temp(:) ! temporary cells array + + if (n_cells == 0) then + ! Allocate cells array + allocate(cells(n)) + else + ! Allocate cells array with increased size + allocate(temp(n_cells + n)) + + ! Copy original cells to temporary array + temp(1:n_cells) = cells + + ! Move allocation from temporary array + call move_alloc(FROM=temp, TO=cells) + end if + + ! Return indices in cells array + if (present(index_start)) index_start = n_cells + 1 + if (present(index_end)) index_end = n_cells + n + n_cells = n_cells + n + + err = 0 + end function openmc_extend_cells + + function openmc_get_cell_index(id, index) result(err) bind(C) ! Return the index in the cells array of a cell with a given ID integer(C_INT32_T), value :: id From 7b56bcead68b27608a8bca3d2e76b0ba04d78485 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Dec 2017 09:52:05 +0700 Subject: [PATCH 21/26] Add unit tests for openmc.capi --- openmc/__init__.py | 1 + openmc/capi/cell.py | 7 + openmc/capi/core.py | 13 +- openmc/capi/error.py | 1 + src/api.F90 | 1 + src/geometry_header.F90 | 17 +++ tests/unit_tests/test_capi.py | 241 ++++++++++++++++++++++++++++++++++ 7 files changed, 275 insertions(+), 6 deletions(-) create mode 100644 tests/unit_tests/test_capi.py diff --git a/openmc/__init__.py b/openmc/__init__.py index 8fb9bcf37c..13dac05e7e 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -27,6 +27,7 @@ from openmc.particle_restart import * from openmc.mixin import * from openmc.plotter import * from openmc.search import * +from . import examples # Import a few convencience functions that used to be here from openmc.model import get_rectangular_prism, get_hexagonal_prism diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index e14df3f4af..f13f64a045 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -25,6 +25,9 @@ _dll.openmc_cell_set_fill.argtypes = [ c_int32, c_int, c_int32, POINTER(c_int32)] _dll.openmc_cell_set_fill.restype = c_int _dll.openmc_cell_set_fill.errcheck = _error_handler +_dll.openmc_cell_set_id.argtypes = [c_int32, c_int32] +_dll.openmc_cell_set_id.restype = c_int +_dll.openmc_cell_set_id.errcheck = _error_handler _dll.openmc_cell_set_temperature.argtypes = [ c_int32, c_double, POINTER(c_int32)] _dll.openmc_cell_set_temperature.restype = c_int @@ -90,6 +93,10 @@ class Cell(_FortranObjectWithID): _dll.openmc_cell_get_id(self._index, cell_id) return cell_id.value + @id.setter + def id(self, cell_id): + _dll.openmc_cell_set_id(self._index, cell_id) + @property def fill(self): fill_type = c_int() diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 1e3b8a5102..c25146fb5b 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -8,6 +8,7 @@ from numpy.ctypeslib import as_array from . import _dll from .error import _error_handler +import openmc.capi class _Bank(Structure): @@ -63,8 +64,8 @@ def find_cell(xyz): Returns ------- - int - ID of the cell. + openmc.capi.Cell + Cell containing the point int If the cell at the given point is repeated in the geometry, this indicates which instance it is, i.e., 0 would be the first instance. @@ -73,7 +74,7 @@ def find_cell(xyz): uid = c_int32() instance = c_int32() _dll.openmc_find((c_double*3)(*xyz), 1, uid, instance) - return uid.value, instance.value + return openmc.capi.cells[uid.value], instance.value def find_material(xyz): @@ -86,14 +87,14 @@ def find_material(xyz): Returns ------- - int or None - ID of the material or None is no material is found + openmc.capi.Material or None + Material containing the point, or None is no material is found """ uid = c_int32() instance = c_int32() _dll.openmc_find((c_double*3)(*xyz), 2, uid, instance) - return uid.value if uid != 0 else None + return openmc.capi.materials[uid.value] if uid != 0 else None def hard_reset(): diff --git a/openmc/capi/error.py b/openmc/capi/error.py index 98d43ae462..c312987bc2 100644 --- a/openmc/capi/error.py +++ b/openmc/capi/error.py @@ -1,4 +1,5 @@ from ctypes import c_int, c_char +from warnings import warn from . import _dll diff --git a/src/api.F90 b/src/api.F90 index d0231c9a8f..4538bbf3ee 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -36,6 +36,7 @@ module openmc_api public :: openmc_cell_get_id public :: openmc_cell_get_fill public :: openmc_cell_set_fill + public :: openmc_cell_set_id public :: openmc_cell_set_temperature public :: openmc_energy_filter_get_bins public :: openmc_energy_filter_set_bins diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 63b1bbe954..b827c055c1 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -570,6 +570,23 @@ contains end function openmc_cell_set_fill + function openmc_cell_set_id(index, id) result(err) bind(C) + ! Set the ID of a cell + integer(C_INT32_T), value, intent(in) :: index + integer(C_INT32_T), value, intent(in) :: id + integer(C_INT) :: err + + if (index >= 1 .and. index <= n_cells) then + cells(index) % id = id + call cell_dict % set(id, index) + err = 0 + else + err = E_OUT_OF_BOUNDS + call set_errmsg("Index in cells array is out of bounds.") + end if + end function openmc_cell_set_id + + function openmc_cell_set_temperature(index, T, instance) result(err) bind(C) ! Set the temperature of a cell integer(C_INT32_T), value, intent(in) :: index ! index in cells diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py new file mode 100644 index 0000000000..08387fe185 --- /dev/null +++ b/tests/unit_tests/test_capi.py @@ -0,0 +1,241 @@ +#!/usr/bin/env python + +from collections.abc import Mapping +import os + +import numpy as np +import pytest +import openmc +import openmc.capi + + +@pytest.fixture(scope='module') +def pincell_model(): + """Set up a model to test with and delete files when done""" + pincell = openmc.examples.pwr_pin_cell() + + # Add a tally + filter1 = openmc.MaterialFilter(pincell.materials) + filter2 = openmc.EnergyFilter([0.0, 1.0, 1.0e3, 20.0e6]) + mat_tally = openmc.Tally() + mat_tally.filters = [filter1, filter2] + mat_tally.nuclides = ['U235', 'U238'] + mat_tally.scores = ['total', 'elastic', '(n,gamma)'] + pincell.tallies.append(mat_tally) + + # Write XML files + pincell.export_to_xml() + + yield + + # Delete generated files + files = ['geometry.xml', 'materials.xml', 'settings.xml', 'tallies.xml', + 'statepoint.10.h5', 'summary.h5', 'test_sp.h5'] + for f in files: + if os.path.exists(f): + os.remove(f) + + +@pytest.fixture(scope='module') +def capi_init(pincell_model): + openmc.capi.init() + yield + openmc.capi.finalize() + + +@pytest.fixture(scope='module') +def capi_run(capi_init): + openmc.capi.run() + + +def test_cell_mapping(capi_init): + cells = openmc.capi.cells + assert isinstance(cells, Mapping) + assert len(cells) == 3 + for cell_id, cell in cells.items(): + assert isinstance(cell, openmc.capi.Cell) + assert cell_id == cell.id + + +def test_cell(capi_init): + cell = openmc.capi.cells[1] + assert isinstance(cell.fill, openmc.capi.Material) + cell.fill = openmc.capi.materials[1] + assert str(cell) == 'Cell[1]' + + +def test_new_cell(capi_init): + with pytest.raises(openmc.capi.AllocationError): + openmc.capi.Cell(1) + new_cell = openmc.capi.Cell() + new_cell_with_id = openmc.capi.Cell(10) + assert len(openmc.capi.cells) == 5 + + +def test_material_mapping(capi_init): + mats = openmc.capi.materials + assert isinstance(mats, Mapping) + assert len(mats) == 3 + for mat_id, mat in mats.items(): + assert isinstance(mat, openmc.capi.Material) + assert mat_id == mat.id + + +def test_material(capi_init): + m = openmc.capi.materials[3] + assert m.nuclides == ['H1', 'O16', 'B10', 'B11'] + + old_dens = m.densities + test_dens = [1.0e-1, 2.0e-1, 2.5e-1, 1.0e-3] + m.set_densities(m.nuclides, test_dens) + assert m.densities == pytest.approx(test_dens) + + rho = 2.25e-2 + m.set_density(rho) + assert sum(m.densities) == pytest.approx(rho) + + +def test_new_material(capi_init): + with pytest.raises(openmc.capi.AllocationError): + openmc.capi.Material(1) + new_mat = openmc.capi.Material() + new_mat_with_id = openmc.capi.Material(10) + assert len(openmc.capi.materials) == 5 + + +def test_nuclide_mapping(capi_init): + nucs = openmc.capi.nuclides + assert isinstance(nucs, Mapping) + assert len(nucs) == 12 + for name, nuc in nucs.items(): + assert isinstance(nuc, openmc.capi.Nuclide) + assert name == nuc.name + + +def test_load_nuclide(capi_init): + openmc.capi.load_nuclide('Pu239') + with pytest.raises(openmc.capi.DataError): + openmc.capi.load_nuclide('Pu3') + + +def test_settings(capi_init): + settings = openmc.capi.settings + assert settings.batches == 10 + settings.batches = 10 + assert settings.inactive == 5 + assert settings.generations_per_batch == 1 + assert settings.particles == 100 + assert settings.seed == 1 + settings.seed = 11 + + assert settings.run_mode == 'eigenvalue' + settings.run_mode = 'volume' + settings.run_mode = 'eigenvalue' + + +def test_tally_mapping(capi_init): + tallies = openmc.capi.tallies + assert isinstance(tallies, Mapping) + assert len(tallies) == 1 + for tally_id, tally in tallies.items(): + assert isinstance(tally, openmc.capi.Tally) + assert tally_id == tally.id + + +def test_tally(capi_init): + t = openmc.capi.tallies[1] + t.id = 1 + assert len(t.filters) == 2 + assert isinstance(t.filters[0], openmc.capi.MaterialFilter) + assert isinstance(t.filters[1], openmc.capi.EnergyFilter) + + # Create new filter and replace existing + with pytest.raises(openmc.capi.AllocationError): + openmc.capi.MaterialFilter(uid=1) + mats = openmc.capi.materials + f = openmc.capi.MaterialFilter([mats[2], mats[1]]) + t.filters = [f] + assert t.filters == [f] + + assert t.nuclides == ['U235', 'U238'] + with pytest.raises(openmc.capi.DataError): + t.nuclides = ['Zr2'] + t.nuclides = ['U234', 'Zr90'] + assert t.nuclides == ['U234', 'Zr90'] + + assert t.scores == ['total', '(n,elastic)', '(n,gamma)'] + new_scores = ['scatter', 'fission', 'nu-fission', '(n,2n)'] + t.scores = new_scores + assert t.scores == new_scores + + +def test_new_tally(capi_init): + with pytest.raises(openmc.capi.AllocationError): + openmc.capi.Material(1) + new_tally = openmc.capi.Tally() + new_tally.scores = ['flux'] + new_tally_with_id = openmc.capi.Tally(10) + new_tally_with_id.scores = ['flux'] + assert len(openmc.capi.tallies) == 3 + + +def test_tally_results(capi_run): + t = openmc.capi.tallies[1] + assert t.num_realizations == 5 + assert np.all(t.mean >= 0) + nonzero = (t.mean > 0.0) + assert np.all(t.std_dev[nonzero] >= 0) + assert np.all(t.ci_width()[nonzero] >= 1.95*t.std_dev[nonzero]) + + +def test_global_tallies(capi_run): + assert openmc.capi.num_realizations() == 5 + gt = openmc.capi.global_tallies() + for mean, std_dev in gt: + assert mean >= 0 + + +def test_statepoint(capi_run): + openmc.capi.statepoint_write('test_sp.h5') + assert os.path.exists('test_sp.h5') + + +def test_source_bank(capi_run): + source = openmc.capi.source_bank() + assert np.all(source['E'] > 0.0) + assert np.all(source['wgt'] == 1.0) + + +def test_keff(capi_run): + mean, std_dev = openmc.capi.keff() + assert 0.0 < mean < 2.5 + assert std_dev > 0.0 + + +def test_by_batch(capi_run): + openmc.capi.hard_reset() + openmc.capi.simulation_init() + for _ in openmc.capi.iter_batches(): + pass + assert openmc.capi.num_realizations() == 5 + + for i in range(3): + openmc.capi.next_batch() + assert openmc.capi.num_realizations() == 8 + openmc.capi.simulation_finalize() + + +def test_find_cell(capi_init): + cell, instance = openmc.capi.find_cell((0., 0., 0.)) + assert cell is openmc.capi.cells[1] + cell, instance = openmc.capi.find_cell((0.4, 0., 0.)) + assert cell is openmc.capi.cells[2] + with pytest.raises(openmc.capi.GeometryError): + openmc.capi.find_cell((100., 100., 100.)) + + +def test_find_material(capi_init): + mat = openmc.capi.find_material((0., 0., 0.)) + assert mat is openmc.capi.materials[1] + mat = openmc.capi.find_material((0.4, 0., 0.)) + assert mat is openmc.capi.materials[2] From 74a28519d01e079128e552ce749def9b5ecec525 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Dec 2017 10:04:14 +0700 Subject: [PATCH 22/26] Fix import for Python 2.7 --- tests/unit_tests/test_capi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 08387fe185..cde3d87245 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -1,6 +1,6 @@ #!/usr/bin/env python -from collections.abc import Mapping +from collections import Mapping import os import numpy as np From 147bbefbe7f99fe5a3aa02cfd598de65cb400a04 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Dec 2017 10:36:28 +0700 Subject: [PATCH 23/26] Run pytest verbosely --- .travis.yml | 2 +- tests/unit_tests/test_element_wo.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.travis.yml b/.travis.yml index 19bb96040b..9ae0b2f5ba 100644 --- a/.travis.yml +++ b/.travis.yml @@ -62,6 +62,6 @@ script: ./check_source.py; else ./run_tests.py -C $OPENMC_CONFIG -j 2 && - pytest --cov=../openmc unit_tests/; + pytest --cov=../openmc -v unit_tests/; fi - cd .. diff --git a/tests/unit_tests/test_element_wo.py b/tests/unit_tests/test_element_wo.py index 7b9487f73d..2431f9e0e7 100644 --- a/tests/unit_tests/test_element_wo.py +++ b/tests/unit_tests/test_element_wo.py @@ -3,7 +3,7 @@ import os import sys -import numpy as np +import pytest from openmc import Material from openmc.data import NATURAL_ABUNDANCE, atomic_mass @@ -31,11 +31,11 @@ def test_element_wo(): if nuc in ('H1', 'H2'): val = 2 * NATURAL_ABUNDANCE[nuc] * atomic_mass(nuc) / water_am - assert np.isclose(densities[nuc][1], val, rtol=1.e-8) + assert densities[nuc][1] == pytest.approx(val) if nuc == 'O16': val = (NATURAL_ABUNDANCE[nuc] + NATURAL_ABUNDANCE['O18']) \ * atomic_mass(nuc) / water_am - assert np.isclose(densities[nuc][1], val, rtol=1.e-8) + assert densities[nuc][1] == pytest.approx(val) if nuc == 'O17': val = NATURAL_ABUNDANCE[nuc] * atomic_mass(nuc) / water_am - assert np.isclose(densities[nuc][1], val, rtol=1.e-8) + assert densities[nuc][1] == pytest.approx(val) From 9fa391aa666e6e0ec28def517a10392c0388e0d4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Dec 2017 11:52:30 +0700 Subject: [PATCH 24/26] Make sure latest version of pytest is used --- .travis.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.travis.yml b/.travis.yml index 9ae0b2f5ba..49b6ed2d12 100644 --- a/.travis.yml +++ b/.travis.yml @@ -41,6 +41,7 @@ before_install: install: - if [[ $OPENMC_CONFIG != "check_source" ]]; then pip install numpy cython; + pip install --upgrade pytest; pip install -e .[test]; fi From 8e1b8d6264266f72573c195e8370dc32cb8a3893 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Dec 2017 11:05:46 +0700 Subject: [PATCH 25/26] Make capi.keff() usable during inactive and active batches --- openmc/capi/core.py | 14 +++++++++++--- openmc/capi/error.py | 20 ++++++++++---------- openmc/capi/tally.py | 13 +++++++++---- src/simulation_header.F90 | 6 ++++-- 4 files changed, 34 insertions(+), 19 deletions(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index c25146fb5b..ced9390435 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -166,9 +166,17 @@ def keff(): Mean k-eigenvalue and standard deviation of the mean """ - k = (c_double*2)() - _dll.openmc_get_keff(k) - return tuple(k) + n = openmc.capi.num_realizations() + if n > 3: + # Use the combined estimator if there are enough realizations + k = (c_double*2)() + _dll.openmc_get_keff(k) + return tuple(k) + else: + # Otherwise, return the tracklength estimator + mean = c_double.in_dll(_dll, 'keff').value + std_dev = c_double.in_dll(_dll, 'keff_std').value if n > 1 else np.inf + return (mean, std_dev) def next_batch(): diff --git a/openmc/capi/error.py b/openmc/capi/error.py index c312987bc2..a11d6ea87d 100644 --- a/openmc/capi/error.py +++ b/openmc/capi/error.py @@ -4,39 +4,39 @@ from warnings import warn from . import _dll -class Error(Exception): +class OpenMCError(Exception): """Root exception class for OpenMC.""" -class GeometryError(Error): +class GeometryError(OpenMCError): """Geometry-related error""" -class InvalidIDError(Error): +class InvalidIDError(OpenMCError): """Use of an ID that is invalid.""" -class AllocationError(Error): +class AllocationError(OpenMCError): """Error related to memory allocation.""" -class OutOfBoundsError(Error): +class OutOfBoundsError(OpenMCError): """Index in array out of bounds.""" -class DataError(Error): +class DataError(OpenMCError): """Error relating to nuclear data.""" -class PhysicsError(Error): +class PhysicsError(OpenMCError): """Error relating to performing physics.""" -class InvalidArgumentError(Error): +class InvalidArgumentError(OpenMCError): """Argument passed was invalid.""" -class InvalidTypeError(Error): +class InvalidTypeError(OpenMCError): """Tried to perform an operation on the wrong type.""" @@ -71,4 +71,4 @@ def _error_handler(err, func, args): elif err == errcode('e_warning'): warn(msg) elif err < 0: - raise Exception("Unknown error encountered (code {}).".format(err)) + raise OpenMCError("Unknown error encountered (code {}).".format(err)) diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 4a76de41d5..751f0e6031 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -92,12 +92,17 @@ def global_tallies(): n = num_realizations() # Determine mean - mean = sum_ / n + if n > 0: + mean = sum_ / n + else: + mean = sum_.copy() # Determine standard deviation nonzero = np.abs(mean) > 0 - stdev = np.zeros_like(mean) - stdev[nonzero] = np.sqrt((sum_sq[nonzero]/n - mean[nonzero]**2)/(n - 1)) + stdev = np.empty_like(mean) + stdev.fill(np.inf) + if n > 1: + stdev[nonzero] = np.sqrt((sum_sq[nonzero]/n - mean[nonzero]**2)/(n - 1)) return list(zip(mean, stdev)) @@ -265,7 +270,7 @@ class Tally(_FortranObjectWithID): def std_dev(self): results = self.results std_dev = np.empty(results.shape[:2]) - std_dev.fill(np.nan) + std_dev.fill(np.inf) n = self.num_realizations if n > 1: diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 9719cc72e3..62c66e5f98 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -1,5 +1,7 @@ module simulation_header + use, intrinsic :: ISO_C_BINDING + use bank_header use constants use settings, only: gen_per_batch @@ -33,8 +35,8 @@ module simulation_header ! Temporary k-effective values type(VectorReal) :: k_generation ! single-generation estimates of k - real(8) :: keff = ONE ! average k over active batches - real(8) :: keff_std ! standard deviation of average k + real(C_DOUBLE), bind(C) :: keff = ONE ! average k over active batches + real(C_DOUBLE), bind(C) :: keff_std ! standard deviation of average k real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength From 08bf32e1ae3ec0b8384dc809a77250bf4d9097b2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Dec 2017 11:38:28 +0700 Subject: [PATCH 26/26] Throw exception if next_batch() is called before simulation is initialized --- openmc/capi/core.py | 10 +++++++--- src/simulation.F90 | 18 ++++++++++++++++++ src/simulation_header.F90 | 6 +++++- tests/unit_tests/test_capi.py | 17 ++++++++++------- 4 files changed, 40 insertions(+), 11 deletions(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index ced9390435..415fa2e936 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -7,7 +7,7 @@ import numpy as np from numpy.ctypeslib import as_array from . import _dll -from .error import _error_handler +from .error import _error_handler, AllocationError import openmc.capi @@ -147,7 +147,7 @@ def iter_batches(): """ while True: # Run next batch - retval = _dll.openmc_next_batch() + retval = next_batch() # Provide opportunity for user to perform action between batches yield @@ -181,7 +181,11 @@ def keff(): def next_batch(): """Run next batch.""" - return _dll.openmc_next_batch() + retval = _dll.openmc_next_batch() + if retval == -3: + raise AllocationError('Simulation has not been initialized. You must call ' + 'openmc.capi.simulation_init() first.') + return retval def plot_geometry(): diff --git a/src/simulation.F90 b/src/simulation.F90 index e3eba0d512..5eacbd40b8 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -73,6 +73,12 @@ contains type(Particle) :: p integer(8) :: i_work + ! Make sure simulation has been initialized + if (.not. simulation_initialized) then + retval = -3 + return + end if + call initialize_batch() ! Handle restart runs @@ -396,6 +402,9 @@ contains subroutine openmc_simulation_init() bind(C) + ! Skip if simulation has already been initialized + if (simulation_initialized) return + ! Set up tally procedure pointers call init_tally_routines() @@ -438,6 +447,9 @@ contains ! Reset current batch current_batch = 0 + ! Set flag indicating initialization is done + simulation_initialized = .true. + end subroutine openmc_simulation_init !=============================================================================== @@ -455,6 +467,9 @@ contains real(8) :: tempr(3) ! temporary array for communication #endif + ! Skip if simulation was never run + if (.not. simulation_initialized) return + ! Stop active batch timer call time_active % stop() @@ -511,6 +526,9 @@ contains if (check_overlaps) call print_overlap_check() end if + ! Reset initialization flag + simulation_initialized = .false. + end subroutine openmc_simulation_finalize !=============================================================================== diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 62c66e5f98..5041f0793c 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -18,11 +18,12 @@ module simulation_header real(8) :: log_spacing ! spacing on logarithmic grid ! ============================================================================ - ! EIGENVALUE SIMULATION VARIABLES + ! SIMULATION VARIABLES integer :: current_batch ! current batch integer :: current_gen ! current generation within a batch integer :: total_gen = 0 ! total number of generations simulated + logical(C_BOOL), bind(C) :: simulation_initialized = .false. ! ============================================================================ ! TALLY PRECISION TRIGGER VARIABLES @@ -33,6 +34,9 @@ module simulation_header integer(8), allocatable :: work_index(:) ! starting index in source bank for each process integer(8) :: current_work ! index in source bank of current history simulated + ! ============================================================================ + ! K-EIGENVALUE SIMULATION VARIABLES + ! Temporary k-effective values type(VectorReal) :: k_generation ! single-generation estimates of k real(C_DOUBLE), bind(C) :: keff = ONE ! average k over active batches diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index cde3d87245..ba9ac5c9e0 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -206,17 +206,20 @@ def test_source_bank(capi_run): assert np.all(source['wgt'] == 1.0) -def test_keff(capi_run): - mean, std_dev = openmc.capi.keff() - assert 0.0 < mean < 2.5 - assert std_dev > 0.0 - - def test_by_batch(capi_run): openmc.capi.hard_reset() + + # Running next batch before simulation is initialized should raise an + # exception + with pytest.raises(openmc.capi.AllocationError): + openmc.capi.next_batch() + openmc.capi.simulation_init() for _ in openmc.capi.iter_batches(): - pass + # Make sure we can get k-effective during inactive/active batches + mean, std_dev = openmc.capi.keff() + assert 0.0 < mean < 2.5 + assert std_dev > 0.0 assert openmc.capi.num_realizations() == 5 for i in range(3):