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