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

Differentiate materials in DAGMC universes (#3056)

Browse source 
Co-authored-by: Baptiste Mouginot <bam.git.not.reply@gmail.com>
Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>
Co-authored-by: azimG <azimgivron@yahoo.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Baptiste Mouginot 2025-01-07 22:50:02 +01:00 • committed by GitHub
parent 393334829d
commit 5ad1a4aa1e
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
19 changed files with 1424 additions and 652 deletions
Download patch file Download diff file Expand all files Collapse all files

34
docs/source/io_formats/geometry.rst
View file

@ -407,13 +407,33 @@ Each ``<dagmc_universe>`` element can have the following attributes or sub-eleme
*Default*: None
:material_overrides:
This element contains information on material overrides to be applied to the
DAGMC universe. It has the following attributes and sub-elements:
.. note:: A geometry.xml file containing only a DAGMC model for a file named `dagmc.h5m` (no CSG)
looks as follows
:cell:
Material override information for a single cell. It contains the following
attributes and sub-elements:
.. code-block:: xml
:id:
The cell ID in the DAGMC geometry for which the material override will
apply.
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<dagmc_universe filename="dagmc.h5m" id="1" />
</geometry>
:materials:
A list of material IDs that will apply to instances of the cell. If the
list contains only one ID, it will replace the original material
assignment of all instances of the DAGMC cell. If the list contains more
than one material, each material ID of the list will be assigned to the
various instances of the DAGMC cell.
*Default*: None
.. note:: A geometry.xml file containing only a DAGMC model for a file named
`dagmc.h5m` (no CSG) looks as follows:
.. code-block:: xml
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<dagmc_universe filename="dagmc.h5m" id="1" />
</geometry>

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

@ -19,6 +19,7 @@ Functions
finalize
find_cell
find_material
dagmc_universe_cell_ids
global_bounding_box
global_tallies
hard_reset

3
include/openmc/capi.h
View file

@ -29,6 +29,9 @@ int openmc_cell_set_temperature(
int32_t index, double T, const int32_t* instance, bool set_contained = false);
int openmc_cell_set_translation(int32_t index, const double xyz[]);
int openmc_cell_set_rotation(int32_t index, const double rot[], size_t rot_len);
int openmc_dagmc_universe_get_cell_ids(
int32_t univ_id, int32_t* ids, size_t* n);
int openmc_dagmc_universe_get_num_cells(int32_t univ_id, size_t* n);
int openmc_energy_filter_get_bins(
int32_t index, const double** energies, size_t* n);
int openmc_energy_filter_set_bins(

8
include/openmc/cell.h
View file

@ -320,7 +320,6 @@ public:
int32_t universe_; //!< Universe # this cell is in
int32_t fill_; //!< Universe # filling this cell
int32_t n_instances_ {0}; //!< Number of instances of this cell
GeometryType geom_type_; //!< Geometric representation type (CSG, DAGMC)
//! \brief Index corresponding to this cell in distribcell arrays
int distribcell_index_ {C_NONE};
@ -350,6 +349,13 @@ public:
vector<double> rotation_;
vector<int32_t> offset_; //!< Distribcell offset table
// Accessors
const GeometryType& geom_type() const { return geom_type_; }
GeometryType& geom_type() { return geom_type_; }
private:
GeometryType geom_type_; //!< Geometric representation type (CSG, DAGMC)
};
struct CellInstanceItem {

15
include/openmc/dagmc.h
View file

@ -29,6 +29,12 @@ void check_dagmc_root_univ();
#include "openmc/particle.h"
#include "openmc/position.h"
#include "openmc/surface.h"
#include "openmc/vector.h"
#include <memory> // for shared_ptr, unique_ptr
#include <string>
#include <unordered_map>
#include <utility> // for pair
class UWUW;
@ -133,6 +139,10 @@ public:
void legacy_assign_material(
std::string mat_string, std::unique_ptr<DAGCell>& c) const;
//! Assign a material overriding normal assignement to a cell
//! \param[in] c The OpenMC cell to which the material is assigned
void override_assign_material(std::unique_ptr<DAGCell>& c) const;
//! Return the index into the model cells vector for a given DAGMC volume
//! handle in the universe
//! \param[in] vol MOAB handle to the DAGMC volume set
@ -184,6 +194,11 @@ private:
//!< generate new material IDs for the universe
bool has_graveyard_; //!< Indicates if the DAGMC geometry has a "graveyard"
//!< volume
std::unordered_map<int32_t, vector<int32_t>>
material_overrides_; //!< Map of material overrides
//!< keys correspond to the DAGMCCell id
//!< values are a list of material ids used
//!< for the override
};
//==============================================================================

8
include/openmc/surface.h
View file

@ -38,7 +38,6 @@ public:
int id_; //!< Unique ID
std::string name_; //!< User-defined name
unique_ptr<BoundaryCondition> bc_; //!< Boundary condition
GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC)
bool surf_source_ {false}; //!< Activate source banking for the surface?
explicit Surface(pugi::xml_node surf_node);
@ -91,6 +90,13 @@ public:
//! Get the BoundingBox for this surface.
virtual BoundingBox bounding_box(bool /*pos_side*/) const { return {}; }
// Accessors
const GeometryType& geom_type() const { return geom_type_; }
GeometryType& geom_type() { return geom_type_; }
private:
GeometryType geom_type_; //!< Geometry type indicator (CSG or DAGMC)
protected:
virtual void to_hdf5_inner(hid_t group_id) const = 0;
};

1
openmc/__init__.py
View file

@ -15,6 +15,7 @@ from openmc.volume import *
from openmc.weight_windows import *
from openmc.surface import *
from openmc.universe import *
from openmc.dagmc import *
from openmc.source import *
from openmc.settings import *
from openmc.lattice import *

625
openmc/dagmc.py Normal file
View file

@ -0,0 +1,625 @@
from collections.abc import Iterable, Mapping
from numbers import Integral
import h5py
import lxml.etree as ET
import numpy as np
import warnings
import openmc
import openmc.checkvalue as cv
from ._xml import get_text
from .checkvalue import check_type, check_value
from .surface import _BOUNDARY_TYPES
from .bounding_box import BoundingBox
from .utility_funcs import input_path
class DAGMCUniverse(openmc.UniverseBase):
"""A reference to a DAGMC file to be used in the model.
.. versionadded:: 0.13.0
Parameters
----------
filename : str
Path to the DAGMC file used to represent this universe.
universe_id : int, optional
Unique identifier of the universe. If not specified, an identifier will
automatically be assigned.
name : str, optional
Name of the universe. If not specified, the name is the empty string.
auto_geom_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between CSG and DAGMC (False)
auto_mat_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between OpenMC and UWUW materials (False)
material_overrides : dict, optional
A dictionary of material overrides. The keys are material name strings
and the values are Iterables of openmc.Material objects. If a material
name is found in the DAGMC file, the material will be replaced with the
openmc.Material object in the value.
Attributes
----------
id : int
Unique identifier of the universe
name : str
Name of the universe
filename : str
Path to the DAGMC file used to represent this universe.
auto_geom_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between CSG and DAGMC (False)
auto_mat_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between OpenMC and UWUW materials (False)
bounding_box : openmc.BoundingBox
Lower-left and upper-right coordinates of an axis-aligned bounding box
of the universe.
.. versionadded:: 0.13.1
material_names : list of str
Return a sorted list of materials names that are contained within the
DAGMC h5m file. This is useful when naming openmc.Material() objects as
each material name present in the DAGMC h5m file must have a matching
openmc.Material() with the same name.
.. versionadded:: 0.13.2
n_cells : int
The number of cells in the DAGMC model. This is the number of cells at
runtime and accounts for the implicit complement whether or not is it
present in the DAGMC file.
.. versionadded:: 0.13.2
n_surfaces : int
The number of surfaces in the model.
.. versionadded:: 0.13.2
material_overrides : dict
A dictionary of material overrides. Keys are cell IDs; values are
iterables of :class:`openmc.Material` objects. The material assignment
of each DAGMC cell ID key will be replaced with the
:class:`~openmc.Material` object in the value. If the value contains
multiple :class:`~openmc.Material` objects, each Material in the list
will be assigned to the corresponding instance of the cell.
.. versionadded:: 0.15.1
"""
def __init__(self,
filename: cv.PathLike,
universe_id=None,
name='',
auto_geom_ids=False,
auto_mat_ids=False,
material_overrides=None):
super().__init__(universe_id, name)
# Initialize class attributes
self.filename = filename
self.auto_geom_ids = auto_geom_ids
self.auto_mat_ids = auto_mat_ids
self._material_overrides = {}
if material_overrides is not None:
self.material_overrides = material_overrides
def __repr__(self):
string = super().__repr__()
string += '{: <16}=\t{}\n'.format('\tGeom', 'DAGMC')
string += '{: <16}=\t{}\n'.format('\tFile', self.filename)
return string
@property
def bounding_box(self):
with h5py.File(self.filename) as dagmc_file:
coords = dagmc_file['tstt']['nodes']['coordinates'][()]
lower_left_corner = coords.min(axis=0)
upper_right_corner = coords.max(axis=0)
return openmc.BoundingBox(lower_left_corner, upper_right_corner)
@property
def filename(self):
return self._filename
@filename.setter
def filename(self, val: cv.PathLike):
cv.check_type('DAGMC filename', val, cv.PathLike)
self._filename = input_path(val)
@property
def material_overrides(self):
return self._material_overrides
@material_overrides.setter
def material_overrides(self, val):
cv.check_type('material overrides', val, Mapping)
for key, value in val.items():
self.add_material_override(key, value)
def replace_material_assignment(self, material_name: str, material: openmc.Material):
"""Replace a material assignment within the DAGMC universe.
Replace the material assignment of all cells filled with a material in
the DAGMC universe. The universe must be synchronized in an initialized
Model (see :meth:`~openmc.DAGMCUniverse.sync_dagmc_cells`) before
calling this method.
.. versionadded:: 0.15.1
Parameters
----------
material_name : str
Material name to replace
material : openmc.Material
Material to replace the material_name with
"""
if material_name not in self.material_names:
raise ValueError(
f"No material with name '{material_name}' found in the DAGMC universe")
if not self.cells:
raise RuntimeError("This DAGMC universe has not been synchronized "
"on an initialized Model.")
for cell in self.cells.values():
if cell.fill is None:
continue
if isinstance(cell.fill, openmc.Iterable):
cell.fill = list(map(lambda x: material if x.name == material_name else x, cell.fill))
else:
cell.fill = material if cell.fill.name == material_name else cell.fill
def add_material_override(self, key, overrides=None):
"""Add a material override to the universe.
.. versionadded:: 0.15
Parameters
----------
key : openmc.DAGMCCell or int
Cell object or ID of the Cell to override
value : openmc.Material or Iterable of openmc.Material
Material(s) to be applied to the Cell passed as the key
"""
# Ensure that they key is a valid type
if not isinstance(key, (int, openmc.DAGMCCell)):
raise ValueError("Unrecognized key type. "
"Must be an integer or openmc.DAGMCCell object")
# Ensure that overrides is an iterable of openmc.Material
overrides = overrides if isinstance(overrides, openmc.Iterable) else [overrides]
cv.check_iterable_type('material objects', overrides, (openmc.Material, type(None)))
# if a DAGMCCell is passed, redcue the key to the ID of the cell
if isinstance(key, openmc.DAGMCCell):
key = key.id
if key not in self.cells:
raise ValueError(f"Cell ID '{key}' not found in DAGMC universe")
self._material_overrides[key] = overrides
@property
def auto_geom_ids(self):
return self._auto_geom_ids
@auto_geom_ids.setter
def auto_geom_ids(self, val):
cv.check_type('DAGMC automatic geometry ids', val, bool)
self._auto_geom_ids = val
@property
def auto_mat_ids(self):
return self._auto_mat_ids
@auto_mat_ids.setter
def auto_mat_ids(self, val):
cv.check_type('DAGMC automatic material ids', val, bool)
self._auto_mat_ids = val
@property
def material_names(self):
dagmc_file_contents = h5py.File(self.filename)
material_tags_hex = dagmc_file_contents['/tstt/tags/NAME'].get(
'values')
material_tags_ascii = []
for tag in material_tags_hex:
candidate_tag = tag.tobytes().decode().replace('\x00', '')
# tags might be for temperature or reflective surfaces
if candidate_tag.startswith('mat:'):
# if name ends with _comp remove it, it is not parsed
if candidate_tag.endswith('_comp'):
candidate_tag = candidate_tag[:-5]
# removes first 4 characters as openmc.Material name should be
# set without the 'mat:' part of the tag
material_tags_ascii.append(candidate_tag[4:])
return sorted(set(material_tags_ascii))
def _n_geom_elements(self, geom_type):
"""
Helper function for retrieving the number geometric entities in a DAGMC
file
Parameters
----------
geom_type : str
The type of geometric entity to count. One of {'Volume', 'Surface'}. Returns
the runtime number of voumes in the DAGMC model (includes implicit complement).
Returns
-------
int
Number of geometry elements of the specified type
"""
cv.check_value('geometry type', geom_type, ('volume', 'surface'))
def decode_str_tag(tag_val):
return tag_val.tobytes().decode().replace('\x00', '')
with h5py.File(self.filename) as dagmc_file:
category_data = dagmc_file['tstt/tags/CATEGORY/values']
category_strs = map(decode_str_tag, category_data)
n = sum([v == geom_type.capitalize() for v in category_strs])
# check for presence of an implicit complement in the file and
# increment the number of cells if it doesn't exist
if geom_type == 'volume':
name_data = dagmc_file['tstt/tags/NAME/values']
name_strs = map(decode_str_tag, name_data)
if not sum(['impl_complement' in n for n in name_strs]):
n += 1
return n
@property
def n_cells(self):
return self._n_geom_elements('volume')
@property
def n_surfaces(self):
return self._n_geom_elements('surface')
def create_xml_subelement(self, xml_element, memo=None):
if memo is None:
memo = set()
if self in memo:
return
memo.add(self)
# Ensure that the material overrides are up-to-date
for cell in self.cells.values():
if cell.fill is None:
continue
self.add_material_override(cell, cell.fill)
# Set xml element values
dagmc_element = ET.Element('dagmc_universe')
dagmc_element.set('id', str(self.id))
if self.auto_geom_ids:
dagmc_element.set('auto_geom_ids', 'true')
if self.auto_mat_ids:
dagmc_element.set('auto_mat_ids', 'true')
dagmc_element.set('filename', str(self.filename))
if self._material_overrides:
mat_element = ET.Element('material_overrides')
for key in self._material_overrides:
cell_overrides = ET.Element('cell_override')
cell_overrides.set("id", str(key))
material_element = ET.Element('material_ids')
material_element.text = ' '.join(
str(t.id) for t in self._material_overrides[key])
cell_overrides.append(material_element)
mat_element.append(cell_overrides)
dagmc_element.append(mat_element)
xml_element.append(dagmc_element)
def bounding_region(
self,
bounded_type: str = 'box',
boundary_type: str = 'vacuum',
starting_id: int = 10000,
padding_distance: float = 0.
):
"""Creates a either a spherical or box shaped bounding region around
the DAGMC geometry.
.. versionadded:: 0.13.1
Parameters
----------
bounded_type : str
The type of bounding surface(s) to use when constructing the region.
Options include a single spherical surface (sphere) or a rectangle
made from six planes (box).
boundary_type : str
Boundary condition that defines the behavior for particles hitting
the surface. Defaults to vacuum boundary condition. Passed into the
surface construction.
starting_id : int
Starting ID of the surface(s) used in the region. For bounded_type
'box', the next 5 IDs will also be used. Defaults to 10000 to reduce
the chance of an overlap of surface IDs with the DAGMC geometry.
padding_distance : float
Distance between the bounding region surfaces and the minimal
bounding box. Allows for the region to be larger than the DAGMC
geometry.
Returns
-------
openmc.Region
Region instance
"""
check_type('boundary type', boundary_type, str)
check_value('boundary type', boundary_type, _BOUNDARY_TYPES)
check_type('starting surface id', starting_id, Integral)
check_type('bounded type', bounded_type, str)
check_value('bounded type', bounded_type, ('box', 'sphere'))
bbox = self.bounding_box.expand(padding_distance, True)
if bounded_type == 'sphere':
radius = np.linalg.norm(bbox.upper_right - bbox.center)
bounding_surface = openmc.Sphere(
surface_id=starting_id,
x0=bbox.center[0],
y0=bbox.center[1],
z0=bbox.center[2],
boundary_type=boundary_type,
r=radius,
)
return -bounding_surface
if bounded_type == 'box':
# defines plane surfaces for all six faces of the bounding box
lower_x = openmc.XPlane(bbox[0][0], surface_id=starting_id)
upper_x = openmc.XPlane(bbox[1][0], surface_id=starting_id+1)
lower_y = openmc.YPlane(bbox[0][1], surface_id=starting_id+2)
upper_y = openmc.YPlane(bbox[1][1], surface_id=starting_id+3)
lower_z = openmc.ZPlane(bbox[0][2], surface_id=starting_id+4)
upper_z = openmc.ZPlane(bbox[1][2], surface_id=starting_id+5)
region = +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z
for surface in region.get_surfaces().values():
surface.boundary_type = boundary_type
return region
def bounded_universe(self, bounding_cell_id=10000, **kwargs):
"""Returns an openmc.Universe filled with this DAGMCUniverse and bounded
with a cell. Defaults to a box cell with a vacuum surface however this
can be changed using the kwargs which are passed directly to
DAGMCUniverse.bounding_region().
Parameters
----------
bounding_cell_id : int
The cell ID number to use for the bounding cell, defaults to 10000 to reduce
the chance of overlapping ID numbers with the DAGMC geometry.
Returns
-------
openmc.Universe
Universe instance
"""
bounding_cell = openmc.Cell(
fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs))
return openmc.Universe(cells=[bounding_cell])
@classmethod
def from_hdf5(cls, group):
"""Create DAGMC universe from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
Returns
-------
openmc.DAGMCUniverse
DAGMCUniverse instance
"""
id = int(group.name.split('/')[-1].lstrip('universe '))
fname = group['filename'][()].decode()
name = group['name'][()].decode() if 'name' in group else None
out = cls(fname, universe_id=id, name=name)
out.auto_geom_ids = bool(group.attrs['auto_geom_ids'])
out.auto_mat_ids = bool(group.attrs['auto_mat_ids'])
return out
@classmethod
def from_xml_element(cls, elem, mats = None):
"""Generate DAGMC universe from XML element
Parameters
----------
elem : lxml.etree._Element
`<dagmc_universe>` element
mats : dict
Dictionary mapping material ID strings to :class:`openmc.Material`
instances (defined in :meth:`openmc.Geometry.from_xml`)
Returns
-------
openmc.DAGMCUniverse
DAGMCUniverse instance
"""
id = int(get_text(elem, 'id'))
fname = get_text(elem, 'filename')
out = cls(fname, universe_id=id)
name = get_text(elem, 'name')
if name is not None:
out.name = name
out.auto_geom_ids = bool(elem.get('auto_geom_ids'))
out.auto_mat_ids = bool(elem.get('auto_mat_ids'))
el_mat_override = elem.find('material_overrides')
if el_mat_override is not None:
if mats is None:
raise ValueError("Material overrides found in DAGMC universe "
"but no materials were provided to populate "
"the mapping.")
out._material_overrides = {}
for elem in el_mat_override.findall('cell_override'):
cell_id = int(get_text(elem, 'id'))
mat_ids = get_text(elem, 'material_ids').split(' ')
mat_objs = [mats[mat_id] for mat_id in mat_ids]
out._material_overrides[cell_id] = mat_objs
return out
def _partial_deepcopy(self):
"""Clone all of the openmc.DAGMCUniverse object's attributes except for
its cells, as they are copied within the clone function. This should
only to be used within the openmc.UniverseBase.clone() context.
"""
clone = openmc.DAGMCUniverse(name=self.name, filename=self.filename)
clone.volume = self.volume
clone.auto_geom_ids = self.auto_geom_ids
clone.auto_mat_ids = self.auto_mat_ids
return clone
def add_cell(self, cell):
"""Add a cell to the universe.
Parameters
----------
cell : openmc.DAGMCCell
Cell to add
"""
if not isinstance(cell, openmc.DAGMCCell):
msg = f'Unable to add a DAGMCCell to DAGMCUniverse ' \
f'ID="{self._id}" since "{cell}" is not a DAGMCCell'
raise TypeError(msg)
cell_id = cell.id
if cell_id not in self._cells:
self._cells[cell_id] = cell
def remove_cell(self, cell):
"""Remove a cell from the universe.
Parameters
----------
cell : openmc.Cell
Cell to remove
"""
if not isinstance(cell, openmc.DAGMCCell):
msg = f'Unable to remove a Cell from Universe ID="{self._id}" ' \
f'since "{cell}" is not a Cell'
raise TypeError(msg)
# If the Cell is in the Universe's list of Cells, delete it
self._cells.pop(cell.id, None)
def sync_dagmc_cells(self, mats: Iterable[openmc.Material]):
"""Synchronize DAGMC cell information between Python and C API
.. versionadded:: 0.15.1
Parameters
----------
mats : iterable of openmc.Material
Iterable of materials to assign to the DAGMC cells
"""
import openmc.lib
if not openmc.lib.is_initialized:
raise RuntimeError("This universe must be part of an openmc.Model "
"initialized via Model.init_lib before calling "
"this method.")
dagmc_cell_ids = openmc.lib.dagmc.dagmc_universe_cell_ids(self.id)
if len(dagmc_cell_ids) != self.n_cells:
raise ValueError(
f"Number of cells in DAGMC universe {self.id} does not match "
f"the number of cells in the Python universe."
)
mats_per_id = {mat.id: mat for mat in mats}
for dag_cell_id in dagmc_cell_ids:
dag_cell = openmc.lib.cells[dag_cell_id]
if isinstance(dag_cell.fill, Iterable):
fill = [mats_per_id[mat.id] for mat in dag_cell.fill if mat]
else:
fill = mats_per_id[dag_cell.fill.id] if dag_cell.fill else None
self.add_cell(openmc.DAGMCCell(cell_id=dag_cell_id, fill=fill))
class DAGMCCell(openmc.Cell):
"""A cell class for DAGMC-based geometries.
.. versionadded:: 0.15.1
Parameters
----------
cell_id : int or None, optional
Unique identifier for the cell. If None, an identifier will be
automatically assigned.
name : str, optional
Name of the cell.
fill : openmc.Material or None, optional
Material filling the cell. If None, the cell is filled with vacuum.
Attributes
----------
DAG_parent_universe : int
The parent universe of the cell.
"""
def __init__(self, cell_id=None, name='', fill=None):
super().__init__(cell_id, name, fill, None)
@property
def DAG_parent_universe(self):
"""Get the parent universe of the cell."""
return self._parent_universe
@DAG_parent_universe.setter
def DAG_parent_universe(self, universe):
"""Set the parent universe of the cell."""
self._parent_universe = universe.id
def bounding_box(self):
return BoundingBox.infinite()
def get_all_cells(self, memo=None):
return {}
def get_all_universes(self, memo=None):
return {}
def clone(self, clone_materials=True, clone_regions=True, memo=None):
warnings.warn("clone is not available for cells in a DAGMC universe")
return self
def plot(self, *args, **kwargs):
raise TypeError("plot is not available for DAGMC cells.")
def create_xml_subelement(self, xml_element, memo=None):
raise TypeError("create_xml_subelement is not available for DAGMC cells.")
@classmethod
def from_xml_element(cls, elem, surfaces, materials, get_universe):
raise TypeError("from_xml_element is not available for DAGMC cells.")

2
openmc/geometry.py
View file

@ -217,7 +217,7 @@ class Geometry:
# Add any DAGMC universes
for e in elem.findall('dagmc_universe'):
dag_univ = openmc.DAGMCUniverse.from_xml_element(e)
dag_univ = openmc.DAGMCUniverse.from_xml_element(e, mats)
universes[dag_univ.id] = dag_univ
# Dictionary that maps each universe to a list of cells/lattices that

1
openmc/lib/__init__.py
View file

@ -68,6 +68,7 @@ from .settings import settings
from .math import *
from .plot import *
from .weight_windows import *
from .dagmc import *
# Flag to denote whether or not openmc.lib.init has been called
# TODO: Establish and use a flag in the C++ code to represent the status of the

43
openmc/lib/dagmc.py Normal file
View file

@ -0,0 +1,43 @@
from ctypes import c_int, c_int32, POINTER, c_size_t
import numpy as np
from . import _dll
from .error import _error_handler
__all__ = [
'dagmc_universe_cell_ids'
]
# DAGMC functions
_dll.openmc_dagmc_universe_get_cell_ids.argtypes = [c_int32, POINTER(c_int32), POINTER(c_size_t)]
_dll.openmc_dagmc_universe_get_cell_ids.restype = c_int
_dll.openmc_dagmc_universe_get_cell_ids.errcheck = _error_handler
_dll.openmc_dagmc_universe_get_num_cells.argtypes = [c_int32, POINTER(c_size_t)]
_dll.openmc_dagmc_universe_get_num_cells.restype = c_int
_dll.openmc_dagmc_universe_get_num_cells.errcheck = _error_handler
def dagmc_universe_cell_ids(universe_id: int) -> np.ndarray:
"""Return an array of cell IDs for a DAGMC universe.
Parameters
----------
dagmc_id : int
ID of the DAGMC universe to get cell IDs from.
Returns
-------
numpy.ndarray
DAGMC cell IDs for the universe.
"""
n = c_size_t()
_dll.openmc_dagmc_universe_get_num_cells(universe_id, n)
cell_ids = np.empty(n.value, dtype=np.int32)
_dll.openmc_dagmc_universe_get_cell_ids(
universe_id, cell_ids.ctypes.data_as(POINTER(c_int32)), n
)
return cell_ids

119
openmc/model/model.py
View file

@ -1,7 +1,6 @@
from __future__ import annotations
from collections.abc import Iterable
from functools import lru_cache
import os
from pathlib import Path
from numbers import Integral
from tempfile import NamedTemporaryFile
@ -15,8 +14,9 @@ import openmc
import openmc._xml as xml
from openmc.dummy_comm import DummyCommunicator
from openmc.executor import _process_CLI_arguments
from openmc.checkvalue import check_type, check_value, PathLike
from openmc.checkvalue import check_type, check_value
from openmc.exceptions import InvalidIDError
import openmc.lib
from openmc.utility_funcs import change_directory
@ -324,6 +324,28 @@ class Model:
# communicator
openmc.lib.init(args=args, intracomm=intracomm, output=output)
def sync_dagmc_universes(self):
"""Synchronize all DAGMC universes in the current geometry.
This method iterates over all DAGMC universes in the geometry and
synchronizes their cells with the current material assignments. Requires
that the model has been initialized via :meth:`Model.init_lib`.
.. versionadded:: 0.15.1
"""
if self.is_initialized:
if self.materials:
materials = self.materials
else:
materials = list(self.geometry.get_all_materials().values())
for univ in self.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
univ.sync_dagmc_cells(materials)
else:
raise ValueError("The model must be initialized before calling "
"this method")
def finalize_lib(self):
"""Finalize simulation and free memory allocated for the C API
@ -1154,51 +1176,86 @@ class Model:
self._change_py_lib_attribs(names_or_ids, volume, 'material', 'volume')
def differentiate_depletable_mats(self, diff_volume_method: str):
def differentiate_depletable_mats(self, diff_volume_method: str = None):
"""Assign distribmats for each depletable material
.. versionadded:: 0.14.0
.. versionchanged:: 0.15.1
diff_volume_method default is None, do not set volumes on the new
material ovjects. Is now a convenience method for
differentiate_mats(diff_volume_method, depletable_only=True)
Parameters
----------
diff_volume_method : str
Specifies how the volumes of the new materials should be found.
Default is to 'divide equally' which divides the original material
volume equally between the new materials, 'match cell' sets the
volume of the material to volume of the cell they fill.
- None: Do not assign volumes to the new materials (Default)
- 'divide_equally': Divide the original material volume equally between the new materials
- 'match cell': Set the volume of the material to the volume of the cell they fill
"""
self.differentiate_mats(diff_volume_method, depletable_only=True)
def differentiate_mats(self, diff_volume_method: str = None, depletable_only: bool = True):
"""Assign distribmats for each material
.. versionadded:: 0.15.1
Parameters
----------
diff_volume_method : str
Specifies how the volumes of the new materials should be found.
- None: Do not assign volumes to the new materials (Default)
- 'divide_equally': Divide the original material volume equally between the new materials
- 'match cell': Set the volume of the material to the volume of the cell they fill
depletable_only : bool
Default is True, only depletable materials will be differentiated. If False, all materials will be
differentiated.
"""
check_value('volume differentiation method', diff_volume_method, ("divide equally", "match cell", None))
# Count the number of instances for each cell and material
self.geometry.determine_paths(instances_only=True)
# Extract all depletable materials which have multiple instances
distribmats = set(
[mat for mat in self.materials
if mat.depletable and mat.num_instances > 1])
if diff_volume_method == 'divide equally':
for mat in distribmats:
if mat.volume is None:
raise RuntimeError("Volume not specified for depletable "
f"material with ID={mat.id}.")
mat.volume /= mat.num_instances
if distribmats:
# Assign distribmats to cells
for cell in self.geometry.get_all_material_cells().values():
if cell.fill in distribmats:
mat = cell.fill
if diff_volume_method == 'divide equally':
cell.fill = [mat.clone() for _ in range(cell.num_instances)]
elif diff_volume_method == 'match cell':
for _ in range(cell.num_instances):
cell.fill = mat.clone()
# Find all or depletable_only materials which have multiple instance
distribmats = set()
for mat in self.materials:
# Differentiate all materials with multiple instances
diff_mat = mat.num_instances > 1
# If depletable_only is True, differentiate only depletable materials
if depletable_only:
diff_mat = diff_mat and mat.depletable
if diff_mat:
# Assign volumes to the materials according to requirements
if diff_volume_method == "divide equally":
if mat.volume is None:
raise RuntimeError(
"Volume not specified for "
f"material with ID={mat.id}.")
else:
mat.volume /= mat.num_instances
elif diff_volume_method == "match cell":
for cell in self.geometry.get_all_material_cells().values():
if cell.fill == mat:
if not cell.volume:
raise ValueError(
f"Volume of cell ID={cell.id} not specified. "
"Set volumes of cells prior to using "
"diff_volume_method='match cell'."
)
cell.fill.volume = cell.volume
"diff_volume_method='match cell'.")
distribmats.add(mat)
if not distribmats:
return
# Assign distribmats to cells
for cell in self.geometry.get_all_material_cells().values():
if cell.fill in distribmats:
mat = cell.fill
if diff_volume_method != 'match cell':
cell.fill = [mat.clone() for _ in range(cell.num_instances)]
elif diff_volume_method == 'match cell':
cell.fill = mat.clone()
cell.fill.volume = cell.volume
if self.materials is not None:
self.materials = openmc.Materials(

837
openmc/universe.py
View file

@ -2,22 +2,16 @@ from __future__ import annotations
import math
from abc import ABC, abstractmethod
from collections.abc import Iterable
from numbers import Integral, Real
from numbers import Real
from pathlib import Path
from tempfile import TemporaryDirectory
import warnings
import h5py
import lxml.etree as ET
import numpy as np
import openmc
import openmc.checkvalue as cv
from ._xml import get_text
from .checkvalue import check_type, check_value
from .mixin import IDManagerMixin
from .surface import _BOUNDARY_TYPES
from .utility_funcs import input_path
class UniverseBase(ABC, IDManagerMixin):
@ -55,6 +49,10 @@ class UniverseBase(ABC, IDManagerMixin):
def name(self):
return self._name
@property
def cells(self):
return self._cells
@name.setter
def name(self, name):
if name is not None:
@ -135,6 +133,130 @@ class UniverseBase(ABC, IDManagerMixin):
"""
def _determine_paths(self, path='', instances_only=False):
"""Count the number of instances for each cell in the universe, and
record the count in the :attr:`Cell.num_instances` properties."""
univ_path = path + f'u{self.id}'
for cell in self.cells.values():
cell_path = f'{univ_path}->c{cell.id}'
fill = cell._fill
fill_type = cell.fill_type
# If universe-filled, recursively count cells in filling universe
if fill_type == 'universe':
fill._determine_paths(cell_path + '->', instances_only)
# If lattice-filled, recursively call for all universes in lattice
elif fill_type == 'lattice':
latt = fill
# Count instances in each universe in the lattice
for index in latt._natural_indices:
latt_path = '{}->l{}({})->'.format(
cell_path, latt.id, ",".join(str(x) for x in index))
univ = latt.get_universe(index)
univ._determine_paths(latt_path, instances_only)
else:
if fill_type == 'material':
mat = fill
elif fill_type == 'distribmat':
mat = fill[cell._num_instances]
else:
mat = None
if mat is not None:
mat._num_instances += 1
if not instances_only:
mat._paths.append(f'{cell_path}->m{mat.id}')
# Append current path
cell._num_instances += 1
if not instances_only:
cell._paths.append(cell_path)
def add_cells(self, cells):
"""Add multiple cells to the universe.
Parameters
----------
cells : Iterable of openmc.Cell
Cells to add
"""
if not isinstance(cells, Iterable):
msg = f'Unable to add Cells to Universe ID="{self._id}" since ' \
f'"{cells}" is not iterable'
raise TypeError(msg)
for cell in cells:
self.add_cell(cell)
@abstractmethod
def add_cell(self, cell):
pass
@abstractmethod
def remove_cell(self, cell):
pass
def clear_cells(self):
"""Remove all cells from the universe."""
self._cells.clear()
def get_all_cells(self, memo=None):
"""Return all cells that are contained within the universe
Returns
-------
cells : dict
Dictionary whose keys are cell IDs and values are :class:`Cell`
instances
"""
if memo is None:
memo = set()
elif self in memo:
return {}
memo.add(self)
# Add this Universe's cells to the dictionary
cells = {}
cells.update(self._cells)
# Append all Cells in each Cell in the Universe to the dictionary
for cell in self._cells.values():
cells.update(cell.get_all_cells(memo))
return cells
def get_all_materials(self, memo=None):
"""Return all materials that are contained within the universe
Returns
-------
materials : dict
Dictionary whose keys are material IDs and values are
:class:`Material` instances
"""
if memo is None:
memo = set()
materials = {}
# Append all Cells in each Cell in the Universe to the dictionary
cells = self.get_all_cells(memo)
for cell in cells.values():
materials.update(cell.get_all_materials(memo))
return materials
@abstractmethod
def _partial_deepcopy(self):
"""Deepcopy all parameters of an openmc.UniverseBase object except its cells.
@ -182,93 +304,6 @@ class UniverseBase(ABC, IDManagerMixin):
return memo[self]
class Universe(UniverseBase):
"""A collection of cells that can be repeated.
Parameters
----------
universe_id : int, optional
Unique identifier of the universe. If not specified, an identifier will
automatically be assigned
name : str, optional
Name of the universe. If not specified, the name is the empty string.
cells : Iterable of openmc.Cell, optional
Cells to add to the universe. By default no cells are added.
Attributes
----------
id : int
Unique identifier of the universe
name : str
Name of the universe
cells : dict
Dictionary whose keys are cell IDs and values are :class:`Cell`
instances
volume : float
Volume of the universe in cm^3. This can either be set manually or
calculated in a stochastic volume calculation and added via the
:meth:`Universe.add_volume_information` method.
bounding_box : openmc.BoundingBox
Lower-left and upper-right coordinates of an axis-aligned bounding box
of the universe.
"""
def __init__(self, universe_id=None, name='', cells=None):
super().__init__(universe_id, name)
if cells is not None:
self.add_cells(cells)
def __repr__(self):
string = super().__repr__()
string += '{: <16}=\t{}\n'.format('\tGeom', 'CSG')
string += '{: <16}=\t{}\n'.format('\tCells', list(self._cells.keys()))
return string
@property
def cells(self):
return self._cells
@property
def bounding_box(self) -> openmc.BoundingBox:
regions = [c.region for c in self.cells.values()
if c.region is not None]
if regions:
return openmc.Union(regions).bounding_box
else:
return openmc.BoundingBox.infinite()
@classmethod
def from_hdf5(cls, group, cells):
"""Create universe from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
cells : dict
Dictionary mapping cell IDs to instances of :class:`openmc.Cell`.
Returns
-------
openmc.Universe
Universe instance
"""
universe_id = int(group.name.split('/')[-1].lstrip('universe '))
cell_ids = group['cells'][()]
# Create this Universe
universe = cls(universe_id)
# Add each Cell to the Universe
for cell_id in cell_ids:
universe.add_cell(cells[cell_id])
return universe
def find(self, point):
"""Find cells/universes/lattices which contain a given point
@ -528,67 +563,6 @@ class Universe(UniverseBase):
axes.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs)
return axes
def add_cell(self, cell):
"""Add a cell to the universe.
Parameters
----------
cell : openmc.Cell
Cell to add
"""
if not isinstance(cell, openmc.Cell):
msg = f'Unable to add a Cell to Universe ID="{self._id}" since ' \
f'"{cell}" is not a Cell'
raise TypeError(msg)
cell_id = cell.id
if cell_id not in self._cells:
self._cells[cell_id] = cell
def add_cells(self, cells):
"""Add multiple cells to the universe.
Parameters
----------
cells : Iterable of openmc.Cell
Cells to add
"""
if not isinstance(cells, Iterable):
msg = f'Unable to add Cells to Universe ID="{self._id}" since ' \
f'"{cells}" is not iterable'
raise TypeError(msg)
for cell in cells:
self.add_cell(cell)
def remove_cell(self, cell):
"""Remove a cell from the universe.
Parameters
----------
cell : openmc.Cell
Cell to remove
"""
if not isinstance(cell, openmc.Cell):
msg = f'Unable to remove a Cell from Universe ID="{self._id}" ' \
f'since "{cell}" is not a Cell'
raise TypeError(msg)
# If the Cell is in the Universe's list of Cells, delete it
self._cells.pop(cell.id, None)
def clear_cells(self):
"""Remove all cells from the universe."""
self._cells.clear()
def get_nuclides(self):
"""Returns all nuclides in the universe
@ -636,55 +610,128 @@ class Universe(UniverseBase):
return nuclides
def get_all_cells(self, memo=None):
"""Return all cells that are contained within the universe
Returns
-------
class Universe(UniverseBase):
"""A collection of cells that can be repeated.
Parameters
----------
universe_id : int, optional
Unique identifier of the universe. If not specified, an identifier will
automatically be assigned
name : str, optional
Name of the universe. If not specified, the name is the empty string.
cells : Iterable of openmc.Cell, optional
Cells to add to the universe. By default no cells are added.
Attributes
----------
id : int
Unique identifier of the universe
name : str
Name of the universe
cells : dict
Dictionary whose keys are cell IDs and values are :class:`Cell`
instances
volume : float
Volume of the universe in cm^3. This can either be set manually or
calculated in a stochastic volume calculation and added via the
:meth:`Universe.add_volume_information` method.
bounding_box : openmc.BoundingBox
Lower-left and upper-right coordinates of an axis-aligned bounding box
of the universe.
"""
def __init__(self, universe_id=None, name='', cells=None):
super().__init__(universe_id, name)
if cells is not None:
self.add_cells(cells)
def __repr__(self):
string = super().__repr__()
string += '{: <16}=\t{}\n'.format('\tGeom', 'CSG')
string += '{: <16}=\t{}\n'.format('\tCells', list(self._cells.keys()))
return string
@property
def bounding_box(self) -> openmc.BoundingBox:
regions = [c.region for c in self.cells.values()
if c.region is not None]
if regions:
return openmc.Union(regions).bounding_box
else:
return openmc.BoundingBox.infinite()
@classmethod
def from_hdf5(cls, group, cells):
"""Create universe from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
cells : dict
Dictionary whose keys are cell IDs and values are :class:`Cell`
instances
"""
if memo is None:
memo = set()
elif self in memo:
return {}
memo.add(self)
# Add this Universe's cells to the dictionary
cells = {}
cells.update(self._cells)
# Append all Cells in each Cell in the Universe to the dictionary
for cell in self._cells.values():
cells.update(cell.get_all_cells(memo))
return cells
def get_all_materials(self, memo=None):
"""Return all materials that are contained within the universe
Dictionary mapping cell IDs to instances of :class:`openmc.Cell`.
Returns
-------
materials : dict
Dictionary whose keys are material IDs and values are
:class:`Material` instances
openmc.Universe
Universe instance
"""
universe_id = int(group.name.split('/')[-1].lstrip('universe '))
cell_ids = group['cells'][()]
# Create this Universe
universe = cls(universe_id)
# Add each Cell to the Universe
for cell_id in cell_ids:
universe.add_cell(cells[cell_id])
return universe
def add_cell(self, cell):
"""Add a cell to the universe.
Parameters
----------
cell : openmc.Cell
Cell to add
"""
if memo is None:
memo = set()
if not isinstance(cell, openmc.Cell):
msg = f'Unable to add a Cell to Universe ID="{self._id}" since ' \
f'"{cell}" is not a Cell'
raise TypeError(msg)
materials = {}
cell_id = cell.id
# Append all Cells in each Cell in the Universe to the dictionary
cells = self.get_all_cells(memo)
for cell in cells.values():
materials.update(cell.get_all_materials(memo))
if cell_id not in self._cells:
self._cells[cell_id] = cell
return materials
def remove_cell(self, cell):
"""Remove a cell from the universe.
Parameters
----------
cell : openmc.Cell
Cell to remove
"""
if not isinstance(cell, openmc.Cell):
msg = f'Unable to remove a Cell from Universe ID="{self._id}" ' \
f'since "{cell}" is not a Cell'
raise TypeError(msg)
# If the Cell is in the Universe's list of Cells, delete it
self._cells.pop(cell.id, None)
def create_xml_subelement(self, xml_element, memo=None):
if memo is None:
@ -706,50 +753,6 @@ class Universe(UniverseBase):
cell_element.set("universe", str(self._id))
xml_element.append(cell_element)
def _determine_paths(self, path='', instances_only=False):
"""Count the number of instances for each cell in the universe, and
record the count in the :attr:`Cell.num_instances` properties."""
univ_path = path + f'u{self.id}'
for cell in self.cells.values():
cell_path = f'{univ_path}->c{cell.id}'
fill = cell._fill
fill_type = cell.fill_type
# If universe-filled, recursively count cells in filling universe
if fill_type == 'universe':
fill._determine_paths(cell_path + '->', instances_only)
# If lattice-filled, recursively call for all universes in lattice
elif fill_type == 'lattice':
latt = fill
# Count instances in each universe in the lattice
for index in latt._natural_indices:
latt_path = '{}->l{}({})->'.format(
cell_path, latt.id, ",".join(str(x) for x in index))
univ = latt.get_universe(index)
univ._determine_paths(latt_path, instances_only)
else:
if fill_type == 'material':
mat = fill
elif fill_type == 'distribmat':
mat = fill[cell._num_instances]
else:
mat = None
if mat is not None:
mat._num_instances += 1
if not instances_only:
mat._paths.append(f'{cell_path}->m{mat.id}')
# Append current path
cell._num_instances += 1
if not instances_only:
cell._paths.append(cell_path)
def _partial_deepcopy(self):
"""Clone all of the openmc.Universe object's attributes except for its cells,
as they are copied within the clone function. This should only to be
@ -758,363 +761,3 @@ class Universe(UniverseBase):
clone = openmc.Universe(name=self.name)
clone.volume = self.volume
return clone
class DAGMCUniverse(UniverseBase):
"""A reference to a DAGMC file to be used in the model.
.. versionadded:: 0.13.0
Parameters
----------
filename : path-like
Path to the DAGMC file used to represent this universe.
universe_id : int, optional
Unique identifier of the universe. If not specified, an identifier will
automatically be assigned.
name : str, optional
Name of the universe. If not specified, the name is the empty string.
auto_geom_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between CSG and DAGMC (False)
auto_mat_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between OpenMC and UWUW materials (False)
Attributes
----------
id : int
Unique identifier of the universe
name : str
Name of the universe
filename : str
Path to the DAGMC file used to represent this universe.
auto_geom_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between CSG and DAGMC (False)
auto_mat_ids : bool
Set IDs automatically on initialization (True) or report overlaps in ID
space between OpenMC and UWUW materials (False)
bounding_box : openmc.BoundingBox
Lower-left and upper-right coordinates of an axis-aligned bounding box
of the universe.
.. versionadded:: 0.13.1
material_names : list of str
Return a sorted list of materials names that are contained within the
DAGMC h5m file. This is useful when naming openmc.Material() objects
as each material name present in the DAGMC h5m file must have a
matching openmc.Material() with the same name.
.. versionadded:: 0.13.2
n_cells : int
The number of cells in the DAGMC model. This is the number of cells at
runtime and accounts for the implicit complement whether or not is it
present in the DAGMC file.
.. versionadded:: 0.13.2
n_surfaces : int
The number of surfaces in the model.
.. versionadded:: 0.13.2
"""
def __init__(self,
filename: cv.PathLike,
universe_id=None,
name='',
auto_geom_ids=False,
auto_mat_ids=False):
super().__init__(universe_id, name)
# Initialize class attributes
self.filename = filename
self.auto_geom_ids = auto_geom_ids
self.auto_mat_ids = auto_mat_ids
def __repr__(self):
string = super().__repr__()
string += '{: <16}=\t{}\n'.format('\tGeom', 'DAGMC')
string += '{: <16}=\t{}\n'.format('\tFile', self.filename)
return string
@property
def bounding_box(self):
with h5py.File(self.filename) as dagmc_file:
coords = dagmc_file['tstt']['nodes']['coordinates'][()]
lower_left_corner = coords.min(axis=0)
upper_right_corner = coords.max(axis=0)
return openmc.BoundingBox(lower_left_corner, upper_right_corner)
@property
def filename(self):
return self._filename
@filename.setter
def filename(self, val: cv.PathLike):
cv.check_type('DAGMC filename', val, cv.PathLike)
self._filename = input_path(val)
@property
def auto_geom_ids(self):
return self._auto_geom_ids
@auto_geom_ids.setter
def auto_geom_ids(self, val):
cv.check_type('DAGMC automatic geometry ids', val, bool)
self._auto_geom_ids = val
@property
def auto_mat_ids(self):
return self._auto_mat_ids
@auto_mat_ids.setter
def auto_mat_ids(self, val):
cv.check_type('DAGMC automatic material ids', val, bool)
self._auto_mat_ids = val
@property
def material_names(self):
dagmc_file_contents = h5py.File(self.filename)
material_tags_hex = dagmc_file_contents['/tstt/tags/NAME'].get(
'values')
material_tags_ascii = []
for tag in material_tags_hex:
candidate_tag = tag.tobytes().decode().replace('\x00', '')
# tags might be for temperature or reflective surfaces
if candidate_tag.startswith('mat:'):
# removes first 4 characters as openmc.Material name should be
# set without the 'mat:' part of the tag
material_tags_ascii.append(candidate_tag[4:])
return sorted(set(material_tags_ascii))
def get_all_cells(self, memo=None):
return {}
def get_all_materials(self, memo=None):
return {}
def _n_geom_elements(self, geom_type):
"""
Helper function for retrieving the number geometric entities in a DAGMC
file
Parameters
----------
geom_type : str
The type of geometric entity to count. One of {'Volume', 'Surface'}. Returns
the runtime number of voumes in the DAGMC model (includes implicit complement).
Returns
-------
int
Number of geometry elements of the specified type
"""
cv.check_value('geometry type', geom_type, ('volume', 'surface'))
def decode_str_tag(tag_val):
return tag_val.tobytes().decode().replace('\x00', '')
with h5py.File(self.filename) as dagmc_file:
category_data = dagmc_file['tstt/tags/CATEGORY/values']
category_strs = map(decode_str_tag, category_data)
n = sum([v == geom_type.capitalize() for v in category_strs])
# check for presence of an implicit complement in the file and
# increment the number of cells if it doesn't exist
if geom_type == 'volume':
name_data = dagmc_file['tstt/tags/NAME/values']
name_strs = map(decode_str_tag, name_data)
if not sum(['impl_complement' in n for n in name_strs]):
n += 1
return n
@property
def n_cells(self):
return self._n_geom_elements('volume')
@property
def n_surfaces(self):
return self._n_geom_elements('surface')
def create_xml_subelement(self, xml_element, memo=None):
if memo is None:
memo = set()
if self in memo:
return
memo.add(self)
# Set xml element values
dagmc_element = ET.Element('dagmc_universe')
dagmc_element.set('id', str(self.id))
if self.auto_geom_ids:
dagmc_element.set('auto_geom_ids', 'true')
if self.auto_mat_ids:
dagmc_element.set('auto_mat_ids', 'true')
dagmc_element.set('filename', str(self.filename))
xml_element.append(dagmc_element)
def bounding_region(
self,
bounded_type: str = 'box',
boundary_type: str = 'vacuum',
starting_id: int = 10000,
padding_distance: float = 0.
):
"""Creates a either a spherical or box shaped bounding region around
the DAGMC geometry.
.. versionadded:: 0.13.1
Parameters
----------
bounded_type : str
The type of bounding surface(s) to use when constructing the region.
Options include a single spherical surface (sphere) or a rectangle
made from six planes (box).
boundary_type : str
Boundary condition that defines the behavior for particles hitting
the surface. Defaults to vacuum boundary condition. Passed into the
surface construction.
starting_id : int
Starting ID of the surface(s) used in the region. For bounded_type
'box', the next 5 IDs will also be used. Defaults to 10000 to reduce
the chance of an overlap of surface IDs with the DAGMC geometry.
padding_distance : float
Distance between the bounding region surfaces and the minimal
bounding box. Allows for the region to be larger than the DAGMC
geometry.
Returns
-------
openmc.Region
Region instance
"""
check_type('boundary type', boundary_type, str)
check_value('boundary type', boundary_type, _BOUNDARY_TYPES)
check_type('starting surface id', starting_id, Integral)
check_type('bounded type', bounded_type, str)
check_value('bounded type', bounded_type, ('box', 'sphere'))
bbox = self.bounding_box.expand(padding_distance, True)
if bounded_type == 'sphere':
radius = np.linalg.norm(bbox.upper_right - bbox.center)
bounding_surface = openmc.Sphere(
surface_id=starting_id,
x0=bbox.center[0],
y0=bbox.center[1],
z0=bbox.center[2],
boundary_type=boundary_type,
r=radius,
)
return -bounding_surface
if bounded_type == 'box':
# defines plane surfaces for all six faces of the bounding box
lower_x = openmc.XPlane(bbox[0][0], surface_id=starting_id)
upper_x = openmc.XPlane(bbox[1][0], surface_id=starting_id+1)
lower_y = openmc.YPlane(bbox[0][1], surface_id=starting_id+2)
upper_y = openmc.YPlane(bbox[1][1], surface_id=starting_id+3)
lower_z = openmc.ZPlane(bbox[0][2], surface_id=starting_id+4)
upper_z = openmc.ZPlane(bbox[1][2], surface_id=starting_id+5)
region = +lower_x & -upper_x & +lower_y & -upper_y & +lower_z & -upper_z
for surface in region.get_surfaces().values():
surface.boundary_type = boundary_type
return region
def bounded_universe(self, bounding_cell_id=10000, **kwargs):
"""Returns an openmc.Universe filled with this DAGMCUniverse and bounded
with a cell. Defaults to a box cell with a vacuum surface however this
can be changed using the kwargs which are passed directly to
DAGMCUniverse.bounding_region().
Parameters
----------
bounding_cell_id : int
The cell ID number to use for the bounding cell, defaults to 10000 to reduce
the chance of overlapping ID numbers with the DAGMC geometry.
Returns
-------
openmc.Universe
Universe instance
"""
bounding_cell = openmc.Cell(
fill=self, cell_id=bounding_cell_id, region=self.bounding_region(**kwargs))
return openmc.Universe(cells=[bounding_cell])
@classmethod
def from_hdf5(cls, group):
"""Create DAGMC universe from HDF5 group
Parameters
----------
group : h5py.Group
Group in HDF5 file
Returns
-------
openmc.DAGMCUniverse
DAGMCUniverse instance
"""
id = int(group.name.split('/')[-1].lstrip('universe '))
fname = group['filename'][()].decode()
name = group['name'][()].decode() if 'name' in group else None
out = cls(fname, universe_id=id, name=name)
out.auto_geom_ids = bool(group.attrs['auto_geom_ids'])
out.auto_mat_ids = bool(group.attrs['auto_mat_ids'])
return out
@classmethod
def from_xml_element(cls, elem):
"""Generate DAGMC universe from XML element
Parameters
----------
elem : lxml.etree._Element
`<dagmc_universe>` element
Returns
-------
openmc.DAGMCUniverse
DAGMCUniverse instance
"""
id = int(get_text(elem, 'id'))
fname = get_text(elem, 'filename')
out = cls(fname, universe_id=id)
name = get_text(elem, 'name')
if name is not None:
out.name = name
out.auto_geom_ids = bool(elem.get('auto_geom_ids'))
out.auto_mat_ids = bool(elem.get('auto_mat_ids'))
return out
def _partial_deepcopy(self):
"""Clone all of the openmc.DAGMCUniverse object's attributes except for
its cells, as they are copied within the clone function. This should
only to be used within the openmc.UniverseBase.clone() context.
"""
clone = openmc.DAGMCUniverse(name=self.name, filename=self.filename)
clone.volume = self.volume
clone.auto_geom_ids = self.auto_geom_ids
clone.auto_mat_ids = self.auto_mat_ids
return clone

4
src/cell.cpp
View file

@ -252,12 +252,12 @@ void Cell::to_hdf5(hid_t cell_group) const
// default constructor
CSGCell::CSGCell()
{
geom_type_ = GeometryType::CSG;
geom_type() = GeometryType::CSG;
}
CSGCell::CSGCell(pugi::xml_node cell_node)
{
geom_type_ = GeometryType::CSG;
geom_type() = GeometryType::CSG;
if (check_for_node(cell_node, "id")) {
id_ = std::stoi(get_node_value(cell_node, "id"));

102
src/dagmc.cpp
View file

@ -72,6 +72,23 @@ DAGUniverse::DAGUniverse(pugi::xml_node node)
adjust_material_ids_ = get_node_value_bool(node, "auto_mat_ids");
}
// get material assignment overloading
if (check_for_node(node, "material_overrides")) {
auto mat_node = node.child("material_overrides");
// loop over all subelements (each subelement corresponds to a material)
for (pugi::xml_node cell_node : mat_node.children("cell_override")) {
// Store assignment reference name
int32_t ref_assignment = std::stoi(get_node_value(cell_node, "id"));
// Get mat name for each assignement instances
vector<int32_t> instance_mats =
get_node_array<int32_t>(cell_node, "material_ids");
// Store mat name for each instances
material_overrides_.emplace(ref_assignment, instance_mats);
}
}
initialize();
}
@ -211,12 +228,13 @@ void DAGUniverse::init_geometry()
if (mat_str == "graveyard") {
graveyard = vol_handle;
}
// material void checks
if (mat_str == "void" || mat_str == "vacuum" || mat_str == "graveyard") {
c->material_.push_back(MATERIAL_VOID);
} else {
if (uses_uwuw()) {
if (material_overrides_.count(c->id_)) {
override_assign_material(c);
} else if (uses_uwuw()) {
uwuw_assign_material(vol_handle, c);
} else {
legacy_assign_material(mat_str, c);
@ -609,6 +627,33 @@ void DAGUniverse::uwuw_assign_material(
fatal_error("DAGMC was not configured with UWUW.");
#endif // OPENMC_UWUW
}
void DAGUniverse::override_assign_material(std::unique_ptr<DAGCell>& c) const
{
// if Cell ID matches an override key, use it to override the material
// assignment else if UWUW is used, get the material assignment from the DAGMC
// metadata
// Notify User that an override is being applied on a DAGMCCell
write_message(fmt::format("Applying override for DAGMCCell {}", c->id_), 8);
if (settings::verbosity >= 10) {
auto msg = fmt::format("Assigning DAGMC cell {} material(s) based on "
"override information (see input XML).",
c->id_);
write_message(msg, 10);
}
// Override the material assignment for each cell instance using the legacy
// assignement
for (auto mat_id : material_overrides_.at(c->id_)) {
if (model::material_map.find(mat_id) == model::material_map.end()) {
fatal_error(fmt::format(
"Material with ID '{}' not found for DAGMC cell {}", mat_id, c->id_));
}
c->material_.push_back(mat_id);
}
}
//==============================================================================
// DAGMC Cell implementation
//==============================================================================
@ -616,7 +661,7 @@ void DAGUniverse::uwuw_assign_material(
DAGCell::DAGCell(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx)
: Cell {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx)
{
geom_type_ = GeometryType::DAG;
geom_type() = GeometryType::DAG;
};
std::pair<double, int32_t> DAGCell::distance(
@ -719,7 +764,7 @@ BoundingBox DAGCell::bounding_box() const
DAGSurface::DAGSurface(std::shared_ptr<moab::DagMC> dag_ptr, int32_t dag_idx)
: Surface {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx)
{
geom_type_ = GeometryType::DAG;
geom_type() = GeometryType::DAG;
} // empty constructor
moab::EntityHandle DAGSurface::mesh_handle() const
@ -818,12 +863,61 @@ int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ)
return univp->cell_index(new_vol);
}
extern "C" int openmc_dagmc_universe_get_cell_ids(
int32_t univ_id, int32_t* ids, size_t* n)
{
// make sure the universe id is a DAGMC Universe
const auto& univ = model::universes[model::universe_map[univ_id]];
if (univ->geom_type() != GeometryType::DAG) {
set_errmsg(fmt::format("Universe {} is not a DAGMC Universe", univ_id));
return OPENMC_E_INVALID_TYPE;
}
std::vector<int32_t> dag_cell_ids;
for (const auto& cell_index : univ->cells_) {
const auto& cell = model::cells[cell_index];
if (cell->geom_type() == GeometryType::CSG) {
set_errmsg(fmt::format("Cell {} is not a DAGMC Cell", cell->id_));
return OPENMC_E_INVALID_TYPE;
}
dag_cell_ids.push_back(cell->id_);
}
std::copy(dag_cell_ids.begin(), dag_cell_ids.end(), ids);
*n = dag_cell_ids.size();
return 0;
}
extern "C" int openmc_dagmc_universe_get_num_cells(int32_t univ_id, size_t* n)
{
// make sure the universe id is a DAGMC Universe
const auto& univ = model::universes[model::universe_map[univ_id]];
if (univ->geom_type() != GeometryType::DAG) {
set_errmsg(fmt::format("Universe {} is not a DAGMC universe", univ_id));
return OPENMC_E_INVALID_TYPE;
}
*n = univ->cells_.size();
return 0;
}
} // namespace openmc
#else
namespace openmc {
extern "C" int openmc_dagmc_universe_get_cell_ids(
int32_t univ_id, int32_t* ids, size_t* n)
{
set_errmsg("OpenMC was not configured with DAGMC");
return OPENMC_E_UNASSIGNED;
};
extern "C" int openmc_dagmc_universe_get_num_cells(int32_t univ_id, size_t* n)
{
set_errmsg("OpenMC was not configured with DAGMC");
return OPENMC_E_UNASSIGNED;
};
void read_dagmc_universes(pugi::xml_node node)
{
if (check_for_node(node, "dagmc_universe")) {

7
src/particle.cpp
View file

@ -533,7 +533,7 @@ void Particle::cross_surface(const Surface& surf)
// if we're crossing a CSG surface, make sure the DAG history is reset
#ifdef DAGMC
if (surf.geom_type_ == GeometryType::CSG)
if (surf.geom_type() == GeometryType::CSG)
history().reset();
#endif
@ -548,7 +548,7 @@ void Particle::cross_surface(const Surface& surf)
#ifdef DAGMC
// in DAGMC, we know what the next cell should be
if (surf.geom_type_ == GeometryType::DAG) {
if (surf.geom_type() == GeometryType::DAG) {
int32_t i_cell = next_cell(std::abs(surface()), cell_last(n_coord() - 1),
lowest_coord().universe) -
1;
@ -668,7 +668,8 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
// the lower universes.
// (unless we're using a dagmc model, which has exactly one universe)
n_coord() = 1;
if (surf.geom_type_ != GeometryType::DAG && !neighbor_list_find_cell(*this)) {
if (surf.geom_type() != GeometryType::DAG &&
!neighbor_list_find_cell(*this)) {
mark_as_lost("Couldn't find particle after reflecting from surface " +
std::to_string(surf.id_) + ".");
return;

2
src/plot.cpp
View file

@ -1303,7 +1303,7 @@ void ProjectionPlot::create_output() const
int32_t i_surface = std::abs(p.surface()) - 1;
if (i_surface > 0 &&
model::surfaces[i_surface]->geom_type_ == GeometryType::DAG) {
model::surfaces[i_surface]->geom_type() == GeometryType::DAG) {
#ifdef DAGMC
int32_t i_cell = next_cell(i_surface,
p.cell_last(p.n_coord() - 1), p.lowest_coord().universe);

8
src/surface.cpp
View file

@ -165,9 +165,9 @@ void Surface::to_hdf5(hid_t group_id) const
{
hid_t surf_group = create_group(group_id, fmt::format("surface {}", id_));
if (geom_type_ == GeometryType::DAG) {
if (geom_type() == GeometryType::DAG) {
write_string(surf_group, "geom_type", "dagmc", false);
} else if (geom_type_ == GeometryType::CSG) {
} else if (geom_type() == GeometryType::CSG) {
write_string(surf_group, "geom_type", "csg", false);
if (bc_) {
@ -189,11 +189,11 @@ void Surface::to_hdf5(hid_t group_id) const
CSGSurface::CSGSurface() : Surface {}
{
geom_type_ = GeometryType::CSG;
geom_type() = GeometryType::CSG;
};
CSGSurface::CSGSurface(pugi::xml_node surf_node) : Surface {surf_node}
{
geom_type_ = GeometryType::CSG;
geom_type() = GeometryType::CSG;
};
//==============================================================================

256
tests/unit_tests/dagmc/test_model.py Normal file
View file

@ -0,0 +1,256 @@
from pathlib import Path
import lxml.etree as ET
import numpy as np
import pytest
import openmc
from openmc.utility_funcs import change_directory
pytestmark = pytest.mark.skipif(
not openmc.lib._dagmc_enabled(),
reason="DAGMC CAD geometry is not enabled.")
@pytest.fixture()
def model(request):
pitch = 1.26
mats = {}
mats["no-void fuel"] = openmc.Material(1, name="no-void fuel")
mats["no-void fuel"].add_nuclide("U235", 0.03)
mats["no-void fuel"].add_nuclide("U238", 0.97)
mats["no-void fuel"].add_nuclide("O16", 2.0)
mats["no-void fuel"].set_density("g/cm3", 10.0)
mats["41"] = openmc.Material(name="41")
mats["41"].add_nuclide("H1", 2.0)
mats["41"].add_element("O", 1.0)
mats["41"].set_density("g/cm3", 1.0)
mats["41"].add_s_alpha_beta("c_H_in_H2O")
p = Path(request.fspath).parent / "dagmc.h5m"
daguniv = openmc.DAGMCUniverse(p, auto_geom_ids=True)
lattice = openmc.RectLattice()
lattice.dimension = [2, 2]
lattice.lower_left = [-pitch, -pitch]
lattice.pitch = [pitch, pitch]
lattice.universes = [
[daguniv, daguniv],
[daguniv, daguniv]]
box = openmc.model.RectangularParallelepiped(-pitch, pitch, -pitch, pitch, -5, 5)
root = openmc.Universe(cells=[openmc.Cell(region=-box, fill=lattice)])
settings = openmc.Settings()
settings.batches = 100
settings.inactive = 10
settings.particles = 1000
ll, ur = root.bounding_box
mat_vol = openmc.VolumeCalculation([mats["no-void fuel"]], 1000000, ll, ur)
cell_vol = openmc.VolumeCalculation(list(root.cells.values()), 1000000, ll, ur)
settings.volume_calculations = [mat_vol, cell_vol]
model = openmc.Model()
model.materials = openmc.Materials(mats.values())
model.geometry = openmc.Geometry(root=root)
model.settings = settings
with change_directory(tmpdir=True):
try:
model.init_lib()
model.sync_dagmc_universes()
yield model
finally:
model.finalize_lib()
openmc.reset_auto_ids()
def test_dagmc_replace_material_assignment(model):
mats = {}
mats["foo"] = openmc.Material(name="foo")
mats["foo"].add_nuclide("H1", 2.0)
mats["foo"].add_element("O", 1.0)
mats["foo"].set_density("g/cm3", 1.0)
mats["foo"].add_s_alpha_beta("c_H_in_H2O")
for univ in model.geometry.get_all_universes().values():
if not isinstance(univ, openmc.DAGMCUniverse):
break
cells_with_41 = []
for cell in univ.cells.values():
if cell.fill is None:
continue
if cell.fill.name == "41":
cells_with_41.append(cell.id)
univ.replace_material_assignment("41", mats["foo"])
for cell_id in cells_with_41:
assert univ.cells[cell_id] == mats["foo"]
def test_dagmc_add_material_override_with_id(model):
mats = {}
mats["foo"] = openmc.Material(name="foo")
mats["foo"].add_nuclide("H1", 2.0)
mats["foo"].add_element("O", 1.0)
mats["foo"].set_density("g/cm3", 1.0)
mats["foo"].add_s_alpha_beta("c_H_in_H2O")
for univ in model.geometry.get_all_universes().values():
if not isinstance(univ, openmc.DAGMCUniverse):
break
cells_with_41 = []
for cell in univ.cells.values():
if cell.fill is None:
continue
if cell.fill.name == "41":
cells_with_41.append(cell.id)
univ.add_material_override(cell.id, mats["foo"])
for cell_id in cells_with_41:
assert univ.cells[cell_id] == mats["foo"]
def test_dagmc_add_material_override_with_cell(model):
mats = {}
mats["foo"] = openmc.Material(name="foo")
mats["foo"].add_nuclide("H1", 2.0)
mats["foo"].add_element("O", 1.0)
mats["foo"].set_density("g/cm3", 1.0)
mats["foo"].add_s_alpha_beta("c_H_in_H2O")
for univ in model.geometry.get_all_universes().values():
if not isinstance(univ, openmc.DAGMCUniverse):
break
cells_with_41 = []
for cell in univ.cells.values():
if cell.fill is None:
continue
if cell.fill.name == "41":
cells_with_41.append(cell.id)
univ.add_material_override(cell, mats["foo"])
for cell_id in cells_with_41:
assert univ.cells[cell_id] == mats["foo"]
def test_model_differentiate_depletable_with_dagmc(model, run_in_tmpdir):
model.calculate_volumes()
# Get the volume of the no-void fuel material before differentiation
volume_before = np.sum([m.volume for m in model.materials if m.name == "no-void fuel"])
# Differentiate the depletable materials
model.differentiate_depletable_mats(diff_volume_method="divide equally")
# Get the volume of the no-void fuel material after differentiation
volume_after = np.sum([m.volume for m in model.materials if "fuel" in m.name])
assert np.isclose(volume_before, volume_after)
assert len(model.materials) == 4*2 +1
def test_model_differentiate_with_dagmc(model):
root = model.geometry.root_universe
ll, ur = root.bounding_box
model.calculate_volumes()
# Get the volume of the no-void fuel material before differentiation
volume_before = np.sum([m.volume for m in model.materials if m.name == "no-void fuel"])
# Differentiate all the materials
model.differentiate_mats(depletable_only=False)
# Get the volume of the no-void fuel material after differentiation
mat_vol = openmc.VolumeCalculation(model.materials, 1000000, ll, ur)
model.settings.volume_calculations = [mat_vol]
model.init_lib() # need to reinitialize the lib after differentiating the materials
model.calculate_volumes()
volume_after = np.sum([m.volume for m in model.materials if "fuel" in m.name])
assert np.isclose(volume_before, volume_after)
assert len(model.materials) == 4*2 + 4
def test_bad_override_cell_id(model):
for univ in model.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
break
with pytest.raises(ValueError, match="Cell ID '1' not found in DAGMC universe"):
univ.material_overrides = {1: model.materials[0]}
def test_bad_override_type(model):
not_a_dag_cell = openmc.Cell()
for univ in model.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
break
with pytest.raises(ValueError, match="Unrecognized key type. Must be an integer or openmc.DAGMCCell object"):
univ.material_overrides = {not_a_dag_cell: model.materials[0]}
def test_bad_replacement_mat_name(model):
for univ in model.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
break
with pytest.raises(ValueError, match="No material with name 'not_a_mat' found in the DAGMC universe"):
univ.replace_material_assignment("not_a_mat", model.materials[0])
def test_dagmc_xml(model):
# Set the environment
mats = {}
mats["no-void fuel"] = openmc.Material(1, name="no-void fuel")
mats["no-void fuel"].add_nuclide("U235", 0.03)
mats["no-void fuel"].add_nuclide("U238", 0.97)
mats["no-void fuel"].add_nuclide("O16", 2.0)
mats["no-void fuel"].set_density("g/cm3", 10.0)
mats[5] = openmc.Material(name="41")
mats[5].add_nuclide("H1", 2.0)
mats[5].add_element("O", 1.0)
mats[5].set_density("g/cm3", 1.0)
mats[5].add_s_alpha_beta("c_H_in_H2O")
for univ in model.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
dag_univ = univ
break
for k, v in mats.items():
if isinstance(k, int):
dag_univ.add_material_override(k, v)
model.materials.append(v)
elif isinstance(k, str):
dag_univ.replace_material_assignment(k, v)
# Tesing the XML subelement generation
root = ET.Element('dagmc_universe')
dag_univ.create_xml_subelement(root)
dagmc_ele = root.find('dagmc_universe')
assert dagmc_ele.get('id') == str(dag_univ.id)
assert dagmc_ele.get('filename') == str(dag_univ.filename)
assert dagmc_ele.get('auto_geom_ids') == str(dag_univ.auto_geom_ids).lower()
override_eles = dagmc_ele.find('material_overrides').findall('cell_override')
assert len(override_eles) == 4
for i, override_ele in enumerate(override_eles):
cell_id = override_ele.get('id')
assert dag_univ.material_overrides[int(cell_id)][0].id == int(override_ele.find('material_ids').text)
model.export_to_model_xml()
xml_model = openmc.Model.from_model_xml()
for univ in xml_model.geometry.get_all_universes().values():
if isinstance(univ, openmc.DAGMCUniverse):
xml_dagmc_univ = univ
break
assert xml_dagmc_univ._material_overrides.keys() == dag_univ._material_overrides.keys()
for xml_mats, model_mats in zip(xml_dagmc_univ._material_overrides.values(), dag_univ._material_overrides.values()):
assert all([xml_mat.id == orig_mat.id for xml_mat, orig_mat in zip(xml_mats, model_mats)])