mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-25 12:35:29 -04:00
db426b66ca
Co-authored-by: GuySten <guyste@post.bgu.ac.il> Co-authored-by: GuySten <62616591+GuySten@users.noreply.github.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
195 lines
6.5 KiB
Python
195 lines
6.5 KiB
Python
from math import pi, log, log10
|
|
from random import uniform, normalvariate
|
|
|
|
import numpy as np
|
|
|
|
import openmc.deplete
|
|
import openmc
|
|
import pytest
|
|
|
|
|
|
@pytest.fixture
|
|
def model():
|
|
"""Sphere of single nuclide"""
|
|
model = openmc.Model()
|
|
|
|
w = openmc.Material(name='tungsten')
|
|
w.add_nuclide('W186', 1.0)
|
|
w.set_density('g/cm3', 19.3)
|
|
w.depletable = True
|
|
|
|
r = uniform(1.0, 10.0)
|
|
w.volume = 4/3 * pi * r**3
|
|
|
|
surf = openmc.Sphere(r=r, boundary_type='vacuum')
|
|
cell = openmc.Cell(fill=w, region=-surf)
|
|
model.geometry = openmc.Geometry([cell])
|
|
|
|
model.settings.batches = 10
|
|
model.settings.particles = 1000
|
|
model.settings.source = openmc.IndependentSource(
|
|
space=openmc.stats.Point(),
|
|
energy=openmc.stats.Discrete([1.0e6], [1.0])
|
|
)
|
|
model.settings.run_mode = 'fixed source'
|
|
|
|
rx_tally = openmc.Tally(name='activation tally')
|
|
rx_tally.scores = ['(n,gamma)']
|
|
model.tallies.append(rx_tally)
|
|
|
|
return model
|
|
|
|
|
|
ENERGIES = np.logspace(log10(1e-5), log10(2e7), 100)
|
|
|
|
|
|
@pytest.mark.parametrize("reaction_rate_mode,reaction_rate_opts,tolerance", [
|
|
("direct", {}, 1e-5),
|
|
("flux", {'energies': ENERGIES}, 0.1),
|
|
("flux", {'energies': ENERGIES, 'reactions': ['(n,gamma)']}, 1e-5),
|
|
("flux", {'energies': ENERGIES, 'reactions': ['(n,gamma)'], 'nuclides': ['W186', 'H3']}, 1e-2),
|
|
])
|
|
@pytest.mark.flaky(reruns=1)
|
|
def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts, tolerance):
|
|
# Determine (n.gamma) reaction rate using initial run
|
|
sp = model.run()
|
|
with openmc.StatePoint(sp) as sp:
|
|
tally = sp.get_tally(name='activation tally')
|
|
capture_rate = tally.mean.flat[0]
|
|
|
|
# Create one-nuclide depletion chain
|
|
chain = openmc.deplete.Chain()
|
|
w186 = openmc.deplete.Nuclide('W186')
|
|
w186.add_reaction('(n,gamma)', None, 0.0, 1.0)
|
|
chain.add_nuclide(w186)
|
|
|
|
# Create transport operator
|
|
op = openmc.deplete.CoupledOperator(
|
|
model, chain,
|
|
normalization_mode="source-rate",
|
|
reaction_rate_mode=reaction_rate_mode,
|
|
reaction_rate_opts=reaction_rate_opts,
|
|
)
|
|
|
|
# To determine the source rate necessary to reduce W186 density in half, we
|
|
# start with the single-nuclide transmutation equation:
|
|
#
|
|
# dn/dt = -f * sigma * phi * n
|
|
# n(t) = n0 * exp(-f * sigma * phi * t)
|
|
#
|
|
# where f is the source rate. The capture rate, r, is sigma * phi * n0,
|
|
# meaning that:
|
|
#
|
|
# n(t) = n0 * exp(-f * r * t / n0)
|
|
#
|
|
# To reduce the density by half, we would need:
|
|
#
|
|
# n(t)/n0 = exp(-f * r * t / n0) = 1/2
|
|
# f = n0 / (r * t) ln(2)
|
|
#
|
|
# So we need to know the initial number of atoms (n0), the capture rate (r),
|
|
# and choose an irradiation time (t)
|
|
|
|
w = model.geometry.get_materials_by_name('tungsten')[0]
|
|
atom_densities = w.get_nuclide_atom_densities()
|
|
atom_per_cc = 1e24 * atom_densities['W186'] # Density in atom/cm^3
|
|
n0 = atom_per_cc * w.volume # Absolute number of atoms
|
|
|
|
# Pick a random irradiation time and then determine necessary source rate to
|
|
# reduce material by half
|
|
t = uniform(1.0, 5.0) * 86400
|
|
source_rates = [n0/(capture_rate*t) * log(2.0)]
|
|
|
|
# Now activate the material
|
|
integrator = openmc.deplete.PredictorIntegrator(
|
|
op, [t], source_rates=source_rates
|
|
)
|
|
integrator.integrate()
|
|
|
|
# Get resulting number of atoms
|
|
results = openmc.deplete.Results('depletion_results.h5')
|
|
_, atoms = results.get_atoms(w, "W186")
|
|
|
|
assert atoms[0] == pytest.approx(n0)
|
|
assert atoms[1] / atoms[0] == pytest.approx(0.5, rel=tolerance)
|
|
|
|
# Check that material name is preserved in depletion results
|
|
step_result = results[0]
|
|
mat_from_results = step_result.get_material(f"{w.id}")
|
|
assert mat_from_results.name == 'tungsten'
|
|
|
|
|
|
def test_decay(run_in_tmpdir):
|
|
"""Test decay-only timesteps where no transport solve is performed"""
|
|
|
|
# Create a model with a single nuclide, Sr89
|
|
mat = openmc.Material()
|
|
mat.add_nuclide('Sr89', 1.0)
|
|
mat.set_density('g/cm3', 1.0)
|
|
mat.depletable = True
|
|
r = 5.0
|
|
mat.volume = 4/3 * pi * r**3
|
|
surf = openmc.Sphere(r=r, boundary_type='vacuum')
|
|
cell = openmc.Cell(fill=mat, region=-surf)
|
|
geometry = openmc.Geometry([cell])
|
|
settings = openmc.Settings()
|
|
settings.batches = 10
|
|
settings.particles = 1000
|
|
settings.run_mode = 'fixed source'
|
|
|
|
# Create depletion chain with only Sr89 and sample its half-life. Note that
|
|
# currently at least one reaction has to exist in the depletion chain
|
|
chain = openmc.deplete.Chain()
|
|
sr89 = openmc.deplete.Nuclide('Sr89')
|
|
sr89.half_life = normalvariate(4365792.0, 6048.0)
|
|
sr89.add_decay_mode('beta-', None, 1.0)
|
|
sr89.add_reaction('(n,gamma)', None, 0.0, 1.0)
|
|
chain.add_nuclide(sr89)
|
|
chain.export_to_xml('test_chain.xml')
|
|
|
|
model = openmc.Model(geometry=geometry, settings=settings)
|
|
# Create transport operator
|
|
op = openmc.deplete.CoupledOperator(
|
|
model, 'test_chain.xml', normalization_mode="source-rate"
|
|
)
|
|
|
|
# Deplete with two decay steps
|
|
integrator = openmc.deplete.PredictorIntegrator(
|
|
op, [sr89.half_life, 2*sr89.half_life], source_rates=[0.0, 0.0]
|
|
)
|
|
integrator.integrate()
|
|
|
|
# Get resulting number of atoms
|
|
results = openmc.deplete.Results('depletion_results.h5')
|
|
_, atoms = results.get_atoms(mat, "Sr89")
|
|
|
|
# Ensure density goes down by a factor of 2 after each half-life
|
|
assert atoms[1] / atoms[0] == pytest.approx(0.5)
|
|
assert atoms[2] / atoms[1] == pytest.approx(0.25)
|
|
|
|
|
|
def test_flux_rr_missing_nuclide(run_in_tmpdir, model):
|
|
# Create two-nuclide depletion chain -- since W184 is not in the model, this
|
|
# test ensures that FluxCollapseHelper loads missing nuclides appropriately
|
|
chain = openmc.deplete.Chain()
|
|
w184 = openmc.deplete.Nuclide('W184')
|
|
w184.add_reaction('(n,gamma)', None, 0.0, 1.0)
|
|
chain.add_nuclide(w184)
|
|
w186 = openmc.deplete.Nuclide('W186')
|
|
w186.add_reaction('(n,gamma)', None, 0.0, 1.0)
|
|
chain.add_nuclide(w186)
|
|
chain.export_to_xml('test_chain.xml')
|
|
|
|
# Create transport operator
|
|
op = openmc.deplete.CoupledOperator(
|
|
model, 'test_chain.xml',
|
|
normalization_mode="source-rate",
|
|
reaction_rate_mode="flux",
|
|
reaction_rate_opts={'energies': [0.0, 20.0e6]},
|
|
)
|
|
|
|
# Deplete with two decay steps
|
|
integrator = openmc.deplete.PredictorIntegrator(
|
|
op, [100.0], source_rates=[10.0]
|
|
)
|
|
integrator.integrate()
|