diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 1016fbaa3d7..86d827b7e36 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -60,13 +60,17 @@ 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. + +This element has the following attributes/sub-elements: :weight: The weight below which particles undergo Russian roulette. @@ -99,6 +103,26 @@ you care. This element has the following attributes/sub-elements: *Default*: 0.0 + :time_neutrons + The time above which neutrons will be killed. + + *Default*: Infinity + + :time_photons + The time above which photons will be killed. + + *Default*: Infinity + + :time_electrons + The time above which electrons will be killed. + + *Default*: Infinity + + :time_positrons + The time above which positorns will be killed. + + *Default*: Infinity + ---------------------------- ```` ---------------------------- diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 376c9abf305..80a561f4d05 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -40,6 +40,10 @@ class Particle : public ParticleData { 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/particle_data.h b/include/openmc/particle_data.h index 6d86d528e70..b74c46a5a91 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -218,9 +218,8 @@ struct BoundaryInfo { * https://doi.org/10.1016/j.anucene.2017.11.032. */ class ParticleData { + public: - //---------------------------------------------------------------------------- - // Constructors ParticleData(); private: @@ -232,7 +231,7 @@ class ParticleData { vector photon_xs_; //!< Microscopic photon cross sections MacroXS macro_xs_; //!< Macroscopic cross sections - int64_t id_; //!< Unique ID + int64_t id_; //!< Unique ID ParticleType type_ {ParticleType::neutron}; //!< Particle type (n, p, e, etc.) int n_coord_ {1}; //!< number of current coordinate levels @@ -303,25 +302,23 @@ class ParticleData { // Secondary particle bank vector secondary_bank_; - int64_t current_work_; // current work index + int64_t current_work_; // current work index - vector flux_derivs_; // for derivatives for this particle + vector flux_derivs_; // for derivatives for this particle vector filter_matches_; // tally filter matches - vector tracks_; // tracks for outputting to file + vector tracks_; // tracks for outputting to file vector nu_bank_; // bank of most recently fissioned particles - vector pht_storage_; // interim pulse-height results - // Global tally accumulators double keff_tally_absorption_ {0.0}; double keff_tally_collision_ {0.0}; double keff_tally_tracklength_ {0.0}; double keff_tally_leakage_ {0.0}; - bool trace_ {false}; //!< flag to show debug information + bool trace_ {false}; //!< flag to show debug information double collision_distance_; // distance to particle's next closest collision @@ -446,7 +443,6 @@ class ParticleData { decltype(tracks_)& tracks() { return tracks_; } decltype(nu_bank_)& nu_bank() { return nu_bank_; } NuBank& nu_bank(int i) { return nu_bank_[i]; } - vector& pht_storage() { return pht_storage_; } double& keff_tally_absorption() { return keff_tally_absorption_; } double& keff_tally_collision() { return keff_tally_collision_; } diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 3ce13c89f90..4cee77d48a7 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -76,9 +76,8 @@ extern std::string weight_windows_file; //!< Location of weight window file to // std::string. Sometimes it is, but it seems libc++ may not be like that // on some computers, like the intel Mac. extern "C" const char* path_statepoint_c; //!< C pointer to statepoint file name - -extern "C" int32_t n_inactive; //!< number of inactive batches -extern "C" int32_t max_lost_particles; //!< maximum number of lost particles +extern "C" int32_t n_inactive; //!< number of inactive batches +extern "C" int32_t max_lost_particles; //!< maximum number of lost particles extern double rel_max_lost_particles; //!< maximum number of lost particles, relative to the //!< total number of particles @@ -92,6 +91,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 7937cd2e7a3..a996ca16c8a 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -50,15 +50,17 @@ class Settings: deviation. 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' + cutoff : dict + 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 @@ -658,6 +660,87 @@ def surf_source_write(self, surf_source_write: dict): self._surf_source_write = surf_source_write + @confidence_intervals.setter + def confidence_intervals(self, confidence_intervals: bool): + cv.check_type('confidence interval', confidence_intervals, bool) + self._confidence_intervals = confidence_intervals + + @electron_treatment.setter + def electron_treatment(self, electron_treatment: str): + cv.check_value('electron treatment', electron_treatment, ['led', 'ttb']) + self._electron_treatment = electron_treatment + + @photon_transport.setter + def photon_transport(self, photon_transport: bool): + cv.check_type('photon transport', photon_transport, bool) + self._photon_transport = photon_transport + + @ptables.setter + def ptables(self, ptables: bool): + cv.check_type('probability tables', ptables, bool) + self._ptables = ptables + + @seed.setter + def seed(self, seed: int): + cv.check_type('random number generator seed', seed, Integral) + cv.check_greater_than('random number generator seed', seed, 0) + self._seed = seed + + @survival_biasing.setter + def survival_biasing(self, survival_biasing: bool): + cv.check_type('survival biasing', survival_biasing, bool) + self._survival_biasing = survival_biasing + + @cutoff.setter + def cutoff(self, cutoff: dict): + if not isinstance(cutoff, Mapping): + msg = f'Unable to set cutoff from "{cutoff}" which is not a '\ + 'Python dictionary' + raise ValueError(msg) + for key in cutoff: + if key == 'weight': + cv.check_type('weight cutoff', cutoff[key], Real) + cv.check_greater_than('weight cutoff', cutoff[key], 0.0) + elif key == 'weight_avg': + cv.check_type('average survival weight', cutoff[key], Real) + cv.check_greater_than('average survival weight', + cutoff[key], 0.0) + elif key in ('energy_neutron', 'energy_photon', 'energy_electron', + 'energy_positron'): + cv.check_type('energy cutoff', cutoff[key], Real) + cv.check_greater_than('energy cutoff', cutoff[key], 0.0) + elif key in ('time_neutron', 'time_photon', 'time_electron', + 'time_positron'): + cv.check_type('time cutoff', cutoff[key], Real) + cv.check_greater_than('time cutoff', cutoff[key], 0.0) + else: + msg = f'Unable to set cutoff to "{key}" which is unsupported ' \ + 'by OpenMC' + + self._cutoff = cutoff + + @entropy_mesh.setter + def entropy_mesh(self, entropy: RegularMesh): + cv.check_type('entropy mesh', entropy, RegularMesh) + self._entropy_mesh = entropy + + @trigger_active.setter + def trigger_active(self, trigger_active: bool): + cv.check_type('trigger active', trigger_active, bool) + self._trigger_active = trigger_active + + @trigger_max_batches.setter + def trigger_max_batches(self, trigger_max_batches: int): + cv.check_type('trigger maximum batches', trigger_max_batches, Integral) + cv.check_greater_than('trigger maximum batches', trigger_max_batches, 0) + self._trigger_max_batches = trigger_max_batches + + @trigger_batch_interval.setter + def trigger_batch_interval(self, trigger_batch_interval: int): + cv.check_type('trigger batch interval', trigger_batch_interval, Integral) + cv.check_greater_than('trigger batch interval', trigger_batch_interval, 0) + self._trigger_batch_interval = trigger_batch_interval + @property def no_reduce(self) -> bool: return self._no_reduce @@ -1511,7 +1594,8 @@ def _cutoff_from_xml_element(self, root): 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 1db5c99d78e..e9490a77f13 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 59294b28c9d..f3fe20c2743 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 8dd4a12e5b8..307410a8646 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 1c039e9ba90..2e790657aa1 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -95,6 +95,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}; @@ -533,6 +534,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 00000000000..e69de29bb2d 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 00000000000..5cd65c05941 --- /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 00000000000..1cafceb772e --- /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 00000000000..8e13d57b73e --- /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 d75a50477a3..ee6d0d981dd 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -26,7 +26,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.) @@ -86,7 +88,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.]