mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 14:35:27 -04:00
Support multiple meshes in R2S calculations (#3860)
This commit is contained in:
Paul Romano
GitHub
parent
806fb4ce77
commit
1116c4bdc0
5 changed files with 308 additions and 163 deletions
|
|
@ -190,6 +190,25 @@ we would run::
|
|||
|
||||
r2s.run(timesteps, source_rates, mat_vol_kwargs={'n_samples': 10_000_000})
|
||||
|
||||
It is also possible to use multiple meshes by passing a list of meshes instead
|
||||
of a single mesh. This can be useful, for example, when different regions of the
|
||||
model require different mesh resolutions. The meshes are assumed to be
|
||||
**non-overlapping**; each element--material combination across all meshes is
|
||||
treated as an independent activation region, and all meshes are handled in a
|
||||
single neutron transport solve. For example::
|
||||
|
||||
# Fine mesh near the activation target
|
||||
mesh_fine = openmc.RegularMesh()
|
||||
mesh_fine.dimension = (10, 10, 10)
|
||||
...
|
||||
|
||||
# Coarse mesh for the surrounding region
|
||||
mesh_coarse = openmc.RegularMesh()
|
||||
mesh_coarse.dimension = (5, 5, 5)
|
||||
...
|
||||
|
||||
r2s = openmc.deplete.R2SManager(model, [mesh_fine, mesh_coarse])
|
||||
|
||||
Direct 1-Step (D1S) Calculations
|
||||
================================
|
||||
|
||||
|
|
|
|||
|
|
@ -36,7 +36,8 @@ DomainTypes: TypeAlias = Union[
|
|||
Sequence[openmc.Cell],
|
||||
Sequence[openmc.Universe],
|
||||
openmc.MeshBase,
|
||||
openmc.Filter
|
||||
openmc.Filter,
|
||||
Sequence[openmc.Filter]
|
||||
]
|
||||
|
||||
|
||||
|
|
@ -69,8 +70,12 @@ def get_microxs_and_flux(
|
|||
----------
|
||||
model : openmc.Model
|
||||
OpenMC model object. Must contain geometry, materials, and settings.
|
||||
domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase, or openmc.Filter
|
||||
domains : list of openmc.Material or openmc.Cell or openmc.Universe, or openmc.MeshBase, or openmc.Filter, or list of openmc.Filter
|
||||
Domains in which to tally reaction rates, or a spatial tally filter.
|
||||
A list of filters can be provided to create one set of tallies per
|
||||
filter (e.g., one :class:`~openmc.MeshMaterialFilter` per mesh) that
|
||||
are all evaluated in a single transport solve. Results are
|
||||
concatenated across all filters in order.
|
||||
nuclides : list of str
|
||||
Nuclides to get cross sections for. If not specified, all burnable
|
||||
nuclides from the depletion chain file are used.
|
||||
|
|
@ -142,26 +147,24 @@ def get_microxs_and_flux(
|
|||
else:
|
||||
energy_filter = openmc.EnergyFilter(energies)
|
||||
|
||||
# Build list of domain filters
|
||||
if isinstance(domains, openmc.Filter):
|
||||
domain_filter = domains
|
||||
domain_filters = [domains]
|
||||
elif isinstance(domains, openmc.MeshBase):
|
||||
domain_filter = openmc.MeshFilter(domains)
|
||||
domain_filters = [openmc.MeshFilter(domains)]
|
||||
elif isinstance(domains, Sequence) and len(domains) > 0 and \
|
||||
isinstance(domains[0], openmc.Filter):
|
||||
domain_filters = list(domains)
|
||||
elif isinstance(domains[0], openmc.Material):
|
||||
domain_filter = openmc.MaterialFilter(domains)
|
||||
domain_filters = [openmc.MaterialFilter(domains)]
|
||||
elif isinstance(domains[0], openmc.Cell):
|
||||
domain_filter = openmc.CellFilter(domains)
|
||||
domain_filters = [openmc.CellFilter(domains)]
|
||||
elif isinstance(domains[0], openmc.Universe):
|
||||
domain_filter = openmc.UniverseFilter(domains)
|
||||
domain_filters = [openmc.UniverseFilter(domains)]
|
||||
else:
|
||||
raise ValueError(f"Unsupported domain type: {type(domains[0])}")
|
||||
|
||||
flux_tally = openmc.Tally(name='MicroXS flux')
|
||||
flux_tally.filters = [domain_filter, energy_filter]
|
||||
flux_tally.scores = ['flux']
|
||||
model.tallies = [flux_tally]
|
||||
|
||||
# Prepare reaction-rate tally for 'direct' or subset for 'flux' with opts
|
||||
rr_tally = None
|
||||
# Prepare reaction-rate nuclides/reactions
|
||||
rr_nuclides: list[str] = []
|
||||
rr_reactions: list[str] = []
|
||||
if reaction_rate_mode == 'direct':
|
||||
|
|
@ -177,20 +180,33 @@ def get_microxs_and_flux(
|
|||
if rr_reactions:
|
||||
rr_reactions = [r for r in rr_reactions if r in set(reactions)]
|
||||
|
||||
# Only construct tally if both lists are non-empty
|
||||
if rr_nuclides and rr_reactions:
|
||||
rr_tally = openmc.Tally(name='MicroXS RR')
|
||||
# Use 1-group energy filter for RR in flux mode
|
||||
if reaction_rate_mode == 'flux':
|
||||
rr_energy_filter = openmc.EnergyFilter(
|
||||
[energy_filter.values[0], energy_filter.values[-1]])
|
||||
else:
|
||||
rr_energy_filter = energy_filter
|
||||
rr_tally.filters = [domain_filter, rr_energy_filter]
|
||||
rr_tally.nuclides = rr_nuclides
|
||||
rr_tally.multiply_density = False
|
||||
rr_tally.scores = rr_reactions
|
||||
model.tallies.append(rr_tally)
|
||||
# Use 1-group energy filter for RR in flux mode
|
||||
has_rr = bool(rr_nuclides and rr_reactions)
|
||||
if has_rr and reaction_rate_mode == 'flux':
|
||||
rr_energy_filter = openmc.EnergyFilter(
|
||||
[energy_filter.values[0], energy_filter.values[-1]])
|
||||
else:
|
||||
rr_energy_filter = energy_filter
|
||||
|
||||
# Create one flux tally (and optionally one RR tally) per domain filter.
|
||||
flux_tallies = []
|
||||
rr_tallies = []
|
||||
model.tallies = []
|
||||
for i, domain_filter in enumerate(domain_filters):
|
||||
flux_tally = openmc.Tally(name=f'MicroXS flux {i}')
|
||||
flux_tally.filters = [domain_filter, energy_filter]
|
||||
flux_tally.scores = ['flux']
|
||||
model.tallies.append(flux_tally)
|
||||
flux_tallies.append(flux_tally)
|
||||
|
||||
if has_rr:
|
||||
rr_tally = openmc.Tally(name=f'MicroXS RR {i}')
|
||||
rr_tally.filters = [domain_filter, rr_energy_filter]
|
||||
rr_tally.nuclides = rr_nuclides
|
||||
rr_tally.multiply_density = False
|
||||
rr_tally.scores = rr_reactions
|
||||
model.tallies.append(rr_tally)
|
||||
rr_tallies.append(rr_tally)
|
||||
|
||||
if openmc.lib.is_initialized:
|
||||
openmc.lib.finalize()
|
||||
|
|
@ -227,40 +243,41 @@ def get_microxs_and_flux(
|
|||
|
||||
# Read in tally results (on all ranks)
|
||||
with StatePoint(statepoint_path) as sp:
|
||||
if rr_tally is not None:
|
||||
rr_tally = sp.tallies[rr_tally.id]
|
||||
rr_tally._read_results()
|
||||
flux_tally = sp.tallies[flux_tally.id]
|
||||
flux_tally._read_results()
|
||||
for i in range(len(flux_tallies)):
|
||||
flux_tallies[i] = sp.tallies[flux_tallies[i].id]
|
||||
flux_tallies[i]._read_results()
|
||||
if rr_tallies:
|
||||
rr_tallies[i] = sp.tallies[rr_tallies[i].id]
|
||||
rr_tallies[i]._read_results()
|
||||
|
||||
# Get flux values and make energy groups last dimension
|
||||
flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1)
|
||||
flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups)
|
||||
# Concatenate results across all domain filters
|
||||
fluxes = []
|
||||
all_flux_arrays = []
|
||||
for flux_tally in flux_tallies:
|
||||
# Get flux values and make energy groups last dimension
|
||||
flux = flux_tally.get_reshaped_data() # (domains, groups, 1, 1)
|
||||
flux = np.moveaxis(flux, 1, -1) # (domains, 1, 1, groups)
|
||||
all_flux_arrays.append(flux)
|
||||
fluxes.extend(flux.squeeze((1, 2)))
|
||||
|
||||
# Create list where each item corresponds to one domain
|
||||
fluxes = list(flux.squeeze((1, 2)))
|
||||
# If we built reaction-rate tallies, compute microscopic cross sections
|
||||
if rr_tallies:
|
||||
direct_micros = []
|
||||
for flux_arr, rr_tally in zip(all_flux_arrays, rr_tallies):
|
||||
flux = flux_arr
|
||||
# Get reaction rates and make energy groups last dimension
|
||||
reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions)
|
||||
reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups)
|
||||
|
||||
# If we built a reaction-rate tally, compute microscopic cross sections
|
||||
if rr_tally is not None:
|
||||
# Get reaction rates
|
||||
reaction_rates = rr_tally.get_reshaped_data() # (domains, groups, nuclides, reactions)
|
||||
# If RR is 1-group, sum flux over groups
|
||||
if reaction_rate_mode == "flux":
|
||||
flux = flux.sum(axis=-1, keepdims=True)
|
||||
|
||||
# Make energy groups last dimension
|
||||
reaction_rates = np.moveaxis(reaction_rates, 1, -1) # (domains, nuclides, reactions, groups)
|
||||
|
||||
# If RR is 1-group, sum flux over groups
|
||||
if reaction_rate_mode == "flux":
|
||||
flux = flux.sum(axis=-1, keepdims=True) # (domains, 1, 1, 1)
|
||||
|
||||
# Divide RR by flux to get microscopic cross sections. The indexing
|
||||
# ensures that only non-zero flux values are used, and broadcasting is
|
||||
# applied to align the shapes of reaction_rates and flux for division.
|
||||
xs = np.zeros_like(reaction_rates) # (domains, nuclides, reactions, groups)
|
||||
d, _, _, g = np.nonzero(flux)
|
||||
xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g]
|
||||
|
||||
# Create lists where each item corresponds to one domain
|
||||
direct_micros = [MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs]
|
||||
xs = np.zeros_like(reaction_rates)
|
||||
d, _, _, g = np.nonzero(flux)
|
||||
xs[d, ..., g] = reaction_rates[d, ..., g] / flux[d, :, :, g]
|
||||
direct_micros.extend(
|
||||
MicroXS(xs_i, rr_nuclides, rr_reactions) for xs_i in xs)
|
||||
|
||||
# If using flux mode, compute flux-collapsed microscopic XS
|
||||
if reaction_rate_mode == 'flux':
|
||||
|
|
@ -273,9 +290,9 @@ def get_microxs_and_flux(
|
|||
) for flux_i in fluxes]
|
||||
|
||||
# Decide which micros to use and merge if needed
|
||||
if reaction_rate_mode == 'flux' and rr_tally is not None:
|
||||
if reaction_rate_mode == 'flux' and rr_tallies:
|
||||
micros = [m1.merge(m2) for m1, m2 in zip(flux_micros, direct_micros)]
|
||||
elif rr_tally is not None:
|
||||
elif rr_tallies:
|
||||
micros = direct_micros
|
||||
else:
|
||||
micros = flux_micros
|
||||
|
|
|
|||
|
|
@ -13,11 +13,11 @@ from .results import Results
|
|||
from ..checkvalue import PathLike
|
||||
from ..mpi import comm
|
||||
from openmc.lib import TemporarySession
|
||||
from openmc.utility_funcs import change_directory
|
||||
|
||||
|
||||
def get_activation_materials(
|
||||
model: openmc.Model, mmv: openmc.MeshMaterialVolumes
|
||||
model: openmc.Model,
|
||||
mmv_list: list[openmc.MeshMaterialVolumes]
|
||||
) -> openmc.Materials:
|
||||
"""Get a list of activation materials for each mesh element/material.
|
||||
|
||||
|
|
@ -31,35 +31,35 @@ def get_activation_materials(
|
|||
----------
|
||||
model : openmc.Model
|
||||
The full model containing the geometry and materials.
|
||||
mmv : openmc.MeshMaterialVolumes
|
||||
The mesh material volumes object containing the materials and their
|
||||
volumes for each mesh element.
|
||||
mmv_list : list of openmc.MeshMaterialVolumes
|
||||
List of mesh material volumes objects, one per mesh, containing the
|
||||
materials and their volumes for each mesh element.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.Materials
|
||||
A list of materials, each corresponding to a unique mesh element and
|
||||
material combination.
|
||||
material combination across all meshes.
|
||||
|
||||
"""
|
||||
# Get the material ID, volume, and element index for each element-material
|
||||
# combination
|
||||
mat_ids = mmv._materials[mmv._materials > -1]
|
||||
volumes = mmv._volumes[mmv._materials > -1]
|
||||
elems, _ = np.where(mmv._materials > -1)
|
||||
|
||||
# Get all materials in the model
|
||||
material_dict = model._get_all_materials()
|
||||
|
||||
# Create a new activation material for each element-material combination
|
||||
# across all meshes
|
||||
materials = openmc.Materials()
|
||||
for elem, mat_id, vol in zip(elems, mat_ids, volumes):
|
||||
mat = material_dict[mat_id]
|
||||
new_mat = mat.clone()
|
||||
new_mat.depletable = True
|
||||
new_mat.name = f'Element {elem}, Material {mat_id}'
|
||||
new_mat.volume = vol
|
||||
materials.append(new_mat)
|
||||
for mesh_idx, mmv in enumerate(mmv_list):
|
||||
mat_ids = mmv._materials[mmv._materials > -1]
|
||||
volumes = mmv._volumes[mmv._materials > -1]
|
||||
elems, _ = np.where(mmv._materials > -1)
|
||||
|
||||
for elem, mat_id, vol in zip(elems, mat_ids, volumes):
|
||||
mat = material_dict[mat_id]
|
||||
new_mat = mat.clone()
|
||||
new_mat.depletable = True
|
||||
new_mat.name = f'Mesh {mesh_idx}, Element {elem}, Material {mat_id}'
|
||||
new_mat.volume = vol
|
||||
materials.append(new_mat)
|
||||
|
||||
return materials
|
||||
|
||||
|
|
@ -70,7 +70,9 @@ class R2SManager:
|
|||
This class is responsible for managing the materials and sources needed for
|
||||
mesh-based or cell-based R2S calculations. It provides methods to get
|
||||
activation materials and decay photon sources based on the mesh/cells and
|
||||
materials in the OpenMC model.
|
||||
materials in the OpenMC model. Multiple meshes can be specified as domains,
|
||||
in which case each element--material combination of each mesh is treated as
|
||||
an activation region (meshes are assumed to be non-overlapping).
|
||||
|
||||
This class supports the use of a different models for the neutron and photon
|
||||
transport calculation. However, for cell-based calculations, it assumes that
|
||||
|
|
@ -83,17 +85,20 @@ class R2SManager:
|
|||
----------
|
||||
neutron_model : openmc.Model
|
||||
The OpenMC model to use for neutron transport.
|
||||
domains : openmc.MeshBase or Sequence[openmc.Cell]
|
||||
The mesh or a sequence of cells that represent the spatial units over
|
||||
which the R2S calculation will be performed.
|
||||
domains : openmc.MeshBase or Sequence[openmc.MeshBase] or Sequence[openmc.Cell]
|
||||
The mesh(es) or a sequence of cells that represent the spatial units
|
||||
over which the R2S calculation will be performed. When a single
|
||||
:class:`~openmc.MeshBase` or a sequence of meshes is given, each
|
||||
element--material combination across all meshes is treated as an
|
||||
activation region.
|
||||
photon_model : openmc.Model, optional
|
||||
The OpenMC model to use for photon transport calculations. If None, a
|
||||
shallow copy of the neutron_model will be created and used.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
domains : openmc.MeshBase or Sequence[openmc.Cell]
|
||||
The mesh or a sequence of cells that represent the spatial units over
|
||||
domains : list of openmc.MeshBase or Sequence[openmc.Cell]
|
||||
The meshes or a sequence of cells that represent the spatial units over
|
||||
which the R2S calculation will be performed.
|
||||
neutron_model : openmc.Model
|
||||
The OpenMC model used for neutron transport.
|
||||
|
|
@ -101,7 +106,7 @@ class R2SManager:
|
|||
The OpenMC model used for photon transport calculations.
|
||||
method : {'mesh-based', 'cell-based'}
|
||||
Indicates whether the R2S calculation uses mesh elements ('mesh-based')
|
||||
as the spatial discetization or a list of a cells ('cell-based').
|
||||
as the spatial discretization or a list of cells ('cell-based').
|
||||
results : dict
|
||||
A dictionary that stores results from the R2S calculation.
|
||||
|
||||
|
|
@ -109,7 +114,7 @@ class R2SManager:
|
|||
def __init__(
|
||||
self,
|
||||
neutron_model: openmc.Model,
|
||||
domains: openmc.MeshBase | Sequence[openmc.Cell],
|
||||
domains: openmc.MeshBase | Sequence[openmc.MeshBase] | Sequence[openmc.Cell],
|
||||
photon_model: openmc.Model | None = None,
|
||||
):
|
||||
self.neutron_model = neutron_model
|
||||
|
|
@ -126,9 +131,14 @@ class R2SManager:
|
|||
self.photon_model = photon_model
|
||||
if isinstance(domains, openmc.MeshBase):
|
||||
self.method = 'mesh-based'
|
||||
self.domains = [domains]
|
||||
elif isinstance(domains, Sequence) and len(domains) > 0 and \
|
||||
isinstance(domains[0], openmc.MeshBase):
|
||||
self.method = 'mesh-based'
|
||||
self.domains = list(domains)
|
||||
else:
|
||||
self.method = 'cell-based'
|
||||
self.domains = domains
|
||||
self.domains = list(domains)
|
||||
self.results = {}
|
||||
|
||||
def run(
|
||||
|
|
@ -243,11 +253,13 @@ class R2SManager:
|
|||
):
|
||||
"""Run the neutron transport step.
|
||||
|
||||
This step computes the material volume fractions on the mesh, creates a
|
||||
mesh-material filter, and retrieves the fluxes and microscopic cross
|
||||
sections for each mesh/material combination. This step will populate the
|
||||
'fluxes' and 'micros' keys in the results dictionary. For a mesh-based
|
||||
calculation, it will also populate the 'mesh_material_volumes' key.
|
||||
This step computes the material volume fractions on each mesh, creates
|
||||
mesh-material filters, and retrieves the fluxes and microscopic cross
|
||||
sections for each mesh/material combination via a single transport
|
||||
solve. This step will populate the 'fluxes' and 'micros' keys in the
|
||||
results dictionary. For a mesh-based calculation, it will also populate
|
||||
the 'mesh_material_volumes' key (a list of
|
||||
:class:`~openmc.MeshMaterialVolumes`, one per mesh).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
|
@ -266,19 +278,28 @@ class R2SManager:
|
|||
output_dir.mkdir(parents=True, exist_ok=True)
|
||||
|
||||
if self.method == 'mesh-based':
|
||||
# Compute material volume fractions on the mesh
|
||||
# Compute material volume fractions on each mesh
|
||||
if mat_vol_kwargs is None:
|
||||
mat_vol_kwargs = {}
|
||||
mat_vol_kwargs.setdefault('bounding_boxes', True)
|
||||
self.results['mesh_material_volumes'] = mmv = comm.bcast(
|
||||
self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs))
|
||||
|
||||
# Save results to file
|
||||
if comm.rank == 0:
|
||||
mmv.save(output_dir / 'mesh_material_volumes.npz')
|
||||
mmv_list = []
|
||||
domain_filters = []
|
||||
for i, mesh in enumerate(self.domains):
|
||||
mmv = comm.bcast(
|
||||
mesh.material_volumes(self.neutron_model, **mat_vol_kwargs))
|
||||
mmv_list.append(mmv)
|
||||
|
||||
# Create mesh-material filter based on what combos were found
|
||||
domains = openmc.MeshMaterialFilter.from_volumes(self.domains, mmv)
|
||||
# Save results to file
|
||||
if comm.rank == 0:
|
||||
mmv.save(output_dir / f'mesh_material_volumes_{i}.npz')
|
||||
|
||||
# Create mesh-material filter for this mesh
|
||||
domain_filters.append(
|
||||
openmc.MeshMaterialFilter.from_volumes(mesh, mmv))
|
||||
|
||||
self.results['mesh_material_volumes'] = mmv_list
|
||||
domains = domain_filters
|
||||
else:
|
||||
domains: Sequence[openmc.Cell] = self.domains
|
||||
|
||||
|
|
@ -357,8 +378,9 @@ class R2SManager:
|
|||
|
||||
if self.method == 'mesh-based':
|
||||
# Get unique material for each (mesh, material) combination
|
||||
mmv = self.results['mesh_material_volumes']
|
||||
self.results['activation_materials'] = get_activation_materials(self.neutron_model, mmv)
|
||||
mmv_list = self.results['mesh_material_volumes']
|
||||
self.results['activation_materials'] = get_activation_materials(
|
||||
self.neutron_model, mmv_list)
|
||||
else:
|
||||
# Create unique material for each cell
|
||||
activation_mats = openmc.Materials()
|
||||
|
|
@ -468,12 +490,20 @@ class R2SManager:
|
|||
# photon model if it is different from the neutron model to account for
|
||||
# potential material changes
|
||||
if self.method == 'mesh-based' and different_photon_model:
|
||||
self.results['mesh_material_volumes_photon'] = photon_mmv = comm.bcast(
|
||||
self.domains.material_volumes(self.photon_model, **mat_vol_kwargs))
|
||||
if mat_vol_kwargs is None:
|
||||
mat_vol_kwargs = {}
|
||||
photon_mmv_list = []
|
||||
for i, mesh in enumerate(self.domains):
|
||||
photon_mmv = comm.bcast(
|
||||
mesh.material_volumes(self.photon_model, **mat_vol_kwargs))
|
||||
photon_mmv_list.append(photon_mmv)
|
||||
|
||||
# Save photon MMV results to file
|
||||
if comm.rank == 0:
|
||||
photon_mmv.save(output_dir / 'mesh_material_volumes.npz')
|
||||
# Save photon MMV results to file
|
||||
if comm.rank == 0:
|
||||
photon_mmv.save(
|
||||
output_dir / f'mesh_material_volumes_{i}.npz')
|
||||
|
||||
self.results['mesh_material_volumes_photon'] = photon_mmv_list
|
||||
|
||||
if comm.rank == 0:
|
||||
tally_ids = [tally.id for tally in self.photon_model.tallies]
|
||||
|
|
@ -543,7 +573,7 @@ class R2SManager:
|
|||
) -> list[openmc.IndependentSource]:
|
||||
"""Create decay photon source for a mesh-based calculation.
|
||||
|
||||
For each mesh element-material combination, an
|
||||
For each mesh element-material combination across all meshes, an
|
||||
:class:`~openmc.IndependentSource` is created with a
|
||||
:class:`~openmc.stats.Box` spatial distribution based on the bounding
|
||||
box of the material within the mesh element. A material constraint is
|
||||
|
|
@ -575,52 +605,56 @@ class R2SManager:
|
|||
index_mat = 0
|
||||
|
||||
# Get various results from previous steps
|
||||
mat_vols = self.results['mesh_material_volumes']
|
||||
mmv_list = self.results['mesh_material_volumes']
|
||||
materials = self.results['activation_materials']
|
||||
results = self.results['depletion_results']
|
||||
photon_mat_vols = self.results.get('mesh_material_volumes_photon')
|
||||
photon_mmv_list = self.results.get('mesh_material_volumes_photon')
|
||||
|
||||
# Total number of mesh elements
|
||||
n_elements = mat_vols.num_elements
|
||||
for mesh_idx, mat_vols in enumerate(mmv_list):
|
||||
photon_mat_vols = photon_mmv_list[mesh_idx] \
|
||||
if photon_mmv_list is not None else None
|
||||
|
||||
for index_elem in range(n_elements):
|
||||
# Determine which materials exist in the photon model for this element
|
||||
if photon_mat_vols is not None:
|
||||
photon_materials = {
|
||||
mat_id
|
||||
for mat_id, _ in photon_mat_vols.by_element(index_elem)
|
||||
if mat_id is not None
|
||||
}
|
||||
# Total number of mesh elements for this mesh
|
||||
n_elements = mat_vols.num_elements
|
||||
|
||||
for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
|
||||
# Skip void volume
|
||||
if mat_id is None:
|
||||
continue
|
||||
for index_elem in range(n_elements):
|
||||
# Determine which materials exist in the photon model for this element
|
||||
if photon_mat_vols is not None:
|
||||
photon_materials = {
|
||||
mat_id
|
||||
for mat_id, _ in photon_mat_vols.by_element(index_elem)
|
||||
if mat_id is not None
|
||||
}
|
||||
|
||||
# Skip if this material doesn't exist in photon model
|
||||
if photon_mat_vols is not None and mat_id not in photon_materials:
|
||||
for mat_id, _, bbox in mat_vols.by_element(index_elem, include_bboxes=True):
|
||||
# Skip void volume
|
||||
if mat_id is None:
|
||||
continue
|
||||
|
||||
# Skip if this material doesn't exist in photon model
|
||||
if photon_mat_vols is not None and mat_id not in photon_materials:
|
||||
index_mat += 1
|
||||
continue
|
||||
|
||||
# Get activated material composition
|
||||
original_mat = materials[index_mat]
|
||||
activated_mat = results[time_index].get_material(str(original_mat.id))
|
||||
|
||||
# Create decay photon source
|
||||
energy = activated_mat.get_decay_photon_energy()
|
||||
if energy is not None:
|
||||
strength = energy.integral()
|
||||
space = openmc.stats.Box(*bbox)
|
||||
sources.append(openmc.IndependentSource(
|
||||
space=space,
|
||||
energy=energy,
|
||||
particle='photon',
|
||||
strength=strength,
|
||||
constraints={'domains': [mat_dict[mat_id]]}
|
||||
))
|
||||
|
||||
# Increment index of activated material
|
||||
index_mat += 1
|
||||
continue
|
||||
|
||||
# Get activated material composition
|
||||
original_mat = materials[index_mat]
|
||||
activated_mat = results[time_index].get_material(str(original_mat.id))
|
||||
|
||||
# Create decay photon source
|
||||
energy = activated_mat.get_decay_photon_energy()
|
||||
if energy is not None:
|
||||
strength = energy.integral()
|
||||
space = openmc.stats.Box(*bbox)
|
||||
sources.append(openmc.IndependentSource(
|
||||
space=space,
|
||||
energy=energy,
|
||||
particle='photon',
|
||||
strength=strength,
|
||||
constraints={'domains': [mat_dict[mat_id]]}
|
||||
))
|
||||
|
||||
# Increment index of activated material
|
||||
index_mat += 1
|
||||
|
||||
return sources
|
||||
|
||||
|
|
@ -638,10 +672,13 @@ class R2SManager:
|
|||
# Load neutron transport results
|
||||
neutron_dir = path / 'neutron_transport'
|
||||
if self.method == 'mesh-based':
|
||||
mmv_file = neutron_dir / 'mesh_material_volumes.npz'
|
||||
if mmv_file.exists():
|
||||
self.results['mesh_material_volumes'] = \
|
||||
openmc.MeshMaterialVolumes.from_npz(mmv_file)
|
||||
mmv_files = sorted(neutron_dir.glob('mesh_material_volumes*.npz'),
|
||||
key=lambda p: int(p.stem.split('_')[-1])
|
||||
if p.stem[-1].isdigit() else 0)
|
||||
if mmv_files:
|
||||
self.results['mesh_material_volumes'] = [
|
||||
openmc.MeshMaterialVolumes.from_npz(f) for f in mmv_files
|
||||
]
|
||||
fluxes_file = neutron_dir / 'fluxes.npy'
|
||||
if fluxes_file.exists():
|
||||
self.results['fluxes'] = list(np.load(fluxes_file, allow_pickle=True))
|
||||
|
|
@ -665,10 +702,15 @@ class R2SManager:
|
|||
|
||||
# Load photon mesh material volumes if they exist (for mesh-based calculations)
|
||||
if self.method == 'mesh-based':
|
||||
photon_mmv_file = photon_dir / 'mesh_material_volumes.npz'
|
||||
if photon_mmv_file.exists():
|
||||
self.results['mesh_material_volumes_photon'] = \
|
||||
openmc.MeshMaterialVolumes.from_npz(photon_mmv_file)
|
||||
photon_mmv_files = sorted(
|
||||
photon_dir.glob('mesh_material_volumes*.npz'),
|
||||
key=lambda p: int(p.stem.split('_')[-1])
|
||||
if p.stem[-1].isdigit() else 0)
|
||||
if photon_mmv_files:
|
||||
self.results['mesh_material_volumes_photon'] = [
|
||||
openmc.MeshMaterialVolumes.from_npz(f)
|
||||
for f in photon_mmv_files
|
||||
]
|
||||
|
||||
# Load tally IDs from JSON file
|
||||
tally_ids_path = photon_dir / 'tally_ids.json'
|
||||
|
|
|
|||
|
|
@ -165,11 +165,11 @@ def test_hybrid_tally_setup():
|
|||
|
||||
# Check that both tallies were created with the expected properties
|
||||
tally_names = [t.name for t in captured['tallies']]
|
||||
assert 'MicroXS flux' in tally_names
|
||||
assert 'MicroXS RR' in tally_names
|
||||
assert 'MicroXS flux 0' in tally_names
|
||||
assert 'MicroXS RR 0' in tally_names
|
||||
|
||||
# Check that the RR tally has the expected nuclides and reactions
|
||||
rr = next(t for t in captured['tallies'] if t.name == 'MicroXS RR')
|
||||
rr = next(t for t in captured['tallies'] if t.name == 'MicroXS RR 0')
|
||||
assert rr.nuclides == ['U235']
|
||||
assert rr.scores == ['fission']
|
||||
|
||||
|
|
|
|||
|
|
@ -6,7 +6,7 @@ from openmc.deplete import Chain, R2SManager
|
|||
|
||||
|
||||
@pytest.fixture
|
||||
def simple_model_and_mesh(tmp_path):
|
||||
def simple_model_and_mesh():
|
||||
# Define two materials: water and Ni
|
||||
h2o = openmc.Material()
|
||||
h2o.add_nuclide("H1", 2.0)
|
||||
|
|
@ -68,7 +68,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
|
|||
nt = Path(outdir) / 'neutron_transport'
|
||||
assert (nt / 'fluxes.npy').exists()
|
||||
assert (nt / 'micros.h5').exists()
|
||||
assert (nt / 'mesh_material_volumes.npz').exists()
|
||||
assert (nt / 'mesh_material_volumes_0.npz').exists()
|
||||
act = Path(outdir) / 'activation'
|
||||
assert (act / 'depletion_results.h5').exists()
|
||||
pt = Path(outdir) / 'photon_transport'
|
||||
|
|
@ -78,7 +78,8 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
|
|||
# Basic results structure checks
|
||||
assert len(r2s.results['fluxes']) == 2
|
||||
assert len(r2s.results['micros']) == 2
|
||||
assert len(r2s.results['mesh_material_volumes']) == 2
|
||||
assert len(r2s.results['mesh_material_volumes']) == 1
|
||||
assert len(r2s.results['mesh_material_volumes'][0]) == 2
|
||||
assert len(r2s.results['activation_materials']) == 2
|
||||
assert len(r2s.results['depletion_results']) == 2
|
||||
|
||||
|
|
@ -93,11 +94,77 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path):
|
|||
r2s_loaded.load_results(outdir)
|
||||
assert len(r2s_loaded.results['fluxes']) == 2
|
||||
assert len(r2s_loaded.results['micros']) == 2
|
||||
assert len(r2s_loaded.results['mesh_material_volumes']) == 2
|
||||
assert len(r2s_loaded.results['mesh_material_volumes']) == 1
|
||||
assert len(r2s_loaded.results['mesh_material_volumes'][0]) == 2
|
||||
assert len(r2s_loaded.results['activation_materials']) == 2
|
||||
assert len(r2s_loaded.results['depletion_results']) == 2
|
||||
|
||||
|
||||
def test_r2s_multi_mesh(simple_model_and_mesh, tmp_path):
|
||||
model, _, _ = simple_model_and_mesh
|
||||
|
||||
# Two 1x1x1 meshes that together cover the full domain, split along y.
|
||||
# Each mesh element spans the full x range [-10, 10], crossing the x=0
|
||||
# material boundary, so both meshes contain both materials within their
|
||||
# single element.
|
||||
mesh1 = openmc.RegularMesh()
|
||||
mesh1.lower_left = (-10.0, -10.0, -10.0)
|
||||
mesh1.upper_right = (10.0, 0.0, 10.0)
|
||||
mesh1.dimension = (1, 1, 1)
|
||||
mesh2 = openmc.RegularMesh()
|
||||
mesh2.lower_left = (-10.0, 0.0, -10.0)
|
||||
mesh2.upper_right = (10.0, 10.0, 10.0)
|
||||
mesh2.dimension = (1, 1, 1)
|
||||
|
||||
r2s = R2SManager(model, [mesh1, mesh2])
|
||||
chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml")
|
||||
|
||||
outdir = r2s.run(
|
||||
timesteps=[(1.0, 'd')],
|
||||
source_rates=[1.0],
|
||||
photon_time_indices=[1],
|
||||
output_dir=tmp_path,
|
||||
chain_file=chain,
|
||||
)
|
||||
|
||||
# Check that per-mesh MMV files were written
|
||||
nt = Path(outdir) / 'neutron_transport'
|
||||
assert (nt / 'fluxes.npy').exists()
|
||||
assert (nt / 'micros.h5').exists()
|
||||
assert (nt / 'mesh_material_volumes_0.npz').exists()
|
||||
assert (nt / 'mesh_material_volumes_1.npz').exists()
|
||||
act = Path(outdir) / 'activation'
|
||||
assert (act / 'depletion_results.h5').exists()
|
||||
pt = Path(outdir) / 'photon_transport'
|
||||
assert (pt / 'tally_ids.json').exists()
|
||||
assert (pt / 'time_1' / 'statepoint.10.h5').exists()
|
||||
|
||||
# Two meshes, each with 1 element containing both materials →
|
||||
# 2 element-material combinations per mesh, 4 total
|
||||
assert len(r2s.results['mesh_material_volumes']) == 2
|
||||
assert len(r2s.results['mesh_material_volumes'][0]) == 2
|
||||
assert len(r2s.results['mesh_material_volumes'][1]) == 2
|
||||
assert len(r2s.results['fluxes']) == 4
|
||||
assert len(r2s.results['micros']) == 4
|
||||
assert len(r2s.results['activation_materials']) == 4
|
||||
assert len(r2s.results['depletion_results']) == 2
|
||||
|
||||
# Activation material names encode mesh index
|
||||
amats = r2s.results['activation_materials']
|
||||
assert all(m.depletable for m in amats)
|
||||
assert any('Mesh 0' in m.name for m in amats)
|
||||
assert any('Mesh 1' in m.name for m in amats)
|
||||
|
||||
# Check loading results
|
||||
r2s_loaded = R2SManager(model, [mesh1, mesh2])
|
||||
r2s_loaded.load_results(outdir)
|
||||
assert len(r2s_loaded.results['mesh_material_volumes']) == 2
|
||||
assert len(r2s_loaded.results['mesh_material_volumes'][0]) == 2
|
||||
assert len(r2s_loaded.results['mesh_material_volumes'][1]) == 2
|
||||
assert len(r2s_loaded.results['activation_materials']) == 4
|
||||
assert len(r2s_loaded.results['depletion_results']) == 2
|
||||
|
||||
|
||||
def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path):
|
||||
model, (c1, c2), _ = simple_model_and_mesh
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue