From 24b741f2b3ce081f4c51fa15b130c50a363a8fc6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 10 Sep 2021 13:57:20 -0500 Subject: [PATCH 1/9] Initial implementation of time filter --- CMakeLists.txt | 1 + docs/source/pythonapi/base.rst | 1 + include/openmc/constants.h | 2 +- include/openmc/particle.h | 2 + include/openmc/particle_data.h | 12 +++- include/openmc/tallies/filter_time.h | 50 +++++++++++++++++ openmc/filter.py | 46 ++++++++++++++- src/particle.cpp | 28 +++++++++- src/tallies/filter.cpp | 3 + src/tallies/filter_time.cpp | 84 ++++++++++++++++++++++++++++ src/tallies/tally_scoring.cpp | 2 +- 11 files changed, 224 insertions(+), 7 deletions(-) create mode 100644 include/openmc/tallies/filter_time.h create mode 100644 src/tallies/filter_time.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index fb5da6f2a..9516971b9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -350,6 +350,7 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_sph_harm.cpp src/tallies/filter_sptl_legendre.cpp src/tallies/filter_surface.cpp + src/tallies/filter_time.cpp src/tallies/filter_universe.cpp src/tallies/filter_zernike.cpp src/tallies/tally.cpp diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 496201ad7..2dea91bb0 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -136,6 +136,7 @@ Constructing Tallies openmc.LegendreFilter openmc.SpatialLegendreFilter openmc.SphericalHarmonicsFilter + openmc.TimeFilter openmc.ZernikeFilter openmc.ZernikeRadialFilter openmc.ParticleFilter diff --git a/include/openmc/constants.h b/include/openmc/constants.h index b01238611..2596a0334 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -90,7 +90,7 @@ constexpr double FINE_STRUCTURE { constexpr double PLANCK_C { 1.2398419839593942e4}; // Planck's constant times c in eV-Angstroms constexpr double AMU {1.66053906660e-27}; // 1 amu in kg -constexpr double C_LIGHT {2.99792458e8}; // speed of light in m/s +constexpr double C_LIGHT {2.99792458e10}; // speed of light in cm/s constexpr double N_AVOGADRO {0.602214076}; // Avogadro's number in 10^24/mol constexpr double K_BOLTZMANN {8.617333262e-5}; // Boltzmann constant in eV/K diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 952c8c4e7..da2f6a61f 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -35,6 +35,8 @@ public: Particle() = default; + double speed() const; + //! 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 67dabe3bf..733ed3b81 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -227,9 +227,11 @@ private: int g_last_; //!< pre-collision energy group (MG only) // Other physical data - double wgt_ {1.0}; //!< particle weight - double mu_; //!< angle of scatter - bool alive_ {true}; //!< is particle alive? + double wgt_ {1.0}; //!< particle weight + double mu_; //!< angle of scatter + double time_ {0.0}; //!< time in [s] + double time_last_ {0.0}; //!< previous time in [s] + bool alive_ {true}; //!< is particle alive? // Other physical data Position r_last_current_; //!< coordinates of the last collision or @@ -349,6 +351,10 @@ public: double& wgt() { return wgt_; } double& mu() { return mu_; } const double& mu() const { return mu_; } + double& time() { return time_; } + const double& time() const { return time_; } + double& time_last() { return time_last_; } + const double& time_last() const { return time_last_; } bool& alive() { return alive_; } Position& r_last_current() { return r_last_current_; } diff --git a/include/openmc/tallies/filter_time.h b/include/openmc/tallies/filter_time.h new file mode 100644 index 000000000..c66481c57 --- /dev/null +++ b/include/openmc/tallies/filter_time.h @@ -0,0 +1,50 @@ +#ifndef OPENMC_TALLIES_FILTER_TIME_H +#define OPENMC_TALLIES_FILTER_TIME_H + +#include + +#include "openmc/tallies/filter.h" +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +//! Bins the incident particle time. +//============================================================================== + +class TimeFilter : public Filter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + ~TimeFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type() const override { return "time"; } + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle& p, TallyEstimator estimator, + FilterMatch& match) const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //---------------------------------------------------------------------------- + // Accessors + + const vector& bins() const { return bins_; } + void set_bins(gsl::span bins); + +protected: + //---------------------------------------------------------------------------- + // Data members + + vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_ENERGY_H diff --git a/openmc/filter.py b/openmc/filter.py index 95be1774d..699fa1aef 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -23,7 +23,7 @@ _FILTER_TYPES = ( 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance', - 'collision' + 'collision', 'time' ) _CURRENT_NAMES = ( @@ -1308,6 +1308,50 @@ class EnergyoutFilter(EnergyFilter): """ +class TimeFilter(RealFilter): + """Bins tally events based on the particle's time. + + Parameters + ---------- + values : iterable of float + A list of values for which each successive pair constitutes a range of + time in [s] for a single bin + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + values : numpy.ndarray + An array of values for which each successive pair constitutes a range of + time in [s] for a single bin + id : int + Unique identifier for the filter + bins : numpy.ndarray + An array of shape (N, 2) where each row is a pair of time in [s] + for a single filter bin + num_bins : int + The number of filter bins + + """ + units = 's' + + def get_bin_index(self, filter_bin): + # Use lower energy bound to find index for RealFilters + deltas = np.abs(self.bins[:, 1] - filter_bin[1]) / filter_bin[1] + min_delta = np.min(deltas) + if min_delta < 1e-3: + return deltas.argmin() + else: + msg = ('Unable to get the bin index for Filter since ' + f'"{filter_bin}" is not one of the bins') + raise ValueError(msg) + + def check_bins(self, bins): + super().check_bins(bins) + for v0, v1 in bins: + cv.check_greater_than('filter value', v0, 0., equality=True) + cv.check_greater_than('filter value', v1, 0., equality=True) + def _path_to_levels(path): """Convert distribcell path to list of levels diff --git a/src/particle.cpp b/src/particle.cpp index fa35f634a..501792e31 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -36,6 +36,30 @@ namespace openmc { +double Particle::speed() const +{ + // Determine mass in eV/c^2 + double mass; + switch (this->type()) { + case ParticleType::neutron: + mass = MASS_NEUTRON_EV; + break; + case ParticleType::photon: + mass = 0.0; + break; + case ParticleType::electron: + case ParticleType::positron: + mass = MASS_ELECTRON_EV; + break; + } + + // Calculate inverse of Lorentz factor + double inv_gamma = mass / (this->E() + mass); + + // Calculate speed via v = c * sqrt(1 - γ^-2) + return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma); +} + void Particle::create_secondary( double wgt, Direction u, double E, ParticleType type) { @@ -99,6 +123,7 @@ void Particle::event_calculate_xs() E_last() = E(); u_last() = u(); r_last() = r(); + time_last() = time(); // Reset event variables event() = TallyEvent::KILL; @@ -170,10 +195,11 @@ void Particle::event_advance() // Select smaller of the two distances double distance = std::min(boundary().distance, collision_distance()); - // Advance particle + // Advance particle in space and time for (int j = 0; j < n_coord(); ++j) { coord(j).r += distance * coord(j).u; } + this->time() += distance / this->speed(); // Score track-length tallies if (!model::active_tracklength_tallies.empty()) { diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index e3934d7ea..a1e5c709f 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -29,6 +29,7 @@ #include "openmc/tallies/filter_sph_harm.h" #include "openmc/tallies/filter_sptl_legendre.h" #include "openmc/tallies/filter_surface.h" +#include "openmc/tallies/filter_time.h" #include "openmc/tallies/filter_universe.h" #include "openmc/tallies/filter_zernike.h" #include "openmc/xml_interface.h" @@ -151,6 +152,8 @@ Filter* Filter::create(const std::string& type, int32_t id) return Filter::create(id); } else if (type == "sphericalharmonics") { return Filter::create(id); + } else if (type == "time") { + return Filter::create(id); } else if (type == "universe") { return Filter::create(id); } else if (type == "zernike") { diff --git a/src/tallies/filter_time.cpp b/src/tallies/filter_time.cpp new file mode 100644 index 000000000..496c8d2aa --- /dev/null +++ b/src/tallies/filter_time.cpp @@ -0,0 +1,84 @@ +#include "openmc/tallies/filter_time.h" + +#include // for min, max + +#include + +#include "openmc/search.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +//============================================================================== +// TimeFilter implementation +//============================================================================== + +void TimeFilter::from_xml(pugi::xml_node node) +{ + auto bins = get_node_array(node, "bins"); + this->set_bins(bins); +} + +void TimeFilter::set_bins(gsl::span bins) +{ + // Clear existing bins + bins_.clear(); + bins_.reserve(bins.size()); + + // Copy bins, ensuring they are valid + for (gsl::index i = 0; i < bins.size(); ++i) { + if (i > 0 && bins[i] <= bins[i - 1]) { + throw std::runtime_error {"Time bins must be monotonically increasing."}; + } + bins_.push_back(bins[i]); + } + + n_bins_ = bins_.size() - 1; +} + +void TimeFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + // Get the start/end time of the particle for this track + auto t_start = p.time_last(); + auto t_end = p.time(); + + // If time interval is entirely out of time bin range, exit + 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; + match.bins_.push_back(i_bin); + match.weights_.push_back(fraction); + + if (t_end < bins_[i_bin + 1]) + break; + } +} + +void TimeFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", bins_); +} + +std::string TimeFilter::text_label(int bin) const +{ + return fmt::format("Time [{}, {})", bins_[bin], bins_[bin + 1]); +} + +} // namespace openmc diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index a0fc79b5b..35d80d0e8 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -628,7 +628,7 @@ void score_general_ce(Particle& p, int i_tally, int start_index, score = flux; } // Score inverse velocity in units of s/cm. - score /= std::sqrt(2. * E / MASS_NEUTRON_EV) * C_LIGHT * 100.; + score /= p.speed(); break; case SCORE_SCATTER: From ef6156307976a325f67d9d26e5289fac7139a8cc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 13 Sep 2021 16:01:48 -0500 Subject: [PATCH 2/9] Ensure time is propagated to secondary particles --- docs/source/io_formats/source.rst | 4 ++-- docs/source/io_formats/statepoint.rst | 10 +++++----- include/openmc/particle_data.h | 7 ++++--- openmc/lib/core.py | 1 + openmc/source.py | 8 ++++++-- src/initialize.cpp | 25 +++++++++++++------------ src/particle.cpp | 4 ++++ src/source.cpp | 8 -------- 8 files changed, 35 insertions(+), 32 deletions(-) diff --git a/docs/source/io_formats/source.rst b/docs/source/io_formats/source.rst index 86ae8a1e5..4ce2229ed 100644 --- a/docs/source/io_formats/source.rst +++ b/docs/source/io_formats/source.rst @@ -20,7 +20,7 @@ following the same format. - **source_bank** (Compound type) -- Source bank information for each particle. The compound type has fields ``r``, ``u``, ``E``, - ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``, - which represent the position, direction, energy, weight, + ``time``, ``wgt``, ``delayed_group``, ``surf_id`` and ``particle``, + which represent the position, direction, energy, time, weight, delayed group, surface ID, and particle type (0=neutron, 1=photon, 2=electron, 3=positron), respectively. diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 46725fa8b..591dc2777 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -52,11 +52,11 @@ The current version of the statepoint file format is 17.0. sum-of-squares for each global tally. - **source_bank** (Compound type) -- Source bank information for each particle. The compound type has fields ``r``, ``u``, ``E``, - ``wgt``, ``delayed_group``, ``surf_id``, and ``particle``, which - represent the position, direction, energy, weight, delayed group, - surface ID, and particle type (0=neutron, 1=photon, 2=electron, - 3=positron), respectively. - Only present when `run_mode` is 'eigenvalue'. + ``time``, ``wgt``, ``delayed_group``, ``surf_id``, and + ``particle``, which represent the position, direction, energy, + time, weight, delayed group, surface ID, and particle type + (0=neutron, 1=photon, 2=electron, 3=positron), respectively. Only + present when `run_mode` is 'eigenvalue'. **/tallies/** diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 733ed3b81..06fccc092 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -44,9 +44,10 @@ struct SourceSite { Position r; Direction u; double E; - double wgt; - int delayed_group; - int surf_id; + double time {0.0}; + double wgt {1.0}; + int delayed_group {0}; + int surf_id {0}; ParticleType particle; int64_t parent_id; int64_t progeny_id; diff --git a/openmc/lib/core.py b/openmc/lib/core.py index e277d527a..e585444fb 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -16,6 +16,7 @@ class _SourceSite(Structure): _fields_ = [('r', c_double*3), ('u', c_double*3), ('E', c_double), + ('time', c_double), ('wgt', c_double), ('delayed_group', c_int), ('surf_id', c_int), diff --git a/openmc/source.py b/openmc/source.py index 415aacd96..8f25bcab0 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -253,6 +253,8 @@ class SourceParticle: Directional cosines E : float Energy of particle in [eV] + time : float + Time of particle in [s] wgt : float Weight of the particle delayed_group : int @@ -263,11 +265,12 @@ class SourceParticle: Type of the particle """ - def __init__(self, r=(0., 0., 0.), u=(0., 0., 1.), E=1.0e6, wgt=1.0, + def __init__(self, r=(0., 0., 0.), u=(0., 0., 1.), E=1.0e6, time=0.0, wgt=1.0, delayed_group=0, surf_id=0, particle=ParticleType.NEUTRON): self.r = tuple(r) self.u = tuple(u) self.E = float(E) + self.time = float(time) self.wgt = float(wgt) self.delayed_group = delayed_group self.surf_id = surf_id @@ -282,7 +285,7 @@ class SourceParticle: Source particle attributes """ - return (self.r, self.u, self.E, self.wgt, + return (self.r, self.u, self.E, self.time, self.wgt, self.delayed_group, self.surf_id, self.particle.value) @@ -309,6 +312,7 @@ def write_source_file(source_particles, filename, **kwargs): ('r', pos_dtype), ('u', pos_dtype), ('E', '= 0; --i) { + MPI_Get_address(&b.time, &disp[3]); + MPI_Get_address(&b.wgt, &disp[4]); + MPI_Get_address(&b.delayed_group, &disp[5]); + MPI_Get_address(&b.surf_id, &disp[6]); + MPI_Get_address(&b.particle, &disp[7]); + MPI_Get_address(&b.parent_id, &disp[8]); + MPI_Get_address(&b.progeny_id, &disp[9]); + for (int i = 9; i >= 0; --i) { disp[i] -= disp[0]; } - int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1}; - MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT, - MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; - MPI_Type_create_struct(9, blocks, disp, types, &mpi::source_site); + int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1}; + MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, + MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; + MPI_Type_create_struct(10, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); } #endif // OPENMC_MPI diff --git a/src/particle.cpp b/src/particle.cpp index 501792e31..7397c318f 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -77,6 +77,7 @@ void Particle::create_secondary( bank.r = r(); bank.u = u; bank.E = settings::run_CE ? E : g(); + bank.time = time(); n_bank_second() += 1; } @@ -111,6 +112,8 @@ void Particle::from_source(const SourceSite* src) E() = data::mg.energy_bin_avg_[g()]; } E_last() = E(); + time() = src->time; + time_last() = src->time; } void Particle::event_calculate_xs() @@ -399,6 +402,7 @@ void Particle::cross_surface() site.r = r(); site.u = u(); site.E = E(); + site.time = time(); site.wgt = wgt(); site.delayed_group = delayed_group(); site.surf_id = surf->id_; diff --git a/src/source.cpp b/src/source.cpp index 9c897dcf3..be4270b51 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -147,9 +147,6 @@ SourceSite IndependentSource::sample(uint64_t* seed) const { SourceSite site; - // Set weight to one by default - site.wgt = 1.0; - // Repeat sampling source location until a good site has been found bool found = false; int n_reject = 0; @@ -227,11 +224,6 @@ SourceSite IndependentSource::sample(uint64_t* seed) const break; } - // Set delayed group - site.delayed_group = 0; - // Set surface ID - site.surf_id = 0; - return site; } From 4831456144e778449c5fd54bfbeeff84caf6d7e5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 14 Sep 2021 06:45:41 -0500 Subject: [PATCH 3/9] Ensure time is propagated for fission neutrons --- src/physics.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/physics.cpp b/src/physics.cpp index 3a52396e1..36161f3ea 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -194,6 +194,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) SourceSite site; site.r = p.r(); site.particle = ParticleType::neutron; + site.time = p.time(); site.wgt = 1. / weight; site.parent_id = p.id(); site.progeny_id = p.n_progeny()++; From 83686f8fdbf7e922cb1049b8f2100f2eb982c275 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 14 Sep 2021 11:29:25 -0500 Subject: [PATCH 4/9] Add a simple test for time filter --- tests/unit_tests/test_time_filter.py | 80 ++++++++++++++++++++++++++++ 1 file changed, 80 insertions(+) create mode 100644 tests/unit_tests/test_time_filter.py diff --git a/tests/unit_tests/test_time_filter.py b/tests/unit_tests/test_time_filter.py new file mode 100644 index 000000000..f9628c7cb --- /dev/null +++ b/tests/unit_tests/test_time_filter.py @@ -0,0 +1,80 @@ +from math import sqrt +from random import uniform + +import numpy as np +import openmc +import pytest + + +def test_time_filter_basics(): + f = openmc.TimeFilter([1.0, 2.0, 5.0, 9.0]) + np.testing.assert_allclose(f.bins[0], [1.0, 2.0]) + np.testing.assert_allclose(f.bins[1], [2.0, 5.0]) + np.testing.assert_allclose(f.bins[-1], [5.0, 9.0]) + assert len(f.bins) == 3 + + # Make sure __repr__ works + repr(f) + + # to_xml_element() + elem = f.to_xml_element() + assert elem.tag == 'filter' + assert elem.attrib['type'] == 'time' + + +@pytest.fixture +def model(): + # Select random sphere radius, source position, and source energy + r = uniform(0., 2.) + x = uniform(2., 10.) + 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) + inner_sphere = openmc.Sphere(r=r) + outer_sphere = openmc.Sphere(r=10.0, boundary_type='vacuum') + center_cell = openmc.Cell(fill=mat, region=-inner_sphere) + outer_void = openmc.Cell(region=+inner_sphere & -outer_sphere) + model.geometry = openmc.Geometry([center_cell, outer_void]) + model.settings = openmc.Settings() + model.settings.run_mode = 'fixed source' + model.settings.particles = 1000 + model.settings.batches = 20 + 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]) + ) + + # 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 + + # Create tally with time filter + tally = openmc.Tally() + tally.filters = [openmc.TimeFilter([0.0, t0, 2*t0])] + tally.scores = ['total'] + model.tallies.append(tally) + return model + + +def test_time_filter_transport(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 + + # After t0, the reaction rate should be positive + assert values[1] > 0.0 From e2921255310bd9da132ef135d945cd84189da30d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Dec 2021 12:02:38 -0500 Subject: [PATCH 5/9] 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 From cad7f6812929093cc9663e62af27e56135266d99 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Dec 2021 13:34:07 -0500 Subject: [PATCH 6/9] Update test results that depend on inverse-velocity score --- .../mgxs_library_condense/results_true.dat | 8 ++++---- .../mgxs_library_distribcell/results_true.dat | 2 +- .../mgxs_library_hdf5/results_true.dat | 4 ++-- .../mgxs_library_mesh/results_true.dat | 8 ++++---- .../mgxs_library_no_nuclides/results_true.dat | 12 ++++++------ .../mgxs_library_nuclides/results_true.dat | 2 +- .../surface_source/surface_source_true.h5 | Bin 98088 -> 106088 bytes .../regression_tests/tallies/results_true.dat | 2 +- 8 files changed, 19 insertions(+), 19 deletions(-) diff --git a/tests/regression_tests/mgxs_library_condense/results_true.dat b/tests/regression_tests/mgxs_library_condense/results_true.dat index 349d9debe..508984e54 100644 --- a/tests/regression_tests/mgxs_library_condense/results_true.dat +++ b/tests/regression_tests/mgxs_library_condense/results_true.dat @@ -162,10 +162,10 @@ 3 2 2 1 1 total 1.0 0.150417 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 4.800872e-07 2.910230e-08 -2 1 2 1 1 total 5.007997e-07 3.613803e-08 -1 2 1 1 1 total 4.872635e-07 4.324235e-08 -3 2 2 1 1 total 4.890451e-07 2.833017e-08 +0 1 1 1 1 total 4.800874e-07 2.910227e-08 +2 1 2 1 1 total 5.008000e-07 3.613802e-08 +1 2 1 1 1 total 4.872638e-07 4.324238e-08 +3 2 2 1 1 total 4.890453e-07 2.833017e-08 mesh 1 group in nuclide mean std. dev. x y z 0 1 1 1 1 total 0.024194 0.001614 diff --git a/tests/regression_tests/mgxs_library_distribcell/results_true.dat b/tests/regression_tests/mgxs_library_distribcell/results_true.dat index 9fff8479d..e6e4b99e4 100644 --- a/tests/regression_tests/mgxs_library_distribcell/results_true.dat +++ b/tests/regression_tests/mgxs_library_distribcell/results_true.dat @@ -49,7 +49,7 @@ sum(distribcell) group out nuclide mean std. dev. 0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),) 1 total 1.0 0.090571 sum(distribcell) group in nuclide mean std. dev. -0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),) 1 total 4.896403e-07 2.047455e-08 +0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),) 1 total 4.896406e-07 2.047455e-08 sum(distribcell) group in nuclide mean std. dev. 0 ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),) 1 total 0.088451 0.003512 sum(distribcell) group in group out nuclide mean std. dev. diff --git a/tests/regression_tests/mgxs_library_hdf5/results_true.dat b/tests/regression_tests/mgxs_library_hdf5/results_true.dat index f49ae97ab..6da73ec9b 100644 --- a/tests/regression_tests/mgxs_library_hdf5/results_true.dat +++ b/tests/regression_tests/mgxs_library_hdf5/results_true.dat @@ -94,8 +94,8 @@ domain=1 type=chi-prompt [1.00000000e+00 0.00000000e+00] [8.07424611e-02 0.00000000e+00] domain=1 type=inverse-velocity -[5.93306434e-08 2.99151552e-06] -[3.31163681e-09 2.75470068e-07] +[5.93309775e-08 2.99151502e-06] +[3.31163250e-09 2.75469881e-07] domain=1 type=prompt-nu-fission [6.07803967e-03 1.32325631e-01] [1.12229103e-04 1.39042100e-02] diff --git a/tests/regression_tests/mgxs_library_mesh/results_true.dat b/tests/regression_tests/mgxs_library_mesh/results_true.dat index 0acc2f407..a8775984f 100644 --- a/tests/regression_tests/mgxs_library_mesh/results_true.dat +++ b/tests/regression_tests/mgxs_library_mesh/results_true.dat @@ -162,10 +162,10 @@ 3 2 2 1 1 total 1.0 0.052287 mesh 1 group in nuclide mean std. dev. x y z -0 1 1 1 1 total 8.523138e-10 3.481466e-11 -2 1 2 1 1 total 9.044552e-10 2.851430e-11 -1 2 1 1 1 total 8.643697e-10 4.046801e-11 -3 2 2 1 1 total 8.667214e-10 1.753703e-11 +0 1 1 1 1 total 8.529639e-10 3.484230e-11 +2 1 2 1 1 total 9.050800e-10 2.850752e-11 +1 2 1 1 1 total 8.650057e-10 4.049297e-11 +3 2 2 1 1 total 8.673543e-10 1.753754e-11 mesh 1 group in nuclide mean std. dev. x y z 0 1 1 1 1 total 0.030356 0.001982 diff --git a/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat b/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat index 5afe3b24c..32ae7f4cf 100644 --- a/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat +++ b/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat @@ -138,8 +138,8 @@ chi-prompt 0 1 2 total 0.0 0.000000 inverse-velocity material group in nuclide mean std. dev. -1 1 1 total 6.148626e-08 3.647428e-09 -0 1 2 total 2.844242e-06 3.098231e-07 +1 1 1 total 6.148661e-08 3.647428e-09 +0 1 2 total 2.844242e-06 3.098230e-07 prompt-nu-fission material group in nuclide mean std. dev. 1 1 1 total 0.019154 0.000623 @@ -440,8 +440,8 @@ chi-prompt 0 2 2 total 0.0 0.0 inverse-velocity material group in nuclide mean std. dev. -1 2 1 total 6.397171e-08 3.117086e-09 -0 2 2 total 2.875630e-06 2.751730e-07 +1 2 1 total 6.397205e-08 3.117085e-09 +0 2 2 total 2.875630e-06 2.751728e-07 prompt-nu-fission material group in nuclide mean std. dev. 1 2 1 total 0.0 0.0 @@ -742,8 +742,8 @@ chi-prompt 0 3 2 total 0.0 0.0 inverse-velocity material group in nuclide mean std. dev. -1 3 1 total 6.447396e-08 2.866174e-09 -0 3 2 total 3.006811e-06 1.808751e-07 +1 3 1 total 6.447429e-08 2.866173e-09 +0 3 2 total 3.006810e-06 1.808750e-07 prompt-nu-fission material group in nuclide mean std. dev. 1 3 1 total 0.0 0.0 diff --git a/tests/regression_tests/mgxs_library_nuclides/results_true.dat b/tests/regression_tests/mgxs_library_nuclides/results_true.dat index fc005b133..c0084547b 100644 --- a/tests/regression_tests/mgxs_library_nuclides/results_true.dat +++ b/tests/regression_tests/mgxs_library_nuclides/results_true.dat @@ -1 +1 @@ -b8706c9586aeeb5ee558829c28772f7f6e18a0d3fe4e30a2059b6e74564d1ebbbd0ab5473965c4e6156dff8b8178ce76c74d8db1d63f032c3d4da5f7eb9eae2c \ No newline at end of file +d36f8abb1131212063470d622dbedeae31602da71006389421bd5dba712c94fe96acbc8ded833cb237687fb1071c1a6c3e0ec67b18ce7cf0f7720451b86d993a \ No newline at end of file diff --git a/tests/regression_tests/surface_source/surface_source_true.h5 b/tests/regression_tests/surface_source/surface_source_true.h5 index 22edb898b57d6a598db85b405b63cf974bd3857c..2ea48b2ac16fa4acb1275d48ab3ef28c2136af37 100644 GIT binary patch delta 21352 zcmYLRcU;Zi{};)-(4cWME6KHIdPgO*Ym^x`kuopJ9=DRrTr1vO~YA>jpFLEUB_B2t6?TPaioS# zbu)iugn4J>mzuU_5~Z0%XR|^0Z);_^g{?|I%gk1}#G+Rd|0J_&|Np;w)wrtvLzr1^ ze#D6q3yGOh!SP(p`lWVLyz2QW!#R;s&)F&!wRx^~Rz&Z2^J_D(wyjcQ$#e62&l=c& z@n2tMN^56jrKOEhjN{#ooou$UT^*tnD(r1oD=) zN-=cQc`j@^CcCZ)ltV`yj!!6kw|sMIH4`WTXaMC|^{8cY@76SdA$;*%ja8gUfXE*DfokZLAr+$*{qj}$_kexg4#)vJ_@Xo}yv;Lf6RF$|(Z%_D_X0h# zWbC}QK3oa#WZaIP!9AV*lrEET)Nj1Gi6<|DMW768H0g5g^yV${WTHPhXgxMh7GjD4Sg07+w%u&rg*bO2HsU$`w(BP%w*E8;5jBQ7jd||I?u43WcM84)$boXf>J;x6nao4d zMo5hYJZI~9;A+-%79V*7cnHTcg2x}|ea4M>s0nmz$aDMh+g_S_pQW3zA*=z<$sdBY z)RQuvTH=Z0@s8UcH2lr7Uk>RK9Pg2y(Bq{U6Qss*1*ETf-(%J@X*)mFlJPLHr~w*E z6P|mP-R;7a>MZ5d5N`l!u-n*}Bu=OvoHw}e+>_B3J?|f`!4!3cqORyNirxoaxHC#4 zKrtw14xaEnX(*$VyGXcer-Pr^>V$~7z(q}YE_~@JkHt~UnbM|+gB#Bs9G7j!yP5f^ zO2$wenxXE2de59<-p`pal)Ax4h+b0tt({pIbB7wDWnel_+-T#9$=_LOfvE;IUh&%P zaUtH!q`46`>YDRh$%nI(Cp7)zyV(O34jK2Gw0_jOnvg4yp#!^hi3j>z^*{o+h?Qaw zp6fY3f4c8*qlvxT1CD9Nb9~(ZMZ-C)EAj5|Hb7C@m0LR#gj#F{mECzRwM*us{Rf#M zYEVjA@LZjm6^Xa)TAQ>}w}2hc?mv00Hf>W_ShdZ;6Vg{6+LSS!D>q5kL%Ie$|1H~d zCL++@)CnrxQJ|1M>|$EC{!TxpoU-$ROh3| zmOb_Nxr%9T`>%N#MMeX5A z9Jk$Fs9D#GMM96`Itcnw@W)b;2n8{<3JDGzq` z;v8lt{0jZd;$iFzMLY7Gm%84BSxFSNj^SLZI!`~wykiuiPEZgOex>KO8X{dxrbs)& zyIpv0&VehFr=$zH;t5KCuFFCO+vN(C-331D#B)=#tM7eqk0lTfsbYYdwjSpD!+E6U;M07>(ykVAh|2jfX8Sa=&>xf*v4d%9!$laJXf&C>gKSqLOXy- z27ndX!i&pp3w^RHBJRy|?_wGr=}--0v}y`_Lp&7CX*fM-b8SmF(hKbcLi*LR_TKF% zWDtZHx?`_wzA*3FS42B^HYTtS)PayC>#rN$s9D?AL$U^%w+@lbczBih&Dd4KeaT4A z8d!_PRO$tO5HhJzYFMk`kV>lX-gZh$UoY<2WbVN)71Gq-7o)BlsKPHdcrb`fNl|x1 z3t(sciXC6xvuQ5vCgBF3D()XNqrTr}5BRT}jZ)^zb15T4)rP!cp3?flejlD2KcGj_ z{6t|Q`A9h9f??xV9k0)(Q6>IePMPFp8IJLUzl~3 zbvPpl0w`MP`yll-3!WMv2cB1sHT$~$G6SXHX@K;G`c z+qZ%sO8CmY0HORzzxk29>qF0yUg4ZB>QAE~$o{JAkoJaOd-GgOTJetJXRw*JP3Jyn z7=3td;GvD}C&e?r@=#jVkLN!4ythiVW!+NS5BdRY&@H;}5CCX$p&!PaB4W5O2A0HpqTWC%L!4u$nv(oAJYF2p-`7{73_2;>$)5nJE9ob4M zf^3*Fb4s+*%VLSEn~Al4Rtv1n-?{CuFpkb z8HUFs1Htd0q&zBENFAud4dS`B9Tuep#!^C{!G{b&n1fIp`|_MNw8El5b#EleI3BR& z#Og8qSjREqxDLntpKe_?zdehWC=l@)?8WUpwB21?uYu{fb}-^K3?{}DXFKxD8QNix zJ{YsLVtu%J?^ia7LI$I)Lj0L!AA;-{6Ks}bX5Vki#;Lx|R-tR193 zD3r$_!~vA6_3fJF%x8|&gHi`6)=#W9xu0N$i-v%52+ujD#$7xY=xXAr2PhuSb9EQy zuDiQ{vKN&QJsj!`#eHFWT#{c$mZ<8Xpv3VT@pXSRScsL9T5U9LK!%Y#7iye$x>IAe zpYfv*M(}hk3OJkbT}b{BPzm06Z*cd0ZbtdC9AyGQ3Q+nF{-_&US7;!^aU6m!PmQbD zzCNom9?~^XNq){OU&65&VH_&qcFxN=l0BC-Bq_Wr8G+cgIHY-H2bV7yAy#t1*!3FR z9?K0Asy`Spkn&ug{pFn=X;_j=rEo4NM~{p9Z{;q*xuDd5-`a?YUo0C78$4Np<58&& zH5&aAswo6uD9W}aHnR#k}I!Mp) z-}s-lP`EselyDxI>(2Pk7j~PG&=HbXr&f0>T*Y!v3&j-(fqnAaDPvk<7|~)@G!l*t z;ko2tacj3V`(YZXdg!Pa!*d-TpSMN4xtR49{m~^HR|Pnt)~z2+_1jp{XhR_clFR1W zO}sH$m>!TUfsnI~bu8RNz;DTDvC;rl8u(>>7SClffrnHBj=C-HnzO1t%S9eiIS857 z_37V=-Kc%D0vlrw(A*t%66KYBb)PvSX`u0sNktE8@`BXRe5r4DBIRW zlW-h@K6N^1)nYA6AN?c==l|$!*W6*Ca+!=0#IY_SryZqUlkZ?b)bRObUDlDteO){D1HuyG+d>VwKxeKwNPe!Yrir{sg+x6Ev7JDOf24r zhrxLpYndqJGXPG<;;Ej<_q!w%7mjNpc`j*&xW1D;+eVF%upN|nJGkHn&zVE?pfrHz z^r(KHwM*En;AcW5m^s9}MgF8?EIuL#Qh;*M&y9{t4zc{v1C-3ct#m>~pU}&!nX2d5 zDf>ryaaU`Yzg6FFY|>K==^S|8>XEu3eKX4+DO_iON)2bX_~%y8l;#hn`%E-atHzRp zu1ts;>Tt7AhFwmWYa-afEu94uAvvK!vsS!{NjA=bGiLLg=lj`<7ycCb=4{;5W_oeW znw)LXPR))RiRWA_n{T$%sHR}fkj{lVIG%CBE_?D%HjsGml!3CqIknTsQ$lz_sat@Z z%D;BVneS~`&=$f&IG$_&U{1r0xOwEx#!Uf=AD=Sk$C#7ssAx8hf?vYV7e4XbnPkyC z2mxi-V*7#>vFs^Q3lkL!d2U7PH;XRJ!aU1eh`nGQ&%JhC`XeL&`9sYE^!sCuWJ{; z5unWS70tV`Naz4^#5$JeHe19VdNPNWtPQBgSacN#IxLRASCPOff`@btJcpD!ryU-_ z@>~j@66k4VobV&(Fk7^z)VNKr%o(50`%{$W3Z3k-+d@LwoQN2Hx zHu=w}R&__Kgh4XMxwYlERop;fZ-X3JJT&qw7`gE%%V=#p0tq!&o&DIha1b~q&aJ-E1lS4%rj#l3vSiW)^OI|H3RK)SzD~J2zzgoj^RA# z{FUruh!H|GE70S6U2I*vn6)1xgoq(IU=~0Bw2?x+9Ft`kV)E>2KiAJJaiz;307o~y zj96aQ-PUjORv5AZ?R_QBef+8}m0V;gDqRUb;Pm-i+p&Fb;|@qGzi2ta0?J$WeEht% z%;8#4$`m->Hol)oxYcSE0C9ZY>xBEYGT6mi3*ic=ICk#+mmh`*!xTJq;2CnS>{^p< zLN84KWdhy-PnmcA)P)TfJ%OunLFD)DSELEM{A#!n(*47inccq2o(!arF2?c86-}Fd z9K|*mIgaZ&-l0_TfX-Jm2N6)Rkh7D3V9S2WCV7AS| zGpv(o*Fwh)c(Z@{J*7)Ai?bfb6`Of(hj-7!QhUZzzZu6NXlmWa#F&S|5&;isHlh{P z?Hj-R3cFSqH$pmiUYxvS`h=ASXe4S*D~Q42pPKv~y(&vTg*i@hF{nq+L|%bjW@)6-Mc z^FuQ9hKjydeb*0Li0>=g#YzPPAD+JTT-GzzI`t520J~pJK4r-IGLXkf11J;i>t5?7 zyeaWJ;KALP{R`&q2;9X4t9K(*knzdzc~<7ujx2uLKnZPTlI@FhM)P z6Dqm)Xk)eV0&6lxsARzL`)@)MmqZFZiu6oD<5FCA^b25BC{2Nbz_Z|b$G1J&3wLe^ zl7UiwL3L+5r6~=<>=eXsAJ6T-7<{RDyrW6FC>4}AZn>epW8G@VAZj1kJCTu)J|m6a zn(>dNl^#45dw9<2OdWq&vGAb02ht(P?aiJBWi74zHt!KDHG6SqZfO4R)+eU7aW8y< zle4ETTrg`5>qv5(6sO{&R!TQ(Q80E?Duw7&Y#pTHi9Uzl{$a8?SD0l14o)Y2N!WHq)@$S9^Cf{54Zang&7p@!T-wRP^{v6t~Mq6yp8> z&e=*)sJ0RM^?oRR058xvPi{WFE)@0wNIr=DFS?(AU@zH9SBa)Y-H0L5FL^x=Za5b0Xqbz zgQe1V?Xp`83#=TcWnlR{y77xU;YgO0G}PE3L^{Y}m82hA%k(%ZgWi#6bYTPk3QGa> z7K8GdZ2yGUn^-x?L8&>6Uv^|X4#?ZiQr~zOzBmGK%+x4%Pxh9lhl*l=9uu;D+=*a1 z$^pt?%vou*mjM;H`=q4fG6a6}GR^0>x(ORI_(?K&?%4S0Zs|W*GOIHH9z#xL_1E4F z7rvn!gBzh^Nz6FwM1N0o5Cruwf(jWI(=!7HjAYp(O-E!h(FPu0wCpw!{-${ml8GS# z6*nGkp#Io{B?=EY8VFa^DsX(>72S{eW@QG7<|vid%~hT53-xgnJdg9-`sJLr@5`E~ z53oEAmLUB4AnVx4!z@>LXsLnf{qp+nuSgU!0H6#SHDCB^?d}OeC(i_Uisv%LdoRo? zWKPqbLa-qxWBUE<*G_CMiB7=hC-4)4Mt#~h5k7Y40qVfhZiUyp0XZyWqN5Tn>W;$o zX%2N2+-IVXq7cEa>-opUZdO=uQ;wtIgIDI^2`N4S!mb5g98^hOe(L#~8feZ~Geb^c zd}QHC{ORNBEwpsyLTN;2!DO8N`!s)Lz(1h|PeMIV-mo!WdSNEIcGcw!lsd4C9nwR7 z--d0%qAY|k8!53hsH0aY%RN2^9y?912z}bk7|I5$@ia2vB+oh2T=rsdJ>epLQo?2Q ztTw^(FdGd_Oor4a&%FbbuY@Z`HUvQG z@CM@+Cs;G@tLf+&d_a2Q(D^_MHY=oO;ArqW_^zR}*t2e}&XI7V z)(`gkJe~0~!Zh7^+@KF!Y-v50T`xozP%LM$0)LeE&8^S!OMVtUgpU5zmcDB+n%!mf z5F&%*QN@;dPlN?ms|8Q!xM@u2-Q}0qw*e6y2Txn`l5KuzLZt(gfpXU}-#e>jx|rsK z7S?c3DSwldb=S+N1OM~z&IMG9{eX?z4=^3&&{1*@FZN!AIWxzg<55#jK8L~prJLK` zgTsfh+|+|o27dJdtSsY4vVNdFkGNdoxpwtBo$XkYRXU%CIOk$XDZX^?QYb3|V=gG6 z_=c%lKfUV(x22?`sv$@N>0?@db)D9mg~50ZieJQqDN3Pg-cxw{yolpasb!N@&Ph+0 z*ZEut_rxl!)H_(%3m`-cbskQ-^KjK3p}SpzI(eA(P39&PE<~Z0q^m zF0X|E38;L}yH$36Vd19=cf_w$lJn4Z^08Y_d(@-mITkHxKBQd6jsDTjW7ca}9?LJo zJy(%D;dQbqX?LYvNP890f<7}c4)ngbj`ch}^wB}Q^kwyPyC(^i4^RP~)~edmtiKBj zH+bqGyLNB`lK~TQ&p!K%>rOYCc4l3bz%`a{@UzD>3rw3@b z&T}u%-^-u)!QL;)4PNrT0WV#{501`^*SuQB5?y)?VFb^n@guACStXPh93p}84?IL^ zU*edT)R3+zz)Xr?KjGv#fF%XE??JN9^vu$(^_biE0`R+u*e?Ec=X)x1oAf65;piWS zrTP?afgdOp&~dbHN|WjrEvn{6(RFkws5q_Y)yg+%^{U#sGcd^w}_1^Ugi5As$fg5_BLGmA1<=blan+`VJ!a#O32u5ed|K{X?Bs- z!VO|5liG4cW4#wESow8mh~r7Y^8&ng33-9zIvCg@a>SxtuUXlM0E+Km7g~5LG^RQ; zP<{ulfJzhZW;`pkVBwNNB^i#7``t7#D3vWATKLBRK|u*$7o_dOHH&(?_bvDW$6X6k zj}QJ#i(DzbY2881Ly*O`m%qM^rMhmj6rr;6Ab#e-6DA5STU_Rs3K5W3-g_z~(H zFL_@tB#5aa0?0kY?d*BM@rX<|8l(@w6URkUF7>x)z}f{5L9!yAn=6kSoX~~69;^;64nIeN8A~y;NU_lMoz|9_HoJy&7vF+4nGBA*6ZDueBCc zux=xVI^svDn9lO2#bT&aQUHdJ5EtmUZosDG*<$7l5twT9c+T&;ZAmksG>v*l0%c)F z{U7?8Y#eG~jRcf)8>;qszYx+HlsZ^5C;5Be_!|-$h6w-^)D)*m07P>*Gh++xnGh83HU|150zi&fUHxXIC7s;C5Ec#ipV(of5MX=Li;q8CsJl>7bice%UMc(6j_ ziw3C-;!|TDFPgnd$d$)-%GFQ2xbCtIg}22lQ}s`9{sn#<>>SvlUtiWtc!d5 zP{$BdWch0S$5yP1=pjf3lXAqll1ER`pk2P<&iEW%`xVcbH?13V_>pj1FXd=|Y%aG;9W>%HdRezCiO z9-sv3)E+ddw?!ks77I~;{w@-0B9wDQQ>veE-WEkh+4lsD#= z$v;KP~TQ!)7LT8bK+Bobqkr#B&{)Kje@jgF(N_!VPZ1 zjaB;|G5f@G7e;jGv;1rAs;LR4 zSi+#r`aPfS?Z6^o+Lu1T8t8Lt>4g7oT@c3VXNdj?50%{+2BxrK_to$@?(!P*R*15w7v`;D#yp?@xRy4zs_cawj9zy z>HgikT6mH$YeA`lo{G=zFGGg1=!&3{0(w4K8Fch~kkCL$&u=KQwkukC-$rgxPfz}a z?*AR%1s~=oS?)$8sRrcV!4odJWwGJ)s$gN#m*MAnFZYaFtVLO1W;3fm&Hasf6wgA3d>l5 z5;WqG{Nss()&qa~s#4n05?6_K@B{BN#ZLC^&Dobc=?{>?CHK{XjyzwCDscIRb^s$Q z=GD$}c+TuGeuLw{uy>I4z{eAr6uA+3@DrQj-s+Yo*9e{cCwv5bUQE zJjRdpS?%Y+yqM26#!Y;U>DxB^O-Hz@m`$N(JwHqUwC>6^?F?H_m? zmYhDO^@_Ot-B(pIgvQZdo2{G~zT?c8nbY|rl*HWJ8e(Rn;%WemulzV-f0?PsRni(} zP~acBYnJ$}vM&pd@t^Hh8A$tQESq-2U0jt|TeYC%pupvGR=GokrvNDwFu-#~!HIo( z?`Fq&=q~}!Z%z;vG`f=i*{BR=U^=Krr>MnZ>(5&ZJ;S9RU$JG%gF0<3ep{I{xVn# zud8$xfZcj%?^~Q@mWXNrHjom3GM7DH+0=BoR`DcJM@)7cANwuif~mW!L}Z|$7T`O_ zx9_XD+_Yh+5K4$7-tFJj)Y>*M^=OsWf|%3+e9C-}`>0z2le&O4V2*P|c1m?Ee$ za%K|FGHmt1QIV$UrjlBMP;3Rl;3gXez4?#US}Qn+BS{awd(|!4VCn)YDM^x$-vZqaeV8&Z* zN4)C-9xE;R@tS=#RLSeY6MByn&x1qi3mg{{?jY*{^OE}WaZp%RVQ~Lggk)V>2ZZhKCGfHI8e`CD0 z#9LwwSaaHWWlNRmyNyb14Om0YJMa5ALXF=MR`qaIWo^hb*aBbqdW3G%Di&;>EYuP6 zI(I`{c4wQ$h)QHjT*x0%w`yieYtzMDWh82u17Pd+WBnV3n|2?S)&a1hKHzp{FWv@~ znEJFzUmq}6kKnpHt@1lDCenJun-qWTRx&BCJN`KwwJ24kHEBrBwu!NubLte!2A&L( z*aOZoj(OQ_F2*0q&{>siPfjF*ZW)@!MqWqH##B`m*#g#)(9f^kw=^oHpF)3gRyo%P zEUq&D%dk$L2mJ9>h1()mB@IAXqr>2*VK-TZs2f0@*b%Vyaf0jGzxXEs7`I}Tj?C`% z!YQ@w1JkuyC2|C%fy`cPCTkT!A6zJ2yggt}1lX&5@Mgz0%t=xaU@|}##=-$d0A^v6BZB$zFosM`b>JB;;F`YFakuzY96ucM{ zzGnmbNxDi(anX?LbcOl8=}VX-sZP9UZ&AseQM?6~Z<{wautu*Z@5xBH_^v2=}h_945rA)ecHJ2Jr_=~B+)hj zIp+$PKYc`hb~Dpi>Izsv%-zci$L)B^oUN~7-Wc%ybhY~|J7&0=Oi_@9kDFxA9xfIN zjNGIn#XpU-|2w6plV4IJ>c6DK&d6UqC5MHPPD;o~iD#+Vf1TL^tJ0DZ5*NT>nPJy9 zy%h$Z3t$Z?@o3^X=Vn(}Gc=MCoIBv2iH)!C&c|Xyb2-GF>}~=0_XSz^ArjVDv@Jj_ zZVtFk?Y_VE#+ddPmApCeChy%nD870)j7=zGGr&4hFmFercNL4j`{MObtTMO(F_bRx z-L%GZlT`6Us3Dbo_K$WSF`uQKkyJL2JJo)Dyn0)+B;^TK5R;udtZ!~f5WMV0)TG3I$+p{`!|Y!Js^nrMIeemdy%C}x|QwAGRPH~V8MoYIZ>SH)CP&aDv-4Uvy< z*xle{=Sp-sJV>Z2J4e~>m-cd)oB{8oV$q%{? zfIr*aEPl;CQ>a880ds8uUmUZq^QSb{*Q9L$%gDxB%eNHtpkMP>&CgcgZDB@qu=OnEZj!&ih<%fDdLgsnpeVDC0Nr5@VPm229(K}4~?Aw zD>?!m-s_%A)_g{-C+o!oZ+#YVsC9^|DQDyaOFRKz+c)*>OQ9jFJqd0PIN3TabjmQH z3z9uD@@!Jvq3f{^ggWd1m?OALxG~;l30rui1S`4#KI!#+Qm>yZsq|d{Yl!z%zN~R< zt}ro(w}i~`I@7Y+z;mqssL32LdG_G5pB>h)4>BrwJHRq>p}1G4Ybw0G(+VP1DZBt_ zvIfswQ_fn7-V2mmSHK@c@2k6=VY{2OD>3i}?0I=|=K=-oFy3IGAzuv0{u6j@8S0+K zdZi~aK}=@MdvEV?a%d;!IWogQaA4G^!=FA2?SMR|>kPQ@Q12^GyPJNWp%M`_N8(%( z<~bb0Z~m()Q5E71<}xCeUa8x3AO1lwb$YQ%MhbVS(YjYeDC3|dg*l@9?@GfFUEEmb zkai|DD6T7Z&E8p-2BmZVV$gO1Uqv^-HDd-eZBjyi4*-7!K~>ZZIOhZ0_o2n2Q?HoT zQXk^t3wT?b?ltK`+HQjZ-KRtrg$_l(UCDrC*1Tc(Yu&DBI*We4G(xs ze9(c;g=||i^2CIAD?3W_M(i?u#!|_BNhF2KxihzQRW76EiCWSfu;}HO3|_}_N!^`b z;vJg(-M-&75o78?WW%O#zWXm$sf92fYtRHPM%i6CWELq@g)nVs-xd0 z7GSBy)mc^K3oRtXL|ko=+UyN$3Tk4)5fhuoKU3=@2)kP^;!Tn1_}#@dD9ZFyq!N*v zWCXV}%8vS67VcFPFd6CmZ%_J&H!keXrtJaPKmofxe&V{g-a^+R?`g=X|9Te39&60@ z1|z{-U%*>8*eOEs!vt#YD5Bngjy;wP8xqVY^<<%#;NI6fi*w)MC+9SAZjk*tD1BF&G`9u3vj;7Z)6<@`hU^BV*l;h)BBFHJ? zI`ZDTwzuoKG-G{CRFw!4z<-8z{G2>YxEx79E$$CEH`8j`k?|}g<^74tK)}3 zTDaB>0;~xDyjn43`G<0r?nd&7f~+*O`Q9w6A=@?dRaO$U!{1YLgHnX`sLF#BnZ6Uw zw$iIuZ`bw(6OKfdMfA&#m7!aps&Q{1%86O?by1B4Q;4sQO!4eBZ|5KOOMR7ybk>mf z(z-kyt9r%68LPZU;k7$mP>^HCs*5LB5(xP4sJ~^gVWy8dDs>=W@gTt6Y`d1f!58|f zzp0>#9)#&1KPLR;OZMAPm7aVd8A5zNejYSfBLr&*X)yrs%?Y7F)u)<%MyxUp04yWx z>g;Wi^>8So)(#{lMD5XF+jf~3#F;K<15jS#!N8-(CT;Ba3N4lX3x9{OJWJ{T!+qJ~}0_k&&lrkFAyH@yqmX+>2i|>}de0Ck^ObQp z$s^knBcCq5R)=MzJ7BN#(9;eD89Q2 zJW|#;vO47{>>A1jeempulPv|21S^68zdiPK+zk552l(Sds-j>dL(MwLiC<2#wVWSG zN)X|_X9fFRKMFf8`NALt{3oK4o6Y{&gNhFUEF;eu_e6JiG*U>05r8>@ZPXX-6BF^C zN8>JJ1l*w@CXz#)tS@b5C90Q#T0&k@9rX9Aah~nwYVwkS-0?bW!=7oQ*s8%3tQiTo z^jyOxpVBRn*CVMN66G=1W$Cel*z%ws21*@8`^>%{o&)*_*I0tZp@4njZ*84x&oWOQ zN_;~A&q=@5t>}^9xDf2+P3)(J-(SztLOPn*js@Wut55U$P7rE2jM$QK-cM`U+Z3=y zBd;=!*m4hN1ZOldeY{br$#I%7fPcxnjmfcCmhl%)uu+TxbU$*=kpX+yTab1X6qqyh z+HC9p8<50OHX4LC)^!+VpTrE%l3pC0ZxHISE8EkmYDt0DN5Sh7a+qYJYviKa!hov^ zbQs`~lg~TV{>Qu{Co?qT05@FPbnN3^%sS&Zz=jCG#WkLlHy_3NJUyTAX(Tr4HnKmke0P%0bkVt(ew0G5y@{+BKe zxW+zftJI{4m^5*|w^Hrdm@$z@5H&f}Gq7{*&udvHP!kgc@t$q5zK8ug;l514){zVG zQEGGwmi;X7wyUZV2esmw#dQ( zkLi`_m_WZor8%w*2dpD4o?Td0H*O-cTr?3dN8(1jbz3KU1@}>DhfIdfI%3kOR?zBs zm$7@BZiNxR4HUA_p+V7-0W9u3xs@aG7e}{wwVT3ZN=c?f2KeSuw<%%QXd{-9(o+EO z5q@25of)Nm3Qe z1XxG5qzwM^ec(x!+9IOnrULdoU@qTPjYUy96|jVuEO5%YJONlt%&9knQakt^D-Voi!libd>Hk>7Z{x65&1*2}0dW5T@&o*_8QE2eigb zeHLP(AdR!G zw^gY%1MpX?$h}*d3QsIEh?P(U``HLeqn0i z&A&n;mXl88GB?lPeeF);L6iPgn5uF%OfXRNrbn!m{aC~9>^!-Wiv@gsieIM;l@O0u zl1a?BHTj`Dc7>(BC>pSg;KE0iwl}QIOqalHIq_Zu_~(`BrK20L#g2~ywV3=cqp4Tk z@8c|ba$=%c2>9)d;z279II$3>dq8Rsx$W>KMQZ z3e`#1ivH_Yu^mTGUXhWrT1=0Bej=QCMN5%VkcOEbv)VMmKZmI5#`;Bo#RU7s70Xa^5CXXpUs1jy#y}9pzMZN~qf^9t4YS z9$6P+CG>xCEWR+n`0hOuoyIe2>0-bJa#iW@Rr!?}!m9uosUh2o6IOkhGl7jbBf&ai za%RP^{_~o##v+O*YErPb^hn3wADHdNDg`O6t|q0WrTMa3pNLAsKo0y^^HNB`PF5N` G-2Fe%6Y_%q delta 12178 zcmYLPdAw8A@i#oB5tZBrA_<$}xm|D}2+GCvg7~<&UKBClc0tjAiUC1`>lGCPD#jdn zH{jA}x7gM;YHed{+i16F7o**R+BS9#T{O~agI4{0=gyqZ&p-Uk%=ylD=FH4F=f215 zyfD1-$vy>hn0fotQ#bZSyl0>4iw%odYT4d1RCgaEVynpT8IgJT|B=cYZUng%5kvXI z=T6$|#>ig(|KF&+qM<^zlV7;ad$v8?i70PhOEsf$RomB4W|U<=YN*V9EZ6rlR6{Vy z5r*ms#@*K_3M36kX#~T&fLQypT-e`GO<+>X6Kcx-F;xg%O;yIaJ0YRAj2)?_jhG#V z&%>{3SV47DB%${0k0~PtP|Q@(1DMP^z)(5ZDA#7pQ1O12v;Bs043^7=q1rJkD4zq8 z*+xud#km3J%E%B@gl|jCsgbCY@a(i=>eP{(ObPbHdJ8@Vbyp3Wo!^Hj3>A_*{`Iz}E#IA7iNGXe(%5 z1^x}AW2#3%v03*ZLuJQV#gh{1<#939jI-HH7%DoBWp5nzBywsz`gIW7lVj101j_~B zOJiB?!k!VHdfZTD!SsMh2_`?*C<-JEh&`TNOQ3I0VA-08tsKm9U7~8dTNw3AKIveB~dPP-`a554YO_%s>Elf~lN|Ojt&s z8U>n`i0q{X4rRjpp_teKM)(8R^&u=LVb2`Ka{VxLWHQUu$qphURAL~e@{(#Bse0m9 zI|LmQyE7S+K7eslNVZSyK6J3SIH8sth9Rbi1Ew;jAjc^f=EIr2cQ|S-ashJm2)27i zVD!vTE>1xO;1-TBj+dOur>%Y}M(-3}?w*2D%?;(^X_!v& zO3g)okyz>$nvAJ7Z?$04VT!AsjVnV4@*rp;0NG90ZIJnMaa2xYTL!Fc$e1xTdmblH zo`<77Nb7Vd%pGL*${>n2pS^PP(V5d(Zk>*UC3$&saGHd<3Ff$9vcTlUw+-LwJob&k zo_1KSI!Eb?{hcvs^eETRDei?;IQhvG-ulxomH`m?x1lklQ3rL-7{IR0FZ7 z&`IQOpNm4DYpRy;rOrhoB&G?ztk@ff<1A#q@5< zhp9g<;^sK8)-S*@FjWPCkwg)>E_=18Z{$IPfk^SsGIw5VY6Yg{a?BmJ3T|0GDEb!0c_K@!XwehV&AgCW+y~ zD+4xSsKZ`L%XQGFK3vLh?;?!;#UZ443F>tb_I#14Y&NOV#i-Av%n-ed<7O{IZ7yS* z?corlSk?|`W%$++Oa`%&%NUbehKhMCSHR!9 zoO(<(-FUfF>2g$E&~2a_S1>xc0{vOeavYePus6YOf`u~hVXq0h0bfsiYgb^R!I@fa z7O!Mz`bunW1TDQC*G!zsgeWhZAbgRdn)k zORwM?pbO3l`xjTCzkrc@;VR7jDwZ>_J6E$@z8WzU9mtW1;Hu&n-h2y0XSzO1E=`*BHExmq_hFoo4 zEBjcBI;~}|JY+M=vY*B9AbBs358G?-T_<*T4N8qDsq0Kt5L*Y>#v1lZt~H9+qtqMm z-Enl|!*a&fJJvuVZ^}9d;`AChVSY*ATn4-crV^bZB0} zPN?K2(Dn`NS-Jt82BNwiZ)A>bWkl(0VA;M2yApeJ19p2W+x=V7^;@{pYc_CYS{v{l z-8>5~2;Qrk7-!s!u5+y7jr}UIfp5z!s9ZqkKY$g9hBqoKbKk^E6_na_?bV1?#;n+~SCOW;3@qz6JXPlIqyb zqb~QdT)P)1OxViEHockc)y+6z!e{{Bg)cqd3pV2iAg~jU56*`0B=6(h)qwQGm$x6Bn&{o0KuHft{ zZblZqIl(p&t0SEDLnzNS=5U5Eu0t%x1C&$b>X1>~!RGW1RKLV>t%U7uXSur_hbds$ zhpQ`TRzP16dne#3yYS0P$o4jrQRXi_img1x za`Q2C2LA1bQ3|nHpG9XS$0WG2;Gz;W9&zY%P`xOe=EIopGOuOJ=mhf@LD>!W#v)_jz7- zzkt#pMmXre6vW;^9PbNkuYJKNeu>5Om$3T^HQH1SSx$Zt6@`r`3KcY0=p3Y%`Xbx& zz?eeqgSsU4F6#I07z|$o? zOs#!tFnmIGK83qN5Sg#wh{z61WLN`EHa^W>ou@IrRYt@i*Wg7}Fo5no&9eIpDk3aS6-|&RQS2Zu_BQZo;i&@S zJ;S~=_@09W`b1e6sLq)#Ha_U<~@f$3rzkx;xqyUI1 zDf&pg35@X^zE%WV2iB;uU#5nU{^o4;T6C;xd=p1O>`mA!9KHy#_8o>BFJLUc&2sD8 z=##*$=P}MNussJ%9GLn$^ss*ToZLV* zVA^2!zR&jB_l@F>{sZQ3`~cPZ5hI#E;s#nj zMx(^52d@Sawcn&N9Q#A|D*X@#;0NsGfuSXG3U=qmY%hbw{|U>jpP*L~CG#VcP%u&0 z>#r~-{!VD?qw8RqS*}$ z_fMJC{Tb>bvX2Pem)Y)suORjg?AEJctJUP6%Z2-M)BrZBaqm?$>t`(2er6PZ!D9Ls z=%>)N5xyaiBwR~B4?h#B&`l|HlZ%^r4WkSO<5m31qsdHO6XkiG<;?3S+pnnM_^IVJ zww2%$hmDfS{tEYzWKc&2Ir!NfYR-<=sQEMoufNVQtzTlEe?tREW!`=m+;yAy>rgySj+lh*cukjlcx5aiB_NwqikyrY+Y_GyMFFbAF>tc_$ za9g~=nEV??@lBd5T+YM-=~$j4}RX6g$*!{Gjm{Hg({Zfi0{;2UQTe54-nQw%7iOqa_p-(4{406=8aU zcY%rio$=n^QNdkw6zQrM_#1CN|2Mhc|AzYgi;?Ajp=RRMfL9q_lu2(F3i3~uEnpgw zNfN#-!KeO#lOUO!z}x>|dJ1pDpV!|X=^d0mHob(B|@}GQ6oBzagy~o)2?tZoK z-C1-7)RXVx6A>KscT8pO=5>EJdioEH)Mvt;j_UHZa$VH-U4|Htq zL>VkeLFYa;wt11zzO5F)lz%Lj_kvGi+ZtG)#|2N`Cxb!vKvXsHDSqoKJW+6_KVi=5 zCm3TyrHxkpX9T4PlLp;L#0Z-74{2Bm%s#~Eeai7%$R^9m|I`@Lp!6`f*!ze}X&-1A zy;i5SimVL{2SedgG<-Z#+x4jtJl+=xf5Z)X2ozX-!cibCg_LMC%%J)iNEkIUXb%^J z^_fV2uygOotiW9qNd$ZPBH{D|-rn#`jsVJjKM`yg5s3wc_Q@c8kf$};T6+XE2h8D> zKP*_jcO)!*&>#vwijeLyQ16U_Bd0l<#1Y40$dQdgGxrzG-yd3iKhff(bu`l67g}v! zSu-ehLHl(G*^J4WX+Y~D%kaVs7VL-c8NFDKVbRqz6(aQx5Xag9(0ZEGB`vL!sqP0Y z+b_sQKQzmLMsHkSuw|rdy*m;M4NaINt$HBU{=kSiq;p794pGhJ5tjo{Sdj98&{}#e zH40jp;!weRv@v^>XeC-pjRbjiBpT};B=*Wdu*b)UmK_5vH(IpjXlNB;pfAb3;Pugx zw?$zK34!4SM_07wlh)CMHqragcy4x%L3Q(E#j!mWi#45zL7JgsXDD`0le*)erN&eC zI{h4 zJ_K6jVAh6+ygZS%i0e7nGZ6&tA%d(BGDmsfz6quuET^wTDRd7NkaHNcnszY`h1S-o zTa%$>bUyV7(5jOW3m>8$AFTnHQl?PlO=x*d+$O5V5wez~HDB+cMSe+*%n`DFm>}ba zLu*ckd*SRSJYEa0+45}Pn zsbJR(n$F{8(VT_FB(mc(x?wt%RWH_$$6{KiP{l@u^irCZvX-zkENIp&(F(-bm?>KF zIB0Rb$1J7RLo(qW-I>rz$BEV@Wc5VRq9;L%>Ifc1aP(UFL})4PmnXk6qSKI|HM6jX z+zCR{JOPXD$)Z(Gh890bv@C^a0~sVwMAJQ_9F%9$oN4b2dAD>DsoBsv@aETQ^4!S$ z@Oce3&qc_}Y^EI|A8!tE%>kFMqqQhn{Zxq^x1kkmpu!`+%a-1`b6~fHDkxFJwDzf< z3N3$%FtktMi$!&Y;hvf=*35ih8i*biI5`(uWstQYTB_+%7^gt9<_WMc51NT|LN}k1 zYn@IT#b9CNKAa~b=bes~GQ5K3QmXYlY&@wm%7Jp|G~q~}24qufCh>RAlr`rpXttyk zlm?;Y^;-K(Xr;M=?9v9Z$Ry0GLDAz2M9VIKRy#v9gZw&r%{~K~ueBB-Gdh<#W#%~I zS9YMaQNAGQpj}OuS`uDa3u$jUtsK$kbd{S_$_zmS*z$0A?=}%Eo0x zsWp$ZI7q@0WU0XAOGWEl3axpGXeP8Ft)wr(QXZMmB+)2rxmYUDy%;O5&cAXYw8lk3 zk)$8^CAa3fV||1!1mb zwxQKfYPyaFEg7q6Q>(G5o=(7}bZl`6O4mZmUL)x=u7OqnlBNe!d$rhmSHo^;tw35g z%mt^#p)E953DU&ckXF*GuvF4j=~7j?YX#_B53Pm3e4A`ot9P^(_Szaj8Kk*dtE_<* z)yaC4te*u#cy;Z}qKR=*>8}}sC98z6PDwb|31RsD~ZLZ=LYE>!Ed!OBltbw5vLUXcn5cp4P@B zg6G%M0lgj?{)Q?@tP?ta9m+@%!xv}CjQ~||5xVHDAV|Z7kKC}}gj?jum2bgH;9E3q-h)#s+BS z9isVnKr?iJ423IezaIG&G_pg;{O#h`z8zZK7cEX}InB}}mZV1J2pQGUJ&IntP5ca6 z>)a}ueH%1iYb~l!bd#)ko1k^U8KivdEv`{nLS;8aW(QAiiX0J??-qRTZrGc5iDr^k zx>L07ozT)c_o@%g)frWg!;p4M--V?V-1rvTd?z~SZ4{Q;M)JFtG>pKL_hPMwT7_3O z9+5Y;NI&g+fy-#hI#D{Cg`>O~n!iP~))r_M?Eqth_O8n7cPNhu$6rSUvo|A3LuZhr zec1Picj+EJ219vRQw6c+3(%@MUX&uVbWKzDKx;g}lM@NXZRN?K(0_* zmi+7YW0&-K6(qJIoix-CTD>2OMqacGrIEi+wDx__l3L4=R(XimhPFW|Y!zTPe{ zg3?3KYGI9WNg3o8Rf6;sw9=!ZbsvT1VJ~4SHOks5idHDf1<@m0|{^u&lVd=al4(S4xN|ubPFj9DvZQMYr5LpxbAUh3HUmdzQKk&j+M~JyT1r Date: Mon, 27 Dec 2021 14:17:10 -0500 Subject: [PATCH 7/9] Address comments on #1934 from @ojschumann --- docs/source/io_formats/particle_restart.rst | 1 + src/particle.cpp | 4 ++- src/particle_restart.cpp | 2 ++ src/tallies/filter_time.cpp | 32 +++++++++++---------- 4 files changed, 23 insertions(+), 16 deletions(-) diff --git a/docs/source/io_formats/particle_restart.rst b/docs/source/io_formats/particle_restart.rst index 6ebf7e7ba..53f86735f 100644 --- a/docs/source/io_formats/particle_restart.rst +++ b/docs/source/io_formats/particle_restart.rst @@ -32,3 +32,4 @@ The current version of the particle restart file format is 2.0. multi-group mode. - **xyz** (*double[3]*) -- Position of the particle. - **uvw** (*double[3]*) -- Direction of the particle. + - **time** (*double*) -- Time of the particle in [s]. diff --git a/src/particle.cpp b/src/particle.cpp index 7397c318f..b3dd72d9a 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -54,7 +54,7 @@ double Particle::speed() const } // Calculate inverse of Lorentz factor - double inv_gamma = mass / (this->E() + mass); + const double inv_gamma = mass / (this->E() + mass); // Calculate speed via v = c * sqrt(1 - γ^-2) return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma); @@ -685,6 +685,7 @@ void Particle::write_restart() const write_dataset(file_id, "energy", simulation::source_bank[i - 1].E); write_dataset(file_id, "xyz", simulation::source_bank[i - 1].r); write_dataset(file_id, "uvw", simulation::source_bank[i - 1].u); + write_dataset(file_id, "time", simulation::source_bank[i - 1].time); } else if (settings::run_mode == RunMode::FIXED_SOURCE) { // re-sample using rng random number seed used to generate source particle int64_t id = (simulation::total_gen + overall_generation() - 1) * @@ -697,6 +698,7 @@ void Particle::write_restart() const write_dataset(file_id, "energy", site.E); write_dataset(file_id, "xyz", site.r); write_dataset(file_id, "uvw", site.u); + write_dataset(file_id, "time", site.time); } // Close file diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index f7159ce44..32d187dd8 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -50,6 +50,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) read_dataset(file_id, "energy", p.E()); read_dataset(file_id, "xyz", p.r()); read_dataset(file_id, "uvw", p.u()); + read_dataset(file_id, "time", p.time()); // Set energy group and average energy in multi-group mode if (!settings::run_CE) { @@ -64,6 +65,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) p.u_last() = p.u(); p.E_last() = p.E(); p.g_last() = p.g(); + p.time_last() = p.time(); // Close hdf5 file file_close(file_id); diff --git a/src/tallies/filter_time.cpp b/src/tallies/filter_time.cpp index 225ef3ad3..6d3b9ab78 100644 --- a/src/tallies/filter_time.cpp +++ b/src/tallies/filter_time.cpp @@ -1,6 +1,7 @@ #include "openmc/tallies/filter_time.h" -#include // for min, max +#include // for min, max, copy +#include // for back_inserter #include @@ -25,14 +26,13 @@ void TimeFilter::set_bins(gsl::span bins) bins_.clear(); bins_.reserve(bins.size()); - // Copy bins, ensuring they are valid - for (gsl::index i = 0; i < bins.size(); ++i) { - if (i > 0 && bins[i] <= bins[i - 1]) { - throw std::runtime_error {"Time bins must be monotonically increasing."}; - } - bins_.push_back(bins[i]); + // Ensure time bins are sorted + if (!std::is_sorted(bins.cbegin(), bins.cend(), std::less_equal())) { + throw std::runtime_error {"Time bins must be monotonically increasing."}; } + // Copy bins + std::copy(bins.cbegin(), bins.cend(), std::back_inserter(bins_)); n_bins_ = bins_.size() - 1; } @@ -40,8 +40,8 @@ void TimeFilter::get_all_bins( const Particle& p, TallyEstimator estimator, FilterMatch& match) const { // Get the start/end time of the particle for this track - auto t_start = p.time_last(); - auto t_end = p.time(); + const auto t_start = p.time_last(); + const auto t_end = p.time(); // If time interval is entirely out of time bin range, exit if (t_end < bins_.front() || t_start >= bins_.back()) @@ -51,7 +51,7 @@ void TimeFilter::get_all_bins( // ------------------------------------------------------------------------- // 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); + const auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_end); match.bins_.push_back(i_bin); match.weights_.push_back(1.0); @@ -60,6 +60,10 @@ void TimeFilter::get_all_bins( // 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 + // Skip if time interval is zero + if (t_start == t_end) + return; + // Determine first bin containing a portion of time interval auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_start); @@ -71,11 +75,9 @@ void TimeFilter::get_all_bins( // 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); - } + const 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; From c0e5ed2971401daaafe19f12e7c3a575cd0d174c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Dec 2021 11:43:44 -0500 Subject: [PATCH 8/9] Remove less_equal from call to is_sorted This causes issues if you compile with -D_GLIBCXX_DEBUG. Underlying reason is described in https://github.com/xtensor-stack/xtensor/issues/2296 --- src/tallies/filter_time.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tallies/filter_time.cpp b/src/tallies/filter_time.cpp index 6d3b9ab78..fc626873b 100644 --- a/src/tallies/filter_time.cpp +++ b/src/tallies/filter_time.cpp @@ -27,7 +27,7 @@ void TimeFilter::set_bins(gsl::span bins) bins_.reserve(bins.size()); // Ensure time bins are sorted - if (!std::is_sorted(bins.cbegin(), bins.cend(), std::less_equal())) { + if (!std::is_sorted(bins.cbegin(), bins.cend())) { throw std::runtime_error {"Time bins must be monotonically increasing."}; } From 5288ffc600722c083d9d45570feb93e296ef2b15 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Jan 2022 08:25:01 -0500 Subject: [PATCH 9/9] Use adjacent_find instead of is_sorted. Handle collision estimators correctly --- src/tallies/filter_time.cpp | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/src/tallies/filter_time.cpp b/src/tallies/filter_time.cpp index fc626873b..3f78d1817 100644 --- a/src/tallies/filter_time.cpp +++ b/src/tallies/filter_time.cpp @@ -1,7 +1,8 @@ #include "openmc/tallies/filter_time.h" -#include // for min, max, copy -#include // for back_inserter +#include // for min, max, copy, adjacent_find +#include // for greater_equal +#include // for back_inserter #include @@ -26,8 +27,9 @@ void TimeFilter::set_bins(gsl::span bins) bins_.clear(); bins_.reserve(bins.size()); - // Ensure time bins are sorted - if (!std::is_sorted(bins.cbegin(), bins.cend())) { + // Ensure time bins are sorted and don't have duplicates + if (std::adjacent_find(bins.cbegin(), bins.cend(), std::greater_equal<>()) != + bins.end()) { throw std::runtime_error {"Time bins must be monotonically increasing."}; } @@ -47,18 +49,11 @@ void TimeFilter::get_all_bins( if (t_end < bins_.front() || t_start >= bins_.back()) return; - if (estimator == TallyEstimator::ANALOG) { + if (estimator == TallyEstimator::TRACKLENGTH) { // ------------------------------------------------------------------------- - // For surface tallies, find a match based on the exact time the particle - // crosses the surface - const auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_end); - match.bins_.push_back(i_bin); - match.weights_.push_back(1.0); - - } 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 + // For tracklength estimator, we have to check the start/end time of + // the current track and find where it overlaps with time bins and score + // accordingly // Skip if time interval is zero if (t_start == t_end) @@ -82,6 +77,13 @@ void TimeFilter::get_all_bins( if (t_end < bins_[i_bin + 1]) break; } + } else { + // ------------------------------------------------------------------------- + // For collision estimator or surface tallies, find a match based on the + // exact time of the particle + const auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_end); + match.bins_.push_back(i_bin); + match.weights_.push_back(1.0); } }