Add the ability to tally microscopic cross sections in void materials with tracklength estimator (#3771)

Co-authored-by: Jonathan Shimwell <drshimwell@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
GuySten 2026-04-04 23:32:51 +03:00 committed by GitHub
parent 60d1dfba7f
commit efc542825e
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
2 changed files with 69 additions and 29 deletions

View file

@ -2436,24 +2436,22 @@ void score_tracklength_tally_general(
if (p.material() != MATERIAL_VOID) {
const auto& mat = model::materials[p.material()];
auto j = mat->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) {
// Determine log union grid index
if (i_log_union == C_NONE) {
int neutron = ParticleType::neutron().transport_index();
i_log_union = std::log(p.E() / data::energy_min[neutron]) /
simulation::log_spacing;
}
// Update micro xs cache
if (!tally.multiply_density()) {
p.update_neutron_xs(i_nuclide, i_log_union);
atom_density = 1.0;
}
} else {
atom_density = tally.multiply_density()
? mat->atom_density(j, p.density_mult())
: 1.0;
if (j != C_NONE)
atom_density = mat->atom_density(j, p.density_mult());
}
if (atom_density > 0) {
if (!tally.multiply_density())
atom_density = 1.0;
} else if (!tally.multiply_density()) {
// Determine log union grid index
if (i_log_union == C_NONE) {
int neutron = ParticleType::neutron().transport_index();
i_log_union = std::log(p.E() / data::energy_min[neutron]) /
simulation::log_spacing;
}
// Update micro xs cache
p.update_neutron_xs(i_nuclide, i_log_union);
atom_density = 1.0;
}
}
@ -2565,25 +2563,25 @@ void score_collision_tally(Particle& p)
double atom_density = 0.;
if (i_nuclide >= 0) {
const auto& mat = model::materials[p.material()];
auto j = mat->mat_nuclide_index_[i_nuclide];
if (j == C_NONE) {
if (p.material() != MATERIAL_VOID) {
const auto& mat = model::materials[p.material()];
auto j = mat->mat_nuclide_index_[i_nuclide];
if (j != C_NONE)
atom_density = mat->atom_density(j, p.density_mult());
}
if (atom_density > 0) {
if (!tally.multiply_density())
atom_density = 1.0;
} else if (!tally.multiply_density()) {
// Determine log union grid index
if (i_log_union == C_NONE) {
int neutron = ParticleType::neutron().transport_index();
i_log_union = std::log(p.E() / data::energy_min[neutron]) /
simulation::log_spacing;
}
// Update micro xs cache
if (!tally.multiply_density()) {
p.update_neutron_xs(i_nuclide, i_log_union);
atom_density = 1.0;
}
} else {
atom_density = tally.multiply_density()
? mat->atom_density(j, p.density_mult())
: 1.0;
p.update_neutron_xs(i_nuclide, i_log_union);
atom_density = 1.0;
}
}

View file

@ -0,0 +1,42 @@
import numpy as np
import openmc
import pytest
@pytest.fixture
def empty_sphere():
openmc.reset_auto_ids()
model = openmc.Model()
surf = openmc.Sphere(r=10, boundary_type='vacuum')
cell = openmc.Cell(region=-surf)
model.geometry = openmc.Geometry([cell])
model.settings.run_mode = 'fixed source'
model.settings.batches = 3
model.settings.particles = 1000
tally = openmc.Tally()
tally.scores = ['total', 'elastic']
tally.nuclides = ['U235']
tally.multiply_density = False
model.tallies.append(tally)
return model
def test_equivalent_microxs(empty_sphere, run_in_tmpdir):
sp_file = empty_sphere.run()
with openmc.StatePoint(sp_file) as sp:
tally1 = sp.tallies[1]
mat = openmc.Material()
mat.add_nuclide('H1', 1e-16)
empty_sphere.geometry.get_all_cells()[1].fill = mat
sp_file = empty_sphere.run()
with openmc.StatePoint(sp_file) as sp:
tally2 = sp.tallies[1]
assert np.isclose(tally1.mean.sum(), tally2.mean.sum(), rtol=1e-10, atol=0)
assert tally1.mean.sum() > 0