mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 14:35:27 -04:00
0169fd9226
Co-authored-by: John Tramm <jtramm@gmail.com> Co-authored-by: Claude Opus 4.6 <noreply@anthropic.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com> Co-authored-by: Copilot <copilot@github.com>
383 lines
13 KiB
Python
383 lines
13 KiB
Python
import os
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
import pytest
|
|
from uncertainties import ufloat
|
|
import openmc
|
|
import openmc.lib
|
|
from openmc.stats import Discrete, Point
|
|
|
|
from tests import cdtemp
|
|
|
|
|
|
@pytest.fixture
|
|
def wws():
|
|
|
|
# weight windows
|
|
ww_files = ('ww_n.txt', 'ww_p.txt')
|
|
cwd = Path(__file__).parent.absolute()
|
|
ww_n_file, ww_p_file = [cwd / Path(f) for f in ww_files]
|
|
|
|
# load pre-generated weight windows
|
|
# (created using the same tally as above)
|
|
ww_n_lower_bnds = np.loadtxt(ww_n_file)
|
|
ww_p_lower_bnds = np.loadtxt(ww_p_file)
|
|
|
|
# create a mesh matching the one used
|
|
# to generate the weight windows
|
|
ww_mesh = openmc.RegularMesh()
|
|
ww_mesh.lower_left = (-240, -240, -240)
|
|
ww_mesh.upper_right = (240, 240, 240)
|
|
ww_mesh.dimension = (5, 6, 7)
|
|
|
|
# energy bounds matching those of the
|
|
# generated weight windows
|
|
e_bnds = [0.0, 0.5, 2E7]
|
|
|
|
ww_n = openmc.WeightWindows(ww_mesh,
|
|
ww_n_lower_bnds,
|
|
None,
|
|
10.0,
|
|
e_bnds,
|
|
survival_ratio=1.01)
|
|
|
|
ww_p = openmc.WeightWindows(ww_mesh,
|
|
ww_p_lower_bnds,
|
|
None,
|
|
10.0,
|
|
e_bnds,
|
|
survival_ratio=1.01)
|
|
|
|
return [ww_n, ww_p]
|
|
|
|
|
|
@pytest.fixture
|
|
def model():
|
|
openmc.reset_auto_ids()
|
|
model = openmc.Model()
|
|
|
|
# materials (M4 steel alloy)
|
|
m4 = openmc.Material()
|
|
m4.set_density('g/cc', 2.3)
|
|
m4.add_nuclide('H1', 0.168018676)
|
|
m4.add_nuclide("H2", 1.93244e-05)
|
|
m4.add_nuclide("O16", 0.561814465)
|
|
m4.add_nuclide("O17", 0.00021401)
|
|
m4.add_nuclide("Na23", 0.021365)
|
|
m4.add_nuclide("Al27", 0.021343)
|
|
m4.add_nuclide("Si28", 0.187439342)
|
|
m4.add_nuclide("Si29", 0.009517714)
|
|
m4.add_nuclide("Si30", 0.006273944)
|
|
m4.add_nuclide("Ca40", 0.018026179)
|
|
m4.add_nuclide("Ca42", 0.00012031)
|
|
m4.add_nuclide("Ca43", 2.51033e-05)
|
|
m4.add_nuclide("Ca44", 0.000387892)
|
|
m4.add_nuclide("Ca46", 7.438e-07)
|
|
m4.add_nuclide("Ca48", 3.47727e-05)
|
|
m4.add_nuclide("Fe54", 0.000248179)
|
|
m4.add_nuclide("Fe56", 0.003895875)
|
|
m4.add_nuclide("Fe57", 8.99727e-05)
|
|
m4.add_nuclide("Fe58", 1.19737e-05)
|
|
|
|
s0 = openmc.Sphere(r=240)
|
|
s1 = openmc.Sphere(r=250, boundary_type='vacuum')
|
|
|
|
c0 = openmc.Cell(fill=m4, region=-s0)
|
|
c1 = openmc.Cell(region=+s0 & -s1)
|
|
|
|
model.geometry = openmc.Geometry([c0, c1])
|
|
|
|
# settings
|
|
settings = model.settings
|
|
settings.run_mode = 'fixed source'
|
|
settings.particles = 500
|
|
settings.batches = 2
|
|
settings.max_history_splits = 100
|
|
settings.photon_transport = True
|
|
space = Point((0.001, 0.001, 0.001))
|
|
energy = Discrete([14E6], [1.0])
|
|
|
|
settings.source = openmc.IndependentSource(space=space, energy=energy)
|
|
|
|
# tally
|
|
mesh = openmc.RegularMesh()
|
|
mesh.lower_left = (-240, -240, -240)
|
|
mesh.upper_right = (240, 240, 240)
|
|
mesh.dimension = (3, 5, 7)
|
|
|
|
mesh_filter = openmc.MeshFilter(mesh)
|
|
|
|
e_bnds = [0.0, 0.5, 2E7]
|
|
energy_filter = openmc.EnergyFilter(e_bnds)
|
|
|
|
particle_filter = openmc.ParticleFilter(['neutron', 'photon'])
|
|
|
|
tally = openmc.Tally()
|
|
tally.filters = [mesh_filter, energy_filter, particle_filter]
|
|
tally.scores = ['flux']
|
|
|
|
model.tallies.append(tally)
|
|
|
|
return model
|
|
|
|
|
|
@pytest.mark.parametrize("shared_secondary", [False, True])
|
|
def test_weightwindows(model, wws, shared_secondary):
|
|
|
|
ww_files = ('ww_n.txt', 'ww_p.txt')
|
|
cwd = Path(__file__).parent.absolute()
|
|
filepaths = [cwd / Path(f) for f in ww_files]
|
|
|
|
with cdtemp(filepaths):
|
|
# run once with variance reduction off
|
|
model.settings.weight_windows_on = False
|
|
model.settings.shared_secondary_bank = shared_secondary
|
|
analog_sp = model.run()
|
|
os.rename(analog_sp, 'statepoint.analog.h5')
|
|
|
|
model.settings.weight_windows = wws
|
|
|
|
# check that string form of the class can be created
|
|
for ww in model.settings.weight_windows:
|
|
str(ww)
|
|
|
|
# run again with variance reduction on
|
|
model.settings.weight_windows_on = True
|
|
ww_sp = model.run()
|
|
os.rename(ww_sp, 'statepoint.ww.h5')
|
|
|
|
# load both statepoints and examine results
|
|
asp = openmc.StatePoint('statepoint.analog.h5')
|
|
wsp = openmc.StatePoint('statepoint.ww.h5')
|
|
|
|
analog_tally = asp.tallies[1]
|
|
ww_tally = wsp.tallies[1]
|
|
|
|
def compare_results(particle, analog_tally, ww_tally):
|
|
# get values from each of the tallies
|
|
an_mean = analog_tally.get_values(filters=[openmc.ParticleFilter],
|
|
filter_bins=[(particle,)])
|
|
ww_mean = ww_tally.get_values(filters=[openmc.ParticleFilter],
|
|
filter_bins=[(particle,)])
|
|
|
|
# expect that more bins were scored with weight windows than
|
|
# the analog run
|
|
assert np.count_nonzero(an_mean) < np.count_nonzero(ww_mean)
|
|
|
|
an_rel_err = analog_tally.get_values(filters=[openmc.ParticleFilter],
|
|
filter_bins=[(particle,)],
|
|
value='rel_err')
|
|
ww_rel_err = ww_tally.get_values(filters=[openmc.ParticleFilter],
|
|
filter_bins=[(particle,)],
|
|
value='rel_err')
|
|
|
|
an_rel_err[an_mean == 0.0] = 1.0
|
|
ww_rel_err[ww_mean == 0.0] = 1.0
|
|
|
|
an_avg_rel_err = np.mean(an_rel_err)
|
|
ww_avg_rel_err = np.mean(ww_rel_err)
|
|
|
|
# expect that the average relative error in the tally
|
|
# decreases
|
|
assert an_avg_rel_err > ww_avg_rel_err
|
|
|
|
# ensure that the value of the mesh bin containing the
|
|
# source is statistically similar in both runs
|
|
an_std_dev = analog_tally.get_values(filters=[openmc.ParticleFilter],
|
|
filter_bins=[(particle,)],
|
|
value='std_dev')
|
|
ww_std_dev = ww_tally.get_values(filters=[openmc.ParticleFilter],
|
|
filter_bins=[(particle,)],
|
|
value='std_dev')
|
|
|
|
# index of the mesh bin containing the source for the higher
|
|
# energy group
|
|
source_bin_idx = (an_mean.shape[0]//2, 0, 0)
|
|
|
|
an_source_bin = ufloat(an_mean[source_bin_idx],
|
|
an_std_dev[source_bin_idx])
|
|
ww_source_bin = ufloat(ww_mean[source_bin_idx],
|
|
ww_std_dev[source_bin_idx])
|
|
|
|
diff = an_source_bin - ww_source_bin
|
|
|
|
# check that values are within two combined standard deviations
|
|
assert abs(diff.nominal_value) / diff.std_dev < 2.0
|
|
|
|
compare_results('neutron', analog_tally, ww_tally)
|
|
compare_results('photon', analog_tally, ww_tally)
|
|
|
|
|
|
def test_lower_ww_bounds_shape():
|
|
"""checks that lower_ww_bounds is reshaped to the mesh dimension when set"""
|
|
ww_mesh = openmc.RegularMesh()
|
|
ww_mesh.lower_left = (-10, -10, -10)
|
|
ww_mesh.upper_right = (10, 10, 10)
|
|
ww_mesh.dimension = (2, 3, 4)
|
|
|
|
ww = openmc.WeightWindows(
|
|
mesh=ww_mesh,
|
|
lower_ww_bounds=[1]*24,
|
|
upper_bound_ratio=5,
|
|
energy_bounds=(1, 1e40)
|
|
)
|
|
assert ww.lower_ww_bounds.shape == (2, 3, 4, 1)
|
|
|
|
|
|
@pytest.mark.parametrize("shared_secondary", [False, True])
|
|
def test_photon_heating(run_in_tmpdir, shared_secondary):
|
|
water = openmc.Material()
|
|
water.add_nuclide('H1', 1.0)
|
|
water.add_nuclide('O16', 2.0)
|
|
water.set_density('g/cm3', 1.0)
|
|
|
|
box = openmc.model.RectangularParallelepiped(
|
|
-300, 300, -300, 300, -300, 300, boundary_type='reflective')
|
|
cell = openmc.Cell(region=-box, fill=water)
|
|
model = openmc.Model()
|
|
model.geometry = openmc.Geometry([cell])
|
|
|
|
mesh = openmc.RegularMesh.from_domain(model.geometry, dimension=(5, 5, 5))
|
|
wwg = openmc.WeightWindowGenerator(mesh, particle_type='photon')
|
|
model.settings.weight_window_generators = [wwg]
|
|
|
|
space = openmc.stats.Point((0, 0, 0))
|
|
energy = openmc.stats.delta_function(5e6)
|
|
model.settings.source = openmc.IndependentSource(
|
|
space=space, energy=energy, particle='photon')
|
|
|
|
model.settings.run_mode = 'fixed source'
|
|
model.settings.batches = 5
|
|
model.settings.particles = 101
|
|
model.settings.shared_secondary_bank = shared_secondary
|
|
|
|
tally = openmc.Tally()
|
|
tally.scores = ['heating']
|
|
tally.filters = [
|
|
openmc.ParticleFilter(['photon']),
|
|
openmc.MeshFilter(mesh)
|
|
]
|
|
model.tallies = [tally]
|
|
|
|
sp_file = model.run()
|
|
with openmc.StatePoint(sp_file) as sp:
|
|
tally_mean = sp.tallies[tally.id].mean
|
|
|
|
# Note: Our current physics model actually does allow this tally to
|
|
# occasionally go slightly negative. However, larger bugs can
|
|
# make this more common. We have selected a particle count for
|
|
# this test that happens to produce no negative tallies for both
|
|
# the shared and non-shared secondary PRNG streams.
|
|
assert np.all(tally_mean >= 0)
|
|
|
|
|
|
def test_roundtrip(run_in_tmpdir, model, wws):
|
|
model.settings.weight_windows = wws
|
|
|
|
# write the model with weight windows to XML
|
|
model.export_to_xml()
|
|
|
|
# ensure that they can be read successfully from XML and that they match the input values
|
|
model_read = openmc.Model.from_xml()
|
|
|
|
zipped_wws = zip(model.settings.weight_windows,
|
|
model_read.settings.weight_windows)
|
|
|
|
# ensure the lower bounds read in from the XML match those of the
|
|
for ww_out, ww_in in zipped_wws:
|
|
assert(ww_out == ww_in)
|
|
|
|
|
|
def test_ww_attrs_python(model):
|
|
mesh = openmc.RegularMesh.from_domain(model.geometry)
|
|
lower_bounds = np.ones(mesh.dimension)
|
|
|
|
# ensure that creation of weight window objects with default arg values
|
|
# is successful
|
|
wws = openmc.WeightWindows(mesh, lower_bounds, upper_bound_ratio=10.0)
|
|
|
|
assert wws.energy_bounds is None
|
|
|
|
wwg = openmc.WeightWindowGenerator(mesh)
|
|
|
|
assert wwg.energy_bounds is None
|
|
|
|
def test_ww_attrs_capi(run_in_tmpdir, model):
|
|
model.export_to_xml()
|
|
|
|
openmc.lib.init()
|
|
|
|
tally = openmc.lib.tallies[model.tallies[0].id]
|
|
|
|
wws = openmc.lib.WeightWindows.from_tally(tally)
|
|
|
|
# this is the first weight window object created
|
|
assert wws.id == 1
|
|
|
|
with pytest.raises(ValueError):
|
|
tally.find_filter(openmc.lib.AzimuthalFilter)
|
|
|
|
mesh_filter = tally.find_filter(openmc.lib.MeshFilter)
|
|
mesh = mesh_filter.mesh
|
|
|
|
assert wws.mesh.id == mesh.id
|
|
|
|
assert wws.particle == openmc.ParticleType.NEUTRON
|
|
|
|
wws.particle = 1
|
|
assert wws.particle == openmc.ParticleType.PHOTON
|
|
wws.particle = 'photon'
|
|
assert wws.particle == openmc.ParticleType.PHOTON
|
|
|
|
with pytest.raises(ValueError):
|
|
wws.particle = '🌠'
|
|
|
|
energy_filter = tally.find_filter(openmc.lib.EnergyFilter)
|
|
np.testing.assert_allclose(np.unique(energy_filter.bins), wws.energy_bounds)
|
|
|
|
# at this point the weight window bounds are uninitialized
|
|
assert all(wws.bounds[0] == -1)
|
|
assert all(wws.bounds[1] == -1)
|
|
|
|
wws = openmc.lib.WeightWindows.from_tally(tally, particle='photon')
|
|
assert wws.id == 2
|
|
assert wws.particle == openmc.ParticleType.PHOTON
|
|
|
|
openmc.lib.finalize()
|
|
|
|
|
|
@pytest.mark.parametrize('library', ('libmesh', 'moab'))
|
|
def test_unstructured_mesh_applied_wws(request, run_in_tmpdir, library):
|
|
"""
|
|
Ensure that weight windows on unstructured mesh work when
|
|
they aren't part of a tally or weight window generator
|
|
"""
|
|
|
|
if library == 'libmesh' and not openmc.lib._libmesh_enabled():
|
|
pytest.skip('LibMesh not enabled in this build.')
|
|
if library == 'moab' and not openmc.lib._dagmc_enabled():
|
|
pytest.skip('DAGMC (and MOAB) mesh not enabled in this build.')
|
|
|
|
water = openmc.Material(name='water')
|
|
water.add_nuclide('H1', 2.0)
|
|
water.add_nuclide('O16', 1.0)
|
|
water.set_density('g/cc', 1.0)
|
|
box = openmc.model.RectangularParallelepiped(*(3*[-10, 10]), boundary_type='vacuum')
|
|
cell = openmc.Cell(region=-box, fill=water)
|
|
|
|
geometry = openmc.Geometry([cell])
|
|
mesh_file = str(request.fspath.dirpath() / 'test_mesh_tets.exo')
|
|
mesh = openmc.UnstructuredMesh(mesh_file, library)
|
|
|
|
dummy_wws = np.ones((12_000,))
|
|
|
|
wws = openmc.WeightWindows(mesh, dummy_wws, upper_bound_ratio=5.0)
|
|
|
|
model = openmc.Model(geometry)
|
|
model.settings.weight_windows = wws
|
|
model.settings.weight_windows_on = True
|
|
model.settings.run_mode = 'fixed source'
|
|
model.settings.particles = 100
|
|
model.settings.batches = 2
|
|
model.run()
|