From e2921255310bd9da132ef135d945cd84189da30d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Dec 2021 12:02:38 -0500 Subject: [PATCH] Support time filters used with surface current tallies --- src/tallies/filter_time.cpp | 49 ++++++++------ src/tallies/tally.cpp | 1 + tests/unit_tests/test_time_filter.py | 97 ++++++++++++++++++++++++---- 3 files changed, 116 insertions(+), 31 deletions(-) diff --git a/src/tallies/filter_time.cpp b/src/tallies/filter_time.cpp index 496c8d2aa..225ef3ad3 100644 --- a/src/tallies/filter_time.cpp +++ b/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; + } } } diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 8e8e3126f..2d673d7bc 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -485,6 +485,7 @@ void Tally::set_scores(const vector& 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 { diff --git a/tests/unit_tests/test_time_filter.py b/tests/unit_tests/test_time_filter.py index f9628c7cb..5b09c6d5f 100644 --- a/tests/unit_tests/test_time_filter.py +++ b/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