diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 51b5d9ffcc5..2f4f9bfb14a 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2437,24 +2437,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; } } @@ -2566,25 +2564,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; } } diff --git a/tests/unit_tests/test_void.py b/tests/unit_tests/test_void.py new file mode 100644 index 00000000000..be7a09d75e0 --- /dev/null +++ b/tests/unit_tests/test_void.py @@ -0,0 +1,42 @@ +import numpy as np +import openmc +import pytest + + +@pytest.fixture +def empty_sphere(): + openmc.reset_auto_ids() + model = openmc.model.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 + model.settings.source = openmc.IndependentSource(space=openmc.stats.Point()) + + tally = openmc.Tally() + tally.scores = ['fission'] + tally.nuclides = ['U235'] + + 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 + empty_sphere.materials.append(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)