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

Merge branch 'c-api' into photon-new

Browse source
This commit is contained in:
Paul Romano 2017-07-20 06:48:54 -05:00
commit 698cb09682
31 changed files with 631 additions and 407 deletions
Download patch file Download diff file Expand all files Collapse all files

58
.travis.yml
View file

@ -9,12 +9,23 @@ addons:
packages:
- gfortran
- g++
- mpich
- libmpich-dev
cache:
directories:
- $HOME/mpich_install
- $HOME/hdf5_install
- $HOME/phdf5_install
- $HOME/nndc_hdf5
env:
- OPENMC_CONFIG="check_source"
- OPENMC_CONFIG="^hdf5-debug$"
- OPENMC_CONFIG="^omp-hdf5-debug$"
- OPENMC_CONFIG="^mpi-hdf5-debug$"
- OPENMC_CONFIG="^phdf5-debug$"
# We aren't testing the check_source script so just run it with Python 3.
matrix:
exclude:
- python: "2.7"
env: OPENMC_CONFIG="check_source"
before_install:
# ============== Handle Python third-party packages ==============
@ -29,30 +40,37 @@ before_install:
- conda config --set always_yes yes --set changeps1 no
- conda update -q conda
- conda info -a
- conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION six numpy scipy h5py=2.5 pandas
- source activate test-environment
# Install GCC, MPICH, HDF5, PHDF5
- ./tests/travis_install.sh
- export FC=gfortran
- export MPI_DIR=$HOME/mpich_install
- export PHDF5_DIR=$HOME/phdf5_install
- export HDF5_DIR=$HOME/hdf5_install
- if [[ $OPENMC_CONFIG != "check_source" ]]; then
conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION six numpy scipy h5py=2.5 pandas;
source activate test-environment;
sudo add-apt-repository ppa:nschloe/hdf5-backports -y;
sudo apt-get update -q;
sudo apt-get install libhdf5-serial-dev libhdf5-mpich-dev -y;
export FC=gfortran;
export MPI_DIR=/usr;
export PHDF5_DIR=/usr;
export HDF5_DIR=/usr;
fi
install: true
before_script:
- if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget https://anl.box.com/shared/static/a6sw2cep34wlz6b9i9jwiotaqoayxcxt.xz -O - | tar -C $HOME -xvJ;
- if [[ $OPENMC_CONFIG != "check_source" ]]; then
if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget https://anl.box.com/shared/static/a6sw2cep34wlz6b9i9jwiotaqoayxcxt.xz -O - | tar -C $HOME -xvJ;
fi;
export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml;
git clone --branch=master git://github.com/smharper/windowed_multipole_library.git wmp_lib;
tar xzvf wmp_lib/multipole_lib.tar.gz;
export OPENMC_MULTIPOLE_LIBRARY=$PWD/multipole_lib;
fi
- export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml
- git clone --branch=master git://github.com/smharper/windowed_multipole_library.git wmp_lib
- tar xzvf wmp_lib/multipole_lib.tar.gz
- export OPENMC_MULTIPOLE_LIBRARY=$PWD/multipole_lib
script:
- cd tests
- export OMP_NUM_THREADS=2
- ./travis.sh
- if [[ $OPENMC_CONFIG == "check_source" ]]; then
./check_source.py;
else
./run_tests.py -C $OPENMC_CONFIG -j 2;
fi
- cd ..

10
CMakeLists.txt
View file

@ -67,10 +67,12 @@ endif()
# this, we check for the environment variable HDF5_ROOT and if it exists, use it
# to check whether its a parallel version.
if(DEFINED ENV{HDF5_ROOT} AND EXISTS $ENV{HDF5_ROOT}/bin/h5pcc)
set(HDF5_PREFER_PARALLEL TRUE)
else()
set(HDF5_PREFER_PARALLEL FALSE)
if(NOT DEFINED HDF5_PREFER_PARALLEL)
if(DEFINED ENV{HDF5_ROOT} AND EXISTS $ENV{HDF5_ROOT}/bin/h5pcc)
set(HDF5_PREFER_PARALLEL TRUE)
else()
set(HDF5_PREFER_PARALLEL FALSE)
endif()
endif()
find_package(HDF5 COMPONENTS Fortran_HL)

93
docs/source/capi/index.rst
View file

@ -4,15 +4,42 @@
C API
=====
.. c:function:: int openmc_find(double* xyz, int rtype)
.. c:function:: void openmc_calculate_voumes()
Return the ID of the cell/material containing a given point
Run a stochastic volume calculation
.. c:function:: int openmc_cell_set_temperature(int id, double T, int* instance)
Set the temperature of a cell.
:param id: ID of the cell
:type id: int
:param T: Temperature in Kelvin
:type T: double
:param instance: Which instance of the cell. To set the temperature for all
instances, pass a null pointer.
:type instance: int*
:return: Return status (negative if an error occurred)
:rtype: int
.. c:function:: void openmc_finalize()
Finalize a simulation
.. c:function:: void openmc_find(double* xyz, int rtype, int* id, int* instance)
Determine the ID of the cell/material containing a given point
:param xyz: Cartesian coordinates
:type xyz: double[3]
:param rtype: Which ID to return (1=cell, 2=material)
:type rtype: int
:rtype: int
:param id: ID of the cell/material found. If a material is requested and the
point is in a void, the ID is 0. If an error occurs, the ID is -1.
:type id: int
:param instance: If a cell is repetaed in the geometry, the instance of the
cell that was found and zero otherwise.
:type instance: int
.. c:function:: void openmc_init(int intracomm)
@ -21,9 +48,54 @@ C API
:param intracomm: MPI intracommunicator
:type intracomm: int
.. c:function:: void openmc_finalize()
.. c:function:: int openmc_load_nuclide(char name[])
Finalize a simulation
Load data for a nuclide from the HDF5 data library.
:param name: Name of the nuclide.
:type name: char[]
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_material_add_nuclide(int id, char name[], double density)
Add a nuclide to an existing material. If the nuclide already exists, the
density is overwritten.
:param id: ID of the material
:type id: int
:param name: Name of the nuclide
:type name: char[]
:param density: Density in atom/b-cm
:type density: double
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_material_get_densities(int id, double* ptr)
Get an array of nuclide densities for a material.
:param id: ID of the material
:type id: int
:param ptr: Pointer to the array of densities
:type ptr: double*
:return: Length of the array
:rtype: int
.. c:function:: int openmc_material_set_density(int id, double density)
Set the density of a material.
:param id: ID of the material
:type id: int
:param density: Density of the material in atom/b-cm
:type density: double
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: void openmc_plot_geometry()
Run plotting mode.
.. c:function:: void openmc_reset()
@ -32,3 +104,14 @@ C API
.. c:function:: void openmc_run()
Run a simulation
.. c:function:: void openmc_tally_results(int id, double** ptr, int shape_[3])
Get a pointer to tally results array.
:param id: ID of the tally
:type id: int
:param ptr: Pointer to the results array
:type ptr: double**
:param shape_: Shape of the results array
:type shape_: int[3]

6
docs/source/io_formats/materials.rst
View file

@ -107,9 +107,13 @@ Each ``material`` element can have the following attributes or sub-elements:
multi-group :ref:`energy_mode`.
:sab:
Associates an S(a,b) table with the material. This element has one
Associates an S(a,b) table with the material. This element has an
attribute/sub-element called ``name``. The ``name`` attribute
is the name of the S(a,b) table that should be associated with the material.
There is also an optional ``fraction`` element which indicates what fraction
of the relevant nuclides will be affected by the S(a,b) table (e.g. which
fraction of a material is crystalline versus amorphous). ``fraction``
defaults to unity.
*Default*: None

17
docs/source/pythonapi/capi.rst Normal file
View file

@ -0,0 +1,17 @@
---------------------------------------------------
:data:`openmc.capi` -- Python bindings to the C API
---------------------------------------------------
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
openmc.capi.lib_context
.. autosummary::
:toctree: generated
:nosignatures:
:template: myclass.rst
openmc.capi.OpenMCLibrary

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

@ -47,5 +47,6 @@ Modules
mgxs
model
data
capi
examples
openmoc

2
openmc/__init__.py
View file

@ -27,6 +27,6 @@ from openmc.particle_restart import *
from openmc.mixin import *
from openmc.plotter import *
from openmc.search import *
from openmc.capi import lib
from openmc.capi import *
__version__ = '0.9.0'

138
openmc/capi.py
View file

@ -1,4 +1,5 @@
from ctypes import CDLL, c_int, POINTER, byref, c_double, c_char_p
from contextlib import contextmanager
from ctypes import CDLL, c_int, c_int32, c_double, c_char_p, POINTER
import sys
from warnings import warn
@ -7,52 +8,63 @@ from numpy.ctypeslib import as_array
import pkg_resources
__all__ = ['OpenMCLibrary', 'lib', 'lib_context']
_int3 = c_int*3
_double3 = c_double*3
_double_array = POINTER(POINTER(c_double))
class _OpenMCLibrary(object):
class OpenMCLibrary(object):
"""Provides bindings to C functions defined by OpenMC shared library.
This class is normally not directly instantiated. Instead, when the
:mod:`openmc` package is imported, an instance is automatically created with
the name :data:`openmc.lib`. Calls to the OpenMC can then be made using that
instance, for example:
.. code-block:: python
openmc.lib.init()
openmc.lib.run()
"""
def __init__(self, filename):
self._dll = CDLL(filename)
# Set argument/return types
self._dll.openmc_calculate_volumes.restype = None
self._dll.openmc_cell_set_temperature.argtypes = [
c_int, c_double]
c_int32, c_double, c_int32]
self._dll.openmc_cell_set_temperature.restype = c_int
self._dll.openmc_finalize.restype = None
self._dll.openmc_find.argtypes = [POINTER(_double3), c_int]
self._dll.openmc_find.restype = c_int
self._dll.openmc_find.argtypes = [
POINTER(_double3), c_int, POINTER(c_int32), POINTER(c_int32)]
self._dll.openmc_find.restype = None
self._dll.openmc_init.argtypes = [POINTER(c_int)]
self._dll.openmc_init.restype = None
self._dll.openmc_load_nuclide.argtypes = [c_char_p]
self._dll.openmc_load_nuclide.restype = c_int
self._dll.openmc_material_add_nuclide.argtypes = [
c_int, c_char_p, c_double]
c_int32, c_char_p, c_double]
self._dll.openmc_material_add_nuclide.restype = c_int
self._dll.openmc_material_get_densities.argtypes = [
c_int, _double_array]
c_int32, _double_array]
self._dll.openmc_material_get_densities.restype = c_int
self._dll.openmc_material_set_density.argtypes = [c_int, c_double]
self._dll.openmc_material_set_density.argtypes = [c_int32, c_double]
self._dll.openmc_material_set_density.restype = c_int
self._dll.openmc_plot_geometry.restype = None
self._dll.openmc_run.restype = None
self._dll.openmc_reset.restype = None
self._dll.openmc_set_density.argtypes = [POINTER(_double3), c_double]
self._dll.openmc_set_density.restype = c_int
self._dll.openmc_set_temperature.argtypes = [
POINTER(_double3), c_double]
self._dll.openmc_set_temperature.restype = c_int
self._dll.openmc_tally_results.argtypes = [
c_int, _double_array, POINTER(_int3)]
c_int32, _double_array, POINTER(_int3)]
self._dll.openmc_tally_results.restype = None
def calculate_volumes(self):
"""Run stochastic volume calculation"""
return self._dll.openmc_calculate_volumes()
def cell_set_temperature(self, cell_id, T):
def cell_set_temperature(self, cell_id, T, instance=None):
"""Set the temperature of a cell
Parameters
@ -61,8 +73,15 @@ class _OpenMCLibrary(object):
ID of the cell
T : float
Temperature in K
instance : int or None
Which instance of the cell
"""
return self._dll.openmc_cell_set_temperature(cell_id, T)
if instance is not None:
return self._dll.openmc_cell_set_temperature(
cell_id, T, instance)
else:
return self._dll.openmc_cell_set_temperature(cell_id, T, None)
def finalize(self):
"""Finalize simulation and free memory"""
@ -83,22 +102,31 @@ class _OpenMCLibrary(object):
int or None
ID of the cell or material. If 'material' is requested and no
material exists at the given coordinate, None is returned.
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.
"""
# Set second argument to openmc_find
if rtype == 'cell':
return self._dll.openmc_find(_double3(*xyz), 1)
r_int = 1
elif rtype == 'material':
uid = self._dll.openmc_find(_double3(*xyz), 2)
return uid if uid != 0 else None
r_int = 2
else:
raise ValueError('Unknown return type: {}'.format(rtype))
# Call openmc_find
uid = c_int32()
instance = c_int32()
self._dll.openmc_find(_double3(*xyz), r_int, uid, instance)
return (uid.value if uid != 0 else None), instance.value
def init(self, intracomm=None):
"""Initialize OpenMC
Parameters
----------
intracomm : int or None
intracomm : mpi4py.MPI.Intracomm or None
MPI intracommunicator
"""
@ -109,7 +137,7 @@ class _OpenMCLibrary(object):
intracomm = intracomm.py2f()
except AttributeError:
pass
return self._dll.openmc_init(byref(c_int(intracomm)))
return self._dll.openmc_init(c_int(intracomm))
else:
return self._dll.openmc_init(None)
@ -165,7 +193,7 @@ class _OpenMCLibrary(object):
"""
data = POINTER(c_double)()
n = self._dll.openmc_material_get_densities(mat_id, byref(data))
n = self._dll.openmc_material_get_densities(mat_id, data)
if data:
return as_array(data, (n,))
else:
@ -201,42 +229,6 @@ class _OpenMCLibrary(object):
"""Run simulation"""
return self._dll.openmc_run()
def set_density(self, xyz, density):
"""Set density at a given point.
Parameters
----------
xyz : iterable of float
Cartesian coordinates at position
density : float
Density in atom/b-cm
Returns
-------
int
Return status (negative if an error occurs)
"""
return self._dll.openmc_set_density(_double3(*xyz), density)
def set_temperature(self, xyz, T):
"""Set temperature at a given point.
Parameters
----------
xyz : iterable of float
Cartesian coordinates at position
T : float
Temperature in K
Returns
-------
int
Return status (negative if an error occurs)
"""
return self._dll.openmc_set_temperature(_double3(*xyz), T)
def tally_results(self, tally_id):
"""Get tally results array
@ -253,7 +245,7 @@ class _OpenMCLibrary(object):
"""
data = POINTER(c_double)()
shape = _int3()
self._dll.openmc_tally_results(tally_id, byref(data), byref(shape))
self._dll.openmc_tally_results(tally_id, data, shape)
if data:
return as_array(data, tuple(shape[::-1]))
else:
@ -268,6 +260,32 @@ class _OpenMCLibrary(object):
.format(key))
@contextmanager
def lib_context(intracomm=None):
"""Provides context manager for calling OpenMC shared library functions.
This function is intended to be used in a 'with' statement and ensures that
OpenMC is properly initialized/finalized. At the completion of the 'with'
block, all memory that was allocated during the block is freed. For
example::
with openmc.lib_context() as lib:
for i in range(n_iters):
lib.reset()
do_stuff()
lib.run()
Parameters
----------
intracomm : mpi4py.MPI.Intracomm or None
MPI intracommunicator
"""
lib.init(comm)
yield lib
lib.finalize()
# Determine shared-library suffix
if sys.platform == 'darwin':
suffix = 'dylib'
@ -278,7 +296,7 @@ else:
filename = pkg_resources.resource_filename(
__name__, '_libopenmc.{}'.format(suffix))
try:
lib = _OpenMCLibrary(filename)
lib = OpenMCLibrary(filename)
except OSError:
warn("OpenMC shared library is not available from the Python API. This "
"means you will not be able to use openmc.lib to make in-memory "

2
openmc/data/njoy.py
View file

@ -358,7 +358,7 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None,
# Determine whether elastic is incoherent (0) or coherent (1)
file_obj = StringIO(ev_thermal.section[7, 2])
elastic_type = endf.get_head_record(file_obj)[2]
elastic_type = endf.get_head_record(file_obj)[2] - 1
else:
elastic = 0
mt_elastic = 0

17
openmc/material.py
View file

@ -618,13 +618,18 @@ class Material(IDManagerMixin):
if element == elm[0]:
self._elements.remove(elm)
def add_s_alpha_beta(self, name):
def add_s_alpha_beta(self, name, fraction=1.0):
r"""Add an :math:`S(\alpha,\beta)` table to the material
Parameters
----------
name : str
Name of the :math:`S(\alpha,\beta)` table
fraction : float
The fraction of relevant nuclei that are affected by the
:math:`S(\alpha,\beta)` table. For example, if the material is a
block of carbon that is 60% graphite and 40% amorphous then add a
graphite :math:`S(\alpha,\beta)` table with fraction=0.6.
"""
@ -638,13 +643,17 @@ class Material(IDManagerMixin):
'non-string table name "{}"'.format(self._id, name)
raise ValueError(msg)
cv.check_type('S(a,b) fraction', fraction, Real)
cv.check_greater_than('S(a,b) fraction', fraction, 0.0, True)
cv.check_less_than('S(a,b) fraction', fraction, 1.0, True)
new_name = openmc.data.get_thermal_name(name)
if new_name != name:
msg = 'OpenMC S(a,b) tables follow the GND naming convention. ' \
'Table "{}" is being renamed as "{}".'.format(name, new_name)
warnings.warn(msg)
self._sab.append(new_name)
self._sab.append((new_name, fraction))
def make_isotropic_in_lab(self):
for nuclide, percent, percent_type in self._nuclides:
@ -972,7 +981,9 @@ class Material(IDManagerMixin):
if len(self._sab) > 0:
for sab in self._sab:
subelement = ET.SubElement(element, "sab")
subelement.set("name", sab)
subelement.set("name", sab[0])
if sab[1] != 1.0:
subelement.set("fraction", str(sab[1]))
return element

110
src/api.F90
View file

@ -24,15 +24,13 @@ module openmc_api
public :: openmc_finalize
public :: openmc_find
public :: openmc_init
public :: openmc_load_nuclide
public :: openmc_material_add_nuclide
public :: openmc_material_get_densities
public :: openmc_material_set_density
public :: openmc_load_nuclide
public :: openmc_plot_geometry
public :: openmc_reset
public :: openmc_run
public :: openmc_set_density
public :: openmc_set_temperature
public :: openmc_tally_results
contains
@ -41,12 +39,13 @@ contains
! OPENMC_CELL_SET_TEMPERATURE sets the temperature of a cell
!===============================================================================
function openmc_cell_set_temperature(id, T) result(err) bind(C)
integer(C_INT), value, intent(in) :: id ! id of cell
function openmc_cell_set_temperature(id, T, instance) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: id ! id of cell
real(C_DOUBLE), value, intent(in) :: T
integer(C_INT32_T), optional, intent(in) :: instance
integer(C_INT) :: err
integer :: i
integer :: i, n
err = -1
if (allocated(cells)) then
@ -54,8 +53,16 @@ contains
i = cell_dict % get_key(id)
associate (c => cells(i))
if (allocated(c % sqrtkT)) then
c % sqrtkT(:) = sqrt(K_BOLTZMANN * T)
err = 0
n = size(c % sqrtkT)
if (present(instance) .and. n > 1) then
if (instance >= 0 .and. instance < n) then
c % sqrtkT(instance + 1) = sqrt(K_BOLTZMANN * T)
err = 0
end if
else
c % sqrtkT(:) = sqrt(K_BOLTZMANN * T)
err = 0
end if
end if
end associate
end if
@ -66,10 +73,11 @@ contains
! OPENMC_FIND determines the ID or a cell or material at a given point in space
!===============================================================================
function openmc_find(xyz, rtype) result(id) bind(C)
subroutine openmc_find(xyz, rtype, id, instance) bind(C)
real(C_DOUBLE), intent(in) :: xyz(3) ! Cartesian point
integer(C_INT), intent(in), value :: rtype ! 1 for cell, 2 for material
integer(C_INT) :: id
integer(C_INT32_T), intent(out) :: id
integer(C_INT32_T), intent(out) :: instance
logical :: found
type(Particle) :: p
@ -80,6 +88,7 @@ contains
call find_cell(p, found)
id = -1
instance = -1
if (found) then
if (rtype == 1) then
id = cells(p % coord(p % n_coord) % cell) % id
@ -90,8 +99,9 @@ contains
id = materials(p % material) % id
end if
end if
instance = p % cell_instance - 1
end if
end function openmc_find
end subroutine openmc_find
!===============================================================================
! OPENMC_LOAD_NUCLIDE loads a nuclide from the cross section library
@ -160,7 +170,7 @@ contains
!===============================================================================
function openmc_material_add_nuclide(id, name, density) result(err) bind(C)
integer(C_INT), value, intent(in) :: id
integer(C_INT32_T), value, intent(in) :: id
character(kind=C_CHAR) :: name(*)
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
@ -232,8 +242,8 @@ contains
!===============================================================================
function openmc_material_get_densities(id, ptr) result(n) bind(C)
integer(C_INT), intent(in), value :: id
type(C_PTR), intent(out) :: ptr
integer(C_INT32_T), intent(in), value :: id
type(C_PTR), intent(out) :: ptr
integer(C_INT) :: n
ptr = C_NULL_PTR
@ -256,7 +266,7 @@ contains
!===============================================================================
function openmc_material_set_density(id, density) result(err) bind(C)
integer(C_INT), value, intent(in) :: id
integer(C_INT32_T), value, intent(in) :: id
real(C_DOUBLE), value, intent(in) :: density
integer(C_INT) :: err
@ -310,77 +320,23 @@ contains
end subroutine openmc_reset
!===============================================================================
! OPENMC_SET_DENSITY sets the density of a material at a given point
!===============================================================================
function openmc_set_density(xyz, density) result(err) bind(C)
real(C_DOUBLE), intent(in) :: xyz(3)
real(C_DOUBLE), intent(in), value :: density
integer(C_INT) :: err
logical :: found
type(Particle) :: p
call p % initialize()
p % coord(1) % xyz(:) = xyz
p % coord(1) % uvw(:) = [ZERO, ZERO, ONE]
call find_cell(p, found)
err = -1
if (found) then
if (p % material /= MATERIAL_VOID) then
associate (m => materials(p % material))
err = m % set_density(density, nuclides)
end associate
end if
end if
end function openmc_set_density
!===============================================================================
! OPENMC_SET_TEMPERATURE sets the temperature of a cell at a given point
!===============================================================================
function openmc_set_temperature(xyz, T) result(err) bind(C)
real(C_DOUBLE), intent(in) :: xyz(3)
real(C_DOUBLE), intent(in), value :: T
integer(C_INT) :: err
logical :: found
type(Particle) :: p
call p % initialize()
p % coord(1) % xyz(:) = xyz
p % coord(1) % uvw(:) = [ZERO, ZERO, ONE]
call find_cell(p, found)
err = -1
if (found) then
associate (c => cells(p % coord(p % n_coord) % cell))
if (size(c % sqrtkT) > 1) then
c % sqrtkT(p % cell_instance) = sqrt(K_BOLTZMANN * T)
else
c % sqrtkT(1) = sqrt(K_BOLTZMANN * T)
end if
err = 0
end associate
end if
end function openmc_set_temperature
!===============================================================================
! OPENMC_TALLY_RESULTS 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.
!===============================================================================
subroutine openmc_tally_results(i, ptr, shape_) bind(C)
integer(C_INT), intent(in), value :: i
type(C_PTR), intent(out) :: ptr
integer(C_INT), intent(out) :: shape_(3)
subroutine openmc_tally_results(id, ptr, shape_) bind(C)
integer(C_INT32_T), intent(in), value :: id
type(C_PTR), intent(out) :: ptr
integer(C_INT), intent(out) :: shape_(3)
integer :: i
ptr = C_NULL_PTR
if (allocated(tallies)) then
if (i >= 1 .and. i <= size(tallies)) then
if (tally_dict % has_key(id)) then
i = tally_dict % get_key(id)
if (allocated(tallies(i) % results)) then
ptr = C_LOC(tallies(i) % results(1,1,1))
shape_(:) = shape(tallies(i) % results)

99
src/cross_section.F90
View file

@ -61,6 +61,7 @@ contains
integer :: i_grid ! index into logarithmic mapping array or material
! union grid
real(8) :: atom_density ! atom density of a nuclide
real(8) :: sab_frac ! fraction of atoms affected by S(a,b)
logical :: check_sab ! should we check for S(a,b) table?
! Exit subroutine if material is void
@ -78,21 +79,25 @@ contains
! Add contribution from each nuclide in material
do i = 1, mat % n_nuclides
! ========================================================================
! ======================================================================
! CHECK FOR S(A,B) TABLE
i_sab = 0
sab_frac = ZERO
! Check if this nuclide matches one of the S(a,b) tables specified -- this
! relies on i_sab_nuclides being in sorted order
! Check if this nuclide matches one of the S(a,b) tables specified.
! This relies on i_sab_nuclides being in sorted order
if (check_sab) then
if (i == mat % i_sab_nuclides(j)) then
! Get index in sab_tables
i_sab = mat % i_sab_tables(j)
sab_frac = mat % sab_fracs(j)
! If particle energy is greater than the highest energy for the S(a,b)
! table, don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) i_sab = 0
! If particle energy is greater than the highest energy for the
! S(a,b) table, then don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then
i_sab = 0
end if
! Increment position in i_sab_nuclides
j = j + 1
@ -102,7 +107,7 @@ contains
end if
end if
! ========================================================================
! ======================================================================
! CALCULATE MICROSCOPIC CROSS SECTION
! Determine microscopic cross sections for this nuclide
@ -110,13 +115,14 @@ contains
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E &
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, p % sqrtkT)
.or. p % sqrtkT /= micro_xs(i_nuclide) % last_sqrtkT &
.or. i_sab /= micro_xs(i_nuclide) % index_sab &
.or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, i_grid, &
p % sqrtkT, sab_frac)
end if
! ========================================================================
! ======================================================================
! ADD TO MACROSCOPIC CROSS SECTION
! Copy atom density of nuclide in material
@ -151,13 +157,15 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_log_union, sqrtkT, &
sab_frac)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer, intent(in) :: i_log_union ! index into logarithmic mapping array or
! material union energy grid
real(8), intent(in) :: sqrtkT ! Square root of kT, material dependent
real(8), intent(in) :: sqrtkT ! square root of kT, material dependent
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
logical :: use_mp ! true if XS can be calculated with windowed multipole
integer :: i_temp ! index for temperature
@ -261,9 +269,11 @@ contains
micro_xs(i_nuclide) % interp_factor = f
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % photon_prod = ZERO
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % thermal = ZERO
micro_xs(i_nuclide) % thermal_elastic = ZERO
micro_xs(i_nuclide) % photon_prod = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * xs % total(i_grid) &
@ -294,20 +304,21 @@ contains
end if
! Initialize sab treatment to false
micro_xs(i_nuclide) % index_sab = NONE
micro_xs(i_nuclide) % elastic_sab = ZERO
micro_xs(i_nuclide) % index_sab = NONE
micro_xs(i_nuclide) % sab_frac = ZERO
! Initialize URR probability table treatment to false
micro_xs(i_nuclide) % use_ptable = .false.
micro_xs(i_nuclide) % use_ptable = .false.
! If there is S(a,b) data for this nuclide, we need to do a few
! things. Since the total cross section was based on non-S(a,b) data, we
! need to correct it by subtracting the non-S(a,b) elastic cross section and
! then add back in the calculated S(a,b) elastic+inelastic cross section.
! If there is S(a,b) data for this nuclide, we need to set the sab_scatter
! and sab_elastic cross sections and correct the total and elastic cross
! sections.
if (i_sab > 0) call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT)
if (i_sab > 0) then
call calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT, sab_frac)
end if
! if the particle is in the unresolved resonance range and there are
! If the particle is in the unresolved resonance range and there are
! probability tables, we need to determine cross sections from the table
if (urr_ptables_on .and. nuc % urr_present .and. .not. use_mp) then
@ -318,7 +329,6 @@ contains
end if
micro_xs(i_nuclide) % last_E = E
micro_xs(i_nuclide) % last_index_sab = i_sab
micro_xs(i_nuclide) % last_sqrtkT = sqrtkT
end associate
@ -326,16 +336,16 @@ contains
!===============================================================================
! CALCULATE_SAB_XS determines the elastic and inelastic scattering
! cross-sections in the thermal energy range. These cross sections replace
! whatever data were taken from the normal Nuclide table.
! cross-sections in the thermal energy range. These cross sections replace a
! fraction of whatever data were taken from the normal Nuclide table.
!===============================================================================
subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT)
subroutine calculate_sab_xs(i_nuclide, i_sab, E, sqrtkT, sab_frac)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
real(8), intent(in) :: sqrtkT ! temperature
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
integer :: i_grid ! index on S(a,b) energy grid
integer :: i_temp ! temperature index
@ -364,7 +374,8 @@ contains
! Randomly sample between temperature i and i+1
f = (kT - sab_tables(i_sab) % kTs(i_temp)) / &
(sab_tables(i_sab) % kTs(i_temp + 1) - sab_tables(i_sab) % kTs(i_temp))
(sab_tables(i_sab) % kTs(i_temp + 1) &
- sab_tables(i_sab) % kTs(i_temp))
if (f > prn()) i_temp = i_temp + 1
end if
@ -425,16 +436,20 @@ contains
end if
end associate
! Store the S(a,b) cross sections.
micro_xs(i_nuclide) % thermal = sab_frac * (elastic + inelastic)
micro_xs(i_nuclide) % thermal_elastic = sab_frac * elastic
! Correct total and elastic cross sections
micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % elastic + inelastic + elastic
micro_xs(i_nuclide) % elastic = inelastic + elastic
micro_xs(i_nuclide) % total = micro_xs(i_nuclide) % total &
+ micro_xs(i_nuclide) % thermal &
- sab_frac * micro_xs(i_nuclide) % elastic
micro_xs(i_nuclide) % elastic = micro_xs(i_nuclide) % thermal &
+ (ONE - sab_frac) * micro_xs(i_nuclide) % elastic
! Store S(a,b) elastic cross section for sampling later
micro_xs(i_nuclide) % elastic_sab = elastic
! Save temperature index
! Save temperature index and thermal fraction
micro_xs(i_nuclide) % index_temp_sab = i_temp
micro_xs(i_nuclide) % sab_frac = sab_frac
end subroutine calculate_sab_xs

2
src/geometry.F90
View file

@ -283,6 +283,8 @@ contains
! Keep track of which instance of the cell the particle is in
p % cell_instance = offset + 1
else
p % cell_instance = 1
end if
! Save the material

15
src/initialize.F90
View file

@ -49,14 +49,23 @@ contains
integer, intent(in), optional :: intracomm ! MPI intracommunicator
! Copy the communicator to a new variable. This is done to avoid changing
! the signature of this subroutine.
! the signature of this subroutine. If MPI is being used but no communicator
! was passed, assume MPI_COMM_WORLD.
#ifdef MPI
#ifdef MPIF08
type(MPI_Comm), intent(in) :: comm ! MPI intracommunicator
comm % MPI_VAL = intracomm
if (present(intracomm)) then
comm % MPI_VAL = intracomm
else
comm = MPI_COMM_WORLD
end if
#else
integer :: comm
comm = intracomm
if (present(intracomm)) then
comm = intracomm
else
comm = MPI_COMM_WORLD
end if
#endif
#endif

44
src/input_xml.F90
View file

@ -2635,6 +2635,7 @@ contains
! table is indeed applied to multiple nuclides.
allocate(mat % sab_names(n_sab))
allocate(mat % i_sab_tables(n_sab))
allocate(mat % sab_fracs(n_sab))
do j = 1, n_sab
! Get pointer to S(a,b) table
@ -2648,6 +2649,13 @@ contains
name = trim(name)
mat % sab_names(j) = name
! Read the fraction of nuclei affected by this S(a,b) table
if (check_for_node(node_sab, "fraction")) then
call get_node_value(node_sab, "fraction", mat % sab_fracs(j))
else
mat % sab_fracs(j) = ONE
end if
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. library_dict % has_key(to_lower(name))) then
call fatal_error("Could not find S(a,b) table " // trim(name) &
@ -3440,6 +3448,15 @@ contains
! =======================================================================
! READ DATA FOR FILTERS
! Check if user is using old XML format and throw an error if so
if (check_for_node(node_tal, "filter")) then
call fatal_error("Tally filters must be specified independently of &
&tallies in a <filter> element. The <tally> element itself should &
&have a list of filters that apply, e.g., <filters>1 2</filters> &
&where 1 and 2 are the IDs of filters specified outside of &
&<tally>.")
end if
! Determine number of filters
if (check_for_node(node_tal, "filters")) then
n_filter = node_word_count(node_tal, "filters")
@ -5271,9 +5288,11 @@ contains
integer :: m ! position for sorting
integer :: temp_nuclide ! temporary value for sorting
integer :: temp_table ! temporary value for sorting
real(8) :: temp_frac ! temporary value for sorting
logical :: found
type(VectorInt) :: i_sab_tables
type(VectorInt) :: i_sab_nuclides
type(VectorInt) :: i_sab_tables
type(VectorInt) :: i_sab_nuclides
type(VectorReal) :: sab_fracs
do i = 1, size(materials)
! Skip materials with no S(a,b) tables
@ -5291,6 +5310,7 @@ contains
if (any(sab % nuclides == nuclides(mat % nuclide(j)) % name)) then
call i_sab_tables % push_back(mat % i_sab_tables(k))
call i_sab_nuclides % push_back(j)
call sab_fracs % push_back(mat % sab_fracs(k))
found = .true.
end if
end do FIND_NUCLIDE
@ -5304,17 +5324,34 @@ contains
end if
end do ASSIGN_SAB
! Make sure each nuclide only appears in one table.
do j = 1, i_sab_nuclides % size()
do k = j+1, i_sab_nuclides % size()
if (i_sab_nuclides % data(j) == i_sab_nuclides % data(k)) then
call fatal_error(trim( &
nuclides(mat % nuclide(i_sab_nuclides % data(j))) % name) &
// " in material " // trim(to_str(mat % id)) // " was found &
&in multiple S(a,b) tables. Each nuclide can only appear in &
&one S(a,b) table per material.")
end if
end do
end do
! Update i_sab_tables and i_sab_nuclides
deallocate(mat % i_sab_tables)
deallocate(mat % sab_fracs)
m = i_sab_tables % size()
allocate(mat % i_sab_tables(m))
allocate(mat % i_sab_nuclides(m))
allocate(mat % sab_fracs(m))
mat % i_sab_tables(:) = i_sab_tables % data(1:m)
mat % i_sab_nuclides(:) = i_sab_nuclides % data(1:m)
mat % sab_fracs(:) = sab_fracs % data(1:m)
! Clear entries in vectors for next material
call i_sab_tables % clear()
call i_sab_nuclides % clear()
call sab_fracs % clear()
! If there are multiple S(a,b) tables, we need to make sure that the
! entries in i_sab_nuclides are sorted or else they won't be applied
@ -5327,6 +5364,7 @@ contains
m = k
temp_nuclide = mat % i_sab_nuclides(k)
temp_table = mat % i_sab_tables(k)
temp_frac = mat % i_sab_tables(k)
MOVE_OVER: do
! Check if insertion value is greater than (m-1)th value
@ -5335,6 +5373,7 @@ contains
! Move values over until hitting one that's not larger
mat % i_sab_nuclides(m) = mat % i_sab_nuclides(m-1)
mat % i_sab_tables(m) = mat % i_sab_tables(m-1)
mat % sab_fracs(m) = mat % sab_fracs(m-1)
m = m - 1
! Exit if we've reached the beginning of the list
@ -5344,6 +5383,7 @@ contains
! Put the original value into its new position
mat % i_sab_nuclides(m) = temp_nuclide
mat % i_sab_tables(m) = temp_table
mat % sab_fracs(m) = temp_frac
end do SORT_SAB
end if

7
src/material_header.F90
View file

@ -26,10 +26,11 @@ module material_header
! Unionized energy grid information
integer, allocatable :: nuclide_grid_index(:,:) ! nuclide e_grid pointers
! S(a,b) data references
! S(a,b) data
integer :: n_sab = 0 ! number of S(a,b) tables
integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide
integer, allocatable :: i_sab_tables(:) ! index in sab_tables
real(8), allocatable :: sab_fracs(:) ! how often to use S(a,b)
! Temporary names read during initialization
character(20), allocatable :: names(:) ! isotope names
@ -48,6 +49,10 @@ module material_header
contains
!===============================================================================
! MATERIAL_SET_DENSITY sets the total density of a material in atom/b-cm.
!===============================================================================
function material_set_density(m, density, nuclides) result(err)
class(Material), intent(inout) :: m
real(8), intent(in) :: density

3
src/mgxs_data.F90
View file

@ -73,9 +73,6 @@ contains
! allocate arrays for MGXS storage and cross section cache
allocate(nuclides_MG(n_nuclides_total))
!$omp parallel
allocate(micro_xs(n_nuclides_total))
!$omp end parallel
! ==========================================================================
! READ ALL MGXS CROSS SECTION TABLES

46
src/nuclide_header.F90
View file

@ -104,34 +104,36 @@ module nuclide_header
end type Nuclide
!===============================================================================
! NUCLIDEMICROXS contains cached microscopic cross sections for a
! particular nuclide at the current energy
! NUCLIDEMICROXS contains cached microscopic cross sections for a particular
! nuclide at the current energy
!===============================================================================
type NuclideMicroXS
integer :: index_grid ! index on nuclide energy grid
integer :: index_temp ! temperature index for nuclide
real(8) :: last_E = ZERO ! last evaluated energy
real(8) :: interp_factor ! interpolation factor on nuc. energy grid
real(8) :: total ! microscropic total xs
real(8) :: elastic ! microscopic elastic scattering xs
real(8) :: absorption ! microscopic absorption xs
real(8) :: fission ! microscopic fission xs
real(8) :: nu_fission ! microscopic production xs
real(8) :: photon_prod ! microscopic photon production xs
! Microscopic cross sections in barns
real(8) :: total
real(8) :: elastic ! If sab_frac is not 1 or 0, then this value is
! averaged over bound and non-bound nuclei
real(8) :: absorption
real(8) :: fission
real(8) :: nu_fission
real(8) :: thermal ! Bound thermal elastic & inelastic scattering
real(8) :: thermal_elastic ! Bound thermal elastic scattering
real(8) :: photon_prod ! microscopic photon production xs
! Information for S(a,b) use
integer :: index_sab ! index in sab_tables (zero means no table)
integer :: last_index_sab = 0 ! index in sab_tables last used by this nuclide
integer :: index_temp_sab ! temperature index for sab_tables
real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table
! Indicies and factors needed to compute cross sections from the data tables
integer :: index_grid ! Index on nuclide energy grid
integer :: index_temp ! Temperature index for nuclide
real(8) :: interp_factor ! Interpolation factor on nuc. energy grid
integer :: index_sab = NONE ! Index in sab_tables
integer :: index_temp_sab ! Temperature index for sab_tables
real(8) :: sab_frac ! Fraction of atoms affected by S(a,b)
logical :: use_ptable ! In URR range with probability tables?
! Information for URR probability table use
logical :: use_ptable ! in URR range with probability tables?
! Information for Doppler broadening
! Energy and temperature last used to evaluate these cross sections. If
! these values have changed, then the cross sections must be re-evaluated.
real(8) :: last_E = ZERO ! Last evaluated energy
real(8) :: last_sqrtkT = ZERO ! Last temperature in sqrt(Boltzmann
! constant * temperature (eV))
! constant * temperature (eV))
end type NuclideMicroXS
!===============================================================================

4
src/particle_restart.F90
View file

@ -32,6 +32,8 @@ contains
! Set verbosity high
verbosity = 10
allocate(micro_xs(n_nuclides_total))
! Initialize the particle to be tracked
call p % initialize()
@ -57,6 +59,8 @@ contains
! Write output if particle made it
call print_particle(p)
deallocate(micro_xs)
end subroutine run_particle_restart
!===============================================================================

62
src/physics.F90
View file

@ -644,6 +644,7 @@ contains
real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering
real(8) :: phi ! azimuthal angle for iso-in-lab scattering
real(8) :: kT ! temperature in eV
logical :: sampled ! whether or not a reaction type has been sampled
type(Nuclide), pointer :: nuc
! copy incoming direction
@ -657,36 +658,43 @@ contains
! For tallying purposes, this routine might be called directly. In that
! case, we need to sample a reaction via the cutoff variable
prob = ZERO
cutoff = prn() * (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption)
sampled = .false.
prob = prob + micro_xs(i_nuclide) % elastic
prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal
if (prob > cutoff) then
! =======================================================================
! ELASTIC SCATTERING
if (micro_xs(i_nuclide) % index_sab /= NONE) then
! S(a,b) scattering
call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, &
p % E, p % coord(1) % uvw, p % mu)
! NON-S(A,B) ELASTIC SCATTERING
! Determine temperature
if (nuc % mp_present) then
kT = p % sqrtkT**2
else
! Determine temperature
if (nuc % mp_present) then
kT = p % sqrtkT**2
else
kT = nuc % kTs(micro_xs(i_nuclide) % index_temp)
end if
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, nuc % reactions(1), kT, &
p % E, p % coord(1) % uvw, p % mu, p % wgt)
kT = nuc % kTs(micro_xs(i_nuclide) % index_temp)
end if
p % event_MT = ELASTIC
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, nuc % reactions(1), kT, p % E, &
p % coord(1) % uvw, p % mu, p % wgt)
else
p % event_MT = ELASTIC
sampled = .true.
end if
prob = micro_xs(i_nuclide) % elastic
if (prob > cutoff .and. .not. sampled) then
! =======================================================================
! S(A,B) SCATTERING
call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, p % E, &
p % coord(1) % uvw, p % mu)
p % event_MT = ELASTIC
sampled = .true.
end if
if (.not. sampled) then
! =======================================================================
! INELASTIC SCATTERING
@ -708,7 +716,7 @@ contains
if (rx % MT == N_FISSION .or. rx % MT == N_F .or. rx % MT == N_NF &
.or. rx % MT == N_2NF .or. rx % MT == N_3NF) cycle
! some materials have gas production cross sections with MT > 200 that
! Some materials have gas production cross sections with MT > 200 that
! are duplicates. Also MT=4 is total level inelastic scattering which
! should be skipped
if (rx % MT >= 200 .or. rx % MT == N_LEVEL) cycle
@ -726,24 +734,24 @@ contains
! Perform collision physics for inelastic scattering
call inelastic_scatter(nuc, nuc%reactions(i), p)
p % event_MT = nuc%reactions(i)%MT
p % event_MT = nuc % reactions(i) % MT
end if
! Set event component
p % event = EVENT_SCATTER
! sample new outgoing angle for isotropic in lab scattering
! Sample new outgoing angle for isotropic in lab scattering
if (materials(p % material) % p0(i_nuc_mat)) then
! sample isotropic-in-lab outgoing direction
! Sample isotropic-in-lab outgoing direction
uvw_new(1) = TWO * prn() - ONE
phi = TWO * PI * prn()
uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1))
uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1))
p % mu = dot_product(uvw_old, uvw_new)
! change direction of particle
! Change direction of particle
p % coord(1) % uvw = uvw_new
end if
@ -879,8 +887,8 @@ contains
associate (sab => sab_tables(i_sab) % data(i_temp))
! Determine whether inelastic or elastic scattering will occur
if (prn() < micro_xs(i_nuclide) % elastic_sab / &
micro_xs(i_nuclide) % elastic) then
if (prn() < micro_xs(i_nuclide) % thermal_elastic / &
micro_xs(i_nuclide) % thermal) then
! elastic scattering
! Get index and interpolation factor for elastic grid

3
src/relaxng/materials.rnc
View file

@ -28,7 +28,8 @@ element materials {
}* &
element sab {
(element name { xsd:string } | attribute name { xsd:string })
(element name { xsd:string } | attribute name { xsd:string }) &
(element fraction { xsd:double } | attribute fraction { xsd:double })?
}*
}+ &

28
src/relaxng/materials.rng
View file

@ -119,14 +119,26 @@
</zeroOrMore>
<zeroOrMore>
<element name="sab">
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<interleave>
<choice>
<element name="name">
<data type="string"/>
</element>
<attribute name="name">
<data type="string"/>
</attribute>
</choice>
<optional>
<choice>
<element name="fraction">
<data type="double"/>
</element>
<attribute name="fraction">
<data type="double"/>
</attribute>
</choice>
</optional>
</interleave>
</element>
</zeroOrMore>
</interleave>

11
tests/run_tests.py
View file

@ -188,6 +188,10 @@ class Test(object):
build_str += "-Dopenmp=OFF "
if self.coverage:
build_str += "-Dcoverage=ON "
if self.phdf5:
build_str += "-DHDF5_PREFER_PARALLEL=ON "
else:
build_str += "-DHDF5_PREFER_PARALLEL=OFF "
self.build_opts = build_str
return self.build_opts
@ -214,6 +218,9 @@ class Test(object):
# Runs cmake when in non-script mode
def run_cmake(self):
build_opts = self.build_opts.split()
self.cmake += build_opts
os.environ['FC'] = self.fc
os.environ['CC'] = self.cc
os.environ['CXX'] = self.cxx
@ -221,10 +228,10 @@ class Test(object):
os.environ['MPI_DIR'] = MPI_DIR
if self.phdf5:
os.environ['HDF5_ROOT'] = PHDF5_DIR
self.cmake.append('-DHDF5_PREFER_PARALLEL=ON')
else:
os.environ['HDF5_ROOT'] = HDF5_DIR
build_opts = self.build_opts.split()
self.cmake += build_opts
self.cmake.append('-DHDF5_PREFER_PARALLEL=OFF')
rc = call(self.cmake)
if rc != 0:
self.success = False

15
tests/test_salphabeta/geometry.xml
View file

@ -1,15 +0,0 @@
<?xml version="1.0"?>
<geometry>
<surface id="1" type="x-plane" coeffs="-10" boundary="vacuum"/>
<surface id="2" type="x-plane" coeffs="-5" />
<surface id="3" type="x-plane" coeffs="0" />
<surface id="4" type="x-plane" coeffs="5" />
<surface id="5" type="x-plane" coeffs="10" boundary="vacuum"/>
<cell id="1" material="1" region="1 -2" />
<cell id="2" material="2" region="2 -3" />
<cell id="3" material="3" region="3 -4" />
<cell id="4" material="4" region="4 -5" />
</geometry>

60
tests/test_salphabeta/inputs_true.dat Normal file
View file

@ -0,0 +1,60 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="1 -2" universe="1" />
<cell id="2" material="2" region="2 -3" universe="1" />
<cell id="3" material="3" region="3 -4" universe="1" />
<cell id="4" material="4" region="4 -5" universe="1" />
<surface boundary="vacuum" coeffs="-10" id="1" type="x-plane" />
<surface coeffs="-5" id="2" type="x-plane" />
<surface coeffs="0" id="3" type="x-plane" />
<surface coeffs="5" id="4" type="x-plane" />
<surface boundary="vacuum" coeffs="10" id="5" type="x-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material id="1">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
<nuclide ao="1.0" name="H1" />
<sab fraction="0.5" name="c_H_in_H2O" />
</material>
<material id="2">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
<nuclide ao="1.0" name="C0" />
<sab name="c_Graphite" />
</material>
<material id="3">
<density units="g/cc" value="4.5" />
<nuclide ao="1.0" name="U235" />
<nuclide ao="1.0" name="Be9" />
<nuclide ao="1.0" name="O16" />
<sab name="c_Be_in_BeO" />
<sab name="c_O_in_BeO" />
</material>
<material id="4">
<density units="g/cm3" value="5.90168" />
<nuclide ao="0.3" name="H1" />
<nuclide ao="0.15" name="Zr90" />
<nuclide ao="0.1" name="Zr91" />
<nuclide ao="0.1" name="Zr92" />
<nuclide ao="0.05" name="Zr94" />
<nuclide ao="0.05" name="Zr96" />
<nuclide ao="0.1" name="U235" />
<nuclide ao="0.15" name="U238" />
<sab name="c_Zr_in_ZrH" />
<sab name="c_H_in_ZrH" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>1000</particles>
<batches>5</batches>
<inactive>0</inactive>
<source strength="1.0">
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

41
tests/test_salphabeta/materials.xml
View file

@ -1,41 +0,0 @@
<?xml version="1.0"?>
<materials>
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
<nuclide name="H1" ao="1.0" />
<sab name="c_H_in_H2O" />
</material>
<material id="2">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
<nuclide name="C0" ao="1.0" />
<sab name="c_Graphite" />
</material>
<material id="3">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
<nuclide name="Be9" ao="1.0" />
<nuclide name="O16" ao="1.0" />
<sab name="c_Be_in_BeO" />
<sab name="c_O_in_BeO" />
</material>
<material id="4">
<density value="5.90168" units="g/cm3" />
<nuclide name="H1" ao="0.3" />
<nuclide name="Zr90" ao="0.15" />
<nuclide name="Zr91" ao="0.1" />
<nuclide name="Zr92" ao="0.1" />
<nuclide name="Zr94" ao="0.05" />
<nuclide name="Zr96" ao="0.05" />
<nuclide name="U235" ao="0.1" />
<nuclide name="U238" ao="0.15" />
<sab name="c_Zr_in_ZrH" />
<sab name="c_H_in_ZrH" />
</material>
</materials>

2
tests/test_salphabeta/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
8.544160E-01 1.274133E-02
8.447580E-01 1.806149E-02

15
tests/test_salphabeta/settings.xml
View file

@ -1,15 +0,0 @@
<?xml version="1.0"?>
<settings>
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

77
tests/test_salphabeta/test_salphabeta.py
View file

@ -3,9 +3,82 @@
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import TestHarness
from testing_harness import PyAPITestHarness
import openmc
import openmc.model
def make_model():
model = openmc.model.Model()
# Materials
m1 = openmc.Material()
m1.set_density('g/cc', 4.5)
m1.add_nuclide('U235', 1.0)
m1.add_nuclide('H1', 1.0)
m1.add_s_alpha_beta('c_H_in_H2O', fraction=0.5)
m2 = openmc.Material()
m2.set_density('g/cc', 4.5)
m2.add_nuclide('U235', 1.0)
m2.add_nuclide('C0', 1.0)
m2.add_s_alpha_beta('c_Graphite')
m3 = openmc.Material()
m3.set_density('g/cc', 4.5)
m3.add_nuclide('U235', 1.0)
m3.add_nuclide('Be9', 1.0)
m3.add_nuclide('O16', 1.0)
m3.add_s_alpha_beta('c_Be_in_BeO')
m3.add_s_alpha_beta('c_O_in_BeO')
m4 = openmc.Material()
m4.set_density('g/cm3', 5.90168)
m4.add_nuclide('H1', 0.3)
m4.add_nuclide('Zr90', 0.15)
m4.add_nuclide('Zr91', 0.1)
m4.add_nuclide('Zr92', 0.1)
m4.add_nuclide('Zr94', 0.05)
m4.add_nuclide('Zr96', 0.05)
m4.add_nuclide('U235', 0.1)
m4.add_nuclide('U238', 0.15)
m4.add_s_alpha_beta('c_Zr_in_ZrH')
m4.add_s_alpha_beta('c_H_in_ZrH')
model.materials += [m1, m2, m3, m4]
# Geometry
x0 = openmc.XPlane(x0=-10, boundary_type='vacuum')
x1 = openmc.XPlane(x0=-5)
x2 = openmc.XPlane(x0=0)
x3 = openmc.XPlane(x0=5)
x4 = openmc.XPlane(x0=10, boundary_type='vacuum')
root_univ = openmc.Universe()
surfs = (x0, x1, x2, x3, x4)
mats = (m1, m2, m3, m4)
cells = []
for i in range(4):
cell = openmc.Cell()
cell.region = +surfs[i] & -surfs[i+1]
cell.fill = mats[i]
root_univ.add_cell(cell)
model.geometry.root_universe = root_univ
# Settings
model.settings.batches = 5
model.settings.inactive = 0
model.settings.particles = 1000
model.settings.source = openmc.Source(space=openmc.stats.Box(
[-4, -4, -4], [4, 4, 4]))
return model
if __name__ == '__main__':
harness = TestHarness('statepoint.10.h5')
model = make_model()
harness = PyAPITestHarness('statepoint.5.h5', model)
harness.main()

11
tests/travis.sh
View file

@ -1,11 +0,0 @@
#!/bin/sh
set -ev
# Run all debug tests
./check_source.py
if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then
./run_tests.py -C "^hdf5-debug$|^omp-hdf5-debug|^mpi-hdf5-debug|^phdf5-debug$|^phdf5-omp-debug$" -j 2 -s
else
./run_tests.py -C "^hdf5-debug$" -j 2
fi

39
tests/travis_install.sh
View file

@ -1,39 +0,0 @@
#!/bin/bash
set -ev
# Build MPICH
if [[ ! -e $HOME/mpich_install/bin/mpiexec ]]; then
wget -q http://www.mpich.org/static/downloads/3.1.3/mpich-3.1.3.tar.gz
tar -xzvf mpich-3.1.3.tar.gz >/dev/null 2>&1
cd mpich-3.1.3
./configure --prefix=$HOME/mpich_install -q
make -j 2 >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
fi
# Build PHDF5
if [[ ! -e $HOME/phdf5_install/bin/h5pfc ]]; then
wget -q http://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8.15/src/hdf5-1.8.15.tar.gz
tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1
mv hdf5-1.8.15 phdf5-1.8.15; cd phdf5-1.8.15
CC=$HOME/mpich_install/bin/mpicc FC=$HOME/mpich_install/bin/mpif90 \
./configure \
--prefix=$HOME/phdf5_install -q --enable-fortran \
--enable-fortran2003 --enable-parallel
make -j 2 >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
fi
# Build HDF5
if [[ ! -e $HOME/hdf5_install/bin/h5fc ]]; then
tar -xzvf hdf5-1.8.15.tar.gz >/dev/null 2>&1
cd hdf5-1.8.15
CC=gcc FC=gfortran ./configure --prefix=$HOME/hdf5_install -q \
--enable-fortran --enable-fortran2003
make -j 2 >/dev/null 2>&1
make install >/dev/null 2>&1
cd ..
fi