crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #944 from paulromano/api-updates

Browse source 
Additions in openmc.capi and elsewhere
This commit is contained in:
Sterling Harper 2017-12-16 17:13:29 -05:00 • committed by GitHub
commit 598b5d4f2b
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
36 changed files with 1323 additions and 568 deletions
Download patch file Download diff file Expand all files Collapse all files

3
.travis.yml
View file

@ -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
@ -62,6 +63,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 ..

2
CMakeLists.txt
View file

@ -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)

5
docs/source/conf.py
View file

@ -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

26
docs/source/index.rst
View file

@ -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 <http://crpg.mit.edu>`_ at the `Massachusetts Institute of Technology
<http://web.mit.edu>`_ 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 <https://groups.google.com/forum/?fromgroups=#!forum/openmc-users>`_.
.. _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 <https://doi.org/10.1016/j.anucene.2014.07.048>`_,"
*Ann. Nucl. Energy*, **82**, 90--97 (2015).
.. only:: html
--------
Contents
--------
--------
Contents
--------
.. toctree::
:maxdepth: 1

92
docs/source/publications.rst
View file

@ -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 <http://dx.doi.org/10.1016/j.anucene.2014.07.048>`_,"
Research and Development <https://doi.org/10.1016/j.anucene.2014.07.048>`_,"
*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 <http://dx.doi.org/10.1016/j.anucene.2012.06.040>`_,"
Code <https://doi.org/10.1016/j.anucene.2012.06.040>`_,"
*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
<http://dx.doi.org/10.1016/j.pnucene.2014.12.018>`_," *Prog. Nucl. Energy*,
<https://doi.org/10.1016/j.pnucene.2014.12.018>`_," *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 <https://doi.org/10.1016/j.anucene.2017.05.002>`_,"
*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 <http://dx.doi.org/10.13182/NSE16-26>`_," *Nucl. Sci. Eng.*,
Simulations <https://doi.org/10.13182/NSE16-26>`_," *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 <http://dx.doi.org/10.1016/j.anucene.2014.10.029>`_,"
acceleration methods <https://doi.org/10.1016/j.anucene.2014.10.029>`_,"
*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 <https://doi.org/10.1016/j.nucengdes.2017.05.020>`_,"
*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
<http://dx.doi.org/10.1016/j.net.2015.12.015>`_," *Nucl. Eng. Tech.*, **48**,
<https://doi.org/10.1016/j.net.2015.12.015>`_," *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
<http://dx.doi.org/10.1016/j.anucene.2015.10.010>`_," *Ann. Nucl. Energy*,
<https://doi.org/10.1016/j.anucene.2015.10.010>`_," *Ann. Nucl. Energy*,
**87**, 601-611 (2016).
- Qicang Shen, William Boyd, Benoit Forget, and Kord Smith, "Tally precision
@ -154,9 +171,25 @@ 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
<https://doi.org/10.1016/j.anucene.2017.10.039>`_," *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).
- 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
<https://doi.org/10.1016/j.anucene.2017.09.052>`_," *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
<http://dx.doi.org/10.11889/j.0253-3219.2017.hjs.40.040502>`_,"
<https://doi.org/10.11889/j.0253-3219.2017.hjs.40.040502>`_,"
*Nucl. Techniques* **40** (4), 040504 (2017). (in Mandarin)
- Nicholas E. Stauff, Changho Lee, Paul K. Romano, and Taek K. Kim,
@ -207,7 +240,7 @@ Doppler Broadening
- Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "`Windowed multipole
for cross section Doppler broadening
<http://dx.doi.org/10.1016/j.jcp.2015.08.013>`_," *J. Comput. Phys.*, **307**,
<https://doi.org/10.1016/j.jcp.2015.08.013>`_," *J. Comput. Phys.*, **307**,
715-727 (2016).
- Jonathan A. Walsh, Benoit Forget, Kord S. Smith, and Forrest B. Brown,
@ -217,12 +250,12 @@ Doppler Broadening
- Colin Josey, Benoit Forget, and Kord Smith, "`Windowed multipole sensitivity
to target accuracy of the optimization procedure
<http://dx.doi.org/10.1080/00223131.2015.1035353>`_,"
<https://doi.org/10.1080/00223131.2015.1035353>`_,"
*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
<http://dx.doi.org/10.1016/j.anucene.2014.08.046>`_," *Ann. Nucl. Energy*,
<https://doi.org/10.1016/j.anucene.2014.08.046>`_," *Ann. Nucl. Energy*,
**75**, 358--364 (2015).
- Tuomas Viitanen, Jaakko Leppanen, and Benoit Forget, "Target motion sampling
@ -231,13 +264,17 @@ Doppler Broadening
- Benoit Forget, Sheng Xu, and Kord Smith, "`Direct Doppler broadening in Monte
Carlo simulations using the multipole representation
<http://dx.doi.org/10.1016/j.anucene.2013.09.043>`_," *Ann. Nucl. Energy*,
<https://doi.org/10.1016/j.anucene.2013.09.043>`_," *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 +292,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
<http://dx.doi.org/10.1051/epjconf/201611106001>`_," *Eur. Phys. J. Web Conf.*
<https://doi.org/10.1051/epjconf/201611106001>`_," *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
<http://dx.doi.org/10.1016/j.cpc.2015.05.025>`_", *Comput. Phys. Commun.*,
<https://doi.org/10.1016/j.cpc.2015.05.025>`_", *Comput. Phys. Commun.*,
**196**, 134-142 (2015).
- Jonathan A. Walsh, Benoit Forget, Kord S. Smith, Brian C. Kiedrowski, and
@ -280,7 +317,7 @@ Nuclear Data
- Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "`Accelerated sampling of
the free gas resonance elastic scattering kernel
<http://dx.doi.org/10.1016/j.anucene.2014.01.017>`_," *Ann. Nucl. Energy*,
<https://doi.org/10.1016/j.anucene.2014.01.017>`_," *Ann. Nucl. Energy*,
**69**, 116--124 (2014).
-----------
@ -325,45 +362,45 @@ Parallelism
- Nicholas Horelik, Andrew Siegel, Benoit Forget, and Kord Smith, "`Monte Carlo
domain decomposition for robust nuclear reactor analysis
<http://dx.doi.org/10.1016/j.parco.2014.10.001>`_," *Parallel Comput.*,
<https://doi.org/10.1016/j.parco.2014.10.001>`_," *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 <http://dx.doi.org/10.1016/j.cpc.2013.10.008>`_,"
using energy banding <https://doi.org/10.1016/j.cpc.2013.10.008>`_,"
*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
<http://dx.doi.org/10.1051/snamc/201404301>`_," *Proc. Joint International
<https://doi.org/10.1051/snamc/201404301>`_," *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 <http://dx.doi.org/10.1051/snamc/201404207>`_," *Proc. Joint
calculations <https://doi.org/10.1051/snamc/201404207>`_," *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
<http://dx.doi.org/10.1051/snamc/201404208>`_," *Proc. Joint International
<https://doi.org/10.1051/snamc/201404208>`_," *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
<http://dx.doi.org/10.1177/1094342013492179>`_," *Int. J. High
<https://doi.org/10.1177/1094342013492179>`_," *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
<http://dx.doi.org/10.1016/j.jcp.2013.06.011>`_," *J. Comput. Phys.*, **252**,
<https://doi.org/10.1016/j.jcp.2013.06.011>`_," *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 <http://dx.doi.org/10.1016/j.jcp.2012.06.012>`_," *J. Comput. Phys.*,
analysis <https://doi.org/10.1016/j.jcp.2012.06.012>`_," *J. Comput. Phys.*,
**235**, 901--911 (2013).
@ -372,13 +409,18 @@ Parallelism
519--522 (2012).
- Paul K. Romano and Benoit Forget, "`Parallel Fission Bank Algorithms in Monte
Carlo Criticality Calculations <http://dx.doi.org/10.13182/NSE10-98>`_,"
Carlo Criticality Calculations <https://doi.org/10.13182/NSE10-98>`_,"
*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 <https://doi.org/10.1016/j.jcp.2017.08.025>`_," *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
<http://hdl.handle.net/1721.1/107880>`_," *Trans. Am. Nucl. Soc.*, **115**,
@ -386,7 +428,7 @@ Depletion
- Anas Gul, K. S. Chaudri, R. Khan, and M. Azeen, "`Development and verification
of LOOP: A Linkage of ORIGEN2.2 and OpenMC
<http://dx.doi.org/10.1016/j.anucene.2016.09.016>`_," *Ann. Nucl. Energy*,
<https://doi.org/10.1016/j.anucene.2016.09.016>`_," *Ann. Nucl. Energy*,
**99**, 321--327 (2017).
- Kai Huang, Hongchun Wu, Yunzhao Li, and Liangzhi Cao, "Generalized depletion

12
docs/source/pythonapi/base.rst
View file

@ -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

7
docs/source/pythonapi/capi.rst
View file

@ -18,12 +18,19 @@ 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.num_realizations
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
-------

8
docs/source/pythonapi/index.rst
View file

@ -40,13 +40,13 @@ Modules
-------
.. toctree::
:maxdepth: 2
:maxdepth: 1
base
stats
mgxs
model
examples
mgxs
stats
data
capi
examples
openmoc

17
docs/source/pythonapi/model.rst
View file

@ -2,6 +2,23 @@
:mod:`openmc.model` -- Model Building
-------------------------------------
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
-------------------

4
openmc/__init__.py
View file

@ -27,5 +27,9 @@ 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
__version__ = '0.9.0'

50
openmc/capi/cell.py
View file

@ -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
@ -22,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
@ -51,11 +57,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):
@ -63,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()
@ -113,11 +147,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

138
openmc/capi/core.py
View file

@ -1,9 +1,22 @@
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, c_char_p,
POINTER, Structure)
from warnings import warn
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
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
@ -18,9 +31,17 @@ _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.argtypes = [POINTER(c_char_p)]
_dll.openmc_statepoint_write.restype = None
def calculate_volumes():
@ -43,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.
@ -53,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):
@ -66,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():
@ -102,6 +123,40 @@ def init(intracomm=None):
_dll.openmc_init(None)
def iter_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 = 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.
@ -111,9 +166,26 @@ 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():
"""Run 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():
@ -131,6 +203,50 @@ 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(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
def run_in_memory(intracomm=None):
"""Provides context manager for calling OpenMC shared library functions.

21
openmc/capi/error.py
View file

@ -1,41 +1,42 @@
from ctypes import c_int, c_char
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."""
@ -70,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))

191
openmc/capi/tally.py
View file

@ -2,23 +2,29 @@ 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
import scipy.stats
from openmc.data.reaction import REACTION_NAME
from . import _dll, Nuclide
from .core import _FortranObjectWithID
from .error import _error_handler, AllocationError, InvalidIDError
from .filter import _get_filter
__all__ = ['Tally', 'tallies']
__all__ = ['Tally', 'tallies', 'global_tallies', 'num_realizations']
# 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
@ -26,10 +32,17 @@ _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
_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 +64,54 @@ _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'
}
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))
# Get sum, sum-of-squares, and number of realizations
sum_ = array[:, 1]
sum_sq = array[:, 2]
n = num_realizations()
# Determine mean
if n > 0:
mean = sum_ / n
else:
mean = sum_.copy()
# Determine standard deviation
nonzero = np.abs(mean) > 0
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))
def num_realizations():
"""Number of realizations of global tallies."""
return c_int32.in_dll(_dll, 'n_realizations').value
class Tally(_FortranObjectWithID):
"""Tally stored internally.
@ -74,10 +135,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()
@ -130,21 +197,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[]
@ -153,15 +205,60 @@ 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):
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):
@ -169,21 +266,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.inf)
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):

67
openmc/clean_xml.py
View file

@ -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

9
openmc/geometry.py
View file

@ -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,13 +81,16 @@ 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
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

9
openmc/material.py
View file

@ -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')

1
openmc/model/__init__.py
View file

@ -1,2 +1,3 @@
from .triso import *
from .model import *
from .funcs import *

279
openmc/model/funcs.py Normal file
View file

@ -0,0 +1,279 @@
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.
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

253
openmc/surface.py
View file

@ -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

18
openmc/universe.py
View file

@ -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

12
readme.rst
View file

@ -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 <https://doi.org/10.1016/j.anucene.2014.07.048>`_,"
*Ann. Nucl. Energy*, **82**, 90--97 (2015).
---------------
Troubleshooting
---------------

3
src/api.F90
View file

@ -36,10 +36,12 @@ 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
public :: openmc_extend_filters
public :: openmc_extend_cells
public :: openmc_extend_materials
public :: openmc_extend_tallies
public :: openmc_filter_get_id
@ -74,6 +76,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

22
src/bank_header.F90
View file

@ -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

3
src/constants.F90
View file

@ -306,7 +306,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

28
src/eigenvalue.F90
View file

@ -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

49
src/geometry_header.F90
View file

@ -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
@ -538,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

17
src/input_xml.F90
View file

@ -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 % reserve(n_max_batches*gen_per_batch)
! Get the trigger information for keff
if (check_for_node(node_base, "keff_trigger")) then
@ -1921,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
@ -2058,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

10
src/output.F90
View file

@ -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
@ -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) * &

145
src/simulation.F90
View file

@ -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
@ -43,6 +43,8 @@ module simulation
implicit none
private
public :: openmc_run
public :: openmc_simulation_init
public :: openmc_simulation_finalize
contains
@ -54,74 +56,81 @@ contains
subroutine openmc_run() bind(C)
call openmc_simulation_init()
do while (openmc_next_batch() == 0)
end do
call openmc_simulation_finalize()
end subroutine openmc_run
!===============================================================================
! OPENMC_NEXT_BATCH
!===============================================================================
function openmc_next_batch() result(retval) bind(C)
integer(C_INT) :: retval
type(Particle) :: p
integer(8) :: i_work
call initialize_simulation()
! Make sure simulation has been initialized
if (.not. simulation_initialized) then
retval = -3
return
end if
! Turn on inactive timer
call time_inactive % start()
call initialize_batch()
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_max_batches
! Handle restart runs
if (restart_run .and. current_batch <= restart_batch) then
call replay_batch_history()
retval = 0
return
end if
call initialize_batch()
! =======================================================================
! LOOP OVER GENERATIONS
GENERATION_LOOP: do current_gen = 1, gen_per_batch
! Handle restart runs
if (restart_run .and. current_batch <= restart_batch) then
call replay_batch_history()
cycle BATCH_LOOP
end if
call initialize_generation()
! =======================================================================
! LOOP OVER GENERATIONS
GENERATION_LOOP: do current_gen = 1, gen_per_batch
! Start timer for transport
call time_transport % start()
call initialize_generation()
! Start timer for transport
call time_transport % start()
! ====================================================================
! LOOP OVER PARTICLES
! ====================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(static) firstprivate(p) copyin(tally_derivs)
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work
! grab source particle from bank
call initialize_history(p, current_work)
! grab source particle from bank
call initialize_history(p, current_work)
! transport particle
call transport(p)
! transport particle
call transport(p)
end do PARTICLE_LOOP
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()
! Accumulate time for transport
call time_transport % stop()
call finalize_generation()
call finalize_generation()
end do GENERATION_LOOP
end do GENERATION_LOOP
call finalize_batch()
call finalize_batch()
if (satisfy_triggers) exit BATCH_LOOP
! 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 do BATCH_LOOP
call time_active % stop()
! ==========================================================================
! END OF RUN WRAPUP
call finalize_simulation()
! Clear particle
call p % clear()
end subroutine openmc_run
end function openmc_next_batch
!===============================================================================
! INITIALIZE_HISTORY
@ -176,6 +185,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)
@ -184,7 +196,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()
@ -336,7 +351,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
@ -385,7 +400,10 @@ contains
! INITIALIZE_SIMULATION
!===============================================================================
subroutine initialize_simulation()
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()
@ -426,14 +444,20 @@ contains
end if
end if
end subroutine initialize_simulation
! Reset current batch
current_batch = 0
! Set flag indicating initialization is done
simulation_initialized = .true.
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
@ -443,12 +467,18 @@ 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()
!$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()
@ -496,7 +526,10 @@ contains
if (check_overlaps) call print_overlap_check()
end if
end subroutine finalize_simulation
! Reset initialization flag
simulation_initialized = .false.
end subroutine openmc_simulation_finalize
!===============================================================================
! CALCULATE_WORK determines how many particles each processor should simulate

24
src/simulation_header.F90
View file

@ -1,8 +1,11 @@
module simulation_header
use, intrinsic :: ISO_C_BINDING
use bank_header
use constants
use settings, only: gen_per_batch
use stl_vector, only: VectorReal
implicit none
@ -15,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
@ -30,16 +34,19 @@ 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
real(8), allocatable :: 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
type(VectorReal) :: k_generation ! single-generation estimates of 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
! 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 +92,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

50
src/state_point.F90
View file

@ -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
@ -40,10 +40,11 @@ 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(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")
@ -121,8 +128,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)
@ -433,7 +443,7 @@ contains
call file_close(file_id)
end if
end subroutine write_state_point
end subroutine openmc_statepoint_write
!===============================================================================
! WRITE_SOURCE_POINT
@ -629,6 +639,7 @@ contains
subroutine load_state_point()
integer :: i
integer :: n
integer :: int_array(3)
integer, allocatable :: array(:)
integer(HID_T) :: file_id
@ -707,10 +718,15 @@ 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), &
file_id, "k_generation")
call read_dataset(entropy(1:restart_batch*gen_per_batch), &
file_id, "entropy")
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 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")
call read_dataset(k_abs_tra, file_id, "k_abs_tra")

5
src/summary.F90
View file

@ -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")

59
src/tallies/tally_header.F90
View file

@ -22,9 +22,11 @@ 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
public :: openmc_tally_get_scores
public :: openmc_tally_results
public :: openmc_tally_set_filters
public :: openmc_tally_set_id
@ -145,7 +147,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
@ -469,6 +471,20 @@ 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
call set_errmsg("Global tallies have not been allocated yet.")
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
@ -510,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
@ -535,6 +567,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.

244
tests/unit_tests/test_capi.py Normal file
View file

@ -0,0 +1,244 @@
#!/usr/bin/env python
from collections 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_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():
# 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):
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]

8
tests/unit_tests/test_element_wo.py
View file

@ -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)