mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
Remove photon heating calculation from Python and use on C++ side
This commit is contained in:
Paul Romano
parent
37cb75a5be
commit
840749070d
3 changed files with 8 additions and 121 deletions
|
|
@ -699,9 +699,6 @@ class IncidentPhoton(EqualityMixin):
|
|||
# Add bremsstrahlung DCS data
|
||||
data._add_bremsstrahlung()
|
||||
|
||||
# Add heating cross sections
|
||||
data._compute_heating()
|
||||
|
||||
return data
|
||||
|
||||
@classmethod
|
||||
|
|
@ -929,83 +926,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
|
||||
|
|
|
|||
|
|
@ -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) {
|
||||
|
|
|
|||
|
|
@ -1218,7 +1218,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
|
|||
if (p->type_ == Particle::Type::neutron) {
|
||||
score = score_neutron_heating(p, tally, flux, HEATING,
|
||||
i_nuclide, atom_density);
|
||||
} else if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
} else {
|
||||
// The energy deposited is the difference between the pre-collision and
|
||||
// post-collision energy...
|
||||
score = E - p->E_;
|
||||
|
|
@ -1231,32 +1231,6 @@ score_general_ce(Particle* p, int i_tally, int start_index,
|
|||
}
|
||||
|
||||
score *= p->wgt_last_ * flux;
|
||||
} else if (p->type_ == Particle::Type::photon) {
|
||||
// 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
|
||||
|
|
@ -2199,11 +2173,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) {
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue