Commit 7b6ec80a authored by Paul Romano's avatar Paul Romano Committed by GitHub

Merge pull request #1203 from liangjg/heat-photon

Photon heating tally
parents 5fe2fbfa 21c693e4
......@@ -165,6 +165,10 @@ Incident Photon Data
- **J** (*double[][]*) -- Compton profile for each subshell in units
of :math:`\hbar / (me^2)`
**/<element>/heating/**
:Datasets: - **xs** (*double[]*) -- Total heating cross section in [b-eV]
**/<element>/incoherent/**
:Datasets: - **xs** (*double[]*) -- Incoherent scattering cross section in [b]
......
......@@ -256,9 +256,12 @@ The following tables show all valid scores:
|inverse-velocity |The flux-weighted inverse velocity where the |
| |velocity is in units of centimeters per second. |
+----------------------+---------------------------------------------------+
|heating |Total neutron heating in units of eV per source |
| |particle. This corresponds to MT=301 produced by |
| |NJOY's HEATR module. |
|heating |Total nuclear heating in units of eV per source |
| |particle. For neutrons, this corresponds to MT=301 |
| |produced by NJOY's HEATR module while for photons, |
| |this is tallied from either direct photon energy |
| |deposition (analog estimator) or pre-generated |
| |photon heating number. |
+----------------------+---------------------------------------------------+
|kappa-fission |The recoverable energy production rate due to |
| |fission. The recoverable energy is defined as the |
......
......@@ -105,9 +105,10 @@ constexpr std::array<const char*, 39> SUBSHELLS = {
"Q1", "Q2", "Q3"
};
// Void material
// Void material and nuclide
// TODO: refactor and remove
constexpr int MATERIAL_VOID {-1};
constexpr int NUCLIDE_NONE {-1};
// ============================================================================
// CROSS SECTION RELATED CONSTANTS
......@@ -127,6 +128,7 @@ constexpr int TEMPERATURE_INTERPOLATION {2};
// Reaction types
// TODO: Convert to enum
constexpr int REACTION_NONE {0};
constexpr int TOTAL_XS {1};
constexpr int ELASTIC {2};
constexpr int N_NONELASTIC {3};
......@@ -230,7 +232,7 @@ constexpr int N_XD {204};
constexpr int N_XT {205};
constexpr int N_X3HE {206};
constexpr int N_XA {207};
constexpr int HEATING {301};
constexpr int NEUTRON_HEATING {301};
constexpr int DAMAGE_ENERGY {444};
constexpr int COHERENT {502};
constexpr int INCOHERENT {504};
......@@ -344,6 +346,7 @@ constexpr int ESTIMATOR_COLLISION {3};
// TODO: Convert to enum
constexpr int EVENT_SURFACE {-2};
constexpr int EVENT_LATTICE {-1};
constexpr int EVENT_KILL {0};
constexpr int EVENT_SCATTER {1};
constexpr int EVENT_ABSORB {2};
......@@ -366,6 +369,7 @@ constexpr int SCORE_INVERSE_VELOCITY {-13}; // flux-weighted inverse velocity
constexpr int SCORE_FISS_Q_PROMPT {-14}; // prompt fission Q-value
constexpr int SCORE_FISS_Q_RECOV {-15}; // recoverable fission Q-value
constexpr int SCORE_DECAY_RATE {-16}; // delayed neutron precursor decay rate
constexpr int SCORE_HEATING {-17}; // nuclear heating (neutron or photon)
// Tally map bin finding
constexpr int NO_BIN_FOUND {-1};
......
......@@ -184,7 +184,7 @@ public:
//! \param u Direction of the secondary particle
//! \param E Energy of the secondary particle in [eV]
//! \param type Particle type
void create_secondary(Direction u, double E, Type type) const;
void create_secondary(Direction u, double E, Type type);
//! initialize from a source site
//
......@@ -261,6 +261,7 @@ public:
// Post-collision physical data
int n_bank_ {0}; //!< number of fission sites banked
int n_bank_second_ {0}; //!< number of secondary particles banked
double wgt_bank_ {0.0}; //!< weight of fission sites banked
int n_delayed_bank_[MAX_DELAYED_GROUPS]; //!< number of delayed fission
//!< sites banked
......
......@@ -67,6 +67,7 @@ public:
xt::xtensor<double, 1> pair_production_total_;
xt::xtensor<double, 1> pair_production_electron_;
xt::xtensor<double, 1> pair_production_nuclear_;
xt::xtensor<double, 1> heating_;
// Form factors
Tabulated1D incoherent_form_factor_;
......
......@@ -10,6 +10,8 @@ std::string& strtrim(std::string& s);
char* strtrim(char* c_str);
std::string to_element(const std::string& name);
void to_lower(std::string& str);
int word_count(std::string const& str);
......
......@@ -89,7 +89,7 @@ _SCORES = {
-5: 'absorption', -6: 'fission', -7: 'nu-fission', -8: 'kappa-fission',
-9: 'current', -10: 'events', -11: 'delayed-nu-fission',
-12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt',
-15: 'fission-q-recoverable', -16: 'decay-rate'
-15: 'fission-q-recoverable', -16: 'decay-rate', -17: 'heating'
}
_ESTIMATORS = {
1: 'analog', 2: 'tracklength', 3: 'collision'
......
......@@ -3,6 +3,7 @@ from collections.abc import Iterable, Callable
from functools import reduce
from itertools import zip_longest
from numbers import Real, Integral
from math import exp, log
import numpy as np
......@@ -153,13 +154,11 @@ class Tabulated1D(Function1D):
self.y = np.asarray(y)
def __call__(self, x):
# Check if input is array or scalar
if isinstance(x, Iterable):
iterable = True
x = np.array(x)
else:
iterable = False
x = np.array([x], dtype=float)
# Check if input is scalar
if not isinstance(x, Iterable):
return self._interpolate_scalar(x)
x = np.array(x)
# Create output array
y = np.zeros_like(x)
......@@ -208,7 +207,46 @@ class Tabulated1D(Function1D):
y[np.isclose(x, self.x[0], atol=1e-14)] = self.y[0]
y[np.isclose(x, self.x[-1], atol=1e-14)] = self.y[-1]
return y if iterable else y[0]
return y
def _interpolate_scalar(self, x):
if x <= self._x[0]:
return self._y[0]
elif x >= self._x[-1]:
return self._y[-1]
# Get the index for interpolation
idx = np.searchsorted(self._x, x, side='right') - 1
# Loop over interpolation regions
for b, p in zip(self.breakpoints, self.interpolation):
if idx < b - 1:
break
xi = self._x[idx] # low edge of the corresponding bin
xi1 = self._x[idx + 1] # high edge of the corresponding bin
yi = self._y[idx]
yi1 = self._y[idx + 1]
if p == 1:
# Histogram
return yi
elif p == 2:
# Linear-linear
return yi + (x - xi)/(xi1 - xi)*(yi1 - yi)
elif p == 3:
# Linear-log
return yi + log(x/xi)/log(xi1/xi)*(yi1 - yi)
elif p == 4:
# Log-linear
return yi*exp((x - xi)/(xi1 - xi)*log(yi1/yi))
elif p == 5:
# Log-log
return yi*exp(log(x/xi)/log(xi1/xi)*log(yi1/yi))
def __len__(self):
return len(self.x)
......
......@@ -3,12 +3,14 @@ from collections.abc import Mapping, Callable
from copy import deepcopy
from io import StringIO
from numbers import Integral, Real
from math import pi, sqrt
import os
import h5py
import numpy as np
import pandas as pd
from scipy.interpolate import CubicSpline
from scipy.integrate import quad
from openmc.mixin import EqualityMixin
import openmc.checkvalue as cv
......@@ -19,6 +21,14 @@ from .endf import Evaluation, get_head_record, get_tab1_record, get_list_record
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
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)
# Electron subshell labels
_SUBSHELLS = [None, 'K', 'L1', 'L2', 'L3', 'M1', 'M2', 'M3', 'M4', 'M5',
'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7', 'O1', 'O2', 'O3',
......@@ -33,6 +43,7 @@ _REACTION_NAME = {
516: ('Total pair production', 'pair_production_total'),
517: ('Pair production, nuclear field', 'pair_production_nuclear'),
522: ('Photoelectric absorption', 'photoelectric'),
525: ('Heating', 'heating'),
526: ('Electro-atomic scattering', 'electro_atomic_scat'),
527: ('Electro-atomic bremsstrahlung', 'electro_atomic_brem'),
528: ('Electro-atomic excitation', 'electro_atomic_excit'),
......@@ -150,6 +161,7 @@ class AtomicRelaxation(EqualityMixin):
self.binding_energy = binding_energy
self.num_electrons = num_electrons
self.transitions = transitions
self._e_fluorescence = {}
@property
def binding_energy(self):
......@@ -379,6 +391,49 @@ 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.
......@@ -499,9 +554,14 @@ class IncidentPhoton(EqualityMixin):
# Read each reaction
data = cls(Z)
for mt in (502, 504, 515, 522):
for mt in (502, 504, 515, 522, 525):
data.reactions[mt] = PhotonReaction.from_ace(ace, mt)
# Get heating cross sections [eV-barn] from factors [eV per collision]
# by multiplying with total xs
data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in
(502, 504, 515, 522)])
# Compton profiles
n_shell = ace.nxs[5]
if n_shell != 0:
......@@ -631,6 +691,9 @@ class IncidentPhoton(EqualityMixin):
# Add bremsstrahlung DCS data
data._add_bremsstrahlung()
# Add heating cross sections
data._compute_heating()
return data
@classmethod
......@@ -751,7 +814,7 @@ class IncidentPhoton(EqualityMixin):
designators = []
for mt, rx in self.reactions.items():
name, key = _REACTION_NAME[mt]
if mt in [502, 504, 515, 517, 522]:
if mt in (502, 504, 515, 517, 522, 525):
sub_group = group.create_group(key)
elif mt >= 534 and mt <= 572:
# Subshell
......@@ -858,6 +921,80 @@ 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::
\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})
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]
e_f = self.atomic_relaxation.energy_fluorescence(shell)
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
......@@ -972,14 +1109,21 @@ class PhotonReaction(EqualityMixin):
elif mt == 522:
# Photoelectric
idx = ace.jxs[1] + 3*n
elif mt == 525:
# Heating
idx = ace.jxs[5]
else:
raise ValueError('ACE photoatomic cross sections do not have '
'data for MT={}.'.format(mt))
# Store cross section
xs = ace.xss[idx : idx+n].copy()
nonzero = (xs != 0.0)
xs[nonzero] = np.exp(xs[nonzero])
if mt == 525:
# Get heating factors in [eV per collision]
xs *= EV_PER_MEV
else:
nonzero = (xs != 0.0)
xs[nonzero] = np.exp(xs[nonzero])
rx.xs = Tabulated1D(energy, xs, [n], [5])
# Get form factors for incoherent/coherent scattering
......
......@@ -239,8 +239,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
already_read.insert(name);
// Check if elemental data has been read, if needed
int pos = name.find_first_of("0123456789");
std::string element = name.substr(0, pos);
std::string element = to_element(name);
if (settings::photon_transport) {
if (already_read.find(element) == already_read.end()) {
// Read photon interaction data from HDF5 photon library
......
......@@ -227,8 +227,7 @@ Material::Material(pugi::xml_node node)
// If the corresponding element hasn't been encountered yet and photon
// transport will be used, we need to add its symbol to the element_dict
if (settings::photon_transport) {
int pos = name.find_first_of("0123456789");
std::string element = name.substr(0, pos);
std::string element = to_element(name);
// Make sure photon cross section data is available
LibraryKey key {Library::Type::photon, element};
......
......@@ -617,6 +617,7 @@ const std::unordered_map<int, const char*> score_names = {
{SCORE_FISS_Q_PROMPT, "Prompt fission power"},
{SCORE_FISS_Q_RECOV, "Recoverable fission power"},
{SCORE_CURRENT, "Current"},
{SCORE_HEATING, "Heating"},
};
//! Create an ASCII output file showing all tally results.
......
......@@ -73,7 +73,7 @@ Particle::clear()
}
void
Particle::create_secondary(Direction u, double E, Type type) const
Particle::create_secondary(Direction u, double E, Type type)
{
simulation::secondary_bank.emplace_back();
......@@ -83,6 +83,8 @@ Particle::create_secondary(Direction u, double E, Type type) const
bank.r = this->r();
bank.u = u;
bank.E = settings::run_CE ? E : g_;
n_bank_second_ += 1;
}
void
......@@ -157,6 +159,11 @@ Particle::transport()
u_last_ = this->u();
r_last_ = this->r();
// Reset event variables
event_ = EVENT_KILL;
event_nuclide_ = NUCLIDE_NONE;
event_mt_ = REACTION_NONE;
// If the cell hasn't been determined based on the particle's location,
// initiate a search for the current cell. This generally happens at the
// beginning of the history and again for any secondary particles
......@@ -309,6 +316,7 @@ Particle::transport()
// Reset banked weight during collision
n_bank_ = 0;
n_bank_second_ = 0;
wgt_bank_ = 0.0;
for (int& v : n_delayed_bank_) v = 0;
......
......@@ -96,6 +96,15 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
read_dataset(rgroup, "xs", photoelectric_total_);
close_group(rgroup);
// Read heating
if (object_exists(group, "heating")) {
rgroup = open_group(group, "heating");
read_dataset(rgroup, "xs", heating_);
close_group(rgroup);
} else {
heating_ = xt::zeros_like(energy_);
}
// Read subshell photoionization cross section and atomic relaxation data
rgroup = open_group(group, "subshells");
std::vector<std::string> designators;
......@@ -280,6 +289,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
xt::log(photoelectric_total_), -500.0);
pair_production_total_ = xt::where(pair_production_total_ > 0.0,
xt::log(pair_production_total_), -500.0);
heating_ = xt::where(heating_ > 0.0, xt::log(heating_), -500.0);
}
void PhotonInteraction::compton_scatter(double alpha, bool doppler,
......
......@@ -17,6 +17,7 @@
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/thermal.h"
#include "openmc/tallies/tally.h"
......@@ -61,12 +62,16 @@ void collision(Particle* p)
// Display information about collision
if (settings::verbosity >= 10 || simulation::trace) {
std::stringstream msg;
if (p->type_ == Particle::Type::neutron) {
if (p->event_ == EVENT_KILL) {
msg << " Killed. Energy = " << p->E_ << " eV.";
} else if (p->type_ == Particle::Type::neutron) {
msg << " " << reaction_name(p->event_mt_) << " with " <<
data::nuclides[p->event_nuclide_]->name_ << ". Energy = " << p->E_ << " eV.";
} else {
} else if (p->type_ == Particle::Type::photon) {
msg << " " << reaction_name(p->event_mt_) << " with " <<
data::elements[p->event_nuclide_].name_ << ". Energy = " << p->E_ << " eV.";
to_element(data::nuclides[p->event_nuclide_]->name_) << ". Energy = " << p->E_ << " eV.";
} else {
msg << " Disappeared. Energy = " << p->E_ << " eV.";
}
write_message(msg, 1);
}
......@@ -209,7 +214,6 @@ void sample_photon_reaction(Particle* p)
// Sample element within material
int i_element = sample_element(p);
p->event_nuclide_ = i_element;
const auto& micro {p->photon_xs_[i_element]};
const auto& element {data::elements[i_element]};
......@@ -226,6 +230,7 @@ void sample_photon_reaction(Particle* p)
if (prob > cutoff) {
double mu = element.rayleigh_scatter(alpha);
p->u() = rotate_angle(p->u(), mu, nullptr);
p->event_ = EVENT_SCATTER;
p->event_mt_ = COHERENT;
return;
}
......@@ -268,6 +273,7 @@ void sample_photon_reaction(Particle* p)
phi += PI;
p->E_ = alpha_out*MASS_ELECTRON_EV;
p->u() = rotate_angle(p->u(), mu, &phi);
p->event_ = EVENT_SCATTER;
p->event_mt_ = INCOHERENT;
return;
}
......@@ -319,6 +325,7 @@ void sample_photon_reaction(Particle* p)
// Allow electrons to fill orbital and produce auger electrons
// and fluorescent photons
element.atomic_relaxation(shell, *p);
p->event_ = EVENT_ABSORB;
p->event_mt_ = 533 + shell.index_subshell;
p->alive_ = false;
p->E_ = 0.0;
......@@ -344,6 +351,7 @@ void sample_photon_reaction(Particle* p)
u = rotate_angle(p->u(), mu_positron, nullptr);
p->create_secondary(u, E_positron, Particle::Type::positron);
p->event_ = EVENT_ABSORB;
p->event_mt_ = PAIR_PROD;
p->alive_ = false;
p->E_ = 0.0;
......@@ -361,6 +369,7 @@ void sample_electron_reaction(Particle* p)
p->E_ = 0.0;
p->alive_ = false;
p->event_ = EVENT_ABSORB;
}
void sample_positron_reaction(Particle* p)
......@@ -386,6 +395,7 @@ void sample_positron_reaction(Particle* p)
p->E_ = 0.0;
p->alive_ = false;
p->event_ = EVENT_ABSORB;
}
int sample_nuclide(const Particle* p)
......@@ -432,7 +442,12 @@ int sample_element(Particle* p)
// Increment probability to compare to cutoff
prob += sigma;
if (prob > cutoff) return i_element;
if (prob > cutoff) {
// Save which nuclide particle had collision with for tally purpose
p->event_nuclide_ = mat->nuclide_[i];
return i_element;
}
}
// If we made it here, no element was sampled
......
......@@ -121,6 +121,8 @@ std::string reaction_name(int mt)
return "fission-q-prompt";
} else if (mt == SCORE_FISS_Q_RECOV) {
return "fission-q-recoverable";
} else if (mt == SCORE_HEATING) {
return "heating";
// Normal ENDF-based reactions
} else if (mt == TOTAL_XS) {
......
......@@ -26,6 +26,12 @@ char* strtrim(char* c_str)
}
std::string to_element(const std::string& name) {
int pos = name.find_first_of("0123456789");
return name.substr(0, pos);
}
void to_lower(std::string& str)
{
for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]);
......
......@@ -9,6 +9,7 @@
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/particle.h"
#include "openmc/reaction.h"
#include "openmc/reaction_product.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
......@@ -109,6 +110,9 @@ score_str_to_int(std::string score_str)
if (score_str == "fission-q-recoverable")
return SCORE_FISS_Q_RECOV;
if (score_str == "heating")
return SCORE_HEATING;
if (score_str == "current")
return SCORE_CURRENT;
......@@ -207,8 +211,6 @@ score_str_to_int(std::string score_str)
return N_X3HE;
if (score_str == "(n,Xa)" || score_str == "He4-production")
return N_XA;
if (score_str == "heating")
return HEATING;
if (score_str == "damage-energy")
return DAMAGE_ENERGY;
......@@ -771,7 +773,7 @@ void read_tallies_xml()
case SCORE_PROMPT_NU_FISSION:
case SCORE_DECAY_RATE:
warning("Particle filter is not used with photon transport"
" on and " + std::to_string(score) + " score.");
" on and " + reaction_name(score) + " score.");
break;
}
}
......
......@@ -7,10 +7,12 @@
#include "openmc/material.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/photon.h"
#include "openmc/reaction_product.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/filter_delayedgroup.h"
......@@ -1149,6 +1151,156 @@ score_general_ce(Particle* p, int i_tally, int start_index,
break;
case SCORE_HEATING:
score = 0.;
if (p->type_ == Particle::Type::neutron) {
if (tally.estimator_ == ESTIMATOR_ANALOG) {
// All events score to a heating tally bin. We actually use a
// collision estimator in place of an analog one since there is no
// reaction-wise heating cross section
if (settings::survival_biasing) {
// We need to account for the fact that some weight was already
// absorbed
score = p->wgt_last_ + p->wgt_absorb_;
} else {
score = p->wgt_last_;
}
if (i_nuclide >= 0) {
// Calculate nuclide heating cross section
double macro_heating = 0.;
const auto& nuc {*data::nuclides[i_nuclide]};
auto m = nuc.reaction_index_[NEUTRON_HEATING];
if (m == C_NONE) continue;
const auto& rxn {*nuc.reactions_[m]};
auto i_temp = p->neutron_xs_[i_nuclide].index_temp;
if (i_temp >= 0) { // Can be false due to multipole
auto i_grid = p->neutron_xs_[i_nuclide].index_grid;
auto f = p->neutron_xs_[i_nuclide].interp_factor;
const auto& xs {rxn.xs_[i_temp]};
if (i_grid >= xs.threshold) {
macro_heating = ((1.0 - f) * xs.value[i_grid-xs.threshold]
+ f * xs.value[i_grid-xs.threshold+1]);
}</