From 7a18108724edbdc0adf004308c3955824d9b102c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Fri, 24 Jan 2025 18:06:15 -0600 Subject: [PATCH] Adjust for secondary particle energy directly in heating scores (#3227) Co-authored-by: Paul Romano --- include/openmc/particle.h | 6 ++++ include/openmc/particle_data.h | 12 ++++--- src/particle.cpp | 19 ++++++++--- src/physics.cpp | 9 +++-- src/physics_mg.cpp | 9 +++-- src/tallies/tally_scoring.cpp | 13 ++------ src/track_output.cpp | 4 +-- src/weight_windows.cpp | 2 +- tests/unit_tests/weightwindows/test.py | 46 ++++++++++++++++++++++++-- 9 files changed, 83 insertions(+), 37 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 6a2e67049fd..b1f9fabd110 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -53,6 +53,12 @@ class Particle : public ParticleData { //! \param type Particle type 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 // //! initializes a particle from data stored in a source site. The source diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 66137c5f688..2255d1a7d88 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -430,7 +430,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]; @@ -544,11 +544,13 @@ class ParticleData : public GeometryState { int& delayed_group() { return delayed_group_; } // delayed group // Post-collision data - int& n_bank() { return n_bank_; } // number of banked fission sites - int& n_bank_second() + double& bank_second_E() { - return n_bank_second_; - } // number of secondaries banked + 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 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 65e8065fded..4dbab3213ac 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -91,17 +91,25 @@ 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(); bank.u = u; bank.E = settings::run_CE ? E : g(); bank.time = time(); + bank_second_E() += bank.E; +} - n_bank_second() += 1; +void Particle::split(double wgt) +{ + auto& bank = secondary_bank().emplace_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) @@ -356,7 +364,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(); @@ -417,6 +425,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() && diff --git a/src/physics.cpp b/src/physics.cpp index 34dc4ce1cbb..c324ce6b9cf 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 e967afd7995..e381b27c9dc 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/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(); } 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/src/weight_windows.cpp b/src/weight_windows.cpp index e798a2f7abd..9579f71c066 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.split(weight / n_split); } // remaining weight is applied to current particle p.wgt() = weight / n_split; diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index 79aadbef0de..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,6 +223,47 @@ def test_lower_ww_bounds_shape(): assert ww.lower_ww_bounds.shape == (2, 3, 4, 1) +def test_photon_heating(run_in_tmpdir): + water = openmc.Material() + water.add_nuclide('H1', 1.0) + water.add_nuclide('O16', 2.0) + water.set_density('g/cm3', 1.0) + + 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.delta_function(5e6) + model.settings.source = openmc.IndependentSource( + space=space, energy=energy, particle='photon') + + model.settings.run_mode = 'fixed source' + model.settings.batches = 5 + model.settings.particles = 100 + + tally = openmc.Tally() + tally.scores = ['heating'] + tally.filters = [ + openmc.ParticleFilter(['photon']), + openmc.MeshFilter(mesh) + ] + model.tallies = [tally] + + sp_file = model.run() + with openmc.StatePoint(sp_file) as sp: + tally_mean = sp.tallies[tally.id].mean + + # these values should be nearly identical + assert np.all(tally_mean >= 0) + + def test_roundtrip(run_in_tmpdir, model, wws): model.settings.weight_windows = wws @@ -249,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()