mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
Merge pull request #1462 from paulromano/ttb-heating-fix
Consistent treatment of energy deposition for photon-electron calculations
This commit is contained in:
Amanda Lund
GitHub
commit
fbbbe8f95b
11 changed files with 177 additions and 267 deletions
|
|
@ -22,9 +22,9 @@ from .function import Tabulated1D
|
|||
|
||||
|
||||
# Constants
|
||||
MASS_ELECTRON_EV = 0.5109989461e6 # Electron mass energy
|
||||
PLANCK_C = 1.2398419739062977e4 # Planck's constant times c in eV-Angstroms
|
||||
FINE_STRUCTURE = 137.035999139 # Inverse fine structure constant
|
||||
MASS_ELECTRON_EV = 0.5109989461e6 # Electron mass energy
|
||||
PLANCK_C = 1.2398419739062977e4 # Planck's constant times c in eV-Angstroms
|
||||
FINE_STRUCTURE = 137.035999139 # Inverse fine structure constant
|
||||
CM_PER_ANGSTROM = 1.0e-8
|
||||
# classical electron radius in cm
|
||||
R0 = CM_PER_ANGSTROM * PLANCK_C / (2.0 * pi * FINE_STRUCTURE * MASS_ELECTRON_EV)
|
||||
|
|
@ -33,7 +33,7 @@ R0 = CM_PER_ANGSTROM * PLANCK_C / (2.0 * pi * FINE_STRUCTURE * MASS_ELECTRON_EV)
|
|||
_SUBSHELLS = [None, 'K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5',
|
||||
'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2', 'O3',
|
||||
'O4', 'O5', 'O6', 'O7', 'O8', 'O9', 'P1', 'P2', 'P3', 'P4',
|
||||
'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11','Q1', 'Q2', 'Q3']
|
||||
'P5', 'P6', 'P7', 'P8', 'P9', 'P10', 'P11', 'Q1', 'Q2', 'Q3']
|
||||
|
||||
_REACTION_NAME = {
|
||||
501: ('Total photon interaction', 'total'),
|
||||
|
|
@ -391,49 +391,6 @@ class AtomicRelaxation(EqualityMixin):
|
|||
_SUBSHELLS, range(len(_SUBSHELLS)))
|
||||
group.create_dataset('transitions', data=df.values.astype(float))
|
||||
|
||||
def energy_fluorescence(self, shell):
|
||||
"""Compute expected energy of fluorescent photons for the shell
|
||||
|
||||
Parameters
|
||||
----------
|
||||
shell : str
|
||||
The subshell to compute
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Energy of fluorescent photons
|
||||
|
||||
"""
|
||||
|
||||
if shell not in self.binding_energy:
|
||||
raise KeyError('Invalid shell {}.'.format(shell))
|
||||
|
||||
if shell in self._e_fluorescence:
|
||||
# Already computed
|
||||
return self._e_fluorescence[shell]
|
||||
e = 0.0
|
||||
if shell not in self.transitions or self.transitions[shell].empty:
|
||||
e = self.binding_energy[shell]
|
||||
else:
|
||||
df = self.transitions[shell]
|
||||
for primary, secondary, energy, prob in df.itertuples(index=False):
|
||||
e_row = 0.0
|
||||
if secondary is None:
|
||||
# Fluorescent photon release in radiative transition
|
||||
e_row += energy
|
||||
else:
|
||||
# Fill the hole left by auger electron
|
||||
e_row += self.energy_fluorescence(secondary)
|
||||
|
||||
# Fill the photoelectron hole
|
||||
e_row += self.energy_fluorescence(primary)
|
||||
|
||||
# Expected fluorescent photon energy
|
||||
e += e_row * prob
|
||||
|
||||
self._e_fluorescence[shell] = e
|
||||
return e
|
||||
|
||||
class IncidentPhoton(EqualityMixin):
|
||||
r"""Photon interaction data.
|
||||
|
|
@ -699,9 +656,6 @@ class IncidentPhoton(EqualityMixin):
|
|||
# Add bremsstrahlung DCS data
|
||||
data._add_bremsstrahlung()
|
||||
|
||||
# Add heating cross sections
|
||||
data._compute_heating()
|
||||
|
||||
return data
|
||||
|
||||
@classmethod
|
||||
|
|
@ -917,10 +871,10 @@ class IncidentPhoton(EqualityMixin):
|
|||
|
||||
for j in range(k):
|
||||
# Cubic spline interpolation in log energy and linear DCS
|
||||
cs = CubicSpline(logx, y[:,j])
|
||||
cs = CubicSpline(logx, y[:, j])
|
||||
|
||||
# Get scaled DCS values (millibarns) on new energy grid
|
||||
dcs[:,j] = cs(log_energy)
|
||||
dcs[:, j] = cs(log_energy)
|
||||
|
||||
_BREMSSTRAHLUNG[i]['dcs'] = dcs
|
||||
|
||||
|
|
@ -929,83 +883,6 @@ class IncidentPhoton(EqualityMixin):
|
|||
self.bremsstrahlung['photon_energy'] = _BREMSSTRAHLUNG['photon_energy']
|
||||
self.bremsstrahlung.update(_BREMSSTRAHLUNG[self.atomic_number])
|
||||
|
||||
def _compute_heating(self):
|
||||
r"""Compute heating cross sections (KERMA)
|
||||
|
||||
Photon energy is deposited as energy loss in three reactions:
|
||||
incoherent scattering, pair production and photoelectric effect.
|
||||
The point-wise heating cross section is calculated as:
|
||||
|
||||
.. math::
|
||||
\begin{aligned}
|
||||
\sigma_{Hx}(E) &= (E - \overline{E}_x(E)) \cdot \sigma_x(E), x \in \left\{I, PP, PE \right\} \\
|
||||
\overline{E}_I(E) &= \frac {\int E' \sigma_I (E,E',\mu) d\mu} {\int \sigma_I (E,E',\mu) d\mu} \\
|
||||
\overline{E}_{PP} &= 2 m_e c^2 = 1.022 \times 10^6 eV \\
|
||||
\overline{E}_{PE} &= E(\text{fluorescent photons})
|
||||
\end{aligned}
|
||||
|
||||
The differential cross section representation for incoherent
|
||||
scattering can be found in the theory manual.
|
||||
|
||||
"""
|
||||
|
||||
# Determine a union energy grid
|
||||
energy = np.array([])
|
||||
for mt in (504, 515, 517, 522):
|
||||
if mt in self:
|
||||
energy = np.union1d(energy, self[mt].xs.x)
|
||||
|
||||
heating_xs = np.zeros_like(energy)
|
||||
|
||||
# Incoherent scattering
|
||||
if 504 in self:
|
||||
rx = self[504]
|
||||
|
||||
def dsigma_dmu(mu, E):
|
||||
k = E / MASS_ELECTRON_EV
|
||||
krat = 1.0 / (1.0 + k * (1.0 - mu))
|
||||
x = E * sqrt(0.5 * (1.0 - mu)) / PLANCK_C
|
||||
return pi * R0*R0 * krat*krat * (krat + 1/krat +
|
||||
mu*mu - 1.0) * rx.scattering_factor(x)
|
||||
|
||||
def eout_dsigma_dmu(mu, E):
|
||||
Eout = E / (1.0 + E / MASS_ELECTRON_EV * (1.0 - mu))
|
||||
return Eout * dsigma_dmu(mu, E)
|
||||
|
||||
def eout_average(E):
|
||||
integral_sigma = quad(dsigma_dmu, -1.0, 1.0,
|
||||
args=(E,), epsabs=0.0, epsrel=1e-3)[0]
|
||||
integral_sigma_e = quad(eout_dsigma_dmu, -1.0, 1.0,
|
||||
args=(E,), epsabs=0.0, epsrel=1e-3)[0]
|
||||
return integral_sigma_e / integral_sigma
|
||||
|
||||
e_out = np.vectorize(eout_average)(energy)
|
||||
heating_xs += (energy - e_out) * rx.xs(energy)
|
||||
|
||||
# Pair production, electron field
|
||||
if 515 in self:
|
||||
heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[515].xs(energy)
|
||||
|
||||
# Pair production, nuclear field
|
||||
if 517 in self:
|
||||
heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[517].xs(energy)
|
||||
|
||||
# Photoelectric effect
|
||||
if 522 in self:
|
||||
# Account for fluorescent photons
|
||||
for mt, rx in self.reactions.items():
|
||||
if mt >= 534 and mt <= 572:
|
||||
shell = _REACTION_NAME[mt][1]
|
||||
relax_data = self.atomic_relaxation
|
||||
if relax_data is not None:
|
||||
e_f = relax_data.energy_fluorescence(shell)
|
||||
else:
|
||||
e_f = 0.0
|
||||
heating_xs += (energy - e_f) * rx.xs(energy)
|
||||
|
||||
heat_rx = PhotonReaction(525)
|
||||
heat_rx.xs = Tabulated1D(energy, heating_xs, [energy.size], [5])
|
||||
self.reactions[525] = heat_rx
|
||||
|
||||
class PhotonReaction(EqualityMixin):
|
||||
"""Photon-induced reaction
|
||||
|
|
|
|||
|
|
@ -1,6 +1,6 @@
|
|||
#include "openmc/output.h"
|
||||
|
||||
#include <algorithm> // for std::transform
|
||||
#include <algorithm> // for transform, max
|
||||
#include <cstring> // for strlen
|
||||
#include <ctime> // for time, localtime
|
||||
#include <iomanip> // for setw, setprecision, put_time
|
||||
|
|
@ -527,8 +527,8 @@ std::pair<double, double>
|
|||
mean_stdev(const double* x, int n)
|
||||
{
|
||||
double mean = x[static_cast<int>(TallyResult::SUM)] / n;
|
||||
double stdev = n > 1 ? std::sqrt((x[static_cast<int>(TallyResult::SUM_SQ)]/n
|
||||
- mean*mean)/(n - 1)) : 0.0;
|
||||
double stdev = n > 1 ? std::sqrt(std::max(0.0, (
|
||||
x[static_cast<int>(TallyResult::SUM_SQ)]/n - mean*mean)/(n - 1))) : 0.0;
|
||||
return {mean, stdev};
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -363,15 +363,6 @@ Tally::Tally(pugi::xml_node node)
|
|||
break;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
const auto& f = model::tally_filters[particle_filter_index].get();
|
||||
auto pf = dynamic_cast<ParticleFilter*>(f);
|
||||
for (auto p : pf->particles()) {
|
||||
if (p == Particle::Type::electron ||
|
||||
p == Particle::Type::positron) {
|
||||
estimator_ = TallyEstimator::ANALOG;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if (particle_filter_index >= 0) {
|
||||
|
|
@ -440,9 +431,9 @@ Tally::Tally(pugi::xml_node node)
|
|||
estimator_ = TallyEstimator::ANALOG;
|
||||
} else if (est == "tracklength" || est == "track-length"
|
||||
|| est == "pathlength" || est == "path-length") {
|
||||
// If the estimator was set to an analog estimator, this means the
|
||||
// tally needs post-collision information
|
||||
if (estimator_ == TallyEstimator::ANALOG) {
|
||||
// If the estimator was set to an analog/collision estimator, this means
|
||||
// the tally needs post-collision information
|
||||
if (estimator_ == TallyEstimator::ANALOG || estimator_ == TallyEstimator::COLLISION) {
|
||||
throw std::runtime_error{"Cannot use track-length estimator for tally "
|
||||
+ std::to_string(id_)};
|
||||
}
|
||||
|
|
@ -654,6 +645,10 @@ Tally::set_scores(const std::vector<std::string>& scores)
|
|||
fatal_error("Cannot tally currents without surface type filters");
|
||||
}
|
||||
break;
|
||||
|
||||
case HEATING:
|
||||
if (settings::photon_transport) estimator_ = TallyEstimator::COLLISION;
|
||||
break;
|
||||
}
|
||||
|
||||
scores_.push_back(score);
|
||||
|
|
|
|||
|
|
@ -158,16 +158,18 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index,
|
|||
|
||||
// Determine the filter scoring index
|
||||
auto filter_index = 0;
|
||||
double filter_weight = 1.;
|
||||
for (auto i = 0; i < tally.filters().size(); ++i) {
|
||||
auto i_filt = tally.filters(i);
|
||||
auto& match {filter_matches[i_filt]};
|
||||
auto i_bin = match.i_bin_;
|
||||
filter_index += match.bins_[i_bin] * tally.strides(i);
|
||||
filter_weight *= match.weights_[i_bin];
|
||||
}
|
||||
|
||||
// Update the tally result
|
||||
#pragma omp atomic
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) += score;
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) += score*filter_weight;
|
||||
|
||||
// Reset the original delayed group bin
|
||||
dg_match.bins_[i_bin] = original_bin;
|
||||
|
|
@ -384,18 +386,19 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin)
|
|||
|| (score_bin == SCORE_PROMPT_NU_FISSION && g == 0)) {
|
||||
|
||||
// Find the filter scoring index for this filter combination
|
||||
//TODO: should this include a weight?
|
||||
int filter_index = 0;
|
||||
double filter_weight = 1.0;
|
||||
for (auto j = 0; j < tally.filters().size(); ++j) {
|
||||
auto i_filt = tally.filters(j);
|
||||
auto& match {p->filter_matches_[i_filt]};
|
||||
auto i_bin = match.i_bin_;
|
||||
filter_index += match.bins_[i_bin] * tally.strides(j);
|
||||
filter_weight *= match.weights_[i_bin];
|
||||
}
|
||||
|
||||
// Update tally results
|
||||
#pragma omp atomic
|
||||
tally.results_(filter_index, i_score, TallyResult::VALUE) += score;
|
||||
tally.results_(filter_index, i_score, TallyResult::VALUE) += score*filter_weight;
|
||||
|
||||
} else if (score_bin == SCORE_DELAYED_NU_FISSION && g != 0) {
|
||||
|
||||
|
|
@ -458,8 +461,8 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin)
|
|||
//! is not used for analog tallies.
|
||||
|
||||
void
|
||||
score_general_ce(Particle* p, int i_tally, int start_index,
|
||||
int filter_index, int i_nuclide, double atom_density, double flux)
|
||||
score_general_ce(Particle* p, int i_tally, int start_index, int filter_index,
|
||||
double filter_weight, int i_nuclide, double atom_density, double flux)
|
||||
{
|
||||
Tally& tally {*model::tallies[i_tally]};
|
||||
|
||||
|
|
@ -469,12 +472,9 @@ score_general_ce(Particle* p, int i_tally, int start_index,
|
|||
for (auto i = 0; i < tally.scores_.size(); ++i) {
|
||||
auto score_bin = tally.scores_[i];
|
||||
auto score_index = start_index + i;
|
||||
|
||||
double score;
|
||||
double score = 0.0;
|
||||
|
||||
switch (score_bin) {
|
||||
|
||||
|
||||
case SCORE_FLUX:
|
||||
if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
// All events score to a flux bin. We actually use a collision estimator
|
||||
|
|
@ -514,7 +514,11 @@ score_general_ce(Particle* p, int i_tally, int start_index,
|
|||
|
||||
} else {
|
||||
if (i_nuclide >= 0) {
|
||||
score = p->neutron_xs_[i_nuclide].total * atom_density * flux;
|
||||
if (p->type_ == Particle::Type::neutron) {
|
||||
score = p->neutron_xs_[i_nuclide].total * atom_density * flux;
|
||||
} else if (p->type_ == Particle::Type::photon) {
|
||||
score = p->photon_xs_[i_nuclide].total * atom_density * flux;
|
||||
}
|
||||
} else {
|
||||
score = p->macro_xs_.total * flux;
|
||||
}
|
||||
|
|
@ -1136,8 +1140,9 @@ score_general_ce(Particle* p, int i_tally, int start_index,
|
|||
|
||||
case SCORE_EVENTS:
|
||||
// Simply count the number of scoring events
|
||||
score = 1.;
|
||||
break;
|
||||
#pragma omp atomic
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) += 1.0;
|
||||
continue;
|
||||
|
||||
|
||||
case ELASTIC:
|
||||
|
|
@ -1214,55 +1219,22 @@ score_general_ce(Particle* p, int i_tally, int start_index,
|
|||
|
||||
|
||||
case HEATING:
|
||||
score = 0.;
|
||||
if (p->type_ == Particle::Type::neutron) {
|
||||
score = score_neutron_heating(p, tally, flux, HEATING,
|
||||
i_nuclide, atom_density);
|
||||
} else if (p->type_ == Particle::Type::photon) {
|
||||
if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
// Score direct energy deposition in the collision
|
||||
score = E - p->E_;
|
||||
// We need to substract the energy of the secondary particles since
|
||||
// they will be transported individually later
|
||||
for (auto i = 0; i < p->n_bank_second_; ++i) {
|
||||
auto i_bank = p->secondary_bank_.size() - p->n_bank_second_ + i;
|
||||
const auto& bank = p->secondary_bank_[i_bank];
|
||||
if (bank.particle == Particle::Type::photon ||
|
||||
bank.particle == Particle::Type::neutron) {
|
||||
score -= bank.E;
|
||||
} else if (bank.particle == Particle::Type::positron) {
|
||||
// Annihilation of the positron will produce two new photons
|
||||
score -= 2*MASS_ELECTRON_EV;
|
||||
}
|
||||
}
|
||||
score *= p->wgt_last_ * flux;
|
||||
} else {
|
||||
// Calculate photon heating cross section on-the-fly
|
||||
if (i_nuclide >= 0) {
|
||||
// Find the element corresponding to the nuclide
|
||||
auto name = data::nuclides[i_nuclide]->name_;
|
||||
std::string element = to_element(name);
|
||||
int i_element = data::element_map[element];
|
||||
auto& heating {data::elements[i_element].heating_};
|
||||
auto i_grid = p->photon_xs_[i_element].index_grid;
|
||||
auto f = p->photon_xs_[i_element].interp_factor;
|
||||
score = std::exp(heating(i_grid) + f * (heating(i_grid+1) -
|
||||
heating(i_grid))) * atom_density * flux;
|
||||
} else {
|
||||
if (p->material_ != MATERIAL_VOID) {
|
||||
const Material& material {*model::materials[p->material_]};
|
||||
for (auto i = 0; i < material.nuclide_.size(); ++i) {
|
||||
auto i_element = material.element_[i];
|
||||
auto atom_density = material.atom_density_(i);
|
||||
auto& heating {data::elements[i_element].heating_};
|
||||
auto i_grid = p->photon_xs_[i_element].index_grid;
|
||||
auto f = p->photon_xs_[i_element].interp_factor;
|
||||
score += std::exp(heating(i_grid) + f * (heating(i_grid+1) -
|
||||
heating(i_grid))) * atom_density * flux;
|
||||
}
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// The energy deposited is the difference between the pre-collision and
|
||||
// post-collision energy...
|
||||
score = E - p->E_;
|
||||
|
||||
// ...less the energy of any secondary particles since they will be
|
||||
// transported individually later
|
||||
const auto& bank = p->secondary_bank_;
|
||||
for (auto it = bank.end() - p->n_bank_second_; it < bank.end(); ++it) {
|
||||
score -= it->E;
|
||||
}
|
||||
|
||||
score *= p->wgt_last_;
|
||||
}
|
||||
break;
|
||||
|
||||
|
|
@ -1332,7 +1304,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
|
|||
|
||||
// Update tally results
|
||||
#pragma omp atomic
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) += score;
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) += score*filter_weight;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -1342,8 +1314,8 @@ score_general_ce(Particle* p, int i_tally, int start_index,
|
|||
//! argument is really just used for filter weights.
|
||||
|
||||
void
|
||||
score_general_mg(Particle* p, int i_tally, int start_index,
|
||||
int filter_index, int i_nuclide, double atom_density, double flux)
|
||||
score_general_mg(Particle* p, int i_tally, int start_index, int filter_index,
|
||||
double filter_weight, int i_nuclide, double atom_density, double flux)
|
||||
{
|
||||
auto& tally {*model::tallies[i_tally]};
|
||||
|
||||
|
|
@ -1984,8 +1956,9 @@ score_general_mg(Particle* p, int i_tally, int start_index,
|
|||
|
||||
case SCORE_EVENTS:
|
||||
// Simply count the number of scoring events
|
||||
score = 1.;
|
||||
break;
|
||||
#pragma omp atomic
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) += 1.0;
|
||||
continue;
|
||||
|
||||
|
||||
default:
|
||||
|
|
@ -1994,7 +1967,7 @@ score_general_mg(Particle* p, int i_tally, int start_index,
|
|||
|
||||
// Update tally results
|
||||
#pragma omp atomic
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) += score;
|
||||
tally.results_(filter_index, score_index, TallyResult::VALUE) += score*filter_weight;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -2002,7 +1975,7 @@ score_general_mg(Particle* p, int i_tally, int start_index,
|
|||
|
||||
void
|
||||
score_all_nuclides(Particle* p, int i_tally, double flux,
|
||||
int filter_index)
|
||||
int filter_index, double filter_weight)
|
||||
{
|
||||
const Tally& tally {*model::tallies[i_tally]};
|
||||
const Material& material {*model::materials[p->material_]};
|
||||
|
|
@ -2015,10 +1988,10 @@ score_all_nuclides(Particle* p, int i_tally, double flux,
|
|||
//TODO: consider replacing this "if" with pointers or templates
|
||||
if (settings::run_CE) {
|
||||
score_general_ce(p, i_tally, i_nuclide*tally.scores_.size(), filter_index,
|
||||
i_nuclide, atom_density, flux);
|
||||
filter_weight, i_nuclide, atom_density, flux);
|
||||
} else {
|
||||
score_general_mg(p, i_tally, i_nuclide*tally.scores_.size(), filter_index,
|
||||
i_nuclide, atom_density, flux);
|
||||
filter_weight, i_nuclide, atom_density, flux);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -2029,15 +2002,22 @@ score_all_nuclides(Particle* p, int i_tally, double flux,
|
|||
//TODO: consider replacing this "if" with pointers or templates
|
||||
if (settings::run_CE) {
|
||||
score_general_ce(p, i_tally, n_nuclides*tally.scores_.size(), filter_index,
|
||||
i_nuclide, atom_density, flux);
|
||||
filter_weight, i_nuclide, atom_density, flux);
|
||||
} else {
|
||||
score_general_mg(p, i_tally, n_nuclides*tally.scores_.size(), filter_index,
|
||||
i_nuclide, atom_density, flux);
|
||||
filter_weight, i_nuclide, atom_density, flux);
|
||||
}
|
||||
}
|
||||
|
||||
void score_analog_tally_ce(Particle* p)
|
||||
{
|
||||
// Since electrons/positrons are not transported, we assign a flux of zero.
|
||||
// Note that the heating score does NOT use the flux and will be non-zero for
|
||||
// electrons/positrons.
|
||||
double flux =
|
||||
(p->type_ == Particle::Type::neutron || p->type_ == Particle::Type::photon) ?
|
||||
1.0 : 0.0;
|
||||
|
||||
for (auto i_tally : model::active_analog_tallies) {
|
||||
const Tally& tally {*model::tallies[i_tally]};
|
||||
|
||||
|
|
@ -2063,7 +2043,7 @@ void score_analog_tally_ce(Particle* p)
|
|||
// density argument for score_general is not used for analog tallies.
|
||||
if (i_nuclide == p->event_nuclide_ || i_nuclide == -1)
|
||||
score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index,
|
||||
i_nuclide, -1., filter_weight);
|
||||
filter_weight, i_nuclide, -1.0, flux);
|
||||
}
|
||||
|
||||
} else {
|
||||
|
|
@ -2072,12 +2052,12 @@ void score_analog_tally_ce(Particle* p)
|
|||
// bins correspond to nuclide indices. First, tally the nuclide.
|
||||
auto i = p->event_nuclide_;
|
||||
score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index,
|
||||
-1, -1., filter_weight);
|
||||
filter_weight, -1, -1.0, flux);
|
||||
|
||||
// Now tally the total material.
|
||||
i = tally.nuclides_.size();
|
||||
score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index,
|
||||
-1, -1., filter_weight);
|
||||
filter_weight, -1, -1.0, flux);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -2122,7 +2102,7 @@ void score_analog_tally_mg(Particle* p)
|
|||
}
|
||||
|
||||
score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index,
|
||||
i_nuclide, atom_density, filter_weight);
|
||||
filter_weight, i_nuclide, atom_density, 1.0);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -2162,7 +2142,8 @@ score_tracklength_tally(Particle* p, double distance)
|
|||
// Loop over nuclide bins.
|
||||
if (tally.all_nuclides_) {
|
||||
if (p->material_ != MATERIAL_VOID)
|
||||
score_all_nuclides(p, i_tally, flux*filter_weight, filter_index);
|
||||
score_all_nuclides(p, i_tally, flux*filter_weight, filter_index,
|
||||
filter_weight);
|
||||
|
||||
} else {
|
||||
for (auto i = 0; i < tally.nuclides_.size(); ++i) {
|
||||
|
|
@ -2180,10 +2161,10 @@ score_tracklength_tally(Particle* p, double distance)
|
|||
//TODO: consider replacing this "if" with pointers or templates
|
||||
if (settings::run_CE) {
|
||||
score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index,
|
||||
i_nuclide, atom_density, flux*filter_weight);
|
||||
filter_weight, i_nuclide, atom_density, flux);
|
||||
} else {
|
||||
score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index,
|
||||
i_nuclide, atom_density, flux*filter_weight);
|
||||
filter_weight, i_nuclide, atom_density, flux);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -2205,11 +2186,13 @@ score_tracklength_tally(Particle* p, double distance)
|
|||
void score_collision_tally(Particle* p)
|
||||
{
|
||||
// Determine the collision estimate of the flux
|
||||
double flux;
|
||||
if (!settings::survival_biasing) {
|
||||
flux = p->wgt_last_ / p->macro_xs_.total;
|
||||
} else {
|
||||
flux = (p->wgt_last_ + p->wgt_absorb_) / p->macro_xs_.total;
|
||||
double flux = 0.0;
|
||||
if (p->type_ == Particle::Type::neutron || p->type_ == Particle::Type::photon) {
|
||||
if (!settings::survival_biasing) {
|
||||
flux = p->wgt_last_ / p->macro_xs_.total;
|
||||
} else {
|
||||
flux = (p->wgt_last_ + p->wgt_absorb_) / p->macro_xs_.total;
|
||||
}
|
||||
}
|
||||
|
||||
for (auto i_tally : model::active_collision_tallies) {
|
||||
|
|
@ -2229,7 +2212,8 @@ void score_collision_tally(Particle* p)
|
|||
|
||||
// Loop over nuclide bins.
|
||||
if (tally.all_nuclides_) {
|
||||
score_all_nuclides(p, i_tally, flux*filter_weight, filter_index);
|
||||
score_all_nuclides(p, i_tally, flux*filter_weight, filter_index,
|
||||
filter_weight);
|
||||
|
||||
} else {
|
||||
for (auto i = 0; i < tally.nuclides_.size(); ++i) {
|
||||
|
|
@ -2245,10 +2229,10 @@ void score_collision_tally(Particle* p)
|
|||
//TODO: consider replacing this "if" with pointers or templates
|
||||
if (settings::run_CE) {
|
||||
score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index,
|
||||
i_nuclide, atom_density, flux*filter_weight);
|
||||
filter_weight, i_nuclide, atom_density, flux);
|
||||
} else {
|
||||
score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index,
|
||||
i_nuclide, atom_density, flux*filter_weight);
|
||||
filter_weight, i_nuclide, atom_density, flux);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -41,7 +41,7 @@
|
|||
<bins>1</bins>
|
||||
</filter>
|
||||
<filter id="2" type="particle">
|
||||
<bins>photon</bins>
|
||||
<bins>neutron photon electron positron</bins>
|
||||
</filter>
|
||||
<tally id="1">
|
||||
<filters>1 2</filters>
|
||||
|
|
@ -50,7 +50,7 @@
|
|||
<tally id="2">
|
||||
<filters>2</filters>
|
||||
<nuclides>Al27 total</nuclides>
|
||||
<scores>total heating</scores>
|
||||
<scores>total</scores>
|
||||
<estimator>tracklength</estimator>
|
||||
</tally>
|
||||
<tally id="3">
|
||||
|
|
|
|||
|
|
@ -1,30 +1,92 @@
|
|||
tally 1:
|
||||
8.601000E-01
|
||||
7.397720E-01
|
||||
9.403000E-01
|
||||
8.841641E-01
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
tally 2:
|
||||
1.819886E-01
|
||||
3.311985E-02
|
||||
1.960159E+05
|
||||
3.842224E+10
|
||||
1.249805E+00
|
||||
1.562013E+00
|
||||
1.249805E+00
|
||||
1.562013E+00
|
||||
8.281718E-01
|
||||
6.858685E-01
|
||||
1.960159E+05
|
||||
3.842224E+10
|
||||
8.281718E-01
|
||||
6.858685E-01
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
tally 3:
|
||||
1.799069E-01
|
||||
3.236650E-02
|
||||
1.945225E+05
|
||||
3.783901E+10
|
||||
1.248100E+00
|
||||
1.557754E+00
|
||||
2.024868E+06
|
||||
4.100089E+12
|
||||
1.248100E+00
|
||||
1.557754E+00
|
||||
2.024868E+06
|
||||
4.100089E+12
|
||||
8.242000E-01
|
||||
6.793056E-01
|
||||
1.945225E+05
|
||||
3.783901E+10
|
||||
5.827819E+03
|
||||
3.396348E+07
|
||||
8.242000E-01
|
||||
6.793056E-01
|
||||
5.827819E+03
|
||||
3.396348E+07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
1.834863E+05
|
||||
3.366722E+10
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
1.834863E+05
|
||||
3.366722E+10
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
5.938849E+03
|
||||
3.526993E+07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
5.938849E+03
|
||||
3.526993E+07
|
||||
tally 4:
|
||||
1.248100E+00
|
||||
1.557754E+00
|
||||
2.024868E+06
|
||||
4.100089E+12
|
||||
1.248100E+00
|
||||
1.557754E+00
|
||||
2.024868E+06
|
||||
4.100089E+12
|
||||
2.308000E-01
|
||||
5.326864E-02
|
||||
1.988563E+05
|
||||
3.954384E+10
|
||||
5.822937E+03
|
||||
3.390660E+07
|
||||
8.242000E-01
|
||||
6.793056E-01
|
||||
1.988612E+05
|
||||
3.954578E+10
|
||||
5.827819E+03
|
||||
3.396348E+07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
1.834863E+05
|
||||
3.366722E+10
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
5.938849E+03
|
||||
3.526993E+07
|
||||
|
|
|
|||
|
|
@ -42,13 +42,13 @@ def model():
|
|||
model.settings.source = source
|
||||
|
||||
surface_filter = openmc.SurfaceFilter(cyl)
|
||||
particle_filter = openmc.ParticleFilter('photon')
|
||||
particle_filter = openmc.ParticleFilter(['neutron', 'photon', 'electron', 'positron'])
|
||||
current_tally = openmc.Tally()
|
||||
current_tally.filters = [surface_filter, particle_filter]
|
||||
current_tally.scores = ['current']
|
||||
tally_tracklength = openmc.Tally()
|
||||
tally_tracklength.filters = [particle_filter]
|
||||
tally_tracklength.scores = ['total', 'heating']
|
||||
tally_tracklength.scores = ['total'] # heating doesn't work with tracklength
|
||||
tally_tracklength.nuclides = ['Al27', 'total']
|
||||
tally_tracklength.estimator = 'tracklength'
|
||||
tally_collision = openmc.Tally()
|
||||
|
|
|
|||
|
|
@ -31,7 +31,7 @@
|
|||
<tally id="1">
|
||||
<filters>1</filters>
|
||||
<nuclides>U235 total</nuclides>
|
||||
<scores>fission heating heating-local</scores>
|
||||
<scores>fission heating-local</scores>
|
||||
<estimator>tracklength</estimator>
|
||||
</tally>
|
||||
<tally id="2">
|
||||
|
|
|
|||
|
|
@ -3,26 +3,18 @@ k-combined:
|
|||
tally 1:
|
||||
2.600054E+00
|
||||
2.253433E+00
|
||||
4.390171E+08
|
||||
6.424554E+16
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
2.600054E+00
|
||||
2.253433E+00
|
||||
4.390171E+08
|
||||
6.424554E+16
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.764667E+07
|
||||
7.568912E+14
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.764667E+07
|
||||
7.568912E+14
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
tally 2:
|
||||
|
|
@ -40,14 +32,14 @@ tally 2:
|
|||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.749088E+07
|
||||
7.519602E+14
|
||||
2.921131E+06
|
||||
2.853657E+12
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.749088E+07
|
||||
7.519602E+14
|
||||
2.921131E+06
|
||||
2.853657E+12
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
tally 3:
|
||||
|
|
@ -65,13 +57,13 @@ tally 3:
|
|||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.729796E+07
|
||||
7.458273E+14
|
||||
2.780823E+06
|
||||
2.586737E+12
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
4.743827E+07
|
||||
7.502589E+14
|
||||
2.921131E+06
|
||||
2.853657E+12
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
|
|
|
|||
|
|
@ -25,7 +25,7 @@ def model():
|
|||
particle_filter = openmc.ParticleFilter(['neutron', 'photon'])
|
||||
tally_tracklength = openmc.Tally()
|
||||
tally_tracklength.filters = [particle_filter]
|
||||
tally_tracklength.scores = ['fission', 'heating', 'heating-local']
|
||||
tally_tracklength.scores = ['fission', 'heating-local']
|
||||
tally_tracklength.nuclides = ['U235', 'total']
|
||||
tally_tracklength.estimator = 'tracklength'
|
||||
tally_collision = openmc.Tally()
|
||||
|
|
@ -43,6 +43,6 @@ def model():
|
|||
return model
|
||||
|
||||
|
||||
def test_photon_production(model):
|
||||
def test_photon_production_fission(model):
|
||||
harness = PyAPITestHarness('statepoint.5.h5', model)
|
||||
harness.main()
|
||||
|
|
|
|||
|
|
@ -1 +1 @@
|
|||
67a1b19f5ee0260a5efd99e0712c5dc8c770b04bcb70aad5c47cdc67e8ea2f9148a9d122598d0f0d2bc2f2ae0222e1a78c8fed10dfb49866123a137b6c1c1264
|
||||
4d9e4a6d4891d02fac2222e0c8665033f0c955ac0082e0c6851c7674a6e967721bad8b77bf0d6c243736128a77a274b4664af54ab10462416ef6c4d18e7218ca
|
||||
Loading…
Add table
Add a link
Reference in a new issue