Support time filters used with surface current tallies

src/tallies/filter_time.cpp

@@ -47,26 +47,39 @@ void TimeFilter::get_all_bins(
   if (t_end < bins_.front() || t_start >= bins_.back())
     return;
 
-  // Determine first bin containing a portion of time interval
-  int i_bin = 0;
-  while (bins_[i_bin + 1] < t_start) {
-    ++i_bin;
-  }
-
-  // Find matching bins
-  double dt_total = t_end - t_start;
-  for (; i_bin < bins_.size() - 1; ++i_bin) {
-    double t_left = std::max(t_start, bins_[i_bin]);
-    double t_right = std::min(t_end, bins_[i_bin + 1]);
-
-    // Add match with weight equal to the fraction of the time interval within
-    // the current time bin
-    double fraction = (t_right - t_left) / dt_total;
+  if (estimator == TallyEstimator::ANALOG) {
+    // -------------------------------------------------------------------------
+    // For surface tallies, find a match based on the exact time the particle
+    // crosses the surface
+    auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_end);
     match.bins_.push_back(i_bin);
-    match.weights_.push_back(fraction);
+    match.weights_.push_back(1.0);
 
-    if (t_end < bins_[i_bin + 1])
-      break;
+  } else {
+    // -------------------------------------------------------------------------
+    // For volume tallies, we have to check the start/end time of the current
+    // track and find where it overlaps with time bins and score accordingly
+
+    // Determine first bin containing a portion of time interval
+    auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_start);
+
+    // Find matching bins
+    double dt_total = t_end - t_start;
+    for (; i_bin < bins_.size() - 1; ++i_bin) {
+      double t_left = std::max(t_start, bins_[i_bin]);
+      double t_right = std::min(t_end, bins_[i_bin + 1]);
+
+      // Add match with weight equal to the fraction of the time interval within
+      // the current time bin
+      if (dt_total > 0.0) {
+        double fraction = (t_right - t_left) / dt_total;
+        match.bins_.push_back(i_bin);
+        match.weights_.push_back(fraction);
+      }
+
+      if (t_end < bins_[i_bin + 1])
+        break;
+    }
   }
 }
 

src/tallies/tally.cpp

@@ -485,6 +485,7 @@ void Tally::set_scores(const vector<std::string>& scores)
           fatal_error("Cannot tally mesh surface currents in the same tally as "
                       "normal surface currents");
         type_ = TallyType::SURFACE;
+        estimator_ = TallyEstimator::ANALOG;
       } else if (meshsurface_present) {
         type_ = TallyType::MESH_SURFACE;
       } else {

tests/unit_tests/test_time_filter.py

@@ -22,8 +22,22 @@ def test_time_filter_basics():
     assert elem.attrib['type'] == 'time'
 
 
-@pytest.fixture
-def model():
+def time(particle, distance, E):
+    """Return the time it takes a particle at a given energy to travel a certain
+    distance"""
+    if particle == 'neutron':
+        mass = 939.56542052e6  # eV/c²
+    elif particle == 'photon':
+        mass = 0.0
+
+    # Calculate speed via v = c * sqrt(1 - γ^-2)
+    inv_gamma = mass / (E + mass)
+    velocity = 2.99792458e10 * sqrt(1 - inv_gamma * inv_gamma)  # cm/s
+    return distance / velocity
+
+
+@pytest.fixture(params=['neutron', 'photon'])
+def model(request):
     # Select random sphere radius, source position, and source energy
     r = uniform(0., 2.)
     x = uniform(2., 10.)
@@ -44,20 +58,16 @@ def model():
     model.settings.run_mode = 'fixed source'
     model.settings.particles = 1000
     model.settings.batches = 20
+    particle = request.param
     model.settings.source = openmc.Source(
         space=openmc.stats.Point((x, 0., 0.)),
         angle=openmc.stats.Monodirectional([-1., 0., 0.]),
-        energy=openmc.stats.Discrete([E], [1.0])
+        energy=openmc.stats.Discrete([E], [1.0]),
+        particle=particle
     )
 
-    # Source particles will take a time of t = x/v to reach the inner sphere,
-    # where x is the distance from the source to the sphere and and v =
-    # sqrt(2E/m). Before this time, there should be no reactions,
-    MASS_NEUTRON_EV = 939.56542052e6  # eV/c²
-    C_LIGHT = 2.99792458e10  # cm/s
-    distance = x - r
-    velocity = sqrt(2*E / MASS_NEUTRON_EV) * C_LIGHT
-    t0 = distance / velocity
+    # Calculate time it will take neutrons to reach sphere
+    t0 = time(particle, x - r, E)
 
     # Create tally with time filter
     tally = openmc.Tally()
@@ -67,14 +77,75 @@ def model():
     return model
 
 
-def test_time_filter_transport(model, run_in_tmpdir):
+@pytest.fixture(params=['neutron', 'photon'])
+def model_surf(request):
+    # Select random distance and source energy
+    x = uniform(50., 100.)
+    E = uniform(0., 20.0e6)
+
+    # Create model
+    model = openmc.Model()
+    mat = openmc.Material()
+    mat.add_nuclide('Zr90', 1.0)
+    mat.set_density('g/cm3', 1.0)
+    model.materials.append(mat)
+    left = openmc.XPlane(-1., boundary_type='vacuum')
+    black_surface = openmc.XPlane(x, boundary_type='vacuum')
+    right = openmc.XPlane(x + 1)
+    void_cell = openmc.Cell(region=+left & -black_surface)
+    black_cell = openmc.Cell(region=+black_surface & -right)
+    model.geometry = openmc.Geometry([void_cell, black_cell])
+    model.settings = openmc.Settings()
+    model.settings.run_mode = 'fixed source'
+    model.settings.particles = 1000
+    model.settings.batches = 20
+    particle = request.param
+    model.settings.source = openmc.Source(
+        space=openmc.stats.Point((0., 0., 0.)),
+        angle=openmc.stats.Monodirectional([1., 0., 0.]),
+        energy=openmc.stats.Discrete([E], [1.0]),
+        particle=particle
+    )
+
+    # Calculate time it will take neutrons to reach purely-absorbing surface
+    t0 = time(particle, x, E)
+
+    # Create tally with surface and time filters
+    tally = openmc.Tally()
+    tally.filters = [
+        openmc.SurfaceFilter([black_surface]),
+        openmc.TimeFilter([0.0, t0*0.999, t0*1.001, 100.0])
+    ]
+    tally.scores = ['current']
+    model.tallies.append(tally)
+    return model
+
+
+def test_time_filter_volume(model, run_in_tmpdir):
     sp_filename = model.run()
     with openmc.StatePoint(sp_filename) as sp:
         t = sp.tallies[model.tallies[0].id]
         values = t.mean.ravel()
 
         # Before t0, the reaction rate should be zero
-        assert values[0] == 0.0
+        assert values[0] == pytest.approx(0.0)
 
         # After t0, the reaction rate should be positive
         assert values[1] > 0.0
+
+
+def test_time_filter_surface(model_surf, run_in_tmpdir):
+    sp_filename = model_surf.run()
+    with openmc.StatePoint(sp_filename) as sp:
+        t = sp.tallies[model_surf.tallies[0].id]
+        values = t.mean.ravel()
+        print(values)
+
+        # Before t0-ε, the current should be zero
+        assert values[0] == 0.0
+
+        # Between t0-ε and t0+ε, the current should be one
+        assert values[1] == 1.0
+
+        # After t0+ε, the current should be zero
+        assert values[2] == 0.0