From 03f98e8f8234294b85406b20b0534623b032e0d9 Mon Sep 17 00:00:00 2001 From: Christopher Fichtlscherer <29277544+cfichtlscherer@users.noreply.github.com> Date: Sat, 19 Aug 2023 22:20:33 +0200 Subject: [PATCH] time-cutoff version 2.0 (#2631) --- docs/source/io_formats/settings.rst | 36 ++++++++++++---- include/openmc/particle.h | 4 ++ include/openmc/settings.h | 2 + openmc/settings.py | 17 ++++---- src/event.cpp | 2 + src/finalize.cpp | 1 + src/particle.cpp | 26 ++++++++++++ src/settings.cpp | 13 ++++++ .../regression_tests/time_cutoff/__init__.py | 0 .../time_cutoff/inputs_true.dat | 34 +++++++++++++++ .../time_cutoff/results_true.dat | 5 +++ tests/regression_tests/time_cutoff/test.py | 42 +++++++++++++++++++ tests/unit_tests/test_settings.py | 8 +++- 13 files changed, 174 insertions(+), 16 deletions(-) create mode 100644 tests/regression_tests/time_cutoff/__init__.py create mode 100644 tests/regression_tests/time_cutoff/inputs_true.dat create mode 100644 tests/regression_tests/time_cutoff/results_true.dat create mode 100755 tests/regression_tests/time_cutoff/test.py diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 36027ef2ff..a52de1b904 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -60,13 +60,15 @@ fission. ```` Element -------------------- -The ```` element indicates two kinds of cutoffs. The first is the weight -cutoff used below which particles undergo Russian roulette. Surviving particles -are assigned a user-determined weight. Note that weight cutoffs and Russian -rouletting are not turned on by default. The second is the energy cutoff which -is used to kill particles under certain energy. The energy cutoff should not be -used unless you know particles under the energy are of no importance to results -you care. This element has the following attributes/sub-elements: +The ```` element indicates three kinds of cutoffs. The first is the +weight cutoff used below which particles undergo Russian roulette. Surviving +particles are assigned a user-determined weight. Note that weight cutoffs and +Russian rouletting are not turned on by default. The second is the energy cutoff +which is used to kill particles under certain energy. The energy cutoff should +not be used unless you know particles under the energy are of no importance to +results you care. The third is the time cutoff used to kill particles whose time +exceeds a specific cutoff. Particles will be killed exactly at the specified +time. :weight: The weight below which particles undergo Russian roulette. @@ -99,6 +101,26 @@ you care. This element has the following attributes/sub-elements: *Default*: 0.0 + :time_neutron + The time above which neutrons will be killed. + + *Default*: Infinity + + :time_photon + The time above which photons will be killed. + + *Default*: Infinity + + :time_electron + The time above which electrons will be killed. + + *Default*: Infinity + + :time_positron + The time above which positorns will be killed. + + *Default*: Infinity + ---------------------------- ```` ---------------------------- diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 376c9abf30..80a561f4d0 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -40,6 +40,10 @@ public: double speed() const; + //! moves the particle by the distance length to its next location + //! \param length the distance the particle is moved + void move_distance(double length); + //! create a secondary particle // //! stores the current phase space attributes of the particle in the diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 3ce13c89f9..9f7882643a 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -92,6 +92,8 @@ extern ElectronTreatment electron_treatment; //!< how to treat secondary electrons extern array energy_cutoff; //!< Energy cutoff in [eV] for each particle type +extern array + time_cutoff; //!< Time cutoff in [s] for each particle type extern int legendre_to_tabular_points; //!< number of points to convert Legendres extern int max_order; //!< Maximum Legendre order for multigroup data diff --git a/openmc/settings.py b/openmc/settings.py index c4be8699fb..ef20d54dbb 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -51,14 +51,16 @@ class Settings: create_fission_neutrons : bool Indicate whether fission neutrons should be created or not. cutoff : dict - Dictionary defining weight cutoff and energy cutoff. The dictionary may - have six keys, 'weight', 'weight_avg', 'energy_neutron', 'energy_photon', - 'energy_electron', and 'energy_positron'. Value for 'weight' + Dictionary defining weight cutoff, energy cutoff and time cutoff. The + dictionary may have ten keys, 'weight', 'weight_avg', 'energy_neutron', + 'energy_photon', 'energy_electron', 'energy_positron', 'time_neutron', + 'time_photon', 'time_electron', and 'time_positron'. Value for 'weight' should be a float indicating weight cutoff below which particle undergo Russian roulette. Value for 'weight_avg' should be a float indicating - weight assigned to particles that are not killed after Russian - roulette. Value of energy should be a float indicating energy in eV - below which particle type will be killed. + weight assigned to particles that are not killed after Russian roulette. + Value of energy should be a float indicating energy in eV below which + particle type will be killed. Value of time should be a float in + seconds. Particles will be killed exactly at the specified time. delayed_photon_scaling : bool Indicate whether to scale the fission photon yield by (EGP + EGD)/EGP where EGP is the energy release of prompt photons and EGD is the energy @@ -1534,7 +1536,8 @@ class Settings: if elem is not None: self.cutoff = {} for key in ('energy_neutron', 'energy_photon', 'energy_electron', - 'energy_positron', 'weight', 'weight_avg'): + 'energy_positron', 'weight', 'weight_avg', 'time_neutron', + 'time_photon', 'time_electron', 'time_positron'): value = get_text(elem, key) if value is not None: self.cutoff[key] = float(value) diff --git a/src/event.cpp b/src/event.cpp index 1db5c99d78..e9490a77f1 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -112,6 +112,8 @@ void process_advance_particle_events() int64_t buffer_idx = simulation::advance_particle_queue[i].idx; Particle& p = simulation::particles[buffer_idx]; p.event_advance(); + if (!p.alive()) + continue; if (p.collision_distance() > p.boundary().distance) { simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx}); } else { diff --git a/src/finalize.cpp b/src/finalize.cpp index 59294b28c9..f3fe20c274 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -78,6 +78,7 @@ int openmc_finalize() settings::electron_treatment = ElectronTreatment::LED; settings::delayed_photon_scaling = true; settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0}; + settings::time_cutoff = {INFTY, INFTY, INFTY, INFTY}; settings::entropy_on = false; settings::event_based = false; settings::gen_per_batch = 1; diff --git a/src/particle.cpp b/src/particle.cpp index 200bcc6d6d..a40ea7024c 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -65,6 +65,13 @@ double Particle::speed() const return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma); } +void Particle::move_distance(double length) +{ + for (int j = 0; j < n_coord(); ++j) { + coord(j).r += length * coord(j).u; + } +} + void Particle::create_secondary( double wgt, Direction u, double E, ParticleType type) { @@ -203,11 +210,25 @@ void Particle::event_advance() double distance = std::min(boundary().distance, collision_distance()); // Advance particle in space and time + // Short-term solution until the surface source is revised and we can use + // this->move_distance(distance) for (int j = 0; j < n_coord(); ++j) { coord(j).r += distance * coord(j).u; } this->time() += distance / this->speed(); + // Kill particle if its time exceeds the cutoff + bool hit_time_boundary = false; + double time_cutoff = settings::time_cutoff[static_cast(type())]; + if (time() > time_cutoff) { + double dt = time() - time_cutoff; + time() = time_cutoff; + + double push_back_distance = speed() * dt; + this->move_distance(-push_back_distance); + hit_time_boundary = true; + } + // Score track-length tallies if (!model::active_tracklength_tallies.empty()) { score_tracklength_tally(*this, distance); @@ -223,6 +244,11 @@ void Particle::event_advance() if (!model::active_tallies.empty()) { score_track_derivative(*this, distance); } + + // Set particle weight to zero if it hit the time boundary + if (hit_time_boundary) { + wgt() = 0.0; + } } void Particle::event_cross_surface() diff --git a/src/settings.cpp b/src/settings.cpp index ace2887a6e..ac56320e8e 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -96,6 +96,7 @@ int64_t max_particles_in_flight {100000}; ElectronTreatment electron_treatment {ElectronTreatment::TTB}; array energy_cutoff {0.0, 1000.0, 0.0, 0.0}; +array time_cutoff {INFTY, INFTY, INFTY, INFTY}; int legendre_to_tabular_points {C_NONE}; int max_order {0}; int n_log_bins {8000}; @@ -540,6 +541,18 @@ void read_settings_xml(pugi::xml_node root) energy_cutoff[3] = std::stod(get_node_value(node_cutoff, "energy_positron")); } + if (check_for_node(node_cutoff, "time_neutron")) { + time_cutoff[0] = std::stod(get_node_value(node_cutoff, "time_neutron")); + } + if (check_for_node(node_cutoff, "time_photon")) { + time_cutoff[1] = std::stod(get_node_value(node_cutoff, "time_photon")); + } + if (check_for_node(node_cutoff, "time_electron")) { + time_cutoff[2] = std::stod(get_node_value(node_cutoff, "time_electron")); + } + if (check_for_node(node_cutoff, "time_positron")) { + time_cutoff[3] = std::stod(get_node_value(node_cutoff, "time_positron")); + } } // Particle trace diff --git a/tests/regression_tests/time_cutoff/__init__.py b/tests/regression_tests/time_cutoff/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/time_cutoff/inputs_true.dat b/tests/regression_tests/time_cutoff/inputs_true.dat new file mode 100644 index 0000000000..7d3e94f50f --- /dev/null +++ b/tests/regression_tests/time_cutoff/inputs_true.dat @@ -0,0 +1,34 @@ + + + + + + + + + + fixed source + 100 + 10 + + + 0.0 0.0 0.0 + + + 10000.0 1.0 + + + + 1e-07 + + + + + 0.0 1e-07 2e-07 + + + 1 + flux + + + diff --git a/tests/regression_tests/time_cutoff/results_true.dat b/tests/regression_tests/time_cutoff/results_true.dat new file mode 100644 index 0000000000..1cafceb772 --- /dev/null +++ b/tests/regression_tests/time_cutoff/results_true.dat @@ -0,0 +1,5 @@ +tally 1: +2.000000E+03 +4.000000E+05 +0.000000E+00 +0.000000E+00 diff --git a/tests/regression_tests/time_cutoff/test.py b/tests/regression_tests/time_cutoff/test.py new file mode 100755 index 0000000000..8e13d57b73 --- /dev/null +++ b/tests/regression_tests/time_cutoff/test.py @@ -0,0 +1,42 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def time_model(): + model = openmc.Model() + time_cutoff = 1e-7 + + # A single sphere + s1 = openmc.Sphere(r=200, boundary_type='vacuum') + sphere = openmc.Cell() + sphere.region = -s1 + model.geometry = openmc.Geometry([sphere]) + + # Set the running parameters + settings_file = openmc.Settings() + settings_file.run_mode = 'fixed source' + settings_file.batches = 10 + settings_file.particles = 100 + settings_file.cutoff = {'time_neutron': time_cutoff} + settings_file.source = openmc.source.Source(space=openmc.stats.Point(), + energy=openmc.stats.Discrete([1e4], [1])) + model.settings = settings_file + + # Tally flux under time cutoff + tallies = openmc.Tallies() + tally = openmc.Tally() + tally.scores = ['flux'] + time_filter = openmc.TimeFilter([0, time_cutoff, 2*time_cutoff]) + tally.filters = [time_filter] + tallies.append(tally) + model.tallies = tallies + + return model + + +def test_time_cutoff(time_model): + harness = PyAPITestHarness('statepoint.10.h5', time_model) + harness.main() diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index 79f2038b21..4b2cc18240 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -27,7 +27,9 @@ def test_export_to_xml(run_in_tmpdir): s.survival_biasing = True s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5, 'energy_photon': 1000.0, 'energy_electron': 1.0e-5, - 'energy_positron': 1.0e-5} + 'energy_positron': 1.0e-5, 'time_neutron': 1.0e-5, + 'time_photon': 1.0e-5, 'time_electron': 1.0e-5, + 'time_positron': 1.0e-5} mesh = openmc.RegularMesh() mesh.lower_left = (-10., -10., -10.) mesh.upper_right = (10., 10., 10.) @@ -88,7 +90,9 @@ def test_export_to_xml(run_in_tmpdir): assert s.survival_biasing assert s.cutoff == {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5, 'energy_photon': 1000.0, - 'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5} + 'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5, + 'time_neutron': 1.0e-5, 'time_photon': 1.0e-5, + 'time_electron': 1.0e-5, 'time_positron': 1.0e-5} assert isinstance(s.entropy_mesh, openmc.RegularMesh) assert s.entropy_mesh.lower_left == [-10., -10., -10.] assert s.entropy_mesh.upper_right == [10., 10., 10.]