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

Compute homogenized materials over mesh elements (#2971)

Browse source
This commit is contained in:
Paul Romano 2024-05-02 12:47:10 -05:00 • committed by GitHub
parent c976653abf
commit 6e57f1dc72
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
7 changed files with 212 additions and 70 deletions
Download patch file Download diff file Expand all files Collapse all files

22
openmc/deplete/abc.py
View file

@ -10,8 +10,6 @@ from collections.abc import Iterable, Callable
from copy import deepcopy
from inspect import signature
from numbers import Real, Integral
from contextlib import contextmanager
import os
from pathlib import Path
import time
from typing import Optional, Union, Sequence
@ -22,6 +20,7 @@ from uncertainties import ufloat
from openmc.checkvalue import check_type, check_greater_than, PathLike
from openmc.mpi import comm
from openmc.utility_funcs import change_directory
from openmc import Material
from .stepresult import StepResult
from .chain import Chain
@ -61,25 +60,6 @@ except AttributeError:
# Can't set __doc__ on properties on Python 3.4
pass
@contextmanager
def change_directory(output_dir):
"""
Helper function for managing the current directory.
Parameters
----------
output_dir : pathlib.Path
Directory to switch to.
"""
orig_dir = os.getcwd()
output_dir = Path(output_dir)
output_dir.mkdir(parents=True, exist_ok=True)
try:
os.chdir(output_dir)
yield
finally:
os.chdir(orig_dir)
class TransportOperator(ABC):
"""Abstract class defining a transport operator

1
openmc/deplete/independent_operator.py
View file

@ -130,7 +130,6 @@ class IndependentOperator(OpenMCOperator):
# Validate micro-xs parameters
check_type('materials', materials, openmc.Materials)
check_type('micros', micros, Iterable, MicroXS)
check_type('fluxes', fluxes, Iterable, float)
if not (len(fluxes) == len(micros) == len(materials)):
msg = (f'The length of fluxes ({len(fluxes)}) should be equal to '

59
openmc/deplete/microxs.py
View file

@ -13,11 +13,11 @@ import numpy as np
from openmc.checkvalue import check_type, check_value, check_iterable_type, PathLike
from openmc.exceptions import DataError
from openmc.utility_funcs import change_directory
from openmc import StatePoint
from openmc.mgxs import GROUP_STRUCTURES
from openmc.data import REACTION_MT
import openmc
from .abc import change_directory
from .chain import Chain, REACTIONS
from .coupled_operator import _find_cross_sections, _get_nuclides_with_data
import openmc.lib
@ -284,38 +284,37 @@ class MicroXS:
# Create 2D array for microscopic cross sections
microxs_arr = np.zeros((len(nuclides), len(mts)))
with TemporaryDirectory() as tmpdir:
# Create a material with all nuclides
mat_all_nucs = openmc.Material()
for nuc in nuclides:
if nuc in nuclides_with_data:
mat_all_nucs.add_nuclide(nuc, 1.0)
mat_all_nucs.set_density("atom/b-cm", 1.0)
# Create a material with all nuclides
mat_all_nucs = openmc.Material()
for nuc in nuclides:
if nuc in nuclides_with_data:
mat_all_nucs.add_nuclide(nuc, 1.0)
mat_all_nucs.set_density("atom/b-cm", 1.0)
# Create simple model containing the above material
surf1 = openmc.Sphere(boundary_type="vacuum")
surf1_cell = openmc.Cell(fill=mat_all_nucs, region=-surf1)
model = openmc.Model()
model.geometry = openmc.Geometry([surf1_cell])
model.settings = openmc.Settings(
particles=1, batches=1, output={'summary': False})
# Create simple model containing the above material
surf1 = openmc.Sphere(boundary_type="vacuum")
surf1_cell = openmc.Cell(fill=mat_all_nucs, region=-surf1)
model = openmc.Model()
model.geometry = openmc.Geometry([surf1_cell])
model.settings = openmc.Settings(
particles=1, batches=1, output={'summary': False})
with change_directory(tmpdir):
# Export model within temporary directory
model.export_to_model_xml()
with change_directory(tmpdir=True):
# Export model within temporary directory
model.export_to_model_xml()
with openmc.lib.run_in_memory(**init_kwargs):
# For each nuclide and reaction, compute the flux-averaged
# cross section
for nuc_index, nuc in enumerate(nuclides):
if nuc not in nuclides_with_data:
continue
lib_nuc = openmc.lib.nuclides[nuc]
for mt_index, mt in enumerate(mts):
xs = lib_nuc.collapse_rate(
mt, temperature, energies, multigroup_flux
)
microxs_arr[nuc_index, mt_index] = xs
with openmc.lib.run_in_memory(**init_kwargs):
# For each nuclide and reaction, compute the flux-averaged
# cross section
for nuc_index, nuc in enumerate(nuclides):
if nuc not in nuclides_with_data:
continue
lib_nuc = openmc.lib.nuclides[nuc]
for mt_index, mt in enumerate(mts):
xs = lib_nuc.collapse_rate(
mt, temperature, energies, multigroup_flux
)
microxs_arr[nuc_index, mt_index] = xs
return cls(microxs_arr, nuclides, reactions)

92
openmc/mesh.py
View file

@ -7,7 +7,8 @@ from functools import wraps
from math import pi, sqrt, atan2
from numbers import Integral, Real
from pathlib import Path
from typing import Optional, Sequence, Tuple
import tempfile
from typing import Optional, Sequence, Tuple, List
import h5py
import lxml.etree as ET
@ -16,6 +17,7 @@ import numpy as np
import openmc
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from openmc.utility_funcs import change_directory
from ._xml import get_text
from .mixin import IDManagerMixin
from .surface import _BOUNDARY_TYPES
@ -145,6 +147,94 @@ class MeshBase(IDManagerMixin, ABC):
else:
raise ValueError(f'Unrecognized mesh type "{mesh_type}" found.')
def get_homogenized_materials(
self,
model: openmc.Model,
n_samples: int = 10_000,
prn_seed: Optional[int] = None,
**kwargs
) -> List[openmc.Material]:
"""Generate homogenized materials over each element in a mesh.
.. versionadded:: 0.14.1
Parameters
----------
model : openmc.Model
Model containing materials to be homogenized and the associated
geometry.
n_samples : int
Number of samples in each mesh element.
prn_seed : int, optional
Pseudorandom number generator (PRNG) seed; if None, one will be
generated randomly.
**kwargs
Keyword-arguments passed to :func:`openmc.lib.init`.
Returns
-------
list of openmc.Material
Homogenized material in each mesh element
"""
import openmc.lib
with change_directory(tmpdir=True):
# In order to get mesh into model, we temporarily replace the
# tallies with a single mesh tally using the current mesh
original_tallies = model.tallies
new_tally = openmc.Tally()
new_tally.filters = [openmc.MeshFilter(self)]
new_tally.scores = ['flux']
model.tallies = [new_tally]
# Export model to XML
model.export_to_model_xml()
# Get material volume fractions
openmc.lib.init(**kwargs)
mesh = openmc.lib.tallies[new_tally.id].filters[0].mesh
mat_volume_by_element = [
[
(mat.id if mat is not None else None, volume)
for mat, volume in mat_volume_list
]
for mat_volume_list in mesh.material_volumes(n_samples, prn_seed)
]
openmc.lib.finalize()
# Restore original tallies
model.tallies = original_tallies
# Create homogenized material for each element
materials = model.geometry.get_all_materials()
homogenized_materials = []
for mat_volume_list in mat_volume_by_element:
material_ids, volumes = [list(x) for x in zip(*mat_volume_list)]
total_volume = sum(volumes)
# Check for void material and remove
try:
index_void = material_ids.index(None)
except ValueError:
pass
else:
material_ids.pop(index_void)
volumes.pop(index_void)
# Compute volume fractions
volume_fracs = np.array(volumes) / total_volume
# Get list of materials and mix 'em up!
mats = [materials[uid] for uid in material_ids]
homogenized_mat = openmc.Material.mix_materials(
mats, volume_fracs, 'vo'
)
homogenized_mat.volume = total_volume
homogenized_materials.append(homogenized_mat)
return homogenized_materials
class StructuredMesh(MeshBase):
"""A base class for structured mesh functionality

22
openmc/model/model.py
View file

@ -1,6 +1,5 @@
from __future__ import annotations
from collections.abc import Iterable
from contextlib import contextmanager
from functools import lru_cache
import os
from pathlib import Path
@ -18,18 +17,7 @@ from openmc.dummy_comm import DummyCommunicator
from openmc.executor import _process_CLI_arguments
from openmc.checkvalue import check_type, check_value, PathLike
from openmc.exceptions import InvalidIDError
@contextmanager
def _change_directory(working_dir):
"""A context manager for executing in a provided working directory"""
start_dir = Path.cwd()
Path.mkdir(working_dir, parents=True, exist_ok=True)
os.chdir(working_dir)
try:
yield
finally:
os.chdir(start_dir)
from openmc.utility_funcs import change_directory
class Model:
@ -395,7 +383,7 @@ class Model:
# Store whether or not the library was initialized when we started
started_initialized = self.is_initialized
with _change_directory(Path(directory)):
with change_directory(directory):
with openmc.lib.quiet_dll(output):
# TODO: Support use of IndependentOperator too
depletion_operator = dep.CoupledOperator(self, **op_kwargs)
@ -677,7 +665,7 @@ class Model:
last_statepoint = None
# Operate in the provided working directory
with _change_directory(Path(cwd)):
with change_directory(cwd):
if self.is_initialized:
# Handle the run options as applicable
# First dont allow ones that must be set via init
@ -767,7 +755,7 @@ class Model:
raise ValueError("The Settings.volume_calculations attribute must"
" be specified before executing this method!")
with _change_directory(Path(cwd)):
with change_directory(cwd):
if self.is_initialized:
if threads is not None:
msg = "Threads must be set via Model.is_initialized(...)"
@ -829,7 +817,7 @@ class Model:
raise ValueError("The Model.plots attribute must be specified "
"before executing this method!")
with _change_directory(Path(cwd)):
with change_directory(cwd):
if self.is_initialized:
# Compute the volumes
openmc.lib.plot_geometry(output)

38
openmc/utility_funcs.py Normal file
View file

@ -0,0 +1,38 @@
from contextlib import contextmanager
import os
from pathlib import Path
from tempfile import TemporaryDirectory
from typing import Optional
from .checkvalue import PathLike
@contextmanager
def change_directory(working_dir: Optional[PathLike] = None, *, tmpdir: bool = False):
"""Context manager for executing in a provided working directory
Parameters
----------
working_dir : path-like
Directory to switch to.
tmpdir : bool
Whether to use a temporary directory instead of a specific working directory
"""
orig_dir = Path.cwd()
# Set up temporary directory if requested
if tmpdir:
tmp = TemporaryDirectory()
working_dir = tmp.name
elif working_dir is None:
raise ValueError('Must pass working_dir argument or specify tmpdir=True.')
working_dir = Path(working_dir)
working_dir.mkdir(parents=True, exist_ok=True)
os.chdir(working_dir)
try:
yield
finally:
os.chdir(orig_dir)
if tmpdir:
tmp.cleanup()

48
tests/unit_tests/test_mesh.py
View file

@ -299,6 +299,7 @@ def test_CylindricalMesh_get_indices_at_coords():
assert mesh.get_indices_at_coords([98, 199.9, 299]) == (0, 2, 0) # third angle quadrant
assert mesh.get_indices_at_coords([102, 199.1, 299]) == (0, 3, 0) # forth angle quadrant
def test_umesh_roundtrip(run_in_tmpdir, request):
umesh = openmc.UnstructuredMesh(request.path.parent / 'test_mesh_tets.e', 'moab')
umesh.output = True
@ -317,3 +318,50 @@ def test_umesh_roundtrip(run_in_tmpdir, request):
xml_mesh = xml_tally.filters[0].mesh
assert umesh.id == xml_mesh.id
def test_mesh_get_homogenized_materials():
"""Test the get_homogenized_materials method"""
# Simple model with 1 cm of Fe56 next to 1 cm of H1
fe = openmc.Material()
fe.add_nuclide('Fe56', 1.0)
fe.set_density('g/cm3', 5.0)
h = openmc.Material()
h.add_nuclide('H1', 1.0)
h.set_density('g/cm3', 1.0)
x0 = openmc.XPlane(-1.0, boundary_type='vacuum')
x1 = openmc.XPlane(0.0)
x2 = openmc.XPlane(1.0)
x3 = openmc.XPlane(2.0, boundary_type='vacuum')
cell1 = openmc.Cell(fill=fe, region=+x0 & -x1)
cell2 = openmc.Cell(fill=h, region=+x1 & -x2)
cell_empty = openmc.Cell(region=+x2 & -x3)
model = openmc.Model(geometry=openmc.Geometry([cell1, cell2, cell_empty]))
model.settings.particles = 1000
model.settings.batches = 10
mesh = openmc.RegularMesh()
mesh.lower_left = (-1., -1., -1.)
mesh.upper_right = (1., 1., 1.)
mesh.dimension = (3, 1, 1)
m1, m2, m3 = mesh.get_homogenized_materials(model, n_samples=1_000_000)
# Left mesh element should be only Fe56
assert m1.get_mass_density('Fe56') == pytest.approx(5.0)
# Middle mesh element should be 50% Fe56 and 50% H1
assert m2.get_mass_density('Fe56') == pytest.approx(2.5, rel=1e-2)
assert m2.get_mass_density('H1') == pytest.approx(0.5, rel=1e-2)
# Right mesh element should be only H1
assert m3.get_mass_density('H1') == pytest.approx(1.0)
mesh_void = openmc.RegularMesh()
mesh_void.lower_left = (0.5, 0.5, -1.)
mesh_void.upper_right = (1.5, 1.5, 1.)
mesh_void.dimension = (1, 1, 1)
m4, = mesh_void.get_homogenized_materials(model, n_samples=1_000_000)
# Mesh element that overlaps void should have half density
assert m4.get_mass_density('H1') == pytest.approx(0.5, rel=1e-2)