From 77d074b247e2c2449128d440cff838a85e9cfbf5 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 16 Dec 2024 13:50:29 -0600 Subject: [PATCH 1/8] Track secondary particle energy instead of number of secondaries --- include/openmc/particle_data.h | 9 ++++----- src/particle.cpp | 4 ++-- src/tallies/tally_scoring.cpp | 13 ++----------- 3 files changed, 8 insertions(+), 18 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 164148cce10..e853706f568 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -418,7 +418,7 @@ class ParticleData : public GeometryState { int delayed_group_ {0}; int n_bank_ {0}; - int n_bank_second_ {0}; + double bank_second_E_ {0.0}; double wgt_bank_ {0.0}; int n_delayed_bank_[MAX_DELAYED_GROUPS]; @@ -532,11 +532,10 @@ class ParticleData : public GeometryState { int& delayed_group() { return delayed_group_; } // delayed group // Post-collision data + double& bank_second_E() { return bank_second_E_; } // energy of last reaction secondaries + const double& bank_second_E() const { return bank_second_E_; } + int& n_bank() { return n_bank_; } // number of banked fission sites - int& n_bank_second() - { - return n_bank_second_; - } // number of secondaries banked double& wgt_bank() { return wgt_bank_; } // weight of banked fission sites int* n_delayed_bank() { diff --git a/src/particle.cpp b/src/particle.cpp index 64c50c9438f..9de792b8306 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -101,7 +101,7 @@ void Particle::create_secondary( bank.E = settings::run_CE ? E : g(); bank.time = time(); - n_bank_second() += 1; + bank_second_E() += bank.E; } void Particle::from_source(const SourceSite* src) @@ -356,7 +356,7 @@ void Particle::event_collide() // Reset banked weight during collision n_bank() = 0; - n_bank_second() = 0; + bank_second_E() = 0.0; wgt_bank() = 0.0; zero_delayed_bank(); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index e798161ec2f..02cb4856719 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -964,14 +964,9 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, // The energy deposited is the difference between the pre-collision // and post-collision energy... score = E - p.E(); - // ...less the energy of any secondary particles since they will be // transported individually later - const auto& bank = p.secondary_bank(); - for (auto it = bank.end() - p.n_bank_second(); it < bank.end(); - ++it) { - score -= it->E; - } + score -= p.bank_second_E(); score *= p.wgt_last(); } else { @@ -1500,13 +1495,9 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, // The energy deposited is the difference between the pre-collision and // post-collision energy... score = E - p.E(); - // ...less the energy of any secondary particles since they will be // transported individually later - const auto& bank = p.secondary_bank(); - for (auto it = bank.end() - p.n_bank_second(); it < bank.end(); ++it) { - score -= it->E; - } + score -= p.bank_second_E(); score *= p.wgt_last(); } From a4ccc4ea01174817fbc66566d8df389d94cec892 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 16 Dec 2024 13:51:49 -0600 Subject: [PATCH 2/8] Enxure that the secondary energy attribute is zeroed-out when reviving a particle from a secondary --- src/particle.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/particle.cpp b/src/particle.cpp index 9de792b8306..48ae33362f1 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -417,6 +417,7 @@ void Particle::event_revive_from_secondary() from_source(&secondary_bank().back()); secondary_bank().pop_back(); n_event() = 0; + bank_second_E() = 0.0; // Subtract secondary particle energy from interim pulse-height results if (!model::active_pulse_height_tallies.empty() && From 88dad1be1df581c6f39a9515312a28501e09542f Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 16 Dec 2024 16:12:47 -0600 Subject: [PATCH 3/8] Style --- include/openmc/particle_data.h | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index e853706f568..1166fa53789 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -532,10 +532,13 @@ class ParticleData : public GeometryState { int& delayed_group() { return delayed_group_; } // delayed group // Post-collision data - double& bank_second_E() { return bank_second_E_; } // energy of last reaction secondaries + double& bank_second_E() + { + return bank_second_E_; + } // energy of last reaction secondaries const double& bank_second_E() const { return bank_second_E_; } - int& n_bank() { return n_bank_; } // number of banked fission sites + int& n_bank() { return n_bank_; } // number of banked fission sites double& wgt_bank() { return wgt_bank_; } // weight of banked fission sites int* n_delayed_bank() { From 32f0da91fb9676bc2b37f2f6ef571f232722e34b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 16 Dec 2024 21:31:57 -0600 Subject: [PATCH 4/8] Do not account for energy from particle splits. --- include/openmc/particle.h | 4 +++- src/particle.cpp | 5 +++-- src/weight_windows.cpp | 2 +- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 6a2e67049fd..ecc4df2d5cf 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -50,8 +50,10 @@ class Particle : public ParticleData { //! \param wgt Weight of the secondary particle //! \param u Direction of the secondary particle //! \param E Energy of the secondary particle in [eV] + //! \param accumulate_E Whether to accumulate the energy of the secondary + //! for adjustment to heating tallies //! \param type Particle type - void create_secondary(double wgt, Direction u, double E, ParticleType type); + void create_secondary(double wgt, Direction u, double E, ParticleType type, bool accumulate_E = true); //! initialize from a source site // diff --git a/src/particle.cpp b/src/particle.cpp index 48ae33362f1..546e2508665 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -83,7 +83,7 @@ void Particle::move_distance(double length) } void Particle::create_secondary( - double wgt, Direction u, double E, ParticleType type) + double wgt, Direction u, double E, ParticleType type, bool accumulate_E) { // If energy is below cutoff for this particle, don't create secondary // particle @@ -101,7 +101,8 @@ void Particle::create_secondary( bank.E = settings::run_CE ? E : g(); bank.time = time(); - bank_second_E() += bank.E; + if (accumulate_E) + bank_second_E() += bank.E; } void Particle::from_source(const SourceSite* src) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 68f7550ae5a..546d5dc7506 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -114,7 +114,7 @@ void apply_weight_windows(Particle& p) // Create secondaries and divide weight among all particles int i_split = std::round(n_split); for (int l = 0; l < i_split - 1; l++) { - p.create_secondary(weight / n_split, p.u(), p.E(), p.type()); + p.create_secondary(weight / n_split, p.u(), p.E(), p.type(), false); } // remaining weight is applied to current particle p.wgt() = weight / n_split; From a36eaf22041ee6eb098befb18ac9f6882f9ffcb2 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 16 Dec 2024 21:33:30 -0600 Subject: [PATCH 5/8] Add test for positive heating when applying weight windows --- tests/unit_tests/weightwindows/test.py | 46 ++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index 79aadbef0de..10f4944b2f3 100644 --- a/tests/unit_tests/weightwindows/test.py +++ b/tests/unit_tests/weightwindows/test.py @@ -224,6 +224,52 @@ def test_lower_ww_bounds_shape(): assert ww.lower_ww_bounds.shape == (2, 3, 4, 1) +def test_photon_heating(): + model = openmc.Model() + + water = openmc.Material() + water.add_element('H', 1.0) + water.add_element('O', 2.0) + water.set_density('g/cm3', 1.0) + + box = openmc.model.RectangularParallelepiped(*(3*[-300, 300]), boundary_type='reflective') + cell = openmc.Cell(region=-box, fill=water) + + 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.Discrete([5E6], [1.0]) + + source = openmc.IndependentSource(space=space, energy=energy, particle='photon') + + model.settings.source = source + + model.settings.run_mode = 'fixed source' + model.settings.batches = 5 + model.settings.particles = 100 + + tally = openmc.Tally() + tally.scores = ['heating'] + particle_filter = openmc.ParticleFilter(['photon']) + mesh_filter = openmc.MeshFilter(mesh) + tally.filters = [particle_filter, mesh_filter] + model.tallies = [tally] + + sp_file = model.run() + with openmc.StatePoint(sp_file) as sp: + tally_sum = sp.tallies[tally.id].mean.sum() + + # these values should be nearly identical + # assert np.all(tally.mean > 0) + assert tally_sum >= 0.0 + + def test_roundtrip(run_in_tmpdir, model, wws): model.settings.weight_windows = wws From ed6cb4ba0056563fbf200e7ecd0f50ba0f063329 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 16 Dec 2024 21:53:55 -0600 Subject: [PATCH 6/8] Creating a new method instead of a flag --- include/openmc/particle.h | 10 +++++++--- src/particle.cpp | 15 ++++++++++++--- src/weight_windows.cpp | 2 +- 3 files changed, 20 insertions(+), 7 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index ecc4df2d5cf..b1f9fabd110 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -50,10 +50,14 @@ class Particle : public ParticleData { //! \param wgt Weight of the secondary particle //! \param u Direction of the secondary particle //! \param E Energy of the secondary particle in [eV] - //! \param accumulate_E Whether to accumulate the energy of the secondary - //! for adjustment to heating tallies //! \param type Particle type - void create_secondary(double wgt, Direction u, double E, ParticleType type, bool accumulate_E = true); + void create_secondary(double wgt, Direction u, double E, ParticleType type); + + //! split a particle + // + //! creates a new particle with weight wgt + //! \param wgt Weight of the new particle + void split(double wgt); //! initialize from a source site // diff --git a/src/particle.cpp b/src/particle.cpp index 546e2508665..b53824b6043 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -83,7 +83,7 @@ void Particle::move_distance(double length) } void Particle::create_secondary( - double wgt, Direction u, double E, ParticleType type, bool accumulate_E) + double wgt, Direction u, double E, ParticleType type) { // If energy is below cutoff for this particle, don't create secondary // particle @@ -100,9 +100,18 @@ void Particle::create_secondary( bank.u = u; bank.E = settings::run_CE ? E : g(); bank.time = time(); + bank_second_E() += bank.E; +} - if (accumulate_E) - bank_second_E() += bank.E; +void Particle::split(double wgt) { + secondary_bank().emplace_back(); + auto& bank {secondary_bank().back()}; + bank.particle = type(); + bank.wgt = wgt; + bank.r = r(); + bank.u = u(); + bank.E = settings::run_CE ? E() : g(); + bank.time = time(); } void Particle::from_source(const SourceSite* src) diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 546d5dc7506..2068b14a42f 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -114,7 +114,7 @@ void apply_weight_windows(Particle& p) // Create secondaries and divide weight among all particles int i_split = std::round(n_split); for (int l = 0; l < i_split - 1; l++) { - p.create_secondary(weight / n_split, p.u(), p.E(), p.type(), false); + p.split(weight / n_split); } // remaining weight is applied to current particle p.wgt() = weight / n_split; From 06ec3f40d5417e591b351235560b2998786d8af0 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 16 Dec 2024 21:54:59 -0600 Subject: [PATCH 7/8] Style --- src/particle.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/particle.cpp b/src/particle.cpp index b53824b6043..2b5871c6c75 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -103,7 +103,8 @@ void Particle::create_secondary( bank_second_E() += bank.E; } -void Particle::split(double wgt) { +void Particle::split(double wgt) +{ secondary_bank().emplace_back(); auto& bank {secondary_bank().back()}; bank.particle = type(); From e2c5787dbd39d3b078ff3a8e72c8e49cb5520d78 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 24 Jan 2025 13:00:11 -0600 Subject: [PATCH 8/8] Simplify usage of emplace_back() --- src/particle.cpp | 7 ++--- src/physics.cpp | 9 +++--- src/physics_mg.cpp | 9 +++--- src/track_output.cpp | 4 +-- tests/unit_tests/weightwindows/test.py | 40 +++++++++++--------------- 5 files changed, 29 insertions(+), 40 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index 2b5871c6c75..2276f33ace5 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -91,9 +91,7 @@ void Particle::create_secondary( return; } - secondary_bank().emplace_back(); - - auto& bank {secondary_bank().back()}; + auto& bank = secondary_bank().emplace_back(); bank.particle = type; bank.wgt = wgt; bank.r = r(); @@ -105,8 +103,7 @@ void Particle::create_secondary( void Particle::split(double wgt) { - secondary_bank().emplace_back(); - auto& bank {secondary_bank().back()}; + auto& bank = secondary_bank().emplace_back(); bank.particle = type(); bank.wgt = wgt; bank.r = r(); diff --git a/src/physics.cpp b/src/physics.cpp index 69e74f2eafd..50c1d5d5f17 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -239,11 +239,10 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) } // Write fission particles to nuBank - p.nu_bank().emplace_back(); - NuBank* nu_bank_entry = &p.nu_bank().back(); - nu_bank_entry->wgt = site.wgt; - nu_bank_entry->E = site.E; - nu_bank_entry->delayed_group = site.delayed_group; + NuBank& nu_bank_entry = p.nu_bank().emplace_back(); + nu_bank_entry.wgt = site.wgt; + nu_bank_entry.E = site.E; + nu_bank_entry.delayed_group = site.delayed_group; } // If shared fission bank was full, and no fissions could be added, diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 3cc0532d2b1..65c9916e2ef 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -188,11 +188,10 @@ void create_fission_sites(Particle& p) } // Write fission particles to nuBank - p.nu_bank().emplace_back(); - NuBank* nu_bank_entry = &p.nu_bank().back(); - nu_bank_entry->wgt = site.wgt; - nu_bank_entry->E = site.E; - nu_bank_entry->delayed_group = site.delayed_group; + NuBank& nu_bank_entry = p.nu_bank().emplace_back(); + nu_bank_entry.wgt = site.wgt; + nu_bank_entry.E = site.E; + nu_bank_entry.delayed_group = site.delayed_group; } // If shared fission bank was full, and no fissions could be added, diff --git a/src/track_output.cpp b/src/track_output.cpp index 5c1436de716..f4344d50f74 100644 --- a/src/track_output.cpp +++ b/src/track_output.cpp @@ -31,8 +31,8 @@ int n_tracks_written; //! Number of tracks written void add_particle_track(Particle& p) { - p.tracks().emplace_back(); - p.tracks().back().particle = p.type(); + auto& track = p.tracks().emplace_back(); + track.particle = p.type(); } void write_particle_track(Particle& p) diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index 10f4944b2f3..d6e509522fa 100644 --- a/tests/unit_tests/weightwindows/test.py +++ b/tests/unit_tests/weightwindows/test.py @@ -4,7 +4,6 @@ import numpy as np import pytest from uncertainties import ufloat - import openmc import openmc.lib from openmc.stats import Discrete, Point @@ -224,31 +223,26 @@ def test_lower_ww_bounds_shape(): assert ww.lower_ww_bounds.shape == (2, 3, 4, 1) -def test_photon_heating(): - model = openmc.Model() - +def test_photon_heating(run_in_tmpdir): water = openmc.Material() - water.add_element('H', 1.0) - water.add_element('O', 2.0) + water.add_nuclide('H1', 1.0) + water.add_nuclide('O16', 2.0) water.set_density('g/cm3', 1.0) - box = openmc.model.RectangularParallelepiped(*(3*[-300, 300]), boundary_type='reflective') + 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.Discrete([5E6], [1.0]) - - source = openmc.IndependentSource(space=space, energy=energy, particle='photon') - - model.settings.source = source + 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 @@ -256,18 +250,18 @@ def test_photon_heating(): tally = openmc.Tally() tally.scores = ['heating'] - particle_filter = openmc.ParticleFilter(['photon']) - mesh_filter = openmc.MeshFilter(mesh) - tally.filters = [particle_filter, mesh_filter] + tally.filters = [ + openmc.ParticleFilter(['photon']), + openmc.MeshFilter(mesh) + ] model.tallies = [tally] sp_file = model.run() with openmc.StatePoint(sp_file) as sp: - tally_sum = sp.tallies[tally.id].mean.sum() + tally_mean = sp.tallies[tally.id].mean # these values should be nearly identical - # assert np.all(tally.mean > 0) - assert tally_sum >= 0.0 + assert np.all(tally_mean >= 0) def test_roundtrip(run_in_tmpdir, model, wws): @@ -295,11 +289,11 @@ def test_ww_attrs_python(model): # is successful wws = openmc.WeightWindows(mesh, lower_bounds, upper_bound_ratio=10.0) - assert wws.energy_bounds == None + assert wws.energy_bounds is None wwg = openmc.WeightWindowGenerator(mesh) - assert wwg.energy_bounds == None + assert wwg.energy_bounds is None def test_ww_attrs_capi(run_in_tmpdir, model): model.export_to_xml()