Time Cutoff

docs/source/io_formats/settings.rst

@@ -60,13 +60,17 @@ fission.
 ``<cutoff>`` Element
 --------------------
 
-The ``<cutoff>`` 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 ``<cutoff>`` 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
+
 ----------------------------
 ``<delayed_photon_scaling>``
 ----------------------------

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

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<ElementMicroXS> 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<SourceSite> secondary_bank_;
 
-  int64_t current_work_; // current work index
+  int64_t current_work_;               // current work index
 
-  vector<double> flux_derivs_; // for derivatives for this particle
+  vector<double> flux_derivs_;         // for derivatives for this particle
 
   vector<FilterMatch> filter_matches_; // tally filter matches
 
-  vector<TrackStateHistory> tracks_; // tracks for outputting to file
+  vector<TrackStateHistory> tracks_;   // tracks for outputting to file
 
   vector<NuBank> nu_bank_; // bank of most recently fissioned particles
 
-  vector<double> 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<double>& pht_storage() { return pht_storage_; }
 
   double& keff_tally_absorption() { return keff_tally_absorption_; }
   double& keff_tally_collision() { return keff_tally_collision_; }

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<double, 4>
   energy_cutoff; //!< Energy cutoff in [eV] for each particle type
+extern array<double, 4>
+  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

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)

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 {

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;

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<int>(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()

src/settings.cpp

@@ -95,6 +95,7 @@ int64_t max_particles_in_flight {100000};
 
 ElectronTreatment electron_treatment {ElectronTreatment::TTB};
 array<double, 4> energy_cutoff {0.0, 1000.0, 0.0, 0.0};
+array<double, 4> 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

tests/regression_tests/time_cutoff/inputs_true.dat

@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+  </materials>
+  <geometry>
+    <cell id="1" material="void" region="-1" universe="1" />
+    <surface boundary="vacuum" coeffs="0.0 0.0 0.0 200" id="1" type="sphere" />
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <source strength="1.0">
+      <space type="point">
+        <parameters>0.0 0.0 0.0</parameters>
+      </space>
+      <energy type="discrete">
+        <parameters>10000.0 1.0</parameters>
+      </energy>
+    </source>
+    <cutoff>
+      <time_neutron>1e-07</time_neutron>
+    </cutoff>
+  </settings>
+  <tallies>
+    <filter id="1" type="time">
+      <bins>0.0 1e-07 2e-07</bins>
+    </filter>
+    <tally id="1">
+      <filters>1</filters>
+      <scores>flux</scores>
+    </tally>
+  </tallies>
+</model>