openmc-dev
TimeFilter reducing flux results
When I try to obtain cell tally flux results on a simple shell of Iron for a time-dependent neutron source, I notice that applying a time filter reduces the flux results. It doesn't matter how large a time bin I use, I always seem to get a reduction in flux results. At the moment I simply have:
time_filter = openmc.TimeFilter(values=[0, 1E+100])
This should easily cover the range of times in the simulation, but when I try:
results = openmc.StatePoint("statepoint.10.h5")
tally = results.get_tally(name='flux_tally', scores=['flux'])
print("flux : ", tally.mean)
I get lower results than without a time filter (~25 vs ~19). My aim is to benchmark a new method for obtaining shutdown dose rates by generating decay photons at certain times, and I noticed that applying a time filter returns 0 for photon flux (versus ~1E-12 for without a time filter). I tried to use OpenMC's prompt photon generator to see if the problem was with the modified code, but even then I got 0 with a time filter. Then I tried simply neutrons, as above.
EDIT: I suspect it may be to do with the 'tracklength' tally estimator, as when I set the estimator to 'analog' or 'collision', the results were unchanged when tallied with/without a time filter.
The full code I used is as follows:
import openmc
import numpy as np
# remove old xml files
!rm *.xml
iron=openmc.Material(name='iron')
iron.set_density('g/cm3', 7.3)
iron.add_nuclide('Fe56', 1.0, percent_type='wo')
mats=openmc.Materials([iron])
sphere_inner_surface=openmc.Sphere(x0=0, y0=0, z0=0, r=30.0)
sphere_outer_surface=openmc.Sphere(x0=0, y0=0, z0=0, r=39.0)
sphere_boundary_surface=openmc.Sphere(x0=0, y0=0, z0=0, r=40.0, boundary_type='vacuum')
inner_sphere_region=-sphere_inner_surface
inner_sphere_cell=openmc.Cell(region=inner_sphere_region)
sphere_shell_region=+sphere_inner_surface & -sphere_outer_surface #& ~inner_sphere_region
sphere_shell_cell=openmc.Cell(region=sphere_shell_region)
sphere_shell_cell.fill=iron
boundary_region=+sphere_outer_surface & -sphere_boundary_surface
boundary_cell=openmc.Cell(region=boundary_region)
universe = openmc.Universe(cells=[inner_sphere_cell, sphere_shell_cell, boundary_cell])
geom = openmc.Geometry(universe)
# creates a 14MeV point source
source = openmc.Source()
source.space = openmc.stats.Point((0, 0, 0))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([14e6], [1])
time_shakes=[0.00000E+00,\
4.22664E+14,\
8.17023E+14,\
1.18497E+15,\
1.52828E+15,\
1.84860E+15,\
2.14747E+15,\
2.42632E+15,\
2.68650E+15,\
2.92926E+15,\
3.15576E+15]
time_seconds=[x*1e-8 for x in time_shakes]
print("time_seconds : ", time_seconds)
source.time=openmc.stats.Discrete(time_seconds, [0.000000000,\
0.133934017,\
0.124964856,\
0.116596334,\
0.108788226,\
0.101503004,\
0.094705652,\
0.088363497,\
0.082446058,\
0.076924892,\
0.071773463])
sett = openmc.Settings()
batches = 10
sett.batches = batches
sett.inactive=0
sett.max_lost_particles=1000
sett.particles = int(1e6)
sett.particle = "neutron"
sett.run_mode = 'fixed source'
sett.photon_transport = True # This line is required to switch on photons tracking
sett.source = source
##############################################################################################
energy_bins = [0, 10000, 15000, 20000, 30000, 40000, 50000, 60000, 70000, 80000, 100000,\
150000, 200000, 300000, 400000, 500000, 600000, 800000, 1000000, 2000000,\
4000000, 6000000, 8000000, 10000000]
time_bins = [0.00000E+00, 1.0E+6,\
4.32664E+6,\
8.27023E+6,\
1.28497E+7,\
1.62828E+7,\
1.94860E+7,\
2.24747E+7,\
2.52632E+7,\
2.78650E+7,\
3.02926E+7,\
3.25576E+7]
time_bins_photon = [0.00000E+0, 1E+100]
energy_filter = openmc.EnergyFilter(energy_bins)
photon_particle_filter = openmc.ParticleFilter(["photon"])
#surface_filter = openmc.SurfaceFilter(sphere_outer_surface)
cell_filter = openmc.CellFilter(sphere_shell_cell)
time_filter = openmc.TimeFilter(values=time_bins_photon)
tallies = openmc.Tallies()
dose_tally = openmc.Tally(name="flux_tally")
dose_tally.scores = ["flux"]
dose_tally.filters = [
cell_filter,
# photon_particle_filter,
# energy_filter,
time_filter
]
tallies.append(dose_tally)
##################################################################################################
######## generate xml files ########
geom.export_to_xml()
mats.export_to_xml()
sett.export_to_xml()
tallies.export_to_xml()
# remove old h5 and xml files
!rm *.h5
# run openmc
openmc.run()
results = openmc.StatePoint("statepoint.10.h5")
tally = results.get_tally(name='flux_tally', scores=['flux'])
print("flux : ", tally.mean)
I have also noticed some underprediction with the time filter when we have been trying to compare with experimental results and MCNP. We have seen better results with surface tallies results than track length estimators, but still some possible issues. I am going to do some digging into the code a bit more.
In src/tallies/filter_time.cpp, the function "get_all_bins" seems to skip particle tracks where t_start == t_end, if the estimator is equal to "TRACKLENGTH":
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;
I'm not sure why the tracklength estimator should just skip tallying a particle track if the interval is zero. Couldn't it just tally in whichever bin that particle fits. It seems that the problem is in dividing by the difference between t_start and t_end to determine weight adjustments for a given 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
const double fraction = (t_right - t_left) / dt_total;
EDIT: This appears to directly affect photon flux tallies as t_start does equal t_end, at least that's what it appears to be if I try to output these values at run-time.
EDIT 2: This also seems to affect neutrons as occasionally t_start-t_end = 0 for neutrons.
I think changing the contents of "get_all_bins" in src/tallies/filter_time.cpp to the following might work:
// 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
// Determine first bin containing a portion of time interval
auto i_bin = lower_bound_index(bins_.begin(), bins_.end(), t_start);
if (t_start != t_end) {
// 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 {
match.bins_.push_back(i_bin);
match.weights_.push_back(1.0);
}
} 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);
}