mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Merge pull request #1603 from paulromano/photon-tally-fix
Bug fix for tallying reaction rates in coupled neutron-photon runs
This commit is contained in:
Amanda Lund
GitHub
commit
80ddc6e379
10 changed files with 174 additions and 49 deletions
|
|
@ -108,9 +108,11 @@ The following tables show all valid scores:
|
|||
+----------------------+---------------------------------------------------+
|
||||
|Score | Description |
|
||||
+======================+===================================================+
|
||||
|absorption |Total absorption rate. This accounts for all |
|
||||
| |reactions which do not produce secondary neutrons |
|
||||
| |as well as fission. |
|
||||
|absorption |Total absorption rate. For incident neutrons, this |
|
||||
| |accounts for all reactions that do not produce |
|
||||
| |secondary neutrons as well as fission. For incident|
|
||||
| |photons, this includes photoelectric and pair |
|
||||
| |production. |
|
||||
+----------------------+---------------------------------------------------+
|
||||
|elastic |Elastic scattering reaction rate. |
|
||||
+----------------------+---------------------------------------------------+
|
||||
|
|
@ -194,6 +196,14 @@ The following tables show all valid scores:
|
|||
+----------------------+---------------------------------------------------+
|
||||
|(n,da) |(n,d\ :math:`\alpha`\ ) reaction rate. |
|
||||
+----------------------+---------------------------------------------------+
|
||||
|coherent-scatter |Coherent (Rayleigh) scattering reaction rate. |
|
||||
+----------------------+---------------------------------------------------+
|
||||
|incoherent-scatter |Incoherent (Compton) scattering reaction rate. |
|
||||
+----------------------+---------------------------------------------------+
|
||||
|photoelectric |Photoelectric absorption reaction rate. |
|
||||
+----------------------+---------------------------------------------------+
|
||||
|pair-production |Pair production reaction rate. |
|
||||
+----------------------+---------------------------------------------------+
|
||||
|*Arbitrary integer* |An arbitrary integer is interpreted to mean the |
|
||||
| |reaction rate for a reaction with a given ENDF MT |
|
||||
| |number. |
|
||||
|
|
|
|||
|
|
@ -159,11 +159,11 @@ const std::unordered_map<int, std::string> REACTION_NAME_MAP {
|
|||
{N_XA, "(n,Xa)"},
|
||||
{HEATING, "heating"},
|
||||
{DAMAGE_ENERGY, "damage-energy"},
|
||||
{COHERENT, "coherent scatter"},
|
||||
{INCOHERENT, "incoherent scatter"},
|
||||
{PAIR_PROD_ELEC, "pair production, electron"},
|
||||
{PAIR_PROD, "pair production"},
|
||||
{PAIR_PROD_NUC, "pair production, nuclear"},
|
||||
{COHERENT, "coherent-scatter"},
|
||||
{INCOHERENT, "incoherent-scatter"},
|
||||
{PAIR_PROD_ELEC, "pair-production-electron"},
|
||||
{PAIR_PROD, "pair-production"},
|
||||
{PAIR_PROD_NUC, "pair-production-nuclear"},
|
||||
{PHOTOELECTRIC, "photoelectric"},
|
||||
{N_PC, "(n,pc)"},
|
||||
{N_DC, "(n,dc)"},
|
||||
|
|
|
|||
|
|
@ -216,6 +216,14 @@ score_str_to_int(std::string score_str)
|
|||
return N_XA;
|
||||
if (score_str == "damage-energy")
|
||||
return DAMAGE_ENERGY;
|
||||
if (score_str == "coherent-scatter")
|
||||
return COHERENT;
|
||||
if (score_str == "incoherent-scatter")
|
||||
return INCOHERENT;
|
||||
if (score_str == "pair-production")
|
||||
return PAIR_PROD;
|
||||
if (score_str == "photoelectric")
|
||||
return PHOTOELECTRIC;
|
||||
|
||||
// So far we have not identified this score string. Check to see if it is a
|
||||
// deprecated score.
|
||||
|
|
|
|||
|
|
@ -530,6 +530,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
// Get the pre-collision energy of the particle.
|
||||
auto E = p.E_last_;
|
||||
|
||||
using Type = Particle::Type;
|
||||
|
||||
for (auto i = 0; i < tally.scores_.size(); ++i) {
|
||||
auto score_bin = tally.scores_[i];
|
||||
auto score_index = start_index + i;
|
||||
|
|
@ -549,8 +551,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
score = p.wgt_last_;
|
||||
}
|
||||
|
||||
if (p.type_ == Particle::Type::neutron ||
|
||||
p.type_ == Particle::Type::photon) {
|
||||
if (p.type_ == Type::neutron || p.type_ == Type::photon) {
|
||||
score *= flux / p.macro_xs_.total;
|
||||
} else {
|
||||
score = 0.;
|
||||
|
|
@ -575,9 +576,9 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
|
||||
} else {
|
||||
if (i_nuclide >= 0) {
|
||||
if (p.type_ == Particle::Type::neutron) {
|
||||
if (p.type_ == Type::neutron) {
|
||||
score = p.neutron_xs_[i_nuclide].total * atom_density * flux;
|
||||
} else if (p.type_ == Particle::Type::photon) {
|
||||
} else if (p.type_ == Type::photon) {
|
||||
score = p.photon_xs_[i_nuclide].total * atom_density * flux;
|
||||
}
|
||||
} else {
|
||||
|
|
@ -588,6 +589,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
|
||||
|
||||
case SCORE_INVERSE_VELOCITY:
|
||||
if (p.type_ != Type::neutron) continue;
|
||||
|
||||
if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
// All events score to an inverse velocity bin. We actually use a
|
||||
// collision estimator in place of an analog one since there is no way
|
||||
|
|
@ -609,6 +612,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
|
||||
|
||||
case SCORE_SCATTER:
|
||||
if (p.type_ != Type::neutron && p.type_ != Type::photon) continue;
|
||||
|
||||
if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
// Skip any event where the particle didn't scatter
|
||||
if (p.event_ != TallyEvent::SCATTER) continue;
|
||||
|
|
@ -617,17 +622,27 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
score = p.wgt_last_ * flux;
|
||||
} else {
|
||||
if (i_nuclide >= 0) {
|
||||
score = (p.neutron_xs_[i_nuclide].total
|
||||
- p.neutron_xs_[i_nuclide].absorption) * atom_density * flux;
|
||||
if (p.type_ == Type::neutron) {
|
||||
const auto& micro = p.neutron_xs_[i_nuclide];
|
||||
score = (micro.total - micro.absorption) * atom_density * flux;
|
||||
} else {
|
||||
const auto& micro = p.photon_xs_[i_nuclide];
|
||||
score = (micro.coherent + micro.incoherent) * atom_density * flux;
|
||||
}
|
||||
} else {
|
||||
score = (p.macro_xs_.total
|
||||
- p.macro_xs_.absorption) * flux;
|
||||
if (p.type_ == Type::neutron) {
|
||||
score = (p.macro_xs_.total - p.macro_xs_.absorption) * flux;
|
||||
} else {
|
||||
score = (p.macro_xs_.coherent + p.macro_xs_.incoherent) * flux;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
|
||||
|
||||
case SCORE_NU_SCATTER:
|
||||
if (p.type_ != Type::neutron) continue;
|
||||
|
||||
// Only analog estimators are available.
|
||||
// Skip any event where the particle didn't scatter
|
||||
if (p.event_ != TallyEvent::SCATTER) continue;
|
||||
|
|
@ -649,6 +664,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
|
||||
|
||||
case SCORE_ABSORPTION:
|
||||
if (p.type_ != Type::neutron && p.type_ != Type::photon) continue;
|
||||
|
||||
if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
if (settings::survival_biasing) {
|
||||
// No absorption events actually occur if survival biasing is on --
|
||||
|
|
@ -663,10 +680,18 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
}
|
||||
} else {
|
||||
if (i_nuclide >= 0) {
|
||||
score = p.neutron_xs_[i_nuclide].absorption * atom_density
|
||||
* flux;
|
||||
if (p.type_ == Type::neutron) {
|
||||
score = p.neutron_xs_[i_nuclide].absorption * atom_density * flux;
|
||||
} else {
|
||||
const auto& xs = p.photon_xs_[i_nuclide];
|
||||
score = (xs.total - xs.coherent - xs.incoherent) * atom_density * flux;
|
||||
}
|
||||
} else {
|
||||
score = p.macro_xs_.absorption * flux;
|
||||
if (p.type_ == Type::neutron) {
|
||||
score = p.macro_xs_.absorption * flux;
|
||||
} else {
|
||||
score = (p.macro_xs_.photoelectric + p.macro_xs_.pair_production) * flux;
|
||||
}
|
||||
}
|
||||
}
|
||||
break;
|
||||
|
|
@ -1205,6 +1230,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
|
||||
|
||||
case ELASTIC:
|
||||
if (p.type_ != Type::neutron) continue;
|
||||
|
||||
if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
// Check if event MT matches
|
||||
if (p.event_mt_ != ELASTIC) continue;
|
||||
|
|
@ -1244,6 +1271,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
case N_GAMMA:
|
||||
case N_P:
|
||||
case N_A:
|
||||
if (p.type_ != Type::neutron) continue;
|
||||
|
||||
if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
// Check if the event MT matches
|
||||
if (p.event_mt_ != score_bin) continue;
|
||||
|
|
@ -1277,8 +1306,46 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
break;
|
||||
|
||||
|
||||
case COHERENT:
|
||||
case INCOHERENT:
|
||||
case PHOTOELECTRIC:
|
||||
case PAIR_PROD:
|
||||
if (p.type_ != Type::photon) continue;
|
||||
|
||||
if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
if (score_bin == PHOTOELECTRIC) {
|
||||
// Photoelectric events are assigned an MT value corresponding to the
|
||||
// shell cross section. Also, photons below the energy cutoff are
|
||||
// assumed to have been absorbed via photoelectric absorption
|
||||
if ((p.event_mt_ < 534 || p.event_mt_ > 572) &&
|
||||
p.event_mt_ != REACTION_NONE) continue;
|
||||
} else {
|
||||
if (p.event_mt_ != score_bin) continue;
|
||||
}
|
||||
score = p.wgt_last_ * flux;
|
||||
} else {
|
||||
if (i_nuclide >= 0) {
|
||||
const auto& micro = p.photon_xs_[i_nuclide];
|
||||
double xs =
|
||||
(score_bin == COHERENT) ? micro.coherent :
|
||||
(score_bin == INCOHERENT) ? micro.incoherent :
|
||||
(score_bin == PHOTOELECTRIC) ? micro.photoelectric :
|
||||
micro.pair_production;
|
||||
score = xs * atom_density * flux;
|
||||
} else {
|
||||
double xs =
|
||||
(score_bin == COHERENT) ? p.macro_xs_.coherent :
|
||||
(score_bin == INCOHERENT) ? p.macro_xs_.incoherent :
|
||||
(score_bin == PHOTOELECTRIC) ? p.macro_xs_.photoelectric :
|
||||
p.macro_xs_.pair_production;
|
||||
score = xs * flux;
|
||||
}
|
||||
}
|
||||
break;
|
||||
|
||||
|
||||
case HEATING:
|
||||
if (p.type_ == Particle::Type::neutron) {
|
||||
if (p.type_ == Type::neutron) {
|
||||
score = score_neutron_heating(p, tally, flux, HEATING,
|
||||
i_nuclide, atom_density);
|
||||
} else {
|
||||
|
|
@ -1301,7 +1368,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index,
|
|||
|
||||
// The default block is really only meant for redundant neutron reactions
|
||||
// (e.g. 444, 901)
|
||||
if (p.type_ != Particle::Type::neutron) continue;
|
||||
if (p.type_ != Type::neutron) continue;
|
||||
|
||||
if (tally.estimator_ == TallyEstimator::ANALOG) {
|
||||
|
||||
|
|
|
|||
|
|
@ -50,19 +50,19 @@
|
|||
<tally id="2">
|
||||
<filters>2</filters>
|
||||
<nuclides>Al27 total</nuclides>
|
||||
<scores>total</scores>
|
||||
<scores>total (n,gamma)</scores>
|
||||
<estimator>tracklength</estimator>
|
||||
</tally>
|
||||
<tally id="3">
|
||||
<filters>2</filters>
|
||||
<nuclides>Al27 total</nuclides>
|
||||
<scores>total heating</scores>
|
||||
<scores>total heating (n,gamma)</scores>
|
||||
<estimator>collision</estimator>
|
||||
</tally>
|
||||
<tally id="4">
|
||||
<filters>2</filters>
|
||||
<nuclides>Al27 total</nuclides>
|
||||
<scores>total heating</scores>
|
||||
<scores>total heating (n,gamma)</scores>
|
||||
<estimator>analog</estimator>
|
||||
</tally>
|
||||
</tallies>
|
||||
|
|
|
|||
|
|
@ -10,10 +10,16 @@ tally 1:
|
|||
tally 2:
|
||||
1.249805E+00
|
||||
1.562013E+00
|
||||
4.284760E-04
|
||||
1.835917E-07
|
||||
1.249805E+00
|
||||
1.562013E+00
|
||||
4.284760E-04
|
||||
1.835917E-07
|
||||
8.281718E-01
|
||||
6.858685E-01
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
8.281718E-01
|
||||
6.858685E-01
|
||||
0.000000E+00
|
||||
|
|
@ -24,52 +30,84 @@ tally 2:
|
|||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
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.248100E+00
|
||||
1.557754E+00
|
||||
2.024868E+06
|
||||
4.100089E+12
|
||||
4.283447E-04
|
||||
1.834792E-07
|
||||
1.248100E+00
|
||||
1.557754E+00
|
||||
2.024868E+06
|
||||
4.100089E+12
|
||||
8.242000E-01
|
||||
6.793056E-01
|
||||
5.827819E+03
|
||||
3.396348E+07
|
||||
4.283447E-04
|
||||
1.834792E-07
|
||||
8.242000E-01
|
||||
6.793056E-01
|
||||
5.827819E+03
|
||||
3.396348E+07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
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
|
||||
1.834863E+05
|
||||
3.366722E+10
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
5.938849E+03
|
||||
3.526993E+07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
5.938849E+03
|
||||
3.526993E+07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
tally 4:
|
||||
1.248100E+00
|
||||
1.557754E+00
|
||||
2.024868E+06
|
||||
4.100089E+12
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
1.248100E+00
|
||||
1.557754E+00
|
||||
2.024868E+06
|
||||
4.100089E+12
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
2.308000E-01
|
||||
5.326864E-02
|
||||
5.822937E+03
|
||||
3.390660E+07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
8.242000E-01
|
||||
6.793056E-01
|
||||
5.827819E+03
|
||||
|
|
@ -80,6 +118,10 @@ tally 4:
|
|||
0.000000E+00
|
||||
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
|
||||
|
|
@ -88,5 +130,11 @@ tally 4:
|
|||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
5.938849E+03
|
||||
3.526993E+07
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
|
|
|
|||
|
|
@ -48,17 +48,17 @@ def model():
|
|||
current_tally.scores = ['current']
|
||||
tally_tracklength = openmc.Tally()
|
||||
tally_tracklength.filters = [particle_filter]
|
||||
tally_tracklength.scores = ['total'] # heating doesn't work with tracklength
|
||||
tally_tracklength.scores = ['total', '(n,gamma)'] # heating doesn't work with tracklength
|
||||
tally_tracklength.nuclides = ['Al27', 'total']
|
||||
tally_tracklength.estimator = 'tracklength'
|
||||
tally_collision = openmc.Tally()
|
||||
tally_collision.filters = [particle_filter]
|
||||
tally_collision.scores = ['total', 'heating']
|
||||
tally_collision.scores = ['total', 'heating', '(n,gamma)']
|
||||
tally_collision.nuclides = ['Al27', 'total']
|
||||
tally_collision.estimator = 'collision'
|
||||
tally_analog = openmc.Tally()
|
||||
tally_analog.filters = [particle_filter]
|
||||
tally_analog.scores = ['total', 'heating']
|
||||
tally_analog.scores = ['total', 'heating', '(n,gamma)']
|
||||
tally_analog.nuclides = ['Al27', 'total']
|
||||
tally_analog.estimator = 'analog'
|
||||
model.tallies.extend([current_tally, tally_tracklength,
|
||||
|
|
|
|||
|
|
@ -40,6 +40,6 @@
|
|||
</filter>
|
||||
<tally id="1">
|
||||
<filters>1</filters>
|
||||
<scores>flux</scores>
|
||||
<scores>flux (n,gamma)</scores>
|
||||
</tally>
|
||||
</tallies>
|
||||
|
|
|
|||
|
|
@ -1,3 +1,5 @@
|
|||
tally 1:
|
||||
sum = 2.275713E+02
|
||||
sum_sq = 5.178870E+04
|
||||
2.275713E+02
|
||||
5.178870E+04
|
||||
0.000000E+00
|
||||
0.000000E+00
|
||||
|
|
|
|||
|
|
@ -41,20 +41,10 @@ class SourceTestHarness(PyAPITestHarness):
|
|||
particle_filter = openmc.ParticleFilter('photon')
|
||||
tally = openmc.Tally()
|
||||
tally.filters = [particle_filter]
|
||||
tally.scores = ['flux']
|
||||
tally.scores = ['flux', '(n,gamma)']
|
||||
tallies = openmc.Tallies([tally])
|
||||
tallies.export_to_xml()
|
||||
|
||||
def _get_results(self):
|
||||
with openmc.StatePoint(self._sp_name) as sp:
|
||||
outstr = ''
|
||||
t = sp.get_tally()
|
||||
outstr += 'tally {}:\n'.format(t.id)
|
||||
outstr += 'sum = {:12.6E}\n'.format(t.sum[0, 0, 0])
|
||||
outstr += 'sum_sq = {:12.6E}\n'.format(t.sum_sq[0, 0, 0])
|
||||
|
||||
return outstr
|
||||
|
||||
|
||||
def test_photon_source():
|
||||
harness = SourceTestHarness('statepoint.1.h5')
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue