time-cutoff version 2.0 (#2631)

docs/source/io_formats/settings.rst

@@ -60,13 +60,15 @@ 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.
 
   :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
+
 ----------------------------
 ``<delayed_photon_scaling>``
 ----------------------------

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

include/openmc/settings.h

@@ -92,6 +92,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

@@ -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)

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

@@ -96,6 +96,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};
@@ -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

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 particle="neutron" strength="1.0" type="independent">
+      <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>

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

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

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.]