mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Merge branch 'develop' into centre_for_cylinder_spherical_meshes
This commit is contained in:
Rémi Delaporte-Mathurin
GitHub
commit
10f2a945ca
6 changed files with 198 additions and 31 deletions
|
|
@ -156,9 +156,9 @@ endif()
|
|||
find_package(HDF5 REQUIRED COMPONENTS C HL)
|
||||
if(HDF5_IS_PARALLEL)
|
||||
if(NOT OPENMC_USE_MPI)
|
||||
message(FATAL_ERROR "Parallel HDF5 was detected, but the detected compiler,\
|
||||
${CMAKE_CXX_COMPILER}, does not support MPI. An MPI-capable compiler must \
|
||||
be used with parallel HDF5.")
|
||||
message(FATAL_ERROR "Parallel HDF5 was detected, but MPI was not enabled.\
|
||||
To use parallel HDF5, OpenMC needs to be built with MPI support by passing\
|
||||
-DOPENMC_USE_MPI=ON when calling cmake.")
|
||||
endif()
|
||||
message(STATUS "Using parallel HDF5")
|
||||
endif()
|
||||
|
|
|
|||
|
|
@ -1029,7 +1029,8 @@ class Material(IDManagerMixin):
|
|||
|
||||
return nuclides
|
||||
|
||||
def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False):
|
||||
def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False,
|
||||
volume: Optional[float] = None):
|
||||
"""Returns the activity of the material or for each nuclide in the
|
||||
material in units of [Bq], [Bq/g] or [Bq/cm3].
|
||||
|
||||
|
|
@ -1044,6 +1045,11 @@ class Material(IDManagerMixin):
|
|||
by_nuclide : bool
|
||||
Specifies if the activity should be returned for the material as a
|
||||
whole or per nuclide. Default is False.
|
||||
volume : float, optional
|
||||
Volume of the material. If not passed, defaults to using the
|
||||
:attr:`Material.volume` attribute.
|
||||
|
||||
.. versionadded:: 0.13.3
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -1057,7 +1063,7 @@ class Material(IDManagerMixin):
|
|||
cv.check_type('by_nuclide', by_nuclide, bool)
|
||||
|
||||
if units == 'Bq':
|
||||
multiplier = self.volume
|
||||
multiplier = volume if volume is not None else self.volume
|
||||
elif units == 'Bq/cm3':
|
||||
multiplier = 1
|
||||
elif units == 'Bq/g':
|
||||
|
|
@ -1070,7 +1076,8 @@ class Material(IDManagerMixin):
|
|||
|
||||
return activity if by_nuclide else sum(activity.values())
|
||||
|
||||
def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False):
|
||||
def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False,
|
||||
volume: Optional[float] = None):
|
||||
"""Returns the decay heat of the material or for each nuclide in the
|
||||
material in units of [W], [W/g] or [W/cm3].
|
||||
|
||||
|
|
@ -1085,6 +1092,11 @@ class Material(IDManagerMixin):
|
|||
by_nuclide : bool
|
||||
Specifies if the decay heat should be returned for the material as a
|
||||
whole or per nuclide. Default is False.
|
||||
volume : float, optional
|
||||
Volume of the material. If not passed, defaults to using the
|
||||
:attr:`Material.volume` attribute.
|
||||
|
||||
.. versionadded:: 0.13.3
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -1098,7 +1110,7 @@ class Material(IDManagerMixin):
|
|||
cv.check_type('by_nuclide', by_nuclide, bool)
|
||||
|
||||
if units == 'W':
|
||||
multiplier = self.volume
|
||||
multiplier = volume if volume is not None else self.volume
|
||||
elif units == 'W/cm3':
|
||||
multiplier = 1
|
||||
elif units == 'W/g':
|
||||
|
|
@ -1113,11 +1125,19 @@ class Material(IDManagerMixin):
|
|||
|
||||
return decayheat if by_nuclide else sum(decayheat.values())
|
||||
|
||||
def get_nuclide_atoms(self):
|
||||
def get_nuclide_atoms(self, volume: Optional[float] = None):
|
||||
"""Return number of atoms of each nuclide in the material
|
||||
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
volume : float, optional
|
||||
Volume of the material. If not passed, defaults to using the
|
||||
:attr:`Material.volume` attribute.
|
||||
|
||||
.. versionadded:: 0.13.3
|
||||
|
||||
Returns
|
||||
-------
|
||||
dict
|
||||
|
|
@ -1125,11 +1145,13 @@ class Material(IDManagerMixin):
|
|||
atoms present in the material.
|
||||
|
||||
"""
|
||||
if self.volume is None:
|
||||
if volume is None:
|
||||
volume = self.volume
|
||||
if volume is None:
|
||||
raise ValueError("Volume must be set in order to determine atoms.")
|
||||
atoms = {}
|
||||
for nuclide, atom_per_bcm in self.get_nuclide_atom_densities().items():
|
||||
atoms[nuclide] = 1.0e24 * atom_per_bcm * self.volume
|
||||
atoms[nuclide] = 1.0e24 * atom_per_bcm * volume
|
||||
return atoms
|
||||
|
||||
def get_mass_density(self, nuclide: Optional[str] = None):
|
||||
|
|
@ -1154,7 +1176,7 @@ class Material(IDManagerMixin):
|
|||
mass_density += density_i
|
||||
return mass_density
|
||||
|
||||
def get_mass(self, nuclide: Optional[str] = None):
|
||||
def get_mass(self, nuclide: Optional[str] = None, volume: Optional[float] = None):
|
||||
"""Return mass of one or all nuclides.
|
||||
|
||||
Note that this method requires that the :attr:`Material.volume` has
|
||||
|
|
@ -1165,6 +1187,12 @@ class Material(IDManagerMixin):
|
|||
nuclides : str, optional
|
||||
Nuclide for which mass is desired. If not specified, the density
|
||||
for the entire material is given.
|
||||
volume : float, optional
|
||||
Volume of the material. If not passed, defaults to using the
|
||||
:attr:`Material.volume` attribute.
|
||||
|
||||
.. versionadded:: 0.13.3
|
||||
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -1172,9 +1200,11 @@ class Material(IDManagerMixin):
|
|||
Mass of the nuclide/material in [g]
|
||||
|
||||
"""
|
||||
if self.volume is None:
|
||||
if volume is None:
|
||||
volume = self.volume
|
||||
if volume is None:
|
||||
raise ValueError("Volume must be set in order to determine mass.")
|
||||
return self.volume*self.get_mass_density(nuclide)
|
||||
return volume*self.get_mass_density(nuclide)
|
||||
|
||||
def clone(self, memo: Optional[dict] = None):
|
||||
"""Create a copy of this material with a new unique ID.
|
||||
|
|
|
|||
|
|
@ -446,7 +446,7 @@ class RegularMesh(StructuredMesh):
|
|||
if self._width is not None:
|
||||
self._width = None
|
||||
warnings.warn("Unsetting width attribute.")
|
||||
|
||||
|
||||
if self.lower_left is not None and any(np.isclose(self.lower_left, upper_right)):
|
||||
raise ValueError("Mesh cannot have zero thickness in any dimension")
|
||||
|
||||
|
|
@ -1578,6 +1578,7 @@ class SphericalMesh(StructuredMesh):
|
|||
pts_cartesian = np.copy(pts_spherical)
|
||||
|
||||
r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2]
|
||||
|
||||
pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) + self.origin[0]
|
||||
pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) + self.origin[1]
|
||||
pts_cartesian[:, 2] = r * np.cos(phi) + self.origin[2]
|
||||
|
|
|
|||
|
|
@ -244,7 +244,7 @@ class Model:
|
|||
def from_model_xml(cls, path='model.xml'):
|
||||
"""Create model from single XML file
|
||||
|
||||
.. vesionadded:: 0.13.3
|
||||
.. versionadded:: 0.13.3
|
||||
|
||||
Parameters
|
||||
----------
|
||||
|
|
@ -740,7 +740,7 @@ class Model:
|
|||
|
||||
if len(self.settings.volume_calculations) == 0:
|
||||
# Then there is no volume calculation specified
|
||||
raise ValueError("The Settings.volume_calculation attribute must"
|
||||
raise ValueError("The Settings.volume_calculations attribute must"
|
||||
" be specified before executing this method!")
|
||||
|
||||
with _change_directory(Path(cwd)):
|
||||
|
|
|
|||
|
|
@ -377,6 +377,9 @@ def test_get_nuclide_atoms():
|
|||
atoms = mat.get_nuclide_atoms()
|
||||
assert atoms['Li6'] == pytest.approx(mat.density * mat.volume)
|
||||
|
||||
atoms = mat.get_nuclide_atoms(volume=10.0)
|
||||
assert atoms['Li6'] == pytest.approx(mat.density * 10.0)
|
||||
|
||||
|
||||
def test_mass():
|
||||
m = openmc.Material()
|
||||
|
|
@ -394,6 +397,9 @@ def test_mass():
|
|||
assert m.get_mass() == pytest.approx(20.0)
|
||||
assert m.fissionable_mass == pytest.approx(10.0)
|
||||
|
||||
# Test with volume specified as argument
|
||||
assert m.get_mass('Zr90', volume=1.0) == pytest.approx(1.0)
|
||||
|
||||
|
||||
def test_materials(run_in_tmpdir):
|
||||
m1 = openmc.Material()
|
||||
|
|
@ -544,11 +550,14 @@ def test_get_activity():
|
|||
m4.volume = 10.
|
||||
assert pytest.approx(m4.get_activity(units='Bq')) == 355978108155965.94*3/2*10 # [Bq]
|
||||
|
||||
# Test with volume specified as argument
|
||||
assert pytest.approx(m4.get_activity(units='Bq', volume=1.0)) == 355978108155965.94*3/2
|
||||
|
||||
|
||||
def test_get_decay_heat():
|
||||
# Set chain file for testing
|
||||
openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml'
|
||||
|
||||
|
||||
"""Tests the decay heat of stable, metastable and active materials"""
|
||||
m1 = openmc.Material()
|
||||
m1.add_nuclide("U235", 0.2)
|
||||
|
|
@ -589,7 +598,10 @@ def test_get_decay_heat():
|
|||
# volume is required to calculate total decay heat
|
||||
m4.volume = 10.
|
||||
assert pytest.approx(m4.get_decay_heat(units='W')) == 40175.15720273193*3/2*10 # [W]
|
||||
|
||||
|
||||
# Test with volume specified as argument
|
||||
assert pytest.approx(m4.get_decay_heat(units='W', volume=1.0)) == 40175.15720273193*3/2
|
||||
|
||||
|
||||
def test_decay_photon_energy():
|
||||
# Set chain file for testing
|
||||
|
|
|
|||
|
|
@ -1,3 +1,5 @@
|
|||
from itertools import product
|
||||
|
||||
import numpy as np
|
||||
from pathlib import Path
|
||||
import pytest
|
||||
|
|
@ -7,29 +9,77 @@ from vtk.util import numpy_support as nps
|
|||
|
||||
import openmc
|
||||
|
||||
@pytest.fixture
|
||||
def model():
|
||||
openmc.reset_auto_ids()
|
||||
|
||||
surf1 = openmc.Sphere(r=10, boundary_type='vacuum')
|
||||
surf2 = openmc.XPlane(x0=-0.001, boundary_type='vacuum')
|
||||
|
||||
cell = openmc.Cell(region=-surf1 & -surf2)
|
||||
|
||||
geometry = openmc.Geometry([cell])
|
||||
|
||||
settings = openmc.Settings()
|
||||
settings.batches = 2
|
||||
settings.particles = 100
|
||||
settings.run_mode = 'fixed source'
|
||||
|
||||
source = openmc.Source()
|
||||
source.angle = openmc.stats.Isotropic()
|
||||
source.energy = openmc.stats.Discrete([1.0e6], [1.0])
|
||||
source.space = openmc.stats.Point((-0.01, -0.01, -0.01))
|
||||
|
||||
settings.source = source
|
||||
|
||||
model = openmc.Model(geometry=geometry, settings=settings)
|
||||
|
||||
return model
|
||||
|
||||
|
||||
regular_mesh = openmc.RegularMesh()
|
||||
regular_mesh.lower_left = (0, 0, 0)
|
||||
regular_mesh.upper_right = (1, 1, 1)
|
||||
regular_mesh.lower_left = (-10, -10, -10)
|
||||
regular_mesh.upper_right = (10, 10, 10)
|
||||
regular_mesh.dimension = [30, 20, 10]
|
||||
|
||||
rectilinear_mesh = openmc.RectilinearMesh()
|
||||
rectilinear_mesh.x_grid = np.linspace(1, 2, num=30)
|
||||
rectilinear_mesh.y_grid = np.linspace(1, 2, num=30)
|
||||
rectilinear_mesh.z_grid = np.linspace(1, 2, num=30)
|
||||
rectilinear_mesh.x_grid = np.linspace(-10, 10, 6)
|
||||
rectilinear_mesh.y_grid = np.logspace(0, 1, 7)
|
||||
rectilinear_mesh.y_grid = \
|
||||
np.concatenate((-rectilinear_mesh.y_grid[::-1], rectilinear_mesh.y_grid))
|
||||
rectilinear_mesh.z_grid = np.linspace(-10, 10, 11)
|
||||
|
||||
cylinder_mesh = openmc.CylindricalMesh()
|
||||
cylinder_mesh.r_grid = np.linspace(1, 2, num=30)
|
||||
cylinder_mesh.phi_grid = np.linspace(0, np.pi, num=50)
|
||||
cylinder_mesh.z_grid = np.linspace(0, 1, num=30)
|
||||
cylinder_mesh.r_grid = np.linspace(0, 10, 23)
|
||||
cylinder_mesh.phi_grid = np.linspace(0, np.pi, 21)
|
||||
cylinder_mesh.z_grid = np.linspace(0, 1, 15)
|
||||
|
||||
spherical_mesh = openmc.SphericalMesh()
|
||||
spherical_mesh.r_grid = np.linspace(1, 2, num=30)
|
||||
spherical_mesh.phi_grid = np.linspace(0, np.pi, num=50)
|
||||
spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30)
|
||||
spherical_mesh.r_grid = np.linspace(1, 10, 30)
|
||||
spherical_mesh.phi_grid = np.linspace(0, 0.8*np.pi, 25)
|
||||
spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, 15)
|
||||
|
||||
MESHES = [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]
|
||||
|
||||
x_plane = openmc.XPlane(x0=-0.001, boundary_type='vacuum')
|
||||
y_plane = openmc.YPlane(y0=-0.001, boundary_type='vacuum')
|
||||
z_plane = openmc.ZPlane(z0=-0.001, boundary_type='vacuum')
|
||||
|
||||
SURFS = [x_plane, y_plane, z_plane]
|
||||
|
||||
|
||||
@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh])
|
||||
def ids(mesh):
|
||||
if isinstance(mesh, openmc.CylindricalMesh):
|
||||
return 'cylindrical_mesh'
|
||||
elif isinstance(mesh, openmc.RegularMesh):
|
||||
return 'regular_mesh'
|
||||
elif isinstance(mesh, openmc.RectilinearMesh):
|
||||
return 'rectilinear_mesh'
|
||||
elif isinstance(mesh, openmc.SphericalMesh):
|
||||
return 'spherical_mesh'
|
||||
|
||||
|
||||
@pytest.mark.parametrize("mesh", MESHES, ids=ids)
|
||||
def test_write_data_to_vtk(mesh, tmpdir):
|
||||
# BUILD
|
||||
filename = Path(tmpdir) / "out.vtk"
|
||||
|
|
@ -69,7 +119,7 @@ def test_write_data_to_vtk(mesh, tmpdir):
|
|||
assert all(data1 == 1.0)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh])
|
||||
@pytest.mark.parametrize("mesh", MESHES, ids=ids)
|
||||
def test_write_data_to_vtk_size_mismatch(mesh):
|
||||
"""Checks that an error is raised when the size of the dataset
|
||||
doesn't match the mesh number of cells
|
||||
|
|
@ -149,3 +199,77 @@ def test_write_data_to_vtk_round_trip(run_in_tmpdir):
|
|||
|
||||
# checks that the vtk cell values are equal to the data
|
||||
assert np.array_equal(vtk_values, data)
|
||||
|
||||
def mesh_surf_id(param):
|
||||
if isinstance(param, openmc.MeshBase):
|
||||
return ids(param)
|
||||
elif isinstance(param, openmc.XPlane):
|
||||
return 'XPlane'
|
||||
elif isinstance(param, openmc.YPlane):
|
||||
return 'YPlane'
|
||||
elif isinstance(param, openmc.ZPlane):
|
||||
return 'ZPlane'
|
||||
|
||||
|
||||
@pytest.mark.parametrize("mesh,surface", product(MESHES, SURFS), ids=mesh_surf_id)
|
||||
def test_vtk_write_ordering(model, mesh, surface):
|
||||
|
||||
tally = openmc.Tally()
|
||||
tally.scores = ['flux']
|
||||
# use the mesh on the specified tally
|
||||
mesh_filter = openmc.MeshFilter(mesh)
|
||||
tally.filters = [mesh_filter]
|
||||
|
||||
model.tallies = openmc.Tallies([tally])
|
||||
|
||||
# run the problem
|
||||
sp_filename = model.run()
|
||||
|
||||
with openmc.StatePoint(sp_filename) as sp:
|
||||
mean = sp.tallies[tally.id].mean
|
||||
|
||||
# write the data to a VTK file
|
||||
vtk_filename = 'test.vtk'
|
||||
mesh.write_data_to_vtk(vtk_filename, datasets={'mean': mean})
|
||||
|
||||
# read file
|
||||
reader = vtk.vtkStructuredGridReader()
|
||||
reader.SetFileName(str(vtk_filename))
|
||||
reader.Update()
|
||||
|
||||
# check name of datasets
|
||||
vtk_grid = reader.GetOutput()
|
||||
array = vtk_grid.GetCellData().GetArray(0)
|
||||
vtk_data = nps.vtk_to_numpy(array)
|
||||
|
||||
# convenience function for determining if a mesh
|
||||
# element has vertices in the geometry. This
|
||||
# particular geometry allows us to assume that tally results
|
||||
# in the element should be zero if none of its vertices lie in the geometry
|
||||
def in_geom(cell):
|
||||
point_ids = cell.GetPointIds()
|
||||
|
||||
for i in range(point_ids.GetNumberOfIds()):
|
||||
p = vtk_grid.GetPoint(point_ids.GetId(i))
|
||||
if model.geometry.find(p):
|
||||
return True
|
||||
|
||||
return False
|
||||
|
||||
# reshape mean according to mesh dimensions
|
||||
mean = mean.reshape(mesh.dimension[::-1]).T
|
||||
centroid = [0.0, 0.0, 0.0]
|
||||
|
||||
# check that tally and vtk array results are zero where expected
|
||||
for ijk in mesh.indices:
|
||||
ijk = tuple(n - 1 for n in ijk)
|
||||
# get the cell from the stuctured mesh object
|
||||
cell = vtk_grid.GetCell(*ijk)
|
||||
if not in_geom(cell):
|
||||
cell.GetCentroid(centroid)
|
||||
err_msg = f'IJK: {ijk} should be zero but is not. Centroid: {centroid}'
|
||||
assert mean[ijk] == 0.0, err_msg
|
||||
|
||||
# need to get flat index with axes reversed due to ordering passed into the VTK file
|
||||
flat_idx = np.ravel_multi_index(tuple(ijk[::-1]), mesh.dimension[::-1])
|
||||
assert vtk_data[flat_idx] == 0.0, err_msg
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue