13703779155

Committed 06 Mar 2025 04:45PM UTC coverage: 85.674% (+0.5%) from 85.129%

Build # 13703779155

Build Type

Pull #3087

github

Committed by

web-flow

Commit Message

Merge 20c292e2a into e360cb467

Pull Request Pull Request #3087: wheel building with scikit build core

Run Details

47261 of 55164 relevant lines covered (85.67%)

29624867.19 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

63.64

/src/reaction.cpp

#include "openmc/reaction.h"

#include <algorithm> // for remove_if
#include <string>
#include <unordered_map>
#include <utility> // for move

#include <fmt/core.h>

#include "openmc/chain.h"
#include "openmc/constants.h"
#include "openmc/endf.h"
#include "openmc/hdf5_interface.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/secondary_uncorrelated.h"
#include "openmc/settings.h"

namespace openmc {

//==============================================================================
// Reaction implementation
//==============================================================================

Reaction::Reaction(
  hid_t group, const vector<int>& temperatures, std::string name)
{
  read_attribute(group, "Q_value", q_value_);
  read_attribute(group, "mt", mt_);
  int tmp;
  read_attribute(group, "center_of_mass", tmp);
  scatter_in_cm_ = (tmp == 1);

  // Checks if redudant attribute exists before loading
  // (for compatibiltiy with legacy .h5 libraries)
  if (attribute_exists(group, "redundant")) {
    read_attribute(group, "redundant", tmp);
    redundant_ = (tmp == 1);
  } else {
    redundant_ = false;
  }

  // Read cross section and threshold_idx data
  for (auto t : temperatures) {
    // Get group corresponding to temperature
    hid_t temp_group = open_group(group, fmt::format("{}K", t).c_str());
    hid_t dset = open_dataset(temp_group, "xs");

    // Get threshold index
    TemperatureXS xs;
    read_attribute(dset, "threshold_idx", xs.threshold);

    // Read cross section values
    read_dataset(dset, xs.value);
    close_dataset(dset);
    close_group(temp_group);

    // create new entry in xs vector
    xs_.push_back(std::move(xs));
  }

  // Read products
  for (const auto& name : group_names(group)) {
    if (name.rfind("product_", 0) == 0) {
      hid_t pgroup = open_group(group, name.c_str());
      products_.emplace_back(pgroup);
      close_group(pgroup);
    }
  }

  if (settings::use_decay_photons) {
    // Remove photon products for D1S method
    products_.erase(
      std::remove_if(products_.begin(), products_.end(),
        [](const auto& p) { return p.particle_ == ParticleType::photon; }),
      products_.end());

    // Determine product for D1S method
    auto nuclide_it = data::chain_nuclide_map.find(name);
    if (nuclide_it != data::chain_nuclide_map.end()) {
      const auto& chain_nuc = data::chain_nuclides[nuclide_it->second];
      const auto& rx_products = chain_nuc->reaction_products();
      auto product_it = rx_products.find(mt_);
      if (product_it != rx_products.end()) {
        auto decay_products = product_it->second;
        for (const auto& decay_product : decay_products) {
          auto product_it = data::chain_nuclide_map.find(decay_product.name);
          if (product_it != data::chain_nuclide_map.end()) {
            const auto& product_nuc = data::chain_nuclides[product_it->second];
            if (product_nuc->photon_energy()) {
              products_.emplace_back(decay_product);
            }
          }
        }
      }
    }
  }
}

double Reaction::xs(int64_t i_temp, int64_t i_grid, double interp_factor) const
{
  // If energy is below threshold, return 0. Otherwise interpolate between
  // nearest grid points
  const auto& x = xs_[i_temp];
  return (i_grid < x.threshold)
           ? 0.0
           : (1.0 - interp_factor) * x.value[i_grid - x.threshold] +
               interp_factor * x.value[i_grid - x.threshold + 1];
}

double Reaction::xs(const NuclideMicroXS& micro) const
{
  return this->xs(micro.index_temp, micro.index_grid, micro.interp_factor);
}

double Reaction::collapse_rate(int64_t i_temp, span<const double> energy,
  span<const double> flux, const vector<double>& grid) const
{
  // Find index corresponding to first energy
  const auto& xs = xs_[i_temp].value;
  int i_low = lower_bound_index(grid.cbegin(), grid.cend(), energy.front());

  // Check for threshold and adjust starting point if necessary
  int j_start = 0;
  int i_threshold = xs_[i_temp].threshold;
  if (i_low < i_threshold) {
    i_low = i_threshold;
    while (energy[j_start + 1] < grid[i_low]) {
      ++j_start;
      if (j_start + 1 == energy.size())
        return 0.0;
    }
  }

  double xs_flux_sum = 0.0;

  for (int j = j_start; j < flux.size(); ++j) {
    double E_group_low = energy[j];
    double E_group_high = energy[j + 1];
    double flux_per_eV = flux[j] / (E_group_high - E_group_low);

    // Determine energy grid index corresponding to group high
    int i_high = i_low;
    while (grid[i_high + 1] < E_group_high && i_high + 1 < grid.size() - 1)
      ++i_high;

    // Loop over energy grid points within [E_group_low, E_group_high]
    for (; i_low <= i_high; ++i_low) {
      // Determine bounding grid energies and cross sections
      double E_l = grid[i_low];
      double E_r = grid[i_low + 1];
      if (E_l == E_r)
        continue;

      double xs_l = xs[i_low - i_threshold];
      double xs_r = xs[i_low + 1 - i_threshold];

      // Determine actual energies
      double E_low = std::max(E_group_low, E_l);
      double E_high = std::min(E_group_high, E_r);

      // Determine average cross section across segment
      double m = (xs_r - xs_l) / (E_r - E_l);
      double xs_low = xs_l + m * (E_low - E_l);
      double xs_high = xs_l + m * (E_high - E_l);
      double xs_avg = 0.5 * (xs_low + xs_high);

      // Add contribution from segment
      double dE = (E_high - E_low);
      xs_flux_sum += flux_per_eV * xs_avg * dE;
    }

    i_low = i_high;

    // Check for end of energy grid
    if (i_low + 1 == grid.size())
      break;
  }

  return xs_flux_sum;
}

//==============================================================================
// Non-member functions
//==============================================================================

std::unordered_map<int, std::string> REACTION_NAME_MAP {
  {SCORE_FLUX, "flux"},
  {SCORE_TOTAL, "total"},
  {SCORE_SCATTER, "scatter"},
  {SCORE_NU_SCATTER, "nu-scatter"},
  {SCORE_ABSORPTION, "absorption"},
  {SCORE_FISSION, "fission"},
  {SCORE_NU_FISSION, "nu-fission"},
  {SCORE_DECAY_RATE, "decay-rate"},
  {SCORE_DELAYED_NU_FISSION, "delayed-nu-fission"},
  {SCORE_PROMPT_NU_FISSION, "prompt-nu-fission"},
  {SCORE_KAPPA_FISSION, "kappa-fission"},
  {SCORE_CURRENT, "current"},
  {SCORE_EVENTS, "events"},
  {SCORE_INVERSE_VELOCITY, "inverse-velocity"},
  {SCORE_FISS_Q_PROMPT, "fission-q-prompt"},
  {SCORE_FISS_Q_RECOV, "fission-q-recoverable"},
  {SCORE_PULSE_HEIGHT, "pulse-height"},
  // Normal ENDF-based reactions
  {TOTAL_XS, "(n,total)"},
  {ELASTIC, "(n,elastic)"},
  {N_NONELASTIC, "(n,nonelastic)"},
  {N_LEVEL, "(n,level)"},
  {N_2ND, "(n,2nd)"},
  {N_2N, "(n,2n)"},
  {N_3N, "(n,3n)"},
  {N_FISSION, "(n,fission)"},
  {N_F, "(n,f)"},
  {N_NF, "(n,nf)"},
  {N_2NF, "(n,2nf)"},
  {N_NA, "(n,na)"},
  {N_N3A, "(n,n3a)"},
  {N_2NA, "(n,2na)"},
  {N_3NA, "(n,3na)"},
  {N_NP, "(n,np)"},
  {N_N2A, "(n,n2a)"},
  {N_2N2A, "(n,2n2a)"},
  {N_ND, "(n,nd)"},
  {N_NT, "(n,nt)"},
  {N_N3HE, "(n,n3He)"},
  {N_ND2A, "(n,nd2a)"},
  {N_NT2A, "(n,nt2a)"},
  {N_4N, "(n,4n)"},
  {N_3NF, "(n,3nf)"},
  {N_2NP, "(n,2np)"},
  {N_3NP, "(n,3np)"},
  {N_N2P, "(n,n2p)"},
  {N_NPA, "(n,npa)"},
  {N_NC, "(n,nc)"},
  {N_DISAPPEAR, "(n,disappear)"},
  {N_GAMMA, "(n,gamma)"},
  {N_P, "(n,p)"},
  {N_D, "(n,d)"},
  {N_T, "(n,t)"},
  {N_3HE, "(n,3He)"},
  {N_A, "(n,a)"},
  {N_2A, "(n,2a)"},
  {N_3A, "(n,3a)"},
  {N_2P, "(n,2p)"},
  {N_PA, "(n,pa)"},
  {N_T2A, "(n,t2a)"},
  {N_D2A, "(n,d2a)"},
  {N_PD, "(n,pd)"},
  {N_PT, "(n,pt)"},
  {N_DA, "(n,da)"},
  {N_5N, "(n,5n)"},
  {N_6N, "(n,6n)"},
  {N_2NT, "(n,2nt)"},
  {N_TA, "(n,ta)"},
  {N_4NP, "(n,4np)"},
  {N_3ND, "(n,3nd)"},
  {N_NDA, "(n,nda)"},
  {N_2NPA, "(n,2npa)"},
  {N_7N, "(n,7n)"},
  {N_8N, "(n,8n)"},
  {N_5NP, "(n,5np)"},
  {N_6NP, "(n,6np)"},
  {N_7NP, "(n,7np)"},
  {N_4NA, "(n,4na)"},
  {N_5NA, "(n,5na)"},
  {N_6NA, "(n,6na)"},
  {N_7NA, "(n,7na)"},
  {N_4ND, "(n,4nd)"},
  {N_5ND, "(n,5nd)"},
  {N_6ND, "(n,6nd)"},
  {N_3NT, "(n,3nt)"},
  {N_4NT, "(n,4nt)"},
  {N_5NT, "(n,5nt)"},
  {N_6NT, "(n,6nt)"},
  {N_2N3HE, "(n,2n3He)"},
  {N_3N3HE, "(n,3n3He)"},
  {N_4N3HE, "(n,4n3He)"},
  {N_3N2P, "(n,3n2p)"},
  {N_3N2A, "(n,3n2a)"},
  {N_3NPA, "(n,3npa)"},
  {N_DT, "(n,dt)"},
  {N_NPD, "(n,npd)"},
  {N_NPT, "(n,npt)"},
  {N_NDT, "(n,ndt)"},
  {N_NP3HE, "(n,np3He)"},
  {N_ND3HE, "(n,nd3He)"},
  {N_NT3HE, "(n,nt3He)"},
  {N_NTA, "(n,nta)"},
  {N_2N2P, "(n,2n2p)"},
  {N_P3HE, "(n,p3He)"},
  {N_D3HE, "(n,d3He)"},
  {N_3HEA, "(n,3Hea)"},
  {N_4N2P, "(n,4n2p)"},
  {N_4N2A, "(n,4n2a)"},
  {N_4NPA, "(n,4npa)"},
  {N_3P, "(n,3p)"},
  {N_N3P, "(n,n3p)"},
  {N_3N2PA, "(n,3n2pa)"},
  {N_5N2P, "(n,5n2p)"},
  {201, "(n,Xn)"},
  {202, "(n,Xgamma)"},
  {N_XP, "(n,Xp)"},
  {N_XD, "(n,Xd)"},
  {N_XT, "(n,Xt)"},
  {N_X3HE, "(n,X3He)"},
  {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"},
  {PHOTOELECTRIC, "photoelectric"},
  {N_PC, "(n,pc)"},
  {N_DC, "(n,dc)"},
  {N_TC, "(n,tc)"},
  {N_3HEC, "(n,3Hec)"},
  {N_AC, "(n,ac)"},
  {N_2NC, "(n,2nc)"},
  {HEATING_LOCAL, "heating-local"},
};

std::unordered_map<std::string, int> REACTION_TYPE_MAP;

void initialize_maps()
{
  // Add level reactions to name map
  for (int level = 0; level <= 48; ++level) {
    if (level >= 1 && level <= 40) {
      REACTION_NAME_MAP[50 + level] = fmt::format("(n,n{})", level);
    }
    REACTION_NAME_MAP[600 + level] = fmt::format("(n,p{})", level);
    REACTION_NAME_MAP[650 + level] = fmt::format("(n,d{})", level);
    REACTION_NAME_MAP[700 + level] = fmt::format("(n,t{})", level);
    REACTION_NAME_MAP[750 + level] = fmt::format("(n,3He{})", level);
    REACTION_NAME_MAP[800 + level] = fmt::format("(n,a{})", level);
    if (level <= 15) {
      REACTION_NAME_MAP[875 + level] = fmt::format("(n,2n{})", level);
    }
  }

  // Create photoelectric subshells
  for (int mt = 534; mt <= 572; ++mt) {
    REACTION_NAME_MAP[mt] =
      fmt::format("photoelectric, {} subshell", SUBSHELLS[mt - 534]);
  }

  // Invert name map to create type map
  for (const auto& kv : REACTION_NAME_MAP) {
    REACTION_TYPE_MAP[kv.second] = kv.first;
  }
}

std::string reaction_name(int mt)
{
  // Initialize remainder of name map and all of type map
  if (REACTION_TYPE_MAP.empty())
    initialize_maps();

  // Get reaction name from map
  auto it = REACTION_NAME_MAP.find(mt);
  if (it != REACTION_NAME_MAP.end()) {
    return it->second;
  } else {
    return fmt::format("MT={}", mt);
  }
}

int reaction_type(std::string name)
{
  // Initialize remainder of name map and all of type map
  if (REACTION_TYPE_MAP.empty())
    initialize_maps();

  // (n,total) exists in REACTION_TYPE_MAP for MT=1, but we need this to return
  // the special SCORE_TOTAL score
  if (name == "(n,total)")
    return SCORE_TOTAL;

  // Check if type map has an entry for this reaction name
  auto it = REACTION_TYPE_MAP.find(name);
  if (it != REACTION_TYPE_MAP.end()) {
    return it->second;
  }

  // Alternate names for several reactions
  if (name == "elastic") {
    return ELASTIC;
  } else if (name == "n2n") {
    return N_2N;
  } else if (name == "n3n") {
    return N_3N;
  } else if (name == "n4n") {
    return N_4N;
  } else if (name == "H1-production") {
    return N_XP;
  } else if (name == "H2-production") {
    return N_XD;
  } else if (name == "H3-production") {
    return N_XT;
  } else if (name == "He3-production") {
    return N_X3HE;
  } else if (name == "He4-production") {
    return N_XA;
  }

  // Assume the given string is a reaction MT number.  Make sure it's a natural
  // number then return.
  int MT = 0;
  try {
    MT = std::stoi(name);
  } catch (const std::invalid_argument& ex) {
    throw std::invalid_argument(
      "Invalid tally score \"" + name +
      "\". See the docs "
      "for details: "
      "https://docs.openmc.org/en/stable/usersguide/tallies.html#scores");
  }
  if (MT < 1)
    throw std::invalid_argument(
      "Invalid tally score \"" + name +
      "\". See the docs "
      "for details: "
      "https://docs.openmc.org/en/stable/usersguide/tallies.html#scores");
  return MT;
}

} // namespace openmc

1	#include "openmc/reaction.h"
2
3	#include <algorithm> // for remove_if
4	#include <string>
5	#include <unordered_map>
6	#include <utility> // for move
7
8	#include <fmt/core.h>
9
10	#include "openmc/chain.h"
11	#include "openmc/constants.h"
12	#include "openmc/endf.h"
13	#include "openmc/hdf5_interface.h"
14	#include "openmc/random_lcg.h"
15	#include "openmc/search.h"
16	#include "openmc/secondary_uncorrelated.h"
17	#include "openmc/settings.h"
18
19	namespace openmc {
20
21	//==============================================================================
22	// Reaction implementation
23	//==============================================================================
24
25	Reaction::Reaction(	728,698✔
26	hid_t group, const vector<int>& temperatures, std::string name)	728,698✔
27	{
28	read_attribute(group, "Q_value", q_value_);	728,698✔
29	read_attribute(group, "mt", mt_);	728,698✔
30	int tmp;
31	read_attribute(group, "center_of_mass", tmp);	728,698✔
32	scatter_in_cm_ = (tmp == 1);	728,698✔
33
34	// Checks if redudant attribute exists before loading
35	// (for compatibiltiy with legacy .h5 libraries)
36	if (attribute_exists(group, "redundant")) {	728,698✔
37	read_attribute(group, "redundant", tmp);	728,698✔
38	redundant_ = (tmp == 1);	728,698✔
39	} else {
40	redundant_ = false;	×
41	}
42
43	// Read cross section and threshold_idx data
44	for (auto t : temperatures) {	1,457,627✔
45	// Get group corresponding to temperature
46	hid_t temp_group = open_group(group, fmt::format("{}K", t).c_str());	728,929✔
47	hid_t dset = open_dataset(temp_group, "xs");	728,929✔
48
49	// Get threshold index
50	TemperatureXS xs;	728,929✔
51	read_attribute(dset, "threshold_idx", xs.threshold);	728,929✔
52
53	// Read cross section values
54	read_dataset(dset, xs.value);	728,929✔
55	close_dataset(dset);	728,929✔
56	close_group(temp_group);	728,929✔
57
58	// create new entry in xs vector
59	xs_.push_back(std::move(xs));	728,929✔
60	}	728,929✔
61
62	// Read products
63	for (const auto& name : group_names(group)) {	3,306,753✔
64	if (name.rfind("product_", 0) == 0) {	2,578,055✔
65	hid_t pgroup = open_group(group, name.c_str());	1,845,113✔
66	products_.emplace_back(pgroup);	1,845,113✔
67	close_group(pgroup);	1,845,113✔
68	}
69	}	728,698✔
70
71	if (settings::use_decay_photons) {	728,698✔
72	// Remove photon products for D1S method
73	products_.erase(	2,486✔
74	std::remove_if(products_.begin(), products_.end(),	1,243✔
75	[](const auto& p) { return p.particle_ == ParticleType::photon; }),	8,866✔
76	products_.end());	1,243✔
77
78	// Determine product for D1S method
79	auto nuclide_it = data::chain_nuclide_map.find(name);	1,243✔
80	if (nuclide_it != data::chain_nuclide_map.end()) {	1,243✔
81	const auto& chain_nuc = data::chain_nuclides[nuclide_it->second];	1,243✔
82	const auto& rx_products = chain_nuc->reaction_products();	1,243✔
83	auto product_it = rx_products.find(mt_);	1,243✔
84	if (product_it != rx_products.end()) {	1,243✔
85	auto decay_products = product_it->second;	220✔
86	for (const auto& decay_product : decay_products) {	440✔
87	auto product_it = data::chain_nuclide_map.find(decay_product.name);	220✔
88	if (product_it != data::chain_nuclide_map.end()) {	220✔
89	const auto& product_nuc = data::chain_nuclides[product_it->second];	220✔
90	if (product_nuc->photon_energy()) {	220✔
91	products_.emplace_back(decay_product);	132✔
92	}
93	}
94	}
95	}	220✔
96	}
97	}
98	}	728,698✔
99
100	double Reaction::xs(int64_t i_temp, int64_t i_grid, double interp_factor) const	537,419,131✔
101	{
102	// If energy is below threshold, return 0. Otherwise interpolate between
103	// nearest grid points
104	const auto& x = xs_[i_temp];	537,419,131✔
105	return (i_grid < x.threshold)	537,419,131✔
106	? 0.0	537,419,131✔
107	: (1.0 - interp_factor) * x.value[i_grid - x.threshold] +	326,280,582✔
108	interp_factor * x.value[i_grid - x.threshold + 1];	537,419,131✔
109	}
110
111	double Reaction::xs(const NuclideMicroXS& micro) const	537,419,131✔
112	{
113	return this->xs(micro.index_temp, micro.index_grid, micro.interp_factor);	537,419,131✔
114	}
115
116	double Reaction::collapse_rate(int64_t i_temp, span<const double> energy,	×
117	span<const double> flux, const vector<double>& grid) const
118	{
119	// Find index corresponding to first energy
120	const auto& xs = xs_[i_temp].value;	×
121	int i_low = lower_bound_index(grid.cbegin(), grid.cend(), energy.front());	×
122
123	// Check for threshold and adjust starting point if necessary
124	int j_start = 0;	×
125	int i_threshold = xs_[i_temp].threshold;	×
126	if (i_low < i_threshold) {	×
127	i_low = i_threshold;	×
128	while (energy[j_start + 1] < grid[i_low]) {	×
129	++j_start;	×
130	if (j_start + 1 == energy.size())	×
131	return 0.0;	×
132	}
133	}
134
135	double xs_flux_sum = 0.0;	×
136
137	for (int j = j_start; j < flux.size(); ++j) {	×
138	double E_group_low = energy[j];	×
139	double E_group_high = energy[j + 1];	×
140	double flux_per_eV = flux[j] / (E_group_high - E_group_low);	×
141
142	// Determine energy grid index corresponding to group high
143	int i_high = i_low;	×
144	while (grid[i_high + 1] < E_group_high && i_high + 1 < grid.size() - 1)	×
145	++i_high;	×
146
147	// Loop over energy grid points within [E_group_low, E_group_high]
148	for (; i_low <= i_high; ++i_low) {	×
149	// Determine bounding grid energies and cross sections
150	double E_l = grid[i_low];	×
151	double E_r = grid[i_low + 1];	×
152	if (E_l == E_r)	×
153	continue;	×
154
155	double xs_l = xs[i_low - i_threshold];	×
156	double xs_r = xs[i_low + 1 - i_threshold];	×
157
158	// Determine actual energies
159	double E_low = std::max(E_group_low, E_l);	×
160	double E_high = std::min(E_group_high, E_r);	×
161
162	// Determine average cross section across segment
163	double m = (xs_r - xs_l) / (E_r - E_l);	×
164	double xs_low = xs_l + m * (E_low - E_l);	×
165	double xs_high = xs_l + m * (E_high - E_l);	×
166	double xs_avg = 0.5 * (xs_low + xs_high);	×
167
168	// Add contribution from segment
169	double dE = (E_high - E_low);	×
170	xs_flux_sum += flux_per_eV * xs_avg * dE;	×
171	}
172
173	i_low = i_high;	×
174
175	// Check for end of energy grid
176	if (i_low + 1 == grid.size())	×
177	break;	×
178	}
179
180	return xs_flux_sum;	×
181	}
182
183	//==============================================================================
184	// Non-member functions
185	//==============================================================================
186
187	std::unordered_map<int, std::string> REACTION_NAME_MAP {
188	{SCORE_FLUX, "flux"},
189	{SCORE_TOTAL, "total"},
190	{SCORE_SCATTER, "scatter"},
191	{SCORE_NU_SCATTER, "nu-scatter"},
192	{SCORE_ABSORPTION, "absorption"},
193	{SCORE_FISSION, "fission"},
194	{SCORE_NU_FISSION, "nu-fission"},
195	{SCORE_DECAY_RATE, "decay-rate"},
196	{SCORE_DELAYED_NU_FISSION, "delayed-nu-fission"},
197	{SCORE_PROMPT_NU_FISSION, "prompt-nu-fission"},
198	{SCORE_KAPPA_FISSION, "kappa-fission"},
199	{SCORE_CURRENT, "current"},
200	{SCORE_EVENTS, "events"},
201	{SCORE_INVERSE_VELOCITY, "inverse-velocity"},
202	{SCORE_FISS_Q_PROMPT, "fission-q-prompt"},
203	{SCORE_FISS_Q_RECOV, "fission-q-recoverable"},
204	{SCORE_PULSE_HEIGHT, "pulse-height"},
205	// Normal ENDF-based reactions
206	{TOTAL_XS, "(n,total)"},
207	{ELASTIC, "(n,elastic)"},
208	{N_NONELASTIC, "(n,nonelastic)"},
209	{N_LEVEL, "(n,level)"},
210	{N_2ND, "(n,2nd)"},
211	{N_2N, "(n,2n)"},
212	{N_3N, "(n,3n)"},
213	{N_FISSION, "(n,fission)"},
214	{N_F, "(n,f)"},
215	{N_NF, "(n,nf)"},
216	{N_2NF, "(n,2nf)"},
217	{N_NA, "(n,na)"},
218	{N_N3A, "(n,n3a)"},
219	{N_2NA, "(n,2na)"},
220	{N_3NA, "(n,3na)"},
221	{N_NP, "(n,np)"},
222	{N_N2A, "(n,n2a)"},
223	{N_2N2A, "(n,2n2a)"},
224	{N_ND, "(n,nd)"},
225	{N_NT, "(n,nt)"},
226	{N_N3HE, "(n,n3He)"},
227	{N_ND2A, "(n,nd2a)"},
228	{N_NT2A, "(n,nt2a)"},
229	{N_4N, "(n,4n)"},
230	{N_3NF, "(n,3nf)"},
231	{N_2NP, "(n,2np)"},
232	{N_3NP, "(n,3np)"},
233	{N_N2P, "(n,n2p)"},
234	{N_NPA, "(n,npa)"},
235	{N_NC, "(n,nc)"},
236	{N_DISAPPEAR, "(n,disappear)"},
237	{N_GAMMA, "(n,gamma)"},
238	{N_P, "(n,p)"},
239	{N_D, "(n,d)"},
240	{N_T, "(n,t)"},
241	{N_3HE, "(n,3He)"},
242	{N_A, "(n,a)"},
243	{N_2A, "(n,2a)"},
244	{N_3A, "(n,3a)"},
245	{N_2P, "(n,2p)"},
246	{N_PA, "(n,pa)"},
247	{N_T2A, "(n,t2a)"},
248	{N_D2A, "(n,d2a)"},
249	{N_PD, "(n,pd)"},
250	{N_PT, "(n,pt)"},
251	{N_DA, "(n,da)"},
252	{N_5N, "(n,5n)"},
253	{N_6N, "(n,6n)"},
254	{N_2NT, "(n,2nt)"},
255	{N_TA, "(n,ta)"},
256	{N_4NP, "(n,4np)"},
257	{N_3ND, "(n,3nd)"},
258	{N_NDA, "(n,nda)"},
259	{N_2NPA, "(n,2npa)"},
260	{N_7N, "(n,7n)"},
261	{N_8N, "(n,8n)"},
262	{N_5NP, "(n,5np)"},
263	{N_6NP, "(n,6np)"},
264	{N_7NP, "(n,7np)"},
265	{N_4NA, "(n,4na)"},
266	{N_5NA, "(n,5na)"},
267	{N_6NA, "(n,6na)"},
268	{N_7NA, "(n,7na)"},
269	{N_4ND, "(n,4nd)"},
270	{N_5ND, "(n,5nd)"},
271	{N_6ND, "(n,6nd)"},
272	{N_3NT, "(n,3nt)"},
273	{N_4NT, "(n,4nt)"},
274	{N_5NT, "(n,5nt)"},
275	{N_6NT, "(n,6nt)"},
276	{N_2N3HE, "(n,2n3He)"},
277	{N_3N3HE, "(n,3n3He)"},
278	{N_4N3HE, "(n,4n3He)"},
279	{N_3N2P, "(n,3n2p)"},
280	{N_3N2A, "(n,3n2a)"},
281	{N_3NPA, "(n,3npa)"},
282	{N_DT, "(n,dt)"},
283	{N_NPD, "(n,npd)"},
284	{N_NPT, "(n,npt)"},
285	{N_NDT, "(n,ndt)"},
286	{N_NP3HE, "(n,np3He)"},
287	{N_ND3HE, "(n,nd3He)"},
288	{N_NT3HE, "(n,nt3He)"},
289	{N_NTA, "(n,nta)"},
290	{N_2N2P, "(n,2n2p)"},
291	{N_P3HE, "(n,p3He)"},
292	{N_D3HE, "(n,d3He)"},
293	{N_3HEA, "(n,3Hea)"},
294	{N_4N2P, "(n,4n2p)"},
295	{N_4N2A, "(n,4n2a)"},
296	{N_4NPA, "(n,4npa)"},
297	{N_3P, "(n,3p)"},
298	{N_N3P, "(n,n3p)"},
299	{N_3N2PA, "(n,3n2pa)"},
300	{N_5N2P, "(n,5n2p)"},
301	{201, "(n,Xn)"},
302	{202, "(n,Xgamma)"},
303	{N_XP, "(n,Xp)"},
304	{N_XD, "(n,Xd)"},
305	{N_XT, "(n,Xt)"},
306	{N_X3HE, "(n,X3He)"},
307	{N_XA, "(n,Xa)"},
308	{HEATING, "heating"},
309	{DAMAGE_ENERGY, "damage-energy"},
310	{COHERENT, "coherent-scatter"},
311	{INCOHERENT, "incoherent-scatter"},
312	{PAIR_PROD_ELEC, "pair-production-electron"},
313	{PAIR_PROD, "pair-production"},
314	{PAIR_PROD_NUC, "pair-production-nuclear"},
315	{PHOTOELECTRIC, "photoelectric"},
316	{N_PC, "(n,pc)"},
317	{N_DC, "(n,dc)"},
318	{N_TC, "(n,tc)"},
319	{N_3HEC, "(n,3Hec)"},
320	{N_AC, "(n,ac)"},
321	{N_2NC, "(n,2nc)"},
322	{HEATING_LOCAL, "heating-local"},
323	};
324
325	std::unordered_map<std::string, int> REACTION_TYPE_MAP;
326
327	void initialize_maps()	3,581✔
328	{
329	// Add level reactions to name map
330	for (int level = 0; level <= 48; ++level) {	179,050✔
331	if (level >= 1 && level <= 40) {	175,469✔
332	REACTION_NAME_MAP[50 + level] = fmt::format("(n,n{})", level);	143,240✔
333	}
334	REACTION_NAME_MAP[600 + level] = fmt::format("(n,p{})", level);	175,469✔
335	REACTION_NAME_MAP[650 + level] = fmt::format("(n,d{})", level);	175,469✔
336	REACTION_NAME_MAP[700 + level] = fmt::format("(n,t{})", level);	175,469✔
337	REACTION_NAME_MAP[750 + level] = fmt::format("(n,3He{})", level);	175,469✔
338	REACTION_NAME_MAP[800 + level] = fmt::format("(n,a{})", level);	175,469✔
339	if (level <= 15) {	175,469✔
340	REACTION_NAME_MAP[875 + level] = fmt::format("(n,2n{})", level);	57,296✔
341	}
342	}
343
344	// Create photoelectric subshells
345	for (int mt = 534; mt <= 572; ++mt) {	143,240✔
346	REACTION_NAME_MAP[mt] =	165,360✔
347	fmt::format("photoelectric, {} subshell", SUBSHELLS[mt - 534]);	393,276✔
348	}
349
350	// Invert name map to create type map
351	for (const auto& kv : REACTION_NAME_MAP) {	1,700,975✔
352	REACTION_TYPE_MAP[kv.second] = kv.first;	1,697,394✔
353	}
354	}	3,581✔
355
356	std::string reaction_name(int mt)	49,948✔
357	{
358	// Initialize remainder of name map and all of type map
359	if (REACTION_TYPE_MAP.empty())	49,948✔
360	initialize_maps();	33✔
361
362	// Get reaction name from map
363	auto it = REACTION_NAME_MAP.find(mt);	49,948✔
364	if (it != REACTION_NAME_MAP.end()) {	49,948✔
365	return it->second;	49,948✔
366	} else {
367	return fmt::format("MT={}", mt);	×
368	}
369	}
370
371	int reaction_type(std::string name)	38,543✔
372	{
373	// Initialize remainder of name map and all of type map
374	if (REACTION_TYPE_MAP.empty())	38,543✔
375	initialize_maps();	3,548✔
376
377	// (n,total) exists in REACTION_TYPE_MAP for MT=1, but we need this to return
378	// the special SCORE_TOTAL score
379	if (name == "(n,total)")	38,543✔
380	return SCORE_TOTAL;	×
381
382	// Check if type map has an entry for this reaction name
383	auto it = REACTION_TYPE_MAP.find(name);	38,543✔
384	if (it != REACTION_TYPE_MAP.end()) {	38,543✔
385	return it->second;	38,307✔
386	}
387
388	// Alternate names for several reactions
389	if (name == "elastic") {	236✔
390	return ELASTIC;	134✔
391	} else if (name == "n2n") {	102✔
392	return N_2N;	×
393	} else if (name == "n3n") {	102✔
394	return N_3N;	×
395	} else if (name == "n4n") {	102✔
396	return N_4N;	×
397	} else if (name == "H1-production") {	102✔
398	return N_XP;	16✔
399	} else if (name == "H2-production") {	86✔
400	return N_XD;	16✔
401	} else if (name == "H3-production") {	70✔
402	return N_XT;	16✔
403	} else if (name == "He3-production") {	54✔
404	return N_X3HE;	38✔
405	} else if (name == "He4-production") {	16✔
406	return N_XA;	16✔
407	}
408
409	// Assume the given string is a reaction MT number. Make sure it's a natural
410	// number then return.
411	int MT = 0;	×
412	try {
413	MT = std::stoi(name);	×
414	} catch (const std::invalid_argument& ex) {	×
415	throw std::invalid_argument(	×
416	"Invalid tally score \"" + name +	×
417	"\". See the docs "
418	"for details: "
419	"https://docs.openmc.org/en/stable/usersguide/tallies.html#scores");	×
420	}	×
421	if (MT < 1)	×
422	throw std::invalid_argument(	×
423	"Invalid tally score \"" + name +	×
424	"\". See the docs "
425	"for details: "
426	"https://docs.openmc.org/en/stable/usersguide/tallies.html#scores");	×
427	return MT;	×
428	}
429
430	} // namespace openmc

openmc-dev / openmc / 13703779155

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous