crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 21:25:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge pull request #1934 from paulromano/time-filter

Browse source 
Track particle time and implement filter for time
This commit is contained in:
Olaf Schumann 2022-01-05 20:24:17 +01:00 • committed by GitHub
commit d12e390f1a
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
30 changed files with 453 additions and 58 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

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

1
docs/source/io_formats/particle_restart.rst
View file

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

4
docs/source/io_formats/source.rst
View file

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

10
docs/source/io_formats/statepoint.rst
View file

@ -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/**

1
docs/source/pythonapi/base.rst
View file

@ -136,6 +136,7 @@ Constructing Tallies
openmc.LegendreFilter
openmc.SpatialLegendreFilter
openmc.SphericalHarmonicsFilter
openmc.TimeFilter
openmc.ZernikeFilter
openmc.ZernikeRadialFilter
openmc.ParticleFilter

2
include/openmc/constants.h
View file

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

2
include/openmc/particle.h
View file

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

19
include/openmc/particle_data.h
View file

@ -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;
@ -227,9 +228,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 +352,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_; }

50
include/openmc/tallies/filter_time.h Normal file
View file

@ -0,0 +1,50 @@
#ifndef OPENMC_TALLIES_FILTER_TIME_H
#define OPENMC_TALLIES_FILTER_TIME_H
#include <gsl/gsl-lite.hpp>
#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<double>& bins() const { return bins_; }
void set_bins(gsl::span<const double> bins);
protected:
//----------------------------------------------------------------------------
// Data members
vector<double> bins_;
};
} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_ENERGY_H

46
openmc/filter.py
View file

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

1
openmc/lib/core.py
View file

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

8
openmc/source.py
View file

@ -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', '<f8'),
('time', '<f8'),
('wgt', '<f8'),
('delayed_group', '<i4'),
('surf_id', '<i4'),

25
src/initialize.cpp
View file

@ -138,24 +138,25 @@ void initialize_mpi(MPI_Comm intracomm)
// Create bank datatype
SourceSite b;
MPI_Aint disp[9];
MPI_Aint disp[10];
MPI_Get_address(&b.r, &disp[0]);
MPI_Get_address(&b.u, &disp[1]);
MPI_Get_address(&b.E, &disp[2]);
MPI_Get_address(&b.wgt, &disp[3]);
MPI_Get_address(&b.delayed_group, &disp[4]);
MPI_Get_address(&b.surf_id, &disp[5]);
MPI_Get_address(&b.particle, &disp[6]);
MPI_Get_address(&b.parent_id, &disp[7]);
MPI_Get_address(&b.progeny_id, &disp[8]);
for (int i = 8; i >= 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

34
src/particle.cpp
View file

@ -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
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);
}
void Particle::create_secondary(
double wgt, Direction u, double E, ParticleType type)
{
@ -53,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;
}
@ -87,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()
@ -99,6 +126,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 +198,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()) {
@ -373,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_;
@ -655,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) *
@ -667,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

2
src/particle_restart.cpp
View file

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

1
src/physics.cpp
View file

@ -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()++;

8
src/source.cpp
View file

@ -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;
}

3
src/tallies/filter.cpp
View file

@ -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<SpatialLegendreFilter>(id);
} else if (type == "sphericalharmonics") {
return Filter::create<SphericalHarmonicsFilter>(id);
} else if (type == "time") {
return Filter::create<TimeFilter>(id);
} else if (type == "universe") {
return Filter::create<UniverseFilter>(id);
} else if (type == "zernike") {

101
src/tallies/filter_time.cpp Normal file
View file

@ -0,0 +1,101 @@
#include "openmc/tallies/filter_time.h"
#include <algorithm> // for min, max, copy, adjacent_find
#include <functional> // for greater_equal
#include <iterator> // for back_inserter
#include <fmt/core.h>
#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<double>(node, "bins");
this->set_bins(bins);
}
void TimeFilter::set_bins(gsl::span<const double> bins)
{
// Clear existing bins
bins_.clear();
bins_.reserve(bins.size());
// 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."};
}
// Copy bins
std::copy(bins.cbegin(), bins.cend(), std::back_inserter(bins_));
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
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())
return;
if (estimator == TallyEstimator::TRACKLENGTH) {
// -------------------------------------------------------------------------
// 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)
return;
// 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
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;
}
} 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);
}
}
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

1
src/tallies/tally.cpp
View file

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

2
src/tallies/tally_scoring.cpp
View file

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

8
tests/regression_tests/mgxs_library_condense/results_true.dat
View file

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

2
tests/regression_tests/mgxs_library_distribcell/results_true.dat
View file

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

4
tests/regression_tests/mgxs_library_hdf5/results_true.dat
View file

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

8
tests/regression_tests/mgxs_library_mesh/results_true.dat
View file

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

12
tests/regression_tests/mgxs_library_no_nuclides/results_true.dat
View file

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

2
tests/regression_tests/mgxs_library_nuclides/results_true.dat
View file

@ -1 +1 @@
b8706c9586aeeb5ee558829c28772f7f6e18a0d3fe4e30a2059b6e74564d1ebbbd0ab5473965c4e6156dff8b8178ce76c74d8db1d63f032c3d4da5f7eb9eae2c
d36f8abb1131212063470d622dbedeae31602da71006389421bd5dba712c94fe96acbc8ded833cb237687fb1071c1a6c3e0ec67b18ce7cf0f7720451b86d993a

BIN
tests/regression_tests/surface_source/surface_source_true.h5
View file

Binary file not shown.

2
tests/regression_tests/tallies/results_true.dat
View file

@ -1 +1 @@
22cd0d5d2764367e2a3ca2b715d0863ac7afd33cbd4f0b2707722e6267d95d3bc07026c9a7fdbc8260c69c40a6617b6e1c831de5f9c4a7cb89c4efa6b74bea76
de09325b517aad9f58940e5fcd53b005c57ea008a15699769ab91aba723de141014101e55be2621e34c369beefb6db8c6cf847e241ee26e73e394c0f155138b0

151
tests/unit_tests/test_time_filter.py Normal file
View file

@ -0,0 +1,151 @@
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'
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.)
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
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]),
particle=particle
)
# Calculate time it will take neutrons to reach sphere
t0 = time(particle, x - r, E)
# 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
@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] == 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