diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 75ca142ba..69e8506d9 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -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. | diff --git a/src/reaction.cpp b/src/reaction.cpp index fbe2c34de..b3559b10c 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -159,11 +159,11 @@ const std::unordered_map 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)"}, diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 924975358..c025755d1 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -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. diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 049c7be6b..3ca005cfa 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -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) { diff --git a/tests/regression_tests/photon_production/inputs_true.dat b/tests/regression_tests/photon_production/inputs_true.dat index 6a6106cbe..09f1fd290 100644 --- a/tests/regression_tests/photon_production/inputs_true.dat +++ b/tests/regression_tests/photon_production/inputs_true.dat @@ -50,19 +50,19 @@ 2 Al27 total - total + total (n,gamma) tracklength 2 Al27 total - total heating + total heating (n,gamma) collision 2 Al27 total - total heating + total heating (n,gamma) analog diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index 35e6e5602..2be99e1e0 100644 --- a/tests/regression_tests/photon_production/results_true.dat +++ b/tests/regression_tests/photon_production/results_true.dat @@ -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 diff --git a/tests/regression_tests/photon_production/test.py b/tests/regression_tests/photon_production/test.py index f56e83535..107b53bbb 100644 --- a/tests/regression_tests/photon_production/test.py +++ b/tests/regression_tests/photon_production/test.py @@ -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, diff --git a/tests/regression_tests/photon_source/inputs_true.dat b/tests/regression_tests/photon_source/inputs_true.dat index 89f4de0e0..f7c9c24b4 100644 --- a/tests/regression_tests/photon_source/inputs_true.dat +++ b/tests/regression_tests/photon_source/inputs_true.dat @@ -40,6 +40,6 @@ 1 - flux + flux (n,gamma) diff --git a/tests/regression_tests/photon_source/results_true.dat b/tests/regression_tests/photon_source/results_true.dat index 1448ad6db..fe069f6cd 100644 --- a/tests/regression_tests/photon_source/results_true.dat +++ b/tests/regression_tests/photon_source/results_true.dat @@ -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 diff --git a/tests/regression_tests/photon_source/test.py b/tests/regression_tests/photon_source/test.py index 94b028003..a2bfb9038 100644 --- a/tests/regression_tests/photon_source/test.py +++ b/tests/regression_tests/photon_source/test.py @@ -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')