diff --git a/.travis.yml b/.travis.yml index c1d3c86b2..64fbdc9df 100644 --- a/.travis.yml +++ b/.travis.yml @@ -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 .. diff --git a/CMakeLists.txt b/CMakeLists.txt index de01e2e54..8364476cd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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) diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst index f3cdbb840..9bb1b2c66 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -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] diff --git a/docs/source/io_formats/materials.rst b/docs/source/io_formats/materials.rst index fc57280f8..5f0aa3a5e 100644 --- a/docs/source/io_formats/materials.rst +++ b/docs/source/io_formats/materials.rst @@ -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 diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst new file mode 100644 index 000000000..af00363cc --- /dev/null +++ b/docs/source/pythonapi/capi.rst @@ -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 diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 5492f362d..bb99dd470 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -47,5 +47,6 @@ Modules mgxs model data + capi examples openmoc diff --git a/openmc/__init__.py b/openmc/__init__.py index 96bcaf6d2..a0fe26491 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -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' diff --git a/openmc/capi.py b/openmc/capi.py index a4ac38c32..bb9241d7d 100644 --- a/openmc/capi.py +++ b/openmc/capi.py @@ -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 " diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index 80817bf71..b6fce5e3e 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -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 diff --git a/openmc/material.py b/openmc/material.py index 762332584..03ef7fbc4 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -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 diff --git a/src/api.F90 b/src/api.F90 index a7a7e3ab4..b00193fcf 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -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) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 37a376845..2a787f788 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -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 diff --git a/src/geometry.F90 b/src/geometry.F90 index fbbe75a18..444456232 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/initialize.F90 b/src/initialize.F90 index 8204307c9..7031da6b4 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 5273a2180..e6930247f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 element. The element itself should & + &have a list of filters that apply, e.g., 1 2 & + &where 1 and 2 are the IDs of filters specified outside of & + &.") + 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 diff --git a/src/material_header.F90 b/src/material_header.F90 index 413296240..4ab59b4da 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -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 diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index a43d45c0e..ad20aca83 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -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 diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index ab937809e..d2a21bfd9 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -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 !=============================================================================== diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index ccf546dd9..4438c48c0 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -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 !=============================================================================== diff --git a/src/physics.F90 b/src/physics.F90 index 6923a5386..46a9dad7d 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 diff --git a/src/relaxng/materials.rnc b/src/relaxng/materials.rnc index 4ad3ea0f0..d55656536 100644 --- a/src/relaxng/materials.rnc +++ b/src/relaxng/materials.rnc @@ -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 })? }* }+ & diff --git a/src/relaxng/materials.rng b/src/relaxng/materials.rng index 0303d361e..0fa81c11b 100644 --- a/src/relaxng/materials.rng +++ b/src/relaxng/materials.rng @@ -119,14 +119,26 @@ - - - - - - - - + + + + + + + + + + + + + + + + + + + + diff --git a/tests/run_tests.py b/tests/run_tests.py index 2d9112e92..cd6c02607 100755 --- a/tests/run_tests.py +++ b/tests/run_tests.py @@ -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 diff --git a/tests/test_salphabeta/geometry.xml b/tests/test_salphabeta/geometry.xml deleted file mode 100644 index 2b978b914..000000000 --- a/tests/test_salphabeta/geometry.xml +++ /dev/null @@ -1,15 +0,0 @@ - - - - - - - - - - - - - - - diff --git a/tests/test_salphabeta/inputs_true.dat b/tests/test_salphabeta/inputs_true.dat new file mode 100644 index 000000000..02e6813f0 --- /dev/null +++ b/tests/test_salphabeta/inputs_true.dat @@ -0,0 +1,60 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 5 + 0 + + + -4 -4 -4 4 4 4 + + + diff --git a/tests/test_salphabeta/materials.xml b/tests/test_salphabeta/materials.xml deleted file mode 100644 index bfe0a6224..000000000 --- a/tests/test_salphabeta/materials.xml +++ /dev/null @@ -1,41 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/tests/test_salphabeta/results_true.dat b/tests/test_salphabeta/results_true.dat index bdfb98990..2a99303f6 100644 --- a/tests/test_salphabeta/results_true.dat +++ b/tests/test_salphabeta/results_true.dat @@ -1,2 +1,2 @@ k-combined: -8.544160E-01 1.274133E-02 +8.447580E-01 1.806149E-02 diff --git a/tests/test_salphabeta/settings.xml b/tests/test_salphabeta/settings.xml deleted file mode 100644 index 70b4e802f..000000000 --- a/tests/test_salphabeta/settings.xml +++ /dev/null @@ -1,15 +0,0 @@ - - - - eigenvalue - 10 - 5 - 1000 - - - - -4 -4 -4 4 4 4 - - - - diff --git a/tests/test_salphabeta/test_salphabeta.py b/tests/test_salphabeta/test_salphabeta.py index b04fcc6eb..600c70332 100644 --- a/tests/test_salphabeta/test_salphabeta.py +++ b/tests/test_salphabeta/test_salphabeta.py @@ -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() diff --git a/tests/travis.sh b/tests/travis.sh deleted file mode 100755 index 443e951cc..000000000 --- a/tests/travis.sh +++ /dev/null @@ -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 diff --git a/tests/travis_install.sh b/tests/travis_install.sh deleted file mode 100755 index 2e35c4f65..000000000 --- a/tests/travis_install.sh +++ /dev/null @@ -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