From c1d213f1361d7834868f8bf4e5651b7ed95357fb Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 4 Mar 2026 05:29:49 +0200 Subject: [PATCH 01/21] fix cell last data --- src/particle.cpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 8d0e5b3f06e..7d6cc39ab19 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -174,11 +174,6 @@ void Particle::event_calculate_xs() r_last() = r(); time_last() = time(); - // Reset event variables - event() = TallyEvent::KILL; - event_nuclide() = NUCLIDE_NONE; - event_mt() = REACTION_NONE; - // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles @@ -192,13 +187,18 @@ void Particle::event_calculate_xs() // Set birth cell attribute if (cell_born() == C_NONE) cell_born() = lowest_coord().cell(); + } - // Initialize last cells from current cell - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } - n_coord_last() = n_coord(); + // Initialize last cells from current cell + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); } + n_coord_last() = n_coord(); + + // Reset event variables + event() = TallyEvent::KILL; + event_nuclide() = NUCLIDE_NONE; + event_mt() = REACTION_NONE; // Write particle track. if (write_track()) From cac4e4dc675cb4c50c0121de780b3864b03a58f7 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 4 Mar 2026 05:34:32 +0200 Subject: [PATCH 02/21] update test result --- .../filter_cellfrom/results_true.dat | 52 +++++++++---------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/tests/regression_tests/filter_cellfrom/results_true.dat b/tests/regression_tests/filter_cellfrom/results_true.dat index 408a4965b1d..c210c480959 100644 --- a/tests/regression_tests/filter_cellfrom/results_true.dat +++ b/tests/regression_tests/filter_cellfrom/results_true.dat @@ -1,53 +1,53 @@ k-combined: 9.035025E-02 2.654309E-03 tally 1: -5.994069E+00 -3.594398E+00 +7.467961E+00 +5.580255E+00 tally 2: -4.559707E-04 -2.080739E-08 +0.000000E+00 +0.000000E+00 tally 3: -5.994525E+00 -3.594945E+00 +7.467961E+00 +5.580255E+00 tally 4: 0.000000E+00 0.000000E+00 tally 5: -1.473892E+00 -2.187509E-01 +0.000000E+00 +0.000000E+00 tally 6: -5.223875E-05 -2.992861E-10 +7.821638E-04 +6.132129E-08 tally 7: -1.473945E+00 -2.187663E-01 +7.821638E-04 +6.132129E-08 tally 8: -1.885798E+01 -3.558423E+01 +0.000000E+00 +0.000000E+00 tally 9: 7.467961E+00 5.580255E+00 tally 10: -5.082094E-04 -2.584432E-08 +7.821638E-04 +6.132129E-08 tally 11: -7.468470E+00 -5.581014E+00 +7.468743E+00 +5.581423E+00 tally 12: -1.885798E+01 -3.558423E+01 +0.000000E+00 +0.000000E+00 tally 13: 0.000000E+00 0.000000E+00 tally 14: -2.739543E-04 -7.600983E-09 +0.000000E+00 +0.000000E+00 tally 15: -2.739543E-04 -7.600983E-09 +0.000000E+00 +0.000000E+00 tally 16: -7.881296E+01 -6.221087E+02 +9.767094E+01 +9.547828E+02 tally 17: 1.051397E+02 1.106292E+03 From ee90dabdcd97e4a6c4c6333aa984d6ecf687523b Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 4 Mar 2026 14:31:33 +0200 Subject: [PATCH 03/21] fix material_last update and added test --- src/particle.cpp | 6 ++ tests/unit_tests/test_energy_balance.py | 105 ++++++++++++++++++++++++ 2 files changed, 111 insertions(+) create mode 100644 tests/unit_tests/test_energy_balance.py diff --git a/src/particle.cpp b/src/particle.cpp index 3be64ff56e5..d62d7973ede 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -232,6 +232,9 @@ void Particle::event_calculate_xs() macro_xs().fission = 0.0; macro_xs().nu_fission = 0.0; } + + // Initialize last material from current material + material_last() = material(); } void Particle::event_advance() @@ -298,6 +301,9 @@ void Particle::event_cross_surface() } n_coord_last() = n_coord(); + // Saving previous material data + material_last() = material(); + // Set surface that particle is on and adjust coordinate levels surface() = boundary().surface(); n_coord() = boundary().coord_level(); diff --git a/tests/unit_tests/test_energy_balance.py b/tests/unit_tests/test_energy_balance.py new file mode 100644 index 00000000000..7913159f692 --- /dev/null +++ b/tests/unit_tests/test_energy_balance.py @@ -0,0 +1,105 @@ +import openmc +import pytest + + +@pytest.fixture +def two_cell_model(): + """Simple two-cell slab model with a fixed source. + Cell1 occupies x in [-10, 0], cell2 x in [0, 10]. + """ + openmc.reset_auto_ids() + model = openmc.Model() + + m1 = openmc.Material() + m1.add_element('Li', 1.0) + m1.set_density('g/cm3', 1.0) + + m2 = openmc.Material() + m2.add_element('Al', 1.0) + m2.set_density('g/cm3', 1.0) + + xmin = openmc.XPlane(-10.0, boundary_type="reflective") + xmid = openmc.XPlane(0.0) + xmax = openmc.XPlane(10.0, boundary_type="reflective") + ymin = openmc.YPlane(-10.0, boundary_type="reflective") + ymax = openmc.YPlane(10.0, boundary_type="reflective") + zmin = openmc.ZPlane(-10.0, boundary_type="reflective") + zmax = openmc.ZPlane(10.0, boundary_type="reflective") + cell1 = openmc.Cell(fill=m1, region=+xmin & -xmid & +ymin & -ymax & +zmin & -zmax) + cell2 = openmc.Cell(fill=m2, region=+xmid & -xmax & +ymin & -ymax & +zmin & -zmax) + model.geometry = openmc.Geometry([cell1, cell2]) + + src = openmc.IndependentSource() + src.space = openmc.stats.Point((-5.0, 0.0, 0.0)) + src.particle = 'photon' + src.energy = openmc.stats.Discrete([1e6],[1.0]) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.particles = 100 + model.settings.source = src + + return model, xmid, cell1, cell2, m1, m2 + + +def test_energy_balance_simple(two_cell_model, run_in_tmpdir): + model, xmid, cell1, cell2, *_ = two_cell_model + + tally1 = openmc.Tally() + tally1.scores = ['heating'] + tally1.filters = [openmc.CellFilter([cell1, cell2])] + + tally2 = openmc.Tally() + tally2.scores = ['current'] + tally2.filters = [openmc.EnergyFunctionFilter([0.0,20e6],[0.0,20e6]), openmc.SurfaceFilter([xmid])] + + + model.tallies = [tally1, tally2] + model.run(apply_tally_results=True) + + assert tally1.mean.sum() == pytest.approx(1e6) + + assert tally2.mean[0] == pytest.approx(tally1.mean[1]) + +def test_current_conservation(two_cell_model, run_in_tmpdir): + model, xmid, cell1, cell2, m1, m2 = two_cell_model + + energyfunc_filter = openmc.EnergyFunctionFilter([0.0,20e6],[0.0,20e6]) + + tally1 = openmc.Tally() + tally1.scores = ['current'] + tally1.filters = [energyfunc_filter, openmc.SurfaceFilter([xmid])] + + tally2 = openmc.Tally() + tally2.scores = ['current'] + tally2.filters = [energyfunc_filter, + openmc.CellFromFilter([cell1]), + openmc.CellFilter([cell2])] + + tally3 = openmc.Tally() + tally3.scores = ['current'] + tally3.filters = [energyfunc_filter, + openmc.CellFromFilter([cell2]), + openmc.CellFilter([cell1])] + + tally4 = openmc.Tally() + tally4.scores = ['current'] + tally4.filters = [energyfunc_filter, + openmc.MaterialFromFilter([m1]), + openmc.MaterialFilter([m2])] + + tally5 = openmc.Tally() + tally5.scores = ['current'] + tally5.filters = [energyfunc_filter, + openmc.MaterialFromFilter([m2]), + openmc.MaterialFilter([m1])] + + + model.tallies = [tally1, tally2, tally3, tally4, tally5] + model.run(apply_tally_results=True) + + assert pytest.approx(tally1.mean.sum()) == tally2.mean[0] + tally3.mean[0] + + assert tally2.mean[0] == pytest.approx(tally4.mean[0]) + + assert tally3.mean[0] == pytest.approx(tally5.mean[0]) From 923a4a107229df981f35abb2260f69656679332b Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 08:38:16 +0200 Subject: [PATCH 04/21] fix cell last update location and added test --- src/particle.cpp | 35 +++++++++++++-------- tests/unit_tests/test_energy_balance.py | 41 +++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 13 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index d62d7973ede..01463b19b60 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -174,6 +174,11 @@ void Particle::event_calculate_xs() r_last() = r(); time_last() = time(); + // Reset event variables + event() = TallyEvent::KILL; + event_nuclide() = NUCLIDE_NONE; + event_mt() = REACTION_NONE; + // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles @@ -187,18 +192,16 @@ void Particle::event_calculate_xs() // Set birth cell attribute if (cell_born() == C_NONE) cell_born() = lowest_coord().cell(); - } - // Initialize last cells from current cell - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } - n_coord_last() = n_coord(); + // Initialize last cells from current cell + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + n_coord_last() = n_coord(); - // Reset event variables - event() = TallyEvent::KILL; - event_nuclide() = NUCLIDE_NONE; - event_mt() = REACTION_NONE; + // Initialize last material from current material + material_last() = material(); + } // Write particle track. if (write_track()) @@ -232,9 +235,6 @@ void Particle::event_calculate_xs() macro_xs().fission = 0.0; macro_xs().nu_fission = 0.0; } - - // Initialize last material from current material - material_last() = material(); } void Particle::event_advance() @@ -425,6 +425,15 @@ void Particle::event_collide() if (!model::active_tallies.empty()) score_collision_derivative(*this); + // Saving previous cell data + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + n_coord_last() = n_coord(); + + // Saving previous material data + material_last() = material(); + #ifdef OPENMC_DAGMC_ENABLED history().reset(); #endif diff --git a/tests/unit_tests/test_energy_balance.py b/tests/unit_tests/test_energy_balance.py index 7913159f692..f8d19a2323b 100644 --- a/tests/unit_tests/test_energy_balance.py +++ b/tests/unit_tests/test_energy_balance.py @@ -1,5 +1,6 @@ import openmc import pytest +import numpy as np @pytest.fixture @@ -41,6 +42,31 @@ def two_cell_model(): return model, xmid, cell1, cell2, m1, m2 +@pytest.fixture +def two_sphere_model(): + openmc.reset_auto_ids() + model = openmc.Model() + + mat = openmc.Material() + mat.add_nuclide('H2', 1.0) + mat.set_density('g/cm3', 1.0) + + surf1 = openmc.Sphere(r=0.01) + surf2 = openmc.Sphere(r=1000, boundary_type='vacuum') + cell1 = openmc.Cell(region=-surf1) + cell2 = openmc.Cell(fill=mat, region=+surf1 & -surf2) + + model.geometry = openmc.Geometry([cell1, cell2]) + + src = openmc.IndependentSource() + src.energy = openmc.stats.Discrete([1e6],[1.0]) + + model.settings.run_mode = 'fixed source' + model.settings.batches = 1 + model.settings.particles = 1000 + model.settings.source = src + + return model, cell1, cell2 def test_energy_balance_simple(two_cell_model, run_in_tmpdir): model, xmid, cell1, cell2, *_ = two_cell_model @@ -103,3 +129,18 @@ def test_current_conservation(two_cell_model, run_in_tmpdir): assert tally2.mean[0] == pytest.approx(tally4.mean[0]) assert tally3.mean[0] == pytest.approx(tally5.mean[0]) + +def test_cellfrom_heating(run_in_tmpdir, two_sphere_model): + + model, cell1, cell2 = two_sphere_model + + tally = openmc.Tally() + tally.scores = ['heating'] + tally.filters = [openmc.CellFromFilter([cell1, cell2]), openmc.CellFilter([cell2])] + + model.tallies = [tally] + + model.run(apply_tally_results=True) + + assert np.all(tally.mean > 0) + From ee97a888a02bfec2c5f139ad7380df6f34659211 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 09:04:26 +0200 Subject: [PATCH 05/21] try fix --- src/particle.cpp | 11 ++-- .../filter_cellfrom/results_true.dat | 52 +++++++++---------- 2 files changed, 32 insertions(+), 31 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 01463b19b60..abbb2c0e0ea 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -182,6 +182,7 @@ void Particle::event_calculate_xs() // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles + bool update_mat = false; if (lowest_coord().cell() == C_NONE) { if (!exhaustive_find_cell(*this)) { mark_as_lost( @@ -199,8 +200,7 @@ void Particle::event_calculate_xs() } n_coord_last() = n_coord(); - // Initialize last material from current material - material_last() = material(); + update_mat = true; } // Write particle track. @@ -235,6 +235,10 @@ void Particle::event_calculate_xs() macro_xs().fission = 0.0; macro_xs().nu_fission = 0.0; } + + // Initialize last material from current material + if (update_mat) + material_last() = material(); } void Particle::event_advance() @@ -431,9 +435,6 @@ void Particle::event_collide() } n_coord_last() = n_coord(); - // Saving previous material data - material_last() = material(); - #ifdef OPENMC_DAGMC_ENABLED history().reset(); #endif diff --git a/tests/regression_tests/filter_cellfrom/results_true.dat b/tests/regression_tests/filter_cellfrom/results_true.dat index c210c480959..1f9d0446104 100644 --- a/tests/regression_tests/filter_cellfrom/results_true.dat +++ b/tests/regression_tests/filter_cellfrom/results_true.dat @@ -1,53 +1,53 @@ k-combined: 9.035025E-02 2.654309E-03 tally 1: -7.467961E+00 -5.580255E+00 +6.586649E+00 +4.340013E+00 tally 2: -0.000000E+00 -0.000000E+00 +4.559707E-04 +2.080739E-08 tally 3: -7.467961E+00 -5.580255E+00 +6.587105E+00 +4.340613E+00 tally 4: 0.000000E+00 0.000000E+00 tally 5: -0.000000E+00 -0.000000E+00 +8.813120E-01 +7.829296E-02 tally 6: -7.821638E-04 -6.132129E-08 +5.223875E-05 +2.992861E-10 tally 7: -7.821638E-04 -6.132129E-08 +8.813643E-01 +7.830214E-02 tally 8: -0.000000E+00 -0.000000E+00 +5.908648E+00 +3.492172E+00 tally 9: 7.467961E+00 5.580255E+00 tally 10: -7.821638E-04 -6.132129E-08 +5.082094E-04 +2.584432E-08 tally 11: -7.468743E+00 -5.581423E+00 +7.468470E+00 +5.581014E+00 tally 12: -0.000000E+00 -0.000000E+00 +5.908648E+00 +3.492172E+00 tally 13: 0.000000E+00 0.000000E+00 tally 14: -0.000000E+00 -0.000000E+00 +2.739543E-04 +7.600983E-09 tally 15: -0.000000E+00 -0.000000E+00 +2.739543E-04 +7.600983E-09 tally 16: -9.767094E+01 -9.547828E+02 +9.176229E+01 +8.428689E+02 tally 17: 1.051397E+02 1.106292E+03 From d2f9a4c2e38ad4fa59e211baa8afd48cec7faf3f Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 09:07:16 +0200 Subject: [PATCH 06/21] try simplification --- src/particle.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index abbb2c0e0ea..93f48a9e28b 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -182,7 +182,6 @@ void Particle::event_calculate_xs() // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles - bool update_mat = false; if (lowest_coord().cell() == C_NONE) { if (!exhaustive_find_cell(*this)) { mark_as_lost( @@ -200,7 +199,8 @@ void Particle::event_calculate_xs() } n_coord_last() = n_coord(); - update_mat = true; + // Initialize last material from current material + material_last() = material(); } // Write particle track. @@ -235,10 +235,6 @@ void Particle::event_calculate_xs() macro_xs().fission = 0.0; macro_xs().nu_fission = 0.0; } - - // Initialize last material from current material - if (update_mat) - material_last() = material(); } void Particle::event_advance() From 527b176499e1447e507699d45533ba941e03a63e Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 09:09:25 +0200 Subject: [PATCH 07/21] Revert "try simplification" This reverts commit d2f9a4c2e38ad4fa59e211baa8afd48cec7faf3f. --- src/particle.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 93f48a9e28b..abbb2c0e0ea 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -182,6 +182,7 @@ void Particle::event_calculate_xs() // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles + bool update_mat = false; if (lowest_coord().cell() == C_NONE) { if (!exhaustive_find_cell(*this)) { mark_as_lost( @@ -199,8 +200,7 @@ void Particle::event_calculate_xs() } n_coord_last() = n_coord(); - // Initialize last material from current material - material_last() = material(); + update_mat = true; } // Write particle track. @@ -235,6 +235,10 @@ void Particle::event_calculate_xs() macro_xs().fission = 0.0; macro_xs().nu_fission = 0.0; } + + // Initialize last material from current material + if (update_mat) + material_last() = material(); } void Particle::event_advance() From dbdb6d53f35244b1793779cefee667ef34198d5a Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 09:45:35 +0200 Subject: [PATCH 08/21] another fix --- include/openmc/constants.h | 2 +- src/particle.cpp | 6 +++--- tests/unit_tests/test_energy_balance.py | 20 +++++++++++++------- 3 files changed, 17 insertions(+), 11 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 0b425a673dc..fa7a48ac618 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -110,7 +110,7 @@ constexpr array SUBSHELLS = {"K", "L1", "L2", "L3", "M1", "M2", // Void material and nuclide // TODO: refactor and remove -constexpr int MATERIAL_VOID {-1}; +constexpr int MATERIAL_VOID {-2}; constexpr int NUCLIDE_NONE {-1}; // ============================================================================ diff --git a/src/particle.cpp b/src/particle.cpp index abbb2c0e0ea..899d236fe9a 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -182,7 +182,6 @@ void Particle::event_calculate_xs() // If the cell hasn't been determined based on the particle's location, // initiate a search for the current cell. This generally happens at the // beginning of the history and again for any secondary particles - bool update_mat = false; if (lowest_coord().cell() == C_NONE) { if (!exhaustive_find_cell(*this)) { mark_as_lost( @@ -200,7 +199,8 @@ void Particle::event_calculate_xs() } n_coord_last() = n_coord(); - update_mat = true; + if (material_last() != C_NONE) + material_last() = material(); } // Write particle track. @@ -237,7 +237,7 @@ void Particle::event_calculate_xs() } // Initialize last material from current material - if (update_mat) + if (material_last() == C_NONE) material_last() = material(); } diff --git a/tests/unit_tests/test_energy_balance.py b/tests/unit_tests/test_energy_balance.py index f8d19a2323b..15360e3be17 100644 --- a/tests/unit_tests/test_energy_balance.py +++ b/tests/unit_tests/test_energy_balance.py @@ -66,7 +66,7 @@ def two_sphere_model(): model.settings.particles = 1000 model.settings.source = src - return model, cell1, cell2 + return model, cell1, cell2, mat def test_energy_balance_simple(two_cell_model, run_in_tmpdir): model, xmid, cell1, cell2, *_ = two_cell_model @@ -132,15 +132,21 @@ def test_current_conservation(two_cell_model, run_in_tmpdir): def test_cellfrom_heating(run_in_tmpdir, two_sphere_model): - model, cell1, cell2 = two_sphere_model + model, cell1, cell2, mat = two_sphere_model - tally = openmc.Tally() - tally.scores = ['heating'] - tally.filters = [openmc.CellFromFilter([cell1, cell2]), openmc.CellFilter([cell2])] + tally1 = openmc.Tally() + tally1.scores = ['heating'] + tally1.filters = [openmc.CellFromFilter([cell1, cell2]), openmc.CellFilter([cell2])] + + tally2 = openmc.Tally() + tally2.scores = ['heating'] + tally2.filters = [openmc.MaterialFromFilter([mat]), openmc.MaterialFilter([mat])] - model.tallies = [tally] + model.tallies = [tally1, tally2] model.run(apply_tally_results=True) - assert np.all(tally.mean > 0) + assert np.all(tally1.mean > 0) + + assert tally1.mean[1] == tally2.mean[0] From 26e472118a1be0374d1c7bf53ff22e8cf23371a1 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 11:29:18 +0200 Subject: [PATCH 09/21] apply fix fix --- include/openmc/constants.h | 3 ++- src/particle.cpp | 12 ++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index fa7a48ac618..44d0b3e9e9b 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -110,7 +110,8 @@ constexpr array SUBSHELLS = {"K", "L1", "L2", "L3", "M1", "M2", // Void material and nuclide // TODO: refactor and remove -constexpr int MATERIAL_VOID {-2}; +constexpr int MATERIAL_VOID {-1}; +constexpr int MATERIAL_UNSET {std::numeric_limits::min()}; constexpr int NUCLIDE_NONE {-1}; // ============================================================================ diff --git a/src/particle.cpp b/src/particle.cpp index 899d236fe9a..ef7b8e30b6f 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -199,7 +199,7 @@ void Particle::event_calculate_xs() } n_coord_last() = n_coord(); - if (material_last() != C_NONE) + if (material_last() != MATERIAL_UNSET) material_last() = material(); } @@ -237,7 +237,7 @@ void Particle::event_calculate_xs() } // Initialize last material from current material - if (material_last() == C_NONE) + if (material_last() == MATERIAL_UNSET) material_last() = material(); } @@ -407,10 +407,6 @@ void Particle::event_collide() // Save coordinates for tallying purposes r_last_current() = r(); - // Set last material to none since cross sections will need to be - // re-evaluated - material_last() = C_NONE; - // Set all directions to base level -- right now, after a collision, only // the base level directions are changed for (int j = 0; j < n_coord() - 1; ++j) { @@ -435,6 +431,10 @@ void Particle::event_collide() } n_coord_last() = n_coord(); + // Set last material to unset since cross sections will need to be + // re-evaluated + material_last() = MATERIAL_UNSET; + #ifdef OPENMC_DAGMC_ENABLED history().reset(); #endif From 05a8c6a6504fcbd40db744448da3f0cfc9218935 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 12:12:32 +0200 Subject: [PATCH 10/21] try another direction --- include/openmc/constants.h | 1 - include/openmc/particle_data.h | 2 ++ src/material.cpp | 1 + src/particle.cpp | 16 +++++----------- 4 files changed, 8 insertions(+), 12 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 44d0b3e9e9b..0b425a673dc 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -111,7 +111,6 @@ constexpr array SUBSHELLS = {"K", "L1", "L2", "L3", "M1", "M2", // Void material and nuclide // TODO: refactor and remove constexpr int MATERIAL_VOID {-1}; -constexpr int MATERIAL_UNSET {std::numeric_limits::min()}; constexpr int NUCLIDE_NONE {-1}; // ============================================================================ diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 5b632bbce53..a8781a0bccf 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -210,6 +210,8 @@ struct MacroXS { double incoherent; //!< macroscopic incoherent xs double photoelectric; //!< macroscopic photoelectric xs double pair_production; //!< macroscopic pair production xs + + double last_E {0.0}; //!< last evaluated energy in [eV] }; //============================================================================== diff --git a/src/material.cpp b/src/material.cpp index 21b11b8b9ce..16ad78cf12c 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -817,6 +817,7 @@ void Material::calculate_xs(Particle& p) const p.macro_xs().absorption = 0.0; p.macro_xs().fission = 0.0; p.macro_xs().nu_fission = 0.0; + p.macro_xs().last_E = p.E(); if (p.type().is_neutron()) { this->calculate_neutron_xs(p); diff --git a/src/particle.cpp b/src/particle.cpp index ef7b8e30b6f..e34c89da0c5 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -199,8 +199,7 @@ void Particle::event_calculate_xs() } n_coord_last() = n_coord(); - if (material_last() != MATERIAL_UNSET) - material_last() = material(); + material_last() = material(); } // Write particle track. @@ -214,9 +213,9 @@ void Particle::event_calculate_xs() if (material() != MATERIAL_VOID) { if (settings::run_CE) { if (material() != material_last() || sqrtkT() != sqrtkT_last() || - density_mult() != density_mult_last()) { + density_mult() != density_mult_last() || E() != macro_xs().last_E) { // If the material is the same as the last material and the - // temperature hasn't changed, we don't need to lookup cross + // energy and temperature hasn't changed, we don't need to lookup cross // sections again. model::materials[material()]->calculate_xs(*this); } @@ -235,10 +234,6 @@ void Particle::event_calculate_xs() macro_xs().fission = 0.0; macro_xs().nu_fission = 0.0; } - - // Initialize last material from current material - if (material_last() == MATERIAL_UNSET) - material_last() = material(); } void Particle::event_advance() @@ -431,9 +426,8 @@ void Particle::event_collide() } n_coord_last() = n_coord(); - // Set last material to unset since cross sections will need to be - // re-evaluated - material_last() = MATERIAL_UNSET; + // Saving previous material data + material_last() = material(); #ifdef OPENMC_DAGMC_ENABLED history().reset(); From 8839a850cb1cc7b758209fc134ef4b24af9e655c Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 12:42:50 +0200 Subject: [PATCH 11/21] wip --- src/particle.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/particle.cpp b/src/particle.cpp index e34c89da0c5..946a6801594 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -233,6 +233,7 @@ void Particle::event_calculate_xs() macro_xs().absorption = 0.0; macro_xs().fission = 0.0; macro_xs().nu_fission = 0.0; + macro_xs().last_E = E(); } } From 6505fb063f26945994b263d5c7792a1834a3382c Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 13:52:26 +0200 Subject: [PATCH 12/21] wip --- openmc/data/decay.py | 30 +++++++++++++++++++++++++----- openmc/deplete/r2s.py | 23 ++++++++++++++++++----- openmc/material.py | 10 +++++++--- tests/unit_tests/test_r2s.py | 8 ++++---- 4 files changed, 54 insertions(+), 17 deletions(-) diff --git a/openmc/data/decay.py b/openmc/data/decay.py index 7cd4bf43d41..ef616384366 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -564,7 +564,7 @@ def sources(self): _DECAY_PHOTON_ENERGY = {} -def decay_photon_energy(nuclide: str) -> Univariate | None: +def decay_photon_energy(nuclide: str, chain = None) -> Univariate | None: """Get photon energy distribution resulting from the decay of a nuclide This function relies on data stored in a depletion chain. Before calling it @@ -577,7 +577,9 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: ---------- nuclide : str Name of nuclide, e.g., 'Co58' - + chain : Chain, optional. + Chain file to get decay energy from + Returns ------- openmc.stats.Univariate or None @@ -585,6 +587,15 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: if no photon source exists. Note that the probabilities represent intensities, given as [Bq/atom] (in other words, decay constants). """ + from openmc.deplete.chain import Chain + if chain is not None: + cv.check_type('chain', chain, Chain) + for nuc in chain.nuclides: + if nuc.name == nuclide: + return nuc.sources.get('photon') + else: + return + if not _DECAY_PHOTON_ENERGY: chain_file = openmc.config.get('chain_file') if chain_file is None: @@ -593,7 +604,6 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: "openmc.config['chain_file'] in order to load decay data." ) - from openmc.deplete import Chain chain = Chain.from_xml(chain_file) for nuc in chain.nuclides: if 'photon' in nuc.sources: @@ -610,7 +620,7 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: _DECAY_ENERGY = {} -def decay_energy(nuclide: str): +def decay_energy(nuclide: str, chain = None): """Get decay energy value resulting from the decay of a nuclide This function relies on data stored in a depletion chain. Before calling it @@ -623,6 +633,8 @@ def decay_energy(nuclide: str): ---------- nuclide : str Name of nuclide, e.g., 'H3' + chain : Chain, optional. + Chain file to get decay energy from Returns ------- @@ -630,6 +642,15 @@ def decay_energy(nuclide: str): Decay energy of nuclide in [eV]. If the nuclide is stable, a value of 0.0 is returned. """ + from openmc.deplete.chain import Chain + if chain is not None: + cv.check_type('chain', chain, Chain) + for nuclide in chain.nuclides: + if nuc.name == nuclide and nuc.decay_energy: + return nuc.decay_energy + else: + return 0.0 + if not _DECAY_ENERGY: chain_file = openmc.config.get('chain_file') if chain_file is None: @@ -638,7 +659,6 @@ def decay_energy(nuclide: str): "openmc.config['chain_file'] in order to load decay data." ) - from openmc.deplete import Chain chain = Chain.from_xml(chain_file) for nuc in chain.nuclides: if nuc.decay_energy: diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 57bbe437ffb..9283a128836 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -230,7 +230,7 @@ def run( ) self.step3_photon_transport( photon_time_indices, bounding_boxes, output_dir / 'photon_transport', - mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs + mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs, chain_file = chain_file ) return output_dir @@ -405,6 +405,7 @@ def step3_photon_transport( output_dir: PathLike = 'photon_transport', mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, + chain_file: PathLike | None = None, ): """Run the photon transport step. @@ -437,12 +438,17 @@ def step3_photon_transport( run_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.Model.run` during the photon transport step. By default, output is disabled. + chain_file : PathLike, optional + Path to the depletion chain XML file to use during activation. If + not provided, the default configured chain file will be used. """ # TODO: Automatically determine bounding box for each cell if bounding_boxes is None and self.method == 'cell-based': raise ValueError("bounding_boxes must be provided for cell-based " "R2S calculations.") + from openmc.deplete import Chain + chain = Chain.from_xml(chain_file) if chain_file is not None else None # Set default run arguments if not provided if run_kwargs is None: @@ -490,7 +496,7 @@ def step3_photon_transport( # Create decay photon source if self.method == 'mesh-based': self.photon_model.settings.source = \ - self.get_decay_photon_source_mesh(time_index) + self.get_decay_photon_source_mesh(time_index, chain_file=chain_file) else: sources = [] results = self.results['depletion_results'] @@ -510,7 +516,7 @@ def step3_photon_transport( # Create decay photon source source space = openmc.stats.Box(*bounding_box) - energy = activated_mat.get_decay_photon_energy() + energy = activated_mat.get_decay_photon_energy(chain=chain) strength = energy.integral() if energy is not None else 0.0 source = openmc.IndependentSource( space=space, @@ -539,7 +545,8 @@ def step3_photon_transport( def get_decay_photon_source_mesh( self, - time_index: int = -1 + time_index: int = -1, + chain_file: PathLike | None = None, ) -> list[openmc.IndependentSource]: """Create decay photon source for a mesh-based calculation. @@ -558,6 +565,9 @@ def get_decay_photon_source_mesh( ---------- time_index : int, optional Time index for the decay photon source. Default is -1 (last time). + chain_file : PathLike, optional + Path to the depletion chain XML file to use during activation. If + not provided, the default configured chain file will be used. Returns ------- @@ -566,6 +576,9 @@ def get_decay_photon_source_mesh( each mesh element-material combination with non-zero source strength. """ + from openmc.deplete import Chain + chain = Chain.from_xml(chain_file) if chain_file is not None else None + mat_dict = self.neutron_model._get_all_materials() # List to hold all sources @@ -607,7 +620,7 @@ def get_decay_photon_source_mesh( activated_mat = results[time_index].get_material(str(original_mat.id)) # Create decay photon source - energy = activated_mat.get_decay_photon_energy() + energy = activated_mat.get_decay_photon_energy(chain=chain) if energy is not None: strength = energy.integral() space = openmc.stats.Box(*bbox) diff --git a/openmc/material.py b/openmc/material.py index 239bc01f771..7b5fb8ceada 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -28,7 +28,6 @@ from openmc.data.function import Tabulated1D from openmc.data import mass_energy_absorption_coefficient, dose_coefficients - # Units for density supported by OpenMC DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum', 'macro') @@ -338,7 +337,8 @@ def get_decay_photon_energy( units: str = 'Bq', volume: float | None = None, exclude_nuclides: list[str] | None = None, - include_nuclides: list[str] | None = None + include_nuclides: list[str] | None = None, + chain: Chain | None = None, ) -> Univariate | None: r"""Return energy distribution of decay photons from unstable nuclides. @@ -359,6 +359,8 @@ def get_decay_photon_energy( include_nuclides : list of str, optional Nuclides to include in the photon source calculation. If specified, only these nuclides are used. + chain : Chain, optional. + Chain file to get decay energy from Returns ------- @@ -367,7 +369,9 @@ def get_decay_photon_energy( the total intensity of the photon source in the requested units. """ + from openmc.deplete import Chain cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'}) + cv.check_type('chain', chain, Chain, none_ok=True) if exclude_nuclides is not None and include_nuclides is not None: raise ValueError("Cannot specify both exclude_nuclides and include_nuclides") @@ -391,7 +395,7 @@ def get_decay_photon_energy( if include_nuclides is not None and nuc not in include_nuclides: continue - source_per_atom = openmc.data.decay_photon_energy(nuc) + source_per_atom = openmc.data.decay_photon_energy(nuc, chain=chain) if source_per_atom is not None and atoms_per_bcm > 0.0: dists.append(source_per_atom) probs.append(1e24 * atoms_per_bcm * multiplier) diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py index a94f85c8c08..f0b325573eb 100644 --- a/tests/unit_tests/test_r2s.py +++ b/tests/unit_tests/test_r2s.py @@ -53,7 +53,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): r2s = R2SManager(model, mesh) # Use custom reduced chain file for Ni - chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml") + chain_file = Path(__file__).parents[1] / "chain_ni.xml" # Run R2S calculation outdir = r2s.run( @@ -61,7 +61,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): source_rates=[1.0], photon_time_indices=[1], output_dir=tmp_path, - chain_file=chain, + chain_file=chain_file, ) # Check directories and files exist @@ -105,7 +105,7 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): r2s = R2SManager(model, [c1, c2]) # Use custom reduced chain file for Ni - chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml") + chain_file = Path(__file__).parents[1] / "chain_ni.xml" # Run R2S calculation bounding_boxes = {c1.id: c1.bounding_box, c2.id: c2.bounding_box} @@ -115,7 +115,7 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): photon_time_indices=[1], output_dir=tmp_path, bounding_boxes=bounding_boxes, - chain_file=chain + chain_file=chain_file ) # Check directories and files exist From 7c59cd6f7a3851fa83339aef997393233301cab1 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 15:02:09 +0200 Subject: [PATCH 13/21] wip --- src/mgxs.cpp | 1 + src/particle.cpp | 5 ----- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/mgxs.cpp b/src/mgxs.cpp index a0fe4060ddb..640c990c634 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -615,6 +615,7 @@ void Mgxs::calculate_xs(Particle& p) p.macro_xs().nu_fission = fissionable ? xs[temperature].nu_fission(angle, p.g()) * p.density_mult() : 0.; + p.macro_xs().last_E = p.E(); } //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index 946a6801594..167d9367b17 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -198,8 +198,6 @@ void Particle::event_calculate_xs() cell_last(j) = coord(j).cell(); } n_coord_last() = n_coord(); - - material_last() = material(); } // Write particle track. @@ -301,9 +299,6 @@ void Particle::event_cross_surface() } n_coord_last() = n_coord(); - // Saving previous material data - material_last() = material(); - // Set surface that particle is on and adjust coordinate levels surface() = boundary().surface(); n_coord() = boundary().coord_level(); From 09649f85bb7879bf3578083234e00ceb2a5e759f Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 16:53:37 +0200 Subject: [PATCH 14/21] wip --- src/geometry.cpp | 21 +++++++++++++++------ src/particle.cpp | 28 +++++++++++++++++++--------- 2 files changed, 34 insertions(+), 15 deletions(-) diff --git a/src/geometry.cpp b/src/geometry.cpp index ddb61385f18..e8e8a065112 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -173,12 +173,21 @@ bool find_cell_inner( } // Set the material, temperature and density multiplier. - p.material_last() = p.material(); - p.material() = c.material(p.cell_instance()); - p.sqrtkT_last() = p.sqrtkT(); - p.sqrtkT() = c.sqrtkT(p.cell_instance()); - p.density_mult_last() = p.density_mult(); - p.density_mult() = c.density_mult(p.cell_instance()); + if (p.sqrtkT() != -1) { + // if current material is set, set last material + p.material_last() = p.material(); + p.material() = c.material(p.cell_instance()); + p.sqrtkT_last() = p.sqrtkT(); + p.sqrtkT() = c.sqrtkT(p.cell_instance()); + p.density_mult_last() = p.density_mult(); + p.density_mult() = c.density_mult(p.cell_instance()); + } else { + // set both current data and last data to current data + p.material_last() = p.material() = c.material(p.cell_instance()); + p.sqrtkT_last() = p.sqrtkT() = c.sqrtkT(p.cell_instance()); + p.density_mult_last() = p.density_mult() = + c.density_mult(p.cell_instance()); + } return true; diff --git a/src/particle.cpp b/src/particle.cpp index 167d9367b17..69929dfc8ea 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -198,6 +198,8 @@ void Particle::event_calculate_xs() cell_last(j) = coord(j).cell(); } n_coord_last() = n_coord(); + + material_last() = material(); } // Write particle track. @@ -293,12 +295,6 @@ void Particle::event_advance() void Particle::event_cross_surface() { - // Saving previous cell data - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } - n_coord_last() = n_coord(); - // Set surface that particle is on and adjust coordinate levels surface() = boundary().surface(); n_coord() = boundary().coord_level(); @@ -310,6 +306,12 @@ void Particle::event_cross_surface() boundary().lattice_translation()[2] != 0) { // Particle crosses lattice boundary + // Saving previous cell data + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + n_coord_last() = n_coord(); + bool verbose = settings::verbosity >= 10 || trace(); cross_lattice(*this, boundary(), verbose); event() = TallyEvent::LATTICE; @@ -478,6 +480,8 @@ void Particle::event_revive_from_secondary() cell_last(j) = coord(j).cell(); } n_coord_last() = n_coord(); + + material_last() = material(); } pht_secondary_particles(); } @@ -574,7 +578,6 @@ void Particle::cross_surface(const Surface& surf) if (surf.geom_type() == GeometryType::CSG) history().reset(); #endif - // Handle any applicable boundary conditions. if (surf.bc_ && settings::run_mode != RunMode::PLOTTING && settings::run_mode != RunMode::VOLUME) { @@ -582,6 +585,12 @@ void Particle::cross_surface(const Surface& surf) return; } + // Saving previous cell data + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + n_coord_last() = n_coord(); + // ========================================================================== // SEARCH NEIGHBOR LISTS FOR NEXT CELL @@ -700,10 +709,11 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u) // Set the new particle direction u() = new_u; - // Reassign particle's cell and surface + // Reassign particle's cell, material and surface coord(0).cell() = cell_last(0); - surface() = -surface(); + material() = material_last(); + surface() = -surface(); // If a reflective surface is coincident with a lattice or universe // boundary, it is necessary to redetermine the particle's coordinates in // the lower universes. From c4ed3e5b7d08e47921a20734127d6288739a1ae6 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 16:55:29 +0200 Subject: [PATCH 15/21] make test less flaky --- tests/unit_tests/test_energy_balance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_energy_balance.py b/tests/unit_tests/test_energy_balance.py index 15360e3be17..d1c07b417b4 100644 --- a/tests/unit_tests/test_energy_balance.py +++ b/tests/unit_tests/test_energy_balance.py @@ -148,5 +148,5 @@ def test_cellfrom_heating(run_in_tmpdir, two_sphere_model): assert np.all(tally1.mean > 0) - assert tally1.mean[1] == tally2.mean[0] + assert tally1.mean[1] == pytest.approx(tally2.mean[0]) From 7bd6f0c7c7f97fad695e34da58cc8c8355bc9b62 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 16:59:30 +0200 Subject: [PATCH 16/21] drop unneeded changes --- openmc/data/decay.py | 30 +++++------------------------- openmc/deplete/r2s.py | 23 +++++------------------ openmc/material.py | 10 +++------- tests/unit_tests/test_r2s.py | 8 ++++---- 4 files changed, 17 insertions(+), 54 deletions(-) diff --git a/openmc/data/decay.py b/openmc/data/decay.py index ef616384366..7cd4bf43d41 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -564,7 +564,7 @@ def sources(self): _DECAY_PHOTON_ENERGY = {} -def decay_photon_energy(nuclide: str, chain = None) -> Univariate | None: +def decay_photon_energy(nuclide: str) -> Univariate | None: """Get photon energy distribution resulting from the decay of a nuclide This function relies on data stored in a depletion chain. Before calling it @@ -577,9 +577,7 @@ def decay_photon_energy(nuclide: str, chain = None) -> Univariate | None: ---------- nuclide : str Name of nuclide, e.g., 'Co58' - chain : Chain, optional. - Chain file to get decay energy from - + Returns ------- openmc.stats.Univariate or None @@ -587,15 +585,6 @@ def decay_photon_energy(nuclide: str, chain = None) -> Univariate | None: if no photon source exists. Note that the probabilities represent intensities, given as [Bq/atom] (in other words, decay constants). """ - from openmc.deplete.chain import Chain - if chain is not None: - cv.check_type('chain', chain, Chain) - for nuc in chain.nuclides: - if nuc.name == nuclide: - return nuc.sources.get('photon') - else: - return - if not _DECAY_PHOTON_ENERGY: chain_file = openmc.config.get('chain_file') if chain_file is None: @@ -604,6 +593,7 @@ def decay_photon_energy(nuclide: str, chain = None) -> Univariate | None: "openmc.config['chain_file'] in order to load decay data." ) + from openmc.deplete import Chain chain = Chain.from_xml(chain_file) for nuc in chain.nuclides: if 'photon' in nuc.sources: @@ -620,7 +610,7 @@ def decay_photon_energy(nuclide: str, chain = None) -> Univariate | None: _DECAY_ENERGY = {} -def decay_energy(nuclide: str, chain = None): +def decay_energy(nuclide: str): """Get decay energy value resulting from the decay of a nuclide This function relies on data stored in a depletion chain. Before calling it @@ -633,8 +623,6 @@ def decay_energy(nuclide: str, chain = None): ---------- nuclide : str Name of nuclide, e.g., 'H3' - chain : Chain, optional. - Chain file to get decay energy from Returns ------- @@ -642,15 +630,6 @@ def decay_energy(nuclide: str, chain = None): Decay energy of nuclide in [eV]. If the nuclide is stable, a value of 0.0 is returned. """ - from openmc.deplete.chain import Chain - if chain is not None: - cv.check_type('chain', chain, Chain) - for nuclide in chain.nuclides: - if nuc.name == nuclide and nuc.decay_energy: - return nuc.decay_energy - else: - return 0.0 - if not _DECAY_ENERGY: chain_file = openmc.config.get('chain_file') if chain_file is None: @@ -659,6 +638,7 @@ def decay_energy(nuclide: str, chain = None): "openmc.config['chain_file'] in order to load decay data." ) + from openmc.deplete import Chain chain = Chain.from_xml(chain_file) for nuc in chain.nuclides: if nuc.decay_energy: diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 9283a128836..57bbe437ffb 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -230,7 +230,7 @@ def run( ) self.step3_photon_transport( photon_time_indices, bounding_boxes, output_dir / 'photon_transport', - mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs, chain_file = chain_file + mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs ) return output_dir @@ -405,7 +405,6 @@ def step3_photon_transport( output_dir: PathLike = 'photon_transport', mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, - chain_file: PathLike | None = None, ): """Run the photon transport step. @@ -438,17 +437,12 @@ def step3_photon_transport( run_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.Model.run` during the photon transport step. By default, output is disabled. - chain_file : PathLike, optional - Path to the depletion chain XML file to use during activation. If - not provided, the default configured chain file will be used. """ # TODO: Automatically determine bounding box for each cell if bounding_boxes is None and self.method == 'cell-based': raise ValueError("bounding_boxes must be provided for cell-based " "R2S calculations.") - from openmc.deplete import Chain - chain = Chain.from_xml(chain_file) if chain_file is not None else None # Set default run arguments if not provided if run_kwargs is None: @@ -496,7 +490,7 @@ def step3_photon_transport( # Create decay photon source if self.method == 'mesh-based': self.photon_model.settings.source = \ - self.get_decay_photon_source_mesh(time_index, chain_file=chain_file) + self.get_decay_photon_source_mesh(time_index) else: sources = [] results = self.results['depletion_results'] @@ -516,7 +510,7 @@ def step3_photon_transport( # Create decay photon source source space = openmc.stats.Box(*bounding_box) - energy = activated_mat.get_decay_photon_energy(chain=chain) + energy = activated_mat.get_decay_photon_energy() strength = energy.integral() if energy is not None else 0.0 source = openmc.IndependentSource( space=space, @@ -545,8 +539,7 @@ def step3_photon_transport( def get_decay_photon_source_mesh( self, - time_index: int = -1, - chain_file: PathLike | None = None, + time_index: int = -1 ) -> list[openmc.IndependentSource]: """Create decay photon source for a mesh-based calculation. @@ -565,9 +558,6 @@ def get_decay_photon_source_mesh( ---------- time_index : int, optional Time index for the decay photon source. Default is -1 (last time). - chain_file : PathLike, optional - Path to the depletion chain XML file to use during activation. If - not provided, the default configured chain file will be used. Returns ------- @@ -576,9 +566,6 @@ def get_decay_photon_source_mesh( each mesh element-material combination with non-zero source strength. """ - from openmc.deplete import Chain - chain = Chain.from_xml(chain_file) if chain_file is not None else None - mat_dict = self.neutron_model._get_all_materials() # List to hold all sources @@ -620,7 +607,7 @@ def get_decay_photon_source_mesh( activated_mat = results[time_index].get_material(str(original_mat.id)) # Create decay photon source - energy = activated_mat.get_decay_photon_energy(chain=chain) + energy = activated_mat.get_decay_photon_energy() if energy is not None: strength = energy.integral() space = openmc.stats.Box(*bbox) diff --git a/openmc/material.py b/openmc/material.py index 7b5fb8ceada..239bc01f771 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -28,6 +28,7 @@ from openmc.data.function import Tabulated1D from openmc.data import mass_energy_absorption_coefficient, dose_coefficients + # Units for density supported by OpenMC DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum', 'macro') @@ -337,8 +338,7 @@ def get_decay_photon_energy( units: str = 'Bq', volume: float | None = None, exclude_nuclides: list[str] | None = None, - include_nuclides: list[str] | None = None, - chain: Chain | None = None, + include_nuclides: list[str] | None = None ) -> Univariate | None: r"""Return energy distribution of decay photons from unstable nuclides. @@ -359,8 +359,6 @@ def get_decay_photon_energy( include_nuclides : list of str, optional Nuclides to include in the photon source calculation. If specified, only these nuclides are used. - chain : Chain, optional. - Chain file to get decay energy from Returns ------- @@ -369,9 +367,7 @@ def get_decay_photon_energy( the total intensity of the photon source in the requested units. """ - from openmc.deplete import Chain cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'}) - cv.check_type('chain', chain, Chain, none_ok=True) if exclude_nuclides is not None and include_nuclides is not None: raise ValueError("Cannot specify both exclude_nuclides and include_nuclides") @@ -395,7 +391,7 @@ def get_decay_photon_energy( if include_nuclides is not None and nuc not in include_nuclides: continue - source_per_atom = openmc.data.decay_photon_energy(nuc, chain=chain) + source_per_atom = openmc.data.decay_photon_energy(nuc) if source_per_atom is not None and atoms_per_bcm > 0.0: dists.append(source_per_atom) probs.append(1e24 * atoms_per_bcm * multiplier) diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py index f0b325573eb..a94f85c8c08 100644 --- a/tests/unit_tests/test_r2s.py +++ b/tests/unit_tests/test_r2s.py @@ -53,7 +53,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): r2s = R2SManager(model, mesh) # Use custom reduced chain file for Ni - chain_file = Path(__file__).parents[1] / "chain_ni.xml" + chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml") # Run R2S calculation outdir = r2s.run( @@ -61,7 +61,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): source_rates=[1.0], photon_time_indices=[1], output_dir=tmp_path, - chain_file=chain_file, + chain_file=chain, ) # Check directories and files exist @@ -105,7 +105,7 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): r2s = R2SManager(model, [c1, c2]) # Use custom reduced chain file for Ni - chain_file = Path(__file__).parents[1] / "chain_ni.xml" + chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml") # Run R2S calculation bounding_boxes = {c1.id: c1.bounding_box, c2.id: c2.bounding_box} @@ -115,7 +115,7 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): photon_time_indices=[1], output_dir=tmp_path, bounding_boxes=bounding_boxes, - chain_file=chain_file + chain_file=chain ) # Check directories and files exist From 3a7fa7768480e159ad470731c23b9cf2ff0425b6 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 17:01:45 +0200 Subject: [PATCH 17/21] update test --- tests/regression_tests/filter_cellfrom/results_true.dat | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/regression_tests/filter_cellfrom/results_true.dat b/tests/regression_tests/filter_cellfrom/results_true.dat index 1f9d0446104..c193c034ab3 100644 --- a/tests/regression_tests/filter_cellfrom/results_true.dat +++ b/tests/regression_tests/filter_cellfrom/results_true.dat @@ -22,8 +22,8 @@ tally 7: 8.813643E-01 7.830214E-02 tally 8: -5.908648E+00 -3.492172E+00 +0.000000E+00 +0.000000E+00 tally 9: 7.467961E+00 5.580255E+00 From 99baa5dacfd6e59a989376b829537cb5b243d0ba Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 18:15:20 +0200 Subject: [PATCH 18/21] wip --- src/particle.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 69929dfc8ea..0ee5d857ab6 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -295,6 +295,13 @@ void Particle::event_advance() void Particle::event_cross_surface() { + + // Saving previous cell data + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + n_coord_last() = n_coord(); + // Set surface that particle is on and adjust coordinate levels surface() = boundary().surface(); n_coord() = boundary().coord_level(); @@ -306,12 +313,6 @@ void Particle::event_cross_surface() boundary().lattice_translation()[2] != 0) { // Particle crosses lattice boundary - // Saving previous cell data - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } - n_coord_last() = n_coord(); - bool verbose = settings::verbosity >= 10 || trace(); cross_lattice(*this, boundary(), verbose); event() = TallyEvent::LATTICE; From 3e13c1074cd31636e75b39e27cee9318d4b144f0 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 18:17:25 +0200 Subject: [PATCH 19/21] wip --- src/particle.cpp | 6 ------ tests/regression_tests/filter_cellfrom/results_true.dat | 4 ++-- 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 0ee5d857ab6..a0d26345aec 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -586,12 +586,6 @@ void Particle::cross_surface(const Surface& surf) return; } - // Saving previous cell data - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } - n_coord_last() = n_coord(); - // ========================================================================== // SEARCH NEIGHBOR LISTS FOR NEXT CELL diff --git a/tests/regression_tests/filter_cellfrom/results_true.dat b/tests/regression_tests/filter_cellfrom/results_true.dat index c193c034ab3..1f9d0446104 100644 --- a/tests/regression_tests/filter_cellfrom/results_true.dat +++ b/tests/regression_tests/filter_cellfrom/results_true.dat @@ -22,8 +22,8 @@ tally 7: 8.813643E-01 7.830214E-02 tally 8: -0.000000E+00 -0.000000E+00 +5.908648E+00 +3.492172E+00 tally 9: 7.467961E+00 5.580255E+00 From 7801f0f1e42be98242fbdb88a7d15132241c6fb3 Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 18:30:33 +0200 Subject: [PATCH 20/21] Revert "drop unneeded changes" This reverts commit 7bd6f0c7c7f97fad695e34da58cc8c8355bc9b62. --- openmc/data/decay.py | 30 +++++++++++++++++++++++++----- openmc/deplete/r2s.py | 23 ++++++++++++++++++----- openmc/material.py | 10 +++++++--- tests/unit_tests/test_r2s.py | 8 ++++---- 4 files changed, 54 insertions(+), 17 deletions(-) diff --git a/openmc/data/decay.py b/openmc/data/decay.py index 7cd4bf43d41..ef616384366 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -564,7 +564,7 @@ def sources(self): _DECAY_PHOTON_ENERGY = {} -def decay_photon_energy(nuclide: str) -> Univariate | None: +def decay_photon_energy(nuclide: str, chain = None) -> Univariate | None: """Get photon energy distribution resulting from the decay of a nuclide This function relies on data stored in a depletion chain. Before calling it @@ -577,7 +577,9 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: ---------- nuclide : str Name of nuclide, e.g., 'Co58' - + chain : Chain, optional. + Chain file to get decay energy from + Returns ------- openmc.stats.Univariate or None @@ -585,6 +587,15 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: if no photon source exists. Note that the probabilities represent intensities, given as [Bq/atom] (in other words, decay constants). """ + from openmc.deplete.chain import Chain + if chain is not None: + cv.check_type('chain', chain, Chain) + for nuc in chain.nuclides: + if nuc.name == nuclide: + return nuc.sources.get('photon') + else: + return + if not _DECAY_PHOTON_ENERGY: chain_file = openmc.config.get('chain_file') if chain_file is None: @@ -593,7 +604,6 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: "openmc.config['chain_file'] in order to load decay data." ) - from openmc.deplete import Chain chain = Chain.from_xml(chain_file) for nuc in chain.nuclides: if 'photon' in nuc.sources: @@ -610,7 +620,7 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: _DECAY_ENERGY = {} -def decay_energy(nuclide: str): +def decay_energy(nuclide: str, chain = None): """Get decay energy value resulting from the decay of a nuclide This function relies on data stored in a depletion chain. Before calling it @@ -623,6 +633,8 @@ def decay_energy(nuclide: str): ---------- nuclide : str Name of nuclide, e.g., 'H3' + chain : Chain, optional. + Chain file to get decay energy from Returns ------- @@ -630,6 +642,15 @@ def decay_energy(nuclide: str): Decay energy of nuclide in [eV]. If the nuclide is stable, a value of 0.0 is returned. """ + from openmc.deplete.chain import Chain + if chain is not None: + cv.check_type('chain', chain, Chain) + for nuclide in chain.nuclides: + if nuc.name == nuclide and nuc.decay_energy: + return nuc.decay_energy + else: + return 0.0 + if not _DECAY_ENERGY: chain_file = openmc.config.get('chain_file') if chain_file is None: @@ -638,7 +659,6 @@ def decay_energy(nuclide: str): "openmc.config['chain_file'] in order to load decay data." ) - from openmc.deplete import Chain chain = Chain.from_xml(chain_file) for nuc in chain.nuclides: if nuc.decay_energy: diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 57bbe437ffb..9283a128836 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -230,7 +230,7 @@ def run( ) self.step3_photon_transport( photon_time_indices, bounding_boxes, output_dir / 'photon_transport', - mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs + mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs, chain_file = chain_file ) return output_dir @@ -405,6 +405,7 @@ def step3_photon_transport( output_dir: PathLike = 'photon_transport', mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, + chain_file: PathLike | None = None, ): """Run the photon transport step. @@ -437,12 +438,17 @@ def step3_photon_transport( run_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.Model.run` during the photon transport step. By default, output is disabled. + chain_file : PathLike, optional + Path to the depletion chain XML file to use during activation. If + not provided, the default configured chain file will be used. """ # TODO: Automatically determine bounding box for each cell if bounding_boxes is None and self.method == 'cell-based': raise ValueError("bounding_boxes must be provided for cell-based " "R2S calculations.") + from openmc.deplete import Chain + chain = Chain.from_xml(chain_file) if chain_file is not None else None # Set default run arguments if not provided if run_kwargs is None: @@ -490,7 +496,7 @@ def step3_photon_transport( # Create decay photon source if self.method == 'mesh-based': self.photon_model.settings.source = \ - self.get_decay_photon_source_mesh(time_index) + self.get_decay_photon_source_mesh(time_index, chain_file=chain_file) else: sources = [] results = self.results['depletion_results'] @@ -510,7 +516,7 @@ def step3_photon_transport( # Create decay photon source source space = openmc.stats.Box(*bounding_box) - energy = activated_mat.get_decay_photon_energy() + energy = activated_mat.get_decay_photon_energy(chain=chain) strength = energy.integral() if energy is not None else 0.0 source = openmc.IndependentSource( space=space, @@ -539,7 +545,8 @@ def step3_photon_transport( def get_decay_photon_source_mesh( self, - time_index: int = -1 + time_index: int = -1, + chain_file: PathLike | None = None, ) -> list[openmc.IndependentSource]: """Create decay photon source for a mesh-based calculation. @@ -558,6 +565,9 @@ def get_decay_photon_source_mesh( ---------- time_index : int, optional Time index for the decay photon source. Default is -1 (last time). + chain_file : PathLike, optional + Path to the depletion chain XML file to use during activation. If + not provided, the default configured chain file will be used. Returns ------- @@ -566,6 +576,9 @@ def get_decay_photon_source_mesh( each mesh element-material combination with non-zero source strength. """ + from openmc.deplete import Chain + chain = Chain.from_xml(chain_file) if chain_file is not None else None + mat_dict = self.neutron_model._get_all_materials() # List to hold all sources @@ -607,7 +620,7 @@ def get_decay_photon_source_mesh( activated_mat = results[time_index].get_material(str(original_mat.id)) # Create decay photon source - energy = activated_mat.get_decay_photon_energy() + energy = activated_mat.get_decay_photon_energy(chain=chain) if energy is not None: strength = energy.integral() space = openmc.stats.Box(*bbox) diff --git a/openmc/material.py b/openmc/material.py index 239bc01f771..7b5fb8ceada 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -28,7 +28,6 @@ from openmc.data.function import Tabulated1D from openmc.data import mass_energy_absorption_coefficient, dose_coefficients - # Units for density supported by OpenMC DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum', 'macro') @@ -338,7 +337,8 @@ def get_decay_photon_energy( units: str = 'Bq', volume: float | None = None, exclude_nuclides: list[str] | None = None, - include_nuclides: list[str] | None = None + include_nuclides: list[str] | None = None, + chain: Chain | None = None, ) -> Univariate | None: r"""Return energy distribution of decay photons from unstable nuclides. @@ -359,6 +359,8 @@ def get_decay_photon_energy( include_nuclides : list of str, optional Nuclides to include in the photon source calculation. If specified, only these nuclides are used. + chain : Chain, optional. + Chain file to get decay energy from Returns ------- @@ -367,7 +369,9 @@ def get_decay_photon_energy( the total intensity of the photon source in the requested units. """ + from openmc.deplete import Chain cv.check_value('units', units, {'Bq', 'Bq/g', 'Bq/kg', 'Bq/cm3'}) + cv.check_type('chain', chain, Chain, none_ok=True) if exclude_nuclides is not None and include_nuclides is not None: raise ValueError("Cannot specify both exclude_nuclides and include_nuclides") @@ -391,7 +395,7 @@ def get_decay_photon_energy( if include_nuclides is not None and nuc not in include_nuclides: continue - source_per_atom = openmc.data.decay_photon_energy(nuc) + source_per_atom = openmc.data.decay_photon_energy(nuc, chain=chain) if source_per_atom is not None and atoms_per_bcm > 0.0: dists.append(source_per_atom) probs.append(1e24 * atoms_per_bcm * multiplier) diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py index a94f85c8c08..f0b325573eb 100644 --- a/tests/unit_tests/test_r2s.py +++ b/tests/unit_tests/test_r2s.py @@ -53,7 +53,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): r2s = R2SManager(model, mesh) # Use custom reduced chain file for Ni - chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml") + chain_file = Path(__file__).parents[1] / "chain_ni.xml" # Run R2S calculation outdir = r2s.run( @@ -61,7 +61,7 @@ def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): source_rates=[1.0], photon_time_indices=[1], output_dir=tmp_path, - chain_file=chain, + chain_file=chain_file, ) # Check directories and files exist @@ -105,7 +105,7 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): r2s = R2SManager(model, [c1, c2]) # Use custom reduced chain file for Ni - chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml") + chain_file = Path(__file__).parents[1] / "chain_ni.xml" # Run R2S calculation bounding_boxes = {c1.id: c1.bounding_box, c2.id: c2.bounding_box} @@ -115,7 +115,7 @@ def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): photon_time_indices=[1], output_dir=tmp_path, bounding_boxes=bounding_boxes, - chain_file=chain + chain_file=chain_file ) # Check directories and files exist From 781a8d9f0b9e61f91b01cfbbc82ff5b5848fa16b Mon Sep 17 00:00:00 2001 From: GuySten Date: Fri, 6 Mar 2026 18:50:58 +0200 Subject: [PATCH 21/21] hard fix --- src/particle.cpp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index a0d26345aec..953497cf376 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -213,7 +213,8 @@ void Particle::event_calculate_xs() if (material() != MATERIAL_VOID) { if (settings::run_CE) { if (material() != material_last() || sqrtkT() != sqrtkT_last() || - density_mult() != density_mult_last() || E() != macro_xs().last_E) { + density_mult() != density_mult_last() || E() != macro_xs().last_E || + true) { // If the material is the same as the last material and the // energy and temperature hasn't changed, we don't need to lookup cross // sections again. @@ -296,10 +297,6 @@ void Particle::event_advance() void Particle::event_cross_surface() { - // Saving previous cell data - for (int j = 0; j < n_coord(); ++j) { - cell_last(j) = coord(j).cell(); - } n_coord_last() = n_coord(); // Set surface that particle is on and adjust coordinate levels @@ -313,6 +310,11 @@ void Particle::event_cross_surface() boundary().lattice_translation()[2] != 0) { // Particle crosses lattice boundary + // Saving previous cell data + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + bool verbose = settings::verbosity >= 10 || trace(); cross_lattice(*this, boundary(), verbose); event() = TallyEvent::LATTICE; @@ -586,6 +588,11 @@ void Particle::cross_surface(const Surface& surf) return; } + // Saving previous cell data + for (int j = 0; j < n_coord(); ++j) { + cell_last(j) = coord(j).cell(); + } + // ========================================================================== // SEARCH NEIGHBOR LISTS FOR NEXT CELL