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

Merge pull request #1 from Mikolaj-A-Kowalski/cam_WP7_task3

Browse source 
Cam wp7 task3
This commit is contained in:
Simon Richards 2020-02-26 20:43:28 +00:00 • committed by GitHub
commit d17b38ba5e
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
39 changed files with 2006 additions and 919 deletions
Download patch file Download diff file Expand all files Collapse all files

22
CMakeLists.txt
View file

@ -9,6 +9,11 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin)
# Set module path
set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules)
# Allow user to specify <project>_ROOT variables
if (NOT (CMAKE_VERSION VERSION_LESS 3.12))
cmake_policy(SET CMP0074 NEW)
endif()
#===============================================================================
# Command line options
#===============================================================================
@ -35,18 +40,18 @@ endif()
#===============================================================================
if(dagmc)
find_package(DAGMC REQUIRED)
add_library(dagmc-imported INTERFACE IMPORTED)
link_directories(${DAGMC_LIBRARY_DIRS})
target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES})
target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS})
else()
set(DAGMC_FOUND false)
endif()
#===============================================================================
# HDF5 for binary output
#===============================================================================
# Allow user to specify HDF5_ROOT
if (NOT (CMAKE_VERSION VERSION_LESS 3.12))
cmake_policy(SET CMP0074 NEW)
endif()
# Unfortunately FindHDF5.cmake will always prefer a serial HDF5 installation
# over a parallel installation if both appear on the user's PATH. To get around
# this, we check for the environment variable HDF5_ROOT and if it exists, use it
@ -357,8 +362,7 @@ target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES
if(dagmc)
target_compile_definitions(libopenmc PRIVATE DAGMC)
target_link_libraries(libopenmc ${DAGMC_LIBRARIES})
target_include_directories(libopenmc PRIVATE ${DAGMC_INCLUDE_DIRS})
target_link_libraries(libopenmc dagmc-imported)
endif()
#===============================================================================
@ -388,6 +392,8 @@ add_custom_command(TARGET libopenmc POST_BUILD
# Install executable, scripts, manpage, license
#===============================================================================
configure_file(cmake/OpenMCConfig.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" @ONLY)
set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC)
install(TARGETS openmc libopenmc faddeeva
EXPORT openmc-targets
@ -401,7 +407,7 @@ install(EXPORT openmc-targets
DESTINATION ${INSTALL_CONFIGDIR})
install(DIRECTORY src/relaxng DESTINATION ${CMAKE_INSTALL_DATADIR}/openmc)
install(FILES cmake/OpenMCConfig.cmake DESTINATION ${INSTALL_CONFIGDIR})
install(FILES "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" DESTINATION ${INSTALL_CONFIGDIR})
install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1)
install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright)
install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})

10
cmake/Modules/FindDAGMC.cmake
View file

@ -12,7 +12,9 @@ find_path(DAGMC_CMAKE_CONFIG NAMES DAGMCConfig.cmake
PATHS ENV LD_LIBRARY_PATH
PATH_SUFFIXES lib Lib cmake lib/cmake
NO_DEFAULT_PATH)
message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}")
include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake)
if(DAGMC_CMAKE_CONFIG)
message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}")
include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake)
else()
message(WARNING "Cound not find DAGMC")
endif()

9
cmake/OpenMCConfig.cmake → cmake/OpenMCConfig.cmake.in
View file

@ -5,6 +5,15 @@ find_package(gsl-lite REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../gsl-lite)
find_package(pugixml REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../pugixml)
find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl)
find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor)
if(@DAGMC_FOUND@)
find_package(DAGMC REQUIRED HINTS @DAGMC_LIBRARY_DIRS@)
add_library(dagmc-imported INTERFACE IMPORTED)
target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES})
foreach(dir ${DAGMC_LIBRARY_DIRS})
target_link_libraries(dagmc-imported INTERFACE "-L${dir}")
endforeach()
target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS})
endif()
if(NOT TARGET OpenMC::libopenmc)
include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake")

72
docs/source/io_formats/settings.rst
View file

@ -142,7 +142,7 @@ materials in the problem and is specified using a :ref:`mesh_element`.
----------------------------
Determines whether to use event-based parallelism instead of the default
history-based parallelism.
history-based parallelism.
*Default*: false
@ -459,6 +459,15 @@ attributes/sub-elements:
*Default*: None
:library:
If this attribute is given, it indicates that the source is to be
instantiated from an externally compiled source function. This source can be
as complex as is required to define the source for your problem. The only
requirement is that there is a function called ``sample_source()``. More
documentation on how to build sources can be found in :ref:`custom_source`.
*Default*: None
:space:
An element specifying the spatial distribution of source sites. This element
has the following attributes:
@ -591,6 +600,67 @@ attributes/sub-elements:
*Default*: false
.. _custom_source:
Custom Sources
++++++++++++++
It is often the case that one may wish to simulate a complex source
distribution, which may include physics not present within OpenMC or to be phase
space complex. It is possible to define a complex source with an externally
defined source function that is loaded at runtime. A simple example source is
shown below.
.. code-block:: c++
#include "openmc/random_lcg.h"
#include "openmc/source.h"
#include "openmc/particle.h"
// you must have external C linkage here
extern "C" openmc::Particle::Bank sample_source(uint64_t* seed) {
openmc::Particle::Bank particle;
// weight
particle.particle = openmc::Particle::Type::neutron;
particle.wgt = 1.0;
// position
double angle = 2.0 * M_PI * openmc::prn(seed);
double radius = 3.0;
particle.r.x = radius * std::cos(angle);
particle.r.y = radius * std::sin(angle);
particle.r.z = 0.0;
// angle
particle.u = {1.0, 0.0, 0.0};
particle.E = 14.08e6;
particle.delayed_group = 0;
return particle;
}
The above source, creates 14.08 MeV neutrons, with an istropic direction
vector but distributed in a ring with a 3 cm radius. This routine is
not particularly complex, but should serve as an example upon which to build
more complicated sources.
.. note:: The function signature must be declared to be extern "C".
.. note:: You should only use the openmc::prn() random number generator
In order to build your external source, you will need to link it against the
OpenMC shared library. This can be done by writing a CMakeLists.txt file:
.. code-block:: cmake
cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
project(openmc_sources CXX)
add_library(source SHARED source_ring.cpp)
find_package(OpenMC REQUIRED HINTS <path to openmc>)
target_link_libraries(source OpenMC::libopenmc)
After running ``cmake`` and ``make``, you will have a libsource.so (or .dylib)
file in your build directory. Setting the :attr:`openmc.Source.library`
attribute to the path of this shared library will indicate that it should be
used for sampling source particles at runtime.
.. _univariate:
Univariate Probability Distributions

2
docs/source/usersguide/install.rst
View file

@ -404,7 +404,7 @@ to install the Python package in :ref:`"editable" mode <devguide_editable>`.
Prerequisites
-------------
The Python API works with Python 3.4+. In addition to Python itself, the API
The Python API works with Python 3.5+. In addition to Python itself, the API
relies on a number of third-party packages. All prerequisites can be installed
using Conda_ (recommended), pip_, or through the package manager in most Linux
distributions. To run simulations in parallel using MPI, it is recommended to

14
docs/source/usersguide/materials.rst
View file

@ -55,6 +55,20 @@ following would add 3.2% enriched uranium to a material::
In addition to U235 and U238, concentrations of U234 and U236 will be present
and are determined through a correlation based on measured data.
It is also possible to perform enrichment of any element that is composed
of two naturally-occurring isotopes (e.g. Li, B) in terms of atomic percent.
To invoke this, provide the additional argument `enrichment_target` to
:meth:`Material.add_element`. For example the following would enrich B10
to 30ao%::
mat.add_element('B', 1.0, enrichment=30.0, enrichment_target='B10')
In order to enrich an isotope in terms of mass percent (wo%), provide the extra
argument `enrichment_type`. For example the following would enrich Li6 to 15wo%::
mat.add_element('Li', 1.0, enrichment=15.0, enrichment_target='Li6',
enrichment_type='wo')
Often, cross section libraries don't actually have all naturally-occurring
isotopes for a given element. For example, in ENDF/B-VII.1, cross section
evaluations are given for O16 and O17 but not for O18. If OpenMC is aware of

8
examples/xml/custom_source/CMakeLists.txt Normal file
View file

@ -0,0 +1,8 @@
cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
project(openmc_sources CXX)
add_library(source SHARED source_ring.cpp)
find_package(OpenMC REQUIRED)
if (OpenMC_FOUND)
message(STATUS "Found OpenMC: ${OpenMC_DIR}")
endif()
target_link_libraries(source OpenMC::libopenmc)

15
examples/xml/custom_source/geometry.xml Normal file
View file

@ -0,0 +1,15 @@
<?xml version="1.0"?>
<geometry>
<!-- Definition of Cells -->
<cell id="1" universe="0" fill="37" region="-2" />
<cell id="100" universe="37" material="40" region="-1" />
<cell id="101" universe="37" material="41" region="1" />
<cell id="2" universe="0" material="41" region="2 -3" />
<!-- Defition of Surfaces -->
<surface id="1" type="z-cylinder" coeffs="0 0 7" />
<surface id="2" type="z-cylinder" coeffs="0 0 9" />
<surface id="3" type="z-cylinder" coeffs="0 0 11" boundary="vacuum" />
</geometry>

16
examples/xml/custom_source/materials.xml Normal file
View file

@ -0,0 +1,16 @@
<?xml version="1.0"?>
<materials>
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
</material>
<material id="41">
<density value="1.0" units="g/cc" />
<nuclide name="H1" ao="2.0" />
<nuclide name="O16" ao="1.0" />
<sab name="c_H_in_H2O"/>
</material>
</materials>

14
examples/xml/custom_source/settings.xml Normal file
View file

@ -0,0 +1,14 @@
<?xml version="1.0"?>
<settings>
<run_mode>fixed source</run_mode>
<batches>10</batches>
<inactive>0</inactive>
<particles>100000</particles>
<!-- Starting source -->
<source>
<library>build/libsource.so</library>
</source>
</settings>

26
examples/xml/custom_source/source_ring.cpp Normal file
View file

@ -0,0 +1,26 @@
#include <cmath> // for M_PI
#include "openmc/random_lcg.h"
#include "openmc/source.h"
#include "openmc/particle.h"
// you must have external C linkage here otherwise
// dlopen will not find the file
extern "C" openmc::Particle::Bank sample_source(uint64_t* seed)
{
openmc::Particle::Bank particle;
// wgt
particle.particle = openmc::Particle::Type::neutron;
particle.wgt = 1.0;
// position
double angle = 2. * M_PI * openmc::prn(seed);
double radius = 3.0;
particle.r.x = radius * std::cos(angle);
particle.r.y = radius * std::sin(angle);
particle.r.z = 0.0;
// angle
particle.u = {1.0, 0.0, 0.0};
particle.E = 14.08e6;
particle.delayed_group = 0;
return particle;
}

17
examples/xml/custom_source/tallies.xml Normal file
View file

@ -0,0 +1,17 @@
<?xml version="1.0"?>
<tallies>
<filter id="1" type="cell">
<bins>100</bins>
</filter>
<filter id="2" type="energy">
<bins>0 20.0e6</bins>
</filter>
<tally id="3">
<filters>1 2 </filters>
<scores>flux</scores>
</tally>
</tallies>

1
include/openmc/settings.h
View file

@ -60,6 +60,7 @@ extern std::string path_input; //!< directory where main .xml files r
extern std::string path_output; //!< directory where output files are written
extern std::string path_particle_restart; //!< path to a particle restart file
extern std::string path_source;
extern std::string path_source_library; //!< path to the source shared object
extern std::string path_sourcepoint; //!< path to a source file
extern "C" std::string path_statepoint; //!< path to a statepoint file

3
include/openmc/source.h
View file

@ -68,6 +68,9 @@ Particle::Bank sample_external_source(uint64_t* seed);
//! Fill source bank at end of generation for fixed source simulations
void fill_source_bank_fixedsource();
//! Fill source bank at the end of a generation for dlopen based source simulation
void fill_source_bank_custom_source();
void free_memory_source();
} // namespace openmc

1
openmc/data/data.py
View file

@ -10,6 +10,7 @@ from warnings import warn
# pp. 293-306 (2013). The "representative isotopic abundance" values from
# column 9 are used except where an interval is given, in which case the
# "best measurement" is used.
# Note that the abundances are given as atomic fractions!
NATURAL_ABUNDANCE = {
'H1': 0.99984426, 'H2': 0.00015574, 'He3': 0.000002,
'He4': 0.999998, 'Li6': 0.07589, 'Li7': 0.92411,

5
openmc/data/photon.py
View file

@ -784,8 +784,9 @@ class IncidentPhoton(EqualityMixin):
sub_group = shell_group.create_group(key)
# Write atomic relaxation
if key in self.atomic_relaxation.subshells:
self.atomic_relaxation.to_hdf5(sub_group, key)
if self.atomic_relaxation is not None:
if key in self.atomic_relaxation.subshells:
self.atomic_relaxation.to_hdf5(sub_group, key)
else:
continue

8
openmc/deplete/nuclide.py
View file

@ -533,6 +533,7 @@ class FissionYield(Mapping):
return zip(self.products, self.yields)
def __add__(self, other):
"""Add one set of fission yields to this set, return new yields"""
if not isinstance(other, FissionYield):
return NotImplemented
new = FissionYield(self.products, self.yields.copy())
@ -550,12 +551,14 @@ class FissionYield(Mapping):
return self + other
def __imul__(self, scalar):
"""Scale these fission yields by a real scalar"""
if not isinstance(scalar, Real):
return NotImplemented
self.yields *= scalar
return self
def __mul__(self, scalar):
"""Return a new set of yields scaled by a real scalar"""
if not isinstance(scalar, Real):
return NotImplemented
new = FissionYield(self.products, self.yields.copy())
@ -568,3 +571,8 @@ class FissionYield(Mapping):
def __repr__(self):
return "<{} containing {} products and yields>".format(
self.__class__.__name__, len(self))
# Avoid greedy numpy operations like np.float64 * fission_yield
# converting this to an array on the fly. Force __rmul__ and
# __radd__. See issue #1492
__array_ufunc__ = None

139
openmc/element.py
View file

@ -4,6 +4,7 @@ import os
from xml.etree import ElementTree as ET
import openmc.checkvalue as cv
from numbers import Real
from openmc.data import NATURAL_ABUNDANCE, atomic_mass
@ -35,6 +36,7 @@ class Element(str):
return self
def expand(self, percent, percent_type, enrichment=None,
enrichment_target=None, enrichment_type='ao',
cross_sections=None):
"""Expand natural element into its naturally-occurring isotopes.
@ -52,9 +54,16 @@ class Element(str):
percent_type : {'ao', 'wo'}
'ao' for atom percent and 'wo' for weight percent
enrichment : float, optional
Enrichment for U235 in weight percent. For example, input 4.95 for
4.95 weight percent enriched U. Default is None
(natural composition).
Enrichment of an enrichment_taget nuclide in percent (ao or wo).
If enrichment_taget is not supplied then it is enrichment for U235
in weight percent. For example, input 4.95 for 4.95 weight percent
enriched U. Default is None (natural composition).
Value must be in <0;100>
enrichment_target: str, optional
Single nuclide name to enrich from a natural composition e.g. O16
enrichment_type: {'ao', 'wo'}, optional
'ao' for enrichment as atom percent and 'wo' for weight percent.
Default is 'ao'
cross_sections : str, optional
Location of cross_sections.xml file. Default is None.
@ -65,16 +74,49 @@ class Element(str):
is a tuple consisting of a nuclide string, the atom/weight percent,
and the string 'ao' or 'wo'.
Raises
------
ValueError
No data is available for any of natural isotopes of the element
ValueError
If only some natural isotopes are avaiable in cross-sections data
library and element is not O, W or Ta
ValueError
Enrichment of isotope which is not present in natural composition
of the element is requested
ValueError
Enrichment is requested of the element that is not composed of
two isotopes.
ValueError
If `enrichment_type` is not 'ao' or 'wo'
ValueError
If `enrichment` is outside <0;100> range
ValueError
If enrichment procedure for Uranium is used when element is not
Uranium.
Notes
-----
When the `enrichment` argument is specified, a correlation from
`ORNL/CSD/TM-244 <https://doi.org/10.2172/5561567>`_ is used to
calculate the weight fractions of U234, U235, U236, and U238. Namely,
the weight fraction of U234 and U236 are taken to be 0.89% and 0.46%,
respectively, of the U235 weight fraction. The remainder of the isotopic
weight is assigned to U238.
respectively, of the U235 weight fraction. The remainder of the
isotopic weight is assigned to U238.
When the `enrichment` argument is specified with `enrichment_target` a
general enrichment procedure is used. It will raise exception unless
element is composed of exactly 2 isotopes. `enrichment` is interpreted
as percent. By default it is atomic. Can be controlled by variable
`enrichment_type`.
"""
# Check input
cv.check_value('enrichment_type', enrichment_type, {'ao', 'wo'})
if enrichment is not None:
cv.check_less_than('enrichment', enrichment, 100.0, equality=True)
cv.check_greater_than('enrichment', enrichment, 0., equality=True)
# Get the nuclides present in nature
natural_nuclides = set()
@ -110,8 +152,8 @@ class Element(str):
mutual_nuclides = sorted(list(mutual_nuclides))
absent_nuclides = sorted(list(absent_nuclides))
# If all natural nuclides are present in the library, expand element
# using all natural nuclides
# If all natural nuclides are present in the library,
# expand element using all natural nuclides
if len(absent_nuclides) == 0:
for nuclide in mutual_nuclides:
abundances[nuclide] = NATURAL_ABUNDANCE[nuclide]
@ -164,7 +206,14 @@ class Element(str):
abundances[nuclide] = NATURAL_ABUNDANCE[nuclide]
# Modify mole fractions if enrichment provided
if enrichment is not None:
# Old treatment for Uranium
if enrichment is not None and enrichment_target is None:
# Check that the element is Uranium
if self.name != 'U':
msg = 'Enrichment procedure for Uranium was requested, '\
'but the isotope is {0} not U'.format(self)
raise ValueError(msg)
# Calculate the mass fractions of isotopes
abundances['U234'] = 0.0089 * enrichment
@ -181,6 +230,80 @@ class Element(str):
for nuclide in abundances.keys():
abundances[nuclide] /= sum_abundances
# Modify mole fractions if enrichment provided
# New treatment for arbitrary element
# Interpret required enrichment as weight
elif enrichment is not None and enrichment_target is not None:
# Provide more informative error message for U235
if enrichment_target == 'U235':
msg = "There is a special procedure for enrichment of U235 "\
"in U. To invoke it, the arguments 'enrichment_taget'"\
"and 'enrichment_type' should be omitted. Provide "\
"only 'enrichment' as 'wo%'. See User Guide for more "\
"details"
raise ValueError(msg)
# Check if it is two-isotope mixture
if len(abundances) != 2:
msg = 'Element {0} does not consist of 2 naturally-occurring '\
'isotopes. Therefore it cannot be enriched with the '\
'in-build procedure. Please enter isotopic abundances '\
'manually.'.format(self)
raise ValueError(msg)
# Check if the target nuclide is present in the mixture
if enrichment_target not in abundances.keys():
msg = 'Could not find the the target nuclide {0} in natural '\
'isotopic composition of element {1}. Following ' \
'isotopes are available: {2} '\
.format(enrichment_target, self, list(abundances.keys()))
raise ValueError(msg)
# If weight percent enrichment is requested convert to mass fractions
if enrichment_type == 'wo':
# Convert the atomic abundances to weight fractions
# Compute the element atomic mass
element_am = 0.0
for nuclide in abundances.keys():
element_am += atomic_mass(nuclide) * abundances[nuclide]
# Convert Molar Fractions to mass fractions
for nuclide in abundances.keys():
abundances[nuclide] *= atomic_mass(nuclide) / element_am
# Normalise to one
sum_abundances = sum(abundances.values())
for nuclide in abundances.keys():
abundances[nuclide] /= sum_abundances
# Enrich the mixture
# The procedure is more generic that it needs to be. It allows
# to enrich mixtures of more then 2 isotopes, keeping the rations
# of non-enriched nuclides the same as in natural composition
# Get fraction of non-enriched isotopes in nat. composition
non_enriched = 1.0 - abundances[enrichment_target]
tail_fraction = 1.0 - enrichment / 100.0
# Enrich all nuclides
# Do bogus operation for enrichment target but overwrite immediatly
# to avoid if statement in the loop
for nuclide, fraction in abundances.items():
abundances[nuclide] = tail_fraction * fraction / non_enriched
abundances[enrichment_target] = enrichment / 100.0
# Convert back to atomic fractions if requested
if enrichment_type == 'wo':
# Convert the mass fractions to mole fractions
for nuclide in abundances.keys():
abundances[nuclide] /= atomic_mass(nuclide)
# Normalize the mole fractions to one
sum_abundances = sum(abundances.values())
for nuclide in abundances.keys():
abundances[nuclide] /= sum_abundances
# Compute the ratio of the nuclide atomic masses to the element
# atomic mass
if percent_type == 'wo':

36
openmc/material.py
View file

@ -499,7 +499,8 @@ class Material(IDManagerMixin):
if macroscopic == self._macroscopic:
self._macroscopic = None
def add_element(self, element, percent, percent_type='ao', enrichment=None):
def add_element(self, element, percent, percent_type='ao', enrichment=None,
enrichment_target=None, enrichment_type='ao'):
"""Add a natural element to the material
Parameters
@ -512,9 +513,24 @@ class Material(IDManagerMixin):
'ao' for atom percent and 'wo' for weight percent. Defaults to atom
percent.
enrichment : float, optional
Enrichment for U235 in weight percent. For example, input 4.95 for
4.95 weight percent enriched U. Default is None
(natural composition).
Enrichment of an enrichment_taget nuclide in percent (ao or wo).
If enrichment_taget is not supplied then it is enrichment for U235
in weight percent. For example, input 4.95 for 4.95 weight percent
enriched U.
Default is None (natural composition).
Value must be in <0;100> for general nuclide
Value must be in <0;100./1.008) for Uranium
enrichment_target: str, optional
Single nuclide name to enrich from a natural composition e.g. O16
enrichment_type: {'ao', 'wo'}, optional
'ao' for enrichment as atom percent and 'wo' for weight percent.
Default is 'ao'
Notes
-----
General enrichment procedure is allowed only for elements composed of
two isotopes. If `enrichment_target` is given without `enrichment`
natural composition is added to the material.
"""
cv.check_type('nuclide', element, str)
@ -526,7 +542,7 @@ class Material(IDManagerMixin):
'macroscopic data-set has already been added'.format(self._id)
raise ValueError(msg)
if enrichment is not None:
if enrichment is not None and enrichment_target is None:
if not isinstance(enrichment, Real):
msg = 'Unable to add an Element to Material ID="{}" with a ' \
'non-floating point enrichment value "{}"'\
@ -558,7 +574,11 @@ class Material(IDManagerMixin):
# Add naturally-occuring isotopes
element = openmc.Element(element)
for nuclide in element.expand(percent, percent_type, enrichment):
for nuclide in element.expand(percent,
percent_type,
enrichment,
enrichment_target,
enrichment_type):
self.add_nuclide(*nuclide)
def add_s_alpha_beta(self, name, fraction=1.0):
@ -927,7 +947,7 @@ class Material(IDManagerMixin):
Fractions of each material to be combined
percent_type : {'ao', 'wo', 'vo'}
Type of percentage, must be one of 'ao', 'wo', or 'vo', to signify atom
percent (molar percent), weight percent, or volume percent,
percent (molar percent), weight percent, or volume percent,
optional. Defaults to 'ao'
name : str
The name for the new material, optional. Defaults to concatenated
@ -996,7 +1016,7 @@ class Material(IDManagerMixin):
zip(materials, fracs)])
new_mat = openmc.Material(name=name)
# Compute atom fractions of nuclides and add them to the new material
# Compute atom fractions of nuclides and add them to the new material
tot_nuclides_per_cc = np.sum([dens for dens in nuclides_per_cc.values()])
for nuc, atom_dens in nuclides_per_cc.items():
new_mat.add_nuclide(nuc, atom_dens/tot_nuclides_per_cc, 'ao')

47
openmc/model/funcs.py
View file

@ -373,52 +373,7 @@ def get_hexagonal_prism(*args, **kwargs):
return hexagonal_prism(*args, **kwargs)
def cylinder_from_points(p1, p2, r, **kwargs):
"""Return cylinder defined by two points passing through its center.
Parameters
----------
p1, p2 : 3-tuples
Coordinates of two points that pass through the center of the cylinder
r : float
Radius of the cylinder
kwargs : dict
Keyword arguments passed to the :class:`openmc.Quadric` constructor
Returns
-------
openmc.Quadric
Quadric surface representing the cylinder.
"""
# Get x, y, z coordinates of two points
x1, y1, z1 = p1
x2, y2, z2 = p2
# Define intermediate terms
dx = x2 - x1
dy = y2 - y1
dz = z2 - z1
cx = y1*z2 - y2*z1
cy = x2*z1 - x1*z2
cz = x1*y2 - x2*y1
# Given p=(x,y,z), p1=(x1, y1, z1), p2=(x2, y2, z2), the equation for the
# cylinder can be derived as r = |(p - p1) ⨯ (p - p2)| / |p2 - p1|.
# Expanding out all terms and grouping according to what Quadric expects
# gives the following coefficients.
kwargs['a'] = dy*dy + dz*dz
kwargs['b'] = dx*dx + dz*dz
kwargs['c'] = dx*dx + dy*dy
kwargs['d'] = -2*dx*dy
kwargs['e'] = -2*dy*dz
kwargs['f'] = -2*dx*dz
kwargs['g'] = 2*(cy*dz - cz*dy)
kwargs['h'] = 2*(cz*dx - cx*dz)
kwargs['j'] = 2*(cx*dy - cy*dx)
kwargs['k'] = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r
return Quadric(**kwargs)
cylinder_from_points = Cylinder.from_points
def subdivide(surfaces):

38
openmc/source.py
View file

@ -8,20 +8,22 @@ from openmc.stats.multivariate import UnitSphere, Spatial
import openmc.checkvalue as cv
class Source(object):
class Source:
"""Distribution of phase space coordinates for source sites.
Parameters
----------
space : openmc.stats.Spatial, optional
space : openmc.stats.Spatial
Spatial distribution of source sites
angle : openmc.stats.UnitSphere, optional
angle : openmc.stats.UnitSphere
Angular distribution of source sites
energy : openmc.stats.Univariate, optional
energy : openmc.stats.Univariate
Energy distribution of source sites
filename : str, optional
filename : str
Source file from which sites should be sampled
strength : Real
library : str
Path to a custom source library
strength : float
Strength of the source
particle : {'neutron', 'photon'}
Source particle type
@ -36,7 +38,9 @@ class Source(object):
Energy distribution of source sites
file : str or None
Source file from which sites should be sampled
strength : Real
library : str or None
Path to a custom source library
strength : float
Strength of the source
particle : {'neutron', 'photon'}
Source particle type
@ -44,11 +48,12 @@ class Source(object):
"""
def __init__(self, space=None, angle=None, energy=None, filename=None,
strength=1.0, particle='neutron'):
library=None, strength=1.0, particle='neutron'):
self._space = None
self._angle = None
self._energy = None
self._file = None
self._library = None
if space is not None:
self.space = space
@ -58,6 +63,8 @@ class Source(object):
self.energy = energy
if filename is not None:
self.file = filename
if library is not None:
self.library = library
self.strength = strength
self.particle = particle
@ -65,6 +72,10 @@ class Source(object):
def file(self):
return self._file
@property
def library(self):
return self._library
@property
def space(self):
return self._space
@ -90,6 +101,11 @@ class Source(object):
cv.check_type('source file', filename, str)
self._file = filename
@library.setter
def library(self, library_name):
cv.check_type('library', library_name, str)
self._library = library_name
@space.setter
def space(self, space):
cv.check_type('spatial distribution', space, Spatial)
@ -131,6 +147,8 @@ class Source(object):
element.set("particle", self.particle)
if self.file is not None:
element.set("file", self.file)
if self.library is not None:
element.set("library", self.library)
if self.space is not None:
element.append(self.space.to_xml_element())
if self.angle is not None:
@ -168,6 +186,10 @@ class Source(object):
if filename is not None:
source.file = filename
library = get_text(elem, 'library')
if library is not None:
source.library = library
space = elem.find('space')
if space is not None:
source.space = Spatial.from_xml_element(space)

1961
openmc/surface.py
View file

File diff suppressed because it is too large Load diff

3
setup.py
View file

@ -56,14 +56,13 @@ kwargs = {
'Topic :: Scientific/Engineering'
'Programming Language :: C++',
'Programming Language :: Python :: 3',
'Programming Language :: Python :: 3.4',
'Programming Language :: Python :: 3.5',
'Programming Language :: Python :: 3.6',
'Programming Language :: Python :: 3.7',
],
# Dependencies
'python_requires': '>=3.4',
'python_requires': '>=3.5',
'install_requires': [
'numpy>=1.9', 'h5py', 'scipy', 'ipython', 'matplotlib',
'pandas', 'lxml', 'uncertainties'

6
src/cell.cpp
View file

@ -377,7 +377,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node)
throw std::runtime_error{"Invalid surface ID " + std::to_string(abs(r))
+ " specified in region for cell " + std::to_string(id_) + "."};
}
r = copysign(it->second + 1, r);
r = (r > 0) ? it->second + 1 : -(it->second + 1);
}
}
@ -536,8 +536,8 @@ CSGCell::to_hdf5(hid_t cell_group) const
region_spec << " |";
} else {
// Note the off-by-one indexing
region_spec << " "
<< copysign(model::surfaces[abs(token)-1]->id_, token);
auto surf_id = model::surfaces[abs(token)-1]->id_;
region_spec << " " << ((token > 0) ? surf_id : -surf_id);
}
}
write_string(group, "region", region_spec.str(), false);

2
src/output.cpp
View file

@ -192,7 +192,7 @@ extern "C" void print_particle(Particle* p)
// Display miscellaneous info.
if (p->surface_ != 0) {
const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]};
fmt::print(" Surface = {}\n", std::copysign(surf.id_, p->surface_));
fmt::print(" Surface = {}\n", (p->surface_ > 0) ? surf.id_ : -surf.id_);
}
fmt::print(" Weight = {}\n", p->wgt_);
if (settings::run_CE) {

4
src/particle.cpp
View file

@ -1,7 +1,7 @@
#include "openmc/particle.h"
#include <algorithm> // copy, min
#include <cmath> // log, abs, copysign
#include <cmath> // log, abs
#include <fmt/core.h>
@ -539,7 +539,7 @@ Particle::cross_surface()
// TODO: off-by-one
surface_ = rotational ?
surf_p->i_periodic_ + 1 :
std::copysign(surf_p->i_periodic_ + 1, surface_);
((surface_ > 0) ? surf_p->i_periodic_ + 1 : -(surf_p->i_periodic_ + 1));
// Figure out what cell particle is in now
n_coord_ = 1;

1
src/settings.cpp
View file

@ -74,6 +74,7 @@ std::string path_input;
std::string path_output;
std::string path_particle_restart;
std::string path_source;
std::string path_source_library;
std::string path_sourcepoint;
std::string path_statepoint;

84
src/source.cpp
View file

@ -1,7 +1,15 @@
#include "openmc/source.h"
#if defined (__unix__) || (defined (__APPLE__) && defined (__MACH__))
#define HAS_DYNAMIC_LINKING
#endif
#include <algorithm> // for move
#ifdef HAS_DYNAMIC_LINKING
#include <dlfcn.h> // for dlopen, dlsym, dlclose, dlerror
#endif
#include <fmt/core.h>
#include "xtensor/xadapt.hpp"
@ -71,7 +79,12 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
fatal_error(fmt::format("Source file '{}' does not exist.",
settings::path_source));
}
} else if (check_for_node(node, "library")) {
settings::path_source_library = get_node_value(node, "library", false, true);
if (!file_exists(settings::path_source_library)) {
fatal_error(fmt::format("Source library '{}' does not exist.",
settings::path_source_library));
}
} else {
// Spatial distribution for external source
@ -237,7 +250,7 @@ void initialize_source()
{
write_message("Initializing source particles...", 5);
if (settings::path_source != "") {
if (!settings::path_source.empty()) {
// Read the source from a binary file instead of sampling from some
// assumed source distribution
@ -261,6 +274,27 @@ void initialize_source()
// Close file
file_close(file_id);
} else if (!settings::path_source_library.empty()) {
#ifdef HAS_DYNAMIC_LINKING
// Get the source from a library object
write_message(fmt::format("Sampling from library source {}...",
settings::path_source), 6);
// Open the library
auto source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY);
if (!source_library) {
fatal_error("Couldn't open source library " + settings::path_source_library);
}
// reset errors
dlerror();
fill_source_bank_custom_source();
#else
fatal_error("Custom source libraries have not yet been implemented for "
"non-POSIX systems");
#endif
} else {
// Generation source sites from specified distribution in user input
@ -321,10 +355,50 @@ void free_memory_source()
model::external_sources.clear();
}
// fill the source bank from the external source
void fill_source_bank_custom_source()
{
#ifdef HAS_DYNAMIC_LINKING
// Open the library
auto source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY);
if (!source_library) {
fatal_error("Couldn't open source library " + settings::path_source_library);
}
// reset errors
dlerror();
// get the function from the library
using sample_t = Particle::Bank (*)(uint64_t* seed);
auto sample_source = reinterpret_cast<sample_t>(dlsym(source_library, "sample_source"));
// check for any dlsym errors
auto dlsym_error = dlerror();
if (dlsym_error) {
dlclose(source_library);
fatal_error(fmt::format("Couldn't open the sample_source symbol: {}", dlsym_error));
}
// Generation source sites from specified distribution in the
// library source
for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
// initialize random number seed
int64_t id = (simulation::total_gen + overall_generation()) *
settings::n_particles + simulation::work_index[mpi::rank] + i + 1;
uint64_t seed = init_seed(id, STREAM_SOURCE);
// sample external source distribution
simulation::source_bank[i] = sample_source(&seed);
}
// release the library
dlclose(source_library);
#endif
}
void fill_source_bank_fixedsource()
{
if (settings::path_source.empty()) {
#pragma omp parallel for
if (settings::path_source.empty() && settings::path_source_library.empty()) {
for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
// initialize random number seed
int64_t id = (simulation::total_gen + overall_generation()) *
@ -334,6 +408,8 @@ void fill_source_bank_fixedsource()
// sample external source distribution
simulation::source_bank[i] = sample_external_source(&seed);
}
} else if (settings::path_source.empty() && !settings::path_source_library.empty()) {
fill_source_bank_custom_source();
}
}

0
tests/regression_tests/source_dlopen/__init__.py Normal file
View file

23
tests/regression_tests/source_dlopen/inputs_true.dat Normal file
View file

@ -0,0 +1,23 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" region="-1" universe="1" />
<surface boundary="vacuum" coeffs="0.0 0.0 0.0 100" id="1" type="sphere" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<material id="1" name="natural_lead">
<density units="g/cm3" value="11.34" />
<nuclide ao="0.014" name="Pb204" />
<nuclide ao="0.241" name="Pb206" />
<nuclide ao="0.221" name="Pb207" />
<nuclide ao="0.524" name="Pb208" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>fixed source</run_mode>
<particles>1000</particles>
<batches>10</batches>
<inactive>0</inactive>
<source library="build/libsource.so" strength="1.0" />
</settings>

0
tests/regression_tests/source_dlopen/results_true.dat Normal file
View file

23
tests/regression_tests/source_dlopen/source_sampling.cpp Normal file
View file

@ -0,0 +1,23 @@
#include <iostream>
#include "openmc/random_lcg.h"
#include "openmc/source.h"
#include "openmc/particle.h"
// you must have external C linkage here otherwise
// dlopen will not find the file
extern "C" openmc::Particle::Bank sample_source(uint64_t *seed) {
openmc::Particle::Bank particle;
// wgt
particle.particle = openmc::Particle::Type::neutron;
particle.wgt = 1.0;
// position
particle.r.x = 0.;
particle.r.y = 0.;
particle.r.z = 0.;
// angle
particle.u = {1.0, 0.0, 0.0};
particle.E = 14.08e6;
particle.delayed_group = 0;
return particle;
}

73
tests/regression_tests/source_dlopen/test.py Normal file
View file

@ -0,0 +1,73 @@
from pathlib import Path
import os
import shutil
import subprocess
import textwrap
import openmc
import pytest
from tests.testing_harness import PyAPITestHarness
@pytest.fixture
def compile_source(request):
"""Compile the external source"""
# Get build directory and write CMakeLists.txt file
openmc_dir = Path(str(request.config.rootdir)) / 'build'
with open('CMakeLists.txt', 'w') as f:
f.write(textwrap.dedent("""
cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
project(openmc_sources CXX)
add_library(source SHARED source_sampling.cpp)
find_package(OpenMC REQUIRED HINTS {})
target_link_libraries(source OpenMC::libopenmc)
""".format(openmc_dir)))
# Create temporary build directory and change to there
local_builddir = Path('build')
local_builddir.mkdir(exist_ok=True)
os.chdir(str(local_builddir))
# Run cmake/make to build the shared libary
subprocess.run(['cmake', os.path.pardir], check=True)
subprocess.run(['make'], check=True)
os.chdir(os.path.pardir)
yield
# Remove local build directory when test is complete
shutil.rmtree('build')
@pytest.fixture
def model():
model = openmc.model.Model()
natural_lead = openmc.Material(name="natural_lead")
natural_lead.add_element('Pb', 1.0)
natural_lead.set_density('g/cm3', 11.34)
model.materials.append(natural_lead)
# geometry
surface_sph1 = openmc.Sphere(r=100, boundary_type='vacuum')
cell_1 = openmc.Cell(fill=natural_lead, region=-surface_sph1)
model.geometry = openmc.Geometry([cell_1])
# settings
model.settings.batches = 10
model.settings.inactive = 0
model.settings.particles = 1000
model.settings.run_mode = 'fixed source'
# custom source from shared library
source = openmc.Source()
source.library = 'build/libsource.so'
model.settings.source = source
return model
def test_dlopen_source(compile_source, model):
harness = PyAPITestHarness('statepoint.10.h5', model)
harness.main()

7
tests/unit_tests/test_data_photon.py
View file

@ -2,6 +2,7 @@
from collections.abc import Mapping, Callable
import os
from pathlib import Path
import numpy as np
import pandas as pd
@ -144,3 +145,9 @@ def test_export_to_hdf5(tmpdir, element):
element2.bremsstrahlung['electron_energy']).all()
# Export to hdf5 again
element2.export_to_hdf5(filename, 'w')
def test_photodat_only(run_in_tmpdir):
endf_dir = Path(os.environ['OPENMC_ENDF_DATA'])
photoatomic_file = endf_dir / 'photoat' / 'photoat-001_H_000.endf'
data = openmc.data.IncidentPhoton.from_endf(photoatomic_file)
data.export_to_hdf5('tmp.h5', 'w')

8
tests/unit_tests/test_deplete_nuclide.py
View file

@ -149,12 +149,18 @@ def test_fission_yield_distribution():
# __getitem__ return yields as a view into yield matrix
assert orig_yields.yields.base is yield_dist.yield_matrix
# Fission yield feature uses scaled and incremented
# Scale and increment fission yields
mod_yields = orig_yields * 2
assert numpy.array_equal(orig_yields.yields * 2, mod_yields.yields)
mod_yields += orig_yields
assert numpy.array_equal(orig_yields.yields * 3, mod_yields.yields)
mod_yields = 2.0 * orig_yields
assert numpy.array_equal(orig_yields.yields * 2, mod_yields.yields)
mod_yields = numpy.float64(2.0) * orig_yields
assert numpy.array_equal(orig_yields.yields * 2, mod_yields.yields)
# Failure modes for adding, multiplying yields
similar = numpy.empty_like(orig_yields.yields)
with pytest.raises(TypeError):

79
tests/unit_tests/test_element.py Normal file
View file

@ -0,0 +1,79 @@
import openmc.element
import pytest as pt
from openmc.data import NATURAL_ABUNDANCE, atomic_mass
# Relative tolerance for float comparison
TOL = 1e-9
def test_expand_no_enrichment():
""" Expand Li in natural compositions"""
lithium = openmc.element.Element('Li')
# Verify the expansion into ATOMIC fraction against natural composition
for isotope in lithium.expand(100.0, 'ao'):
assert isotope[1] == pt.approx(NATURAL_ABUNDANCE[isotope[0]] * 100.0, rel=TOL)
# Verify the expansion into WEIGHT fraction against natural composition
natural = {'Li6': NATURAL_ABUNDANCE['Li6'] * atomic_mass('Li6'),
'Li7': NATURAL_ABUNDANCE['Li7'] * atomic_mass('Li7')}
li_am = sum(natural.values())
for key in natural.keys():
natural[key] /= li_am
for isotope in lithium.expand(100.0, 'wo'):
assert isotope[1] == pt.approx(natural[isotope[0]] * 100.0, rel=TOL)
def test_expand_enrichment():
""" Expand and verify enrichment of Li """
lithium = openmc.element.Element('Li')
# Verify the enrichment by atoms
ref = {'Li6': 75.0, 'Li7': 25.0}
for isotope in lithium.expand(100.0, 'ao', 25.0, 'Li7', 'ao'):
assert isotope[1] == pt.approx(ref[isotope[0]], rel=TOL)
# Verify the enrichment by weight
for isotope in lithium.expand(100.0, 'wo', 25.0, 'Li7', 'wo'):
assert isotope[1] == pt.approx(ref[isotope[0]], rel=TOL)
def test_expand_exceptions():
""" Test that correct exceptions are raised for invalid input """
# 1 Isotope Element
with pt.raises(ValueError):
element = openmc.element.Element('Be')
fail = element.expand(70.0, 'ao', 4.0, 'Be9')
# 3 Isotope Element
with pt.raises(ValueError):
element = openmc.element.Element('Cr')
fail = element.expand(70.0, 'ao', 4.0, 'Cr52')
# Non-present Enrichment Target
with pt.raises(ValueError):
element = openmc.element.Element('H')
fail = element.expand(70.0, 'ao', 4.0, 'H4')
# Enrichment Procedure for Uranium if not Uranium
with pt.raises(ValueError):
element = openmc.element.Element('Li')
fail = element.expand(70.0, 'ao', 4.0)
# Missing Enrichment Target
with pt.raises(ValueError):
element = openmc.element.Element('Li')
fail = element.expand(70.0, 'ao', 4.0, enrichment_type='ao')
# Invalid Enrichment Type Entry
with pt.raises(ValueError):
element = openmc.element.Element('Li')
fail = element.expand(70.0, 'ao', 4.0, 'Li7', 'Grand Moff Tarkin')
# Trying to enrich Uranium
with pt.raises(ValueError):
element = openmc.element.Element('U')
fail = element.expand(70.0, 'ao', 4.0, 'U235', 'wo')

7
tests/unit_tests/test_material.py
View file

@ -29,10 +29,17 @@ def test_elements():
m = openmc.Material()
m.add_element('Zr', 1.0)
m.add_element('U', 1.0, enrichment=4.5)
m.add_element('Li', 1.0, enrichment=60.0, enrichment_target='Li7')
m.add_element('H', 1.0, enrichment=50.0, enrichment_target='H2',
enrichment_type='wo')
with pytest.raises(ValueError):
m.add_element('U', 1.0, enrichment=100.0)
with pytest.raises(ValueError):
m.add_element('Pu', 1.0, enrichment=3.0)
with pytest.raises(ValueError):
m.add_element('U', 1.0, enrichment=70.0, enrichment_target='U235')
with pytest.raises(ValueError):
m.add_element('He', 1.0, enrichment=17.0, enrichment_target='He6')
def test_density():

8
tests/unit_tests/test_source.py
View file

@ -34,3 +34,11 @@ def test_source_file():
elem = src.to_xml_element()
assert 'strength' in elem.attrib
assert 'file' in elem.attrib
def test_source_dlopen():
library = './libsource.so'
src = openmc.Source(library=library)
assert src.library == library
elem = src.to_xml_element()
assert 'library' in elem.attrib

133
tests/unit_tests/test_surface.py
View file

@ -2,6 +2,7 @@ from functools import partial
from random import uniform, seed
import numpy as np
import math
import openmc
import pytest
@ -139,6 +140,51 @@ def test_zplane():
repr(s)
def test_cylinder():
x0, y0, z0, r = 2, 3, 4, 2
dx, dy, dz = 1, -1, 1
s = openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=r)
assert s.x0 == 2
assert s.y0 == 3
assert s.z0 == 4
assert s.dx == 1
assert s.dy == -1
assert s.dz == 1
assert s.r == 2
# Check bounding box
assert_infinite_bb(s)
# evaluate method
# |(p - p1) ⨯ (p - p2)|^2 / |p2 - p1|^2 - r^2
p1 = s._origin
p2 = p1 + s._axis
perp = np.array((1, -2, 1))*(1 / s._axis)
divisor = np.linalg.norm(p2 - p1)
pin = p1 + 5*s._axis # point inside cylinder
pout = np.array((4., 0., 2.5)) # point outside the cylinder
pon = p1 + s.r*perp / np.linalg.norm(perp) # point on cylinder
for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)):
c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / divisor
val = c1*c1 - s.r*s.r
p_eval = s.evaluate(p)
assert fn(p_eval, 0.)
assert p_eval == pytest.approx(val)
# translate method
st = s.translate((1.0, 1.0, 1.0))
assert st.x0 == s.x0 + 1
assert st.y0 == s.y0 + 1
assert st.z0 == s.z0 + 1
assert st.dx == s.dx
assert st.dy == s.dy
assert st.dz == s.dz
assert st.r == s.r
# Make sure repr works
repr(s)
def test_xcylinder():
y, z, r = 3, 5, 2
s = openmc.XCylinder(y0=y, z0=z, r=r)
@ -284,6 +330,58 @@ def cone_common(apex, r2, cls):
repr(s)
def test_cone():
x0, y0, z0, r2 = 2, 3, 4, 4
dx, dy, dz = 1, -1, 1
s = openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=r2)
assert s.x0 == 2
assert s.y0 == 3
assert s.z0 == 4
assert s.dx == 1
assert s.dy == -1
assert s.dz == 1
assert s.r2 == 4
# Check bounding box
assert_infinite_bb(s)
# evaluate method
# cos^2(theta) * ((p - p1))**2 - (d @ (p - p1))^2
# The argument r2 for cones is actually tan^2(theta) so that
# cos^2(theta) = 1 / (1 + r2)
#
# This makes the evaluation equation shown below where p is the evaluation
# point (x, y, z) p1 is the apex (origin) of the cone and r2 is related to
# the aperature of the cone as described above
# (p - p1) @ (p - p1) / (1 + r2) - (d @ (p - p1))^2
# point inside
p1 = s._origin
d = s._axis
perp = np.array((1, -2, 1))*(1 / d)
perp /= np.linalg.norm(perp)
pin = p1 + 5*d # point inside cone
pout = p1 + 3.2*perp # point outside cone
pon = p1 + 3.2*d + 3.2*math.sqrt(s.r2)*perp # point on cone
for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)):
val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2)
p_eval = s.evaluate(p)
assert fn(p_eval, 0.)
assert p_eval == pytest.approx(val)
# translate method
st = s.translate((1.0, 1.0, 1.0))
assert st.x0 == s.x0 + 1
assert st.y0 == s.y0 + 1
assert st.z0 == s.z0 + 1
assert st.dx == s.dx
assert st.dy == s.dy
assert st.dz == s.dz
assert st.r2 == s.r2
# Make sure repr works
repr(s)
def test_xcone():
apex = (10, 0, 0)
r2 = 4
@ -339,7 +437,7 @@ def test_cylinder_from_points():
p1 = np.array([xi(), xi(), xi()])
p2 = np.array([xi(), xi(), xi()])
r = uniform(1.0, 100.0)
s = openmc.model.cylinder_from_points(p1, p2, r)
s = openmc.Cylinder.from_points(p1, p2, r)
# Points p1 and p2 need to be inside cylinder
assert p1 in -s
@ -369,24 +467,27 @@ def test_cylinder_from_points_axis():
# (x - 3)^2 + (y - 4)^2 = 2^2
# x^2 + y^2 - 6x - 8y + 21 = 0
s = openmc.model.cylinder_from_points((3., 4., 0.), (3., 4., 1.), 2.)
assert (s.a, s.b, s.c) == pytest.approx((1., 1., 0.))
assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
assert (s.g, s.h, s.j) == pytest.approx((-6., -8., 0.))
assert s.k == pytest.approx(21.)
s = openmc.Cylinder.from_points((3., 4., 0.), (3., 4., 1.), 2.)
a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs()
assert (a, b, c) == pytest.approx((1., 1., 0.))
assert (d, e, f) == pytest.approx((0., 0., 0.))
assert (g, h, j) == pytest.approx((-6., -8., 0.))
assert k == pytest.approx(21.)
# (y + 7)^2 + (z - 1)^2 = 3^2
# y^2 + z^2 + 14y - 2z + 41 = 0
s = openmc.model.cylinder_from_points((0., -7, 1.), (1., -7., 1.), 3.)
assert (s.a, s.b, s.c) == pytest.approx((0., 1., 1.))
assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
assert (s.g, s.h, s.j) == pytest.approx((0., 14., -2.))
assert s.k == 41.
s = openmc.Cylinder.from_points((0., -7, 1.), (1., -7., 1.), 3.)
a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs()
assert (a, b, c) == pytest.approx((0., 1., 1.))
assert (d, e, f) == pytest.approx((0., 0., 0.))
assert (g, h, j) == pytest.approx((0., 14., -2.))
assert k == 41.
# (x - 2)^2 + (z - 5)^2 = 4^2
# x^2 + z^2 - 4x - 10z + 13 = 0
s = openmc.model.cylinder_from_points((2., 0., 5.), (2., 1., 5.), 4.)
assert (s.a, s.b, s.c) == pytest.approx((1., 0., 1.))
assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
assert (s.g, s.h, s.j) == pytest.approx((-4., 0., -10.))
assert s.k == pytest.approx(13.)
s = openmc.Cylinder.from_points((2., 0., 5.), (2., 1., 5.), 4.)
a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs()
assert (a, b, c) == pytest.approx((1., 0., 1.))
assert (d, e, f) == pytest.approx((0., 0., 0.))
assert (g, h, j) == pytest.approx((-4., 0., -10.))
assert k == pytest.approx(13.)