1 #include "openmc/material.h"
2 
3 #include <algorithm> // for min, max, sort, fill
4 #include <cmath>
5 #include <iterator>
6 #include <string>
7 #include <sstream>
8 #include <unordered_set>
9 
10 #include "xtensor/xbuilder.hpp"
11 #include "xtensor/xoperation.hpp"
12 #include "xtensor/xview.hpp"
13 
14 #include "openmc/capi.h"
15 #include "openmc/cross_sections.h"
16 #include "openmc/container_util.h"
17 #include "openmc/dagmc.h"
18 #include "openmc/error.h"
19 #include "openmc/file_utils.h"
20 #include "openmc/hdf5_interface.h"
21 #include "openmc/math_functions.h"
22 #include "openmc/message_passing.h"
23 #include "openmc/mgxs_interface.h"
24 #include "openmc/nuclide.h"
25 #include "openmc/photon.h"
26 #include "openmc/search.h"
27 #include "openmc/settings.h"
28 #include "openmc/simulation.h"
29 #include "openmc/string_utils.h"
30 #include "openmc/thermal.h"
31 #include "openmc/xml_interface.h"
32 
33 namespace openmc {
34 
35 //==============================================================================
36 // Global variables
37 //==============================================================================
38 
39 namespace model {
40 
41 std::unordered_map<int32_t, int32_t> material_map;
42 vector<unique_ptr<Material>> materials;
43 
44 } // namespace model
45 
46 //==============================================================================
47 // Material implementation
48 //==============================================================================
49 
Material(pugi::xml_node node)50 Material::Material(pugi::xml_node node)
51 {
52   index_ = model::materials.size(); // Avoids warning about narrowing
53 
54   if (check_for_node(node, "id")) {
55     this->set_id(std::stoi(get_node_value(node, "id")));
56   } else {
57     fatal_error("Must specify id of material in materials XML file.");
58   }
59 
60   if (check_for_node(node, "name")) {
61     name_ = get_node_value(node, "name");
62   }
63 
64   if (check_for_node(node, "depletable")) {
65     depletable_ = get_node_value_bool(node, "depletable");
66   }
67 
68   bool sum_density {false};
69   pugi::xml_node density_node = node.child("density");
70   std::string units;
71   if (density_node) {
72     units = get_node_value(density_node, "units");
73     if (units == "sum") {
74       sum_density = true;
75     } else if (units == "macro") {
76       if (check_for_node(density_node, "value")) {
77         density_ = std::stod(get_node_value(density_node, "value"));
78       } else {
79         density_ = 1.0;
80       }
81     } else {
82       double val = std::stod(get_node_value(density_node, "value"));
83       if (val <= 0.0) {
84         fatal_error("Need to specify a positive density on material "
85           + std::to_string(id_) + ".");
86       }
87 
88       if (units == "g/cc" || units == "g/cm3") {
89         density_ = -val;
90       } else if (units == "kg/m3") {
91         density_ = -1.0e-3 * val;
92       } else if (units == "atom/b-cm") {
93         density_ = val;
94       } else if (units == "atom/cc" || units == "atom/cm3") {
95         density_ = 1.0e-24 * val;
96       } else {
97         fatal_error("Unknown units '" + units + "' specified on material "
98           + std::to_string(id_) + ".");
99       }
100     }
101   } else {
102     fatal_error("Must specify <density> element in material "
103       + std::to_string(id_) + ".");
104   }
105 
106   if (node.child("element")) {
107     fatal_error("Unable to add an element to material " + std::to_string(id_) +
108       " since the element option has been removed from the xml input. "
109       "Elements can only be added via the Python API, which will expand "
110       "elements into their natural nuclides.");
111   }
112 
113   // =======================================================================
114   // READ AND PARSE <nuclide> TAGS
115 
116   // Check to ensure material has at least one nuclide
117   if (!check_for_node(node, "nuclide") && !check_for_node(node, "macroscopic")) {
118     fatal_error("No macroscopic data or nuclides specified on material "
119       + std::to_string(id_));
120   }
121 
122   // Create list of macroscopic x/s based on those specified, just treat
123   // them as nuclides. This is all really a facade so the user thinks they
124   // are entering in macroscopic data but the code treats them the same
125   // as nuclides internally.
126   // Get pointer list of XML <macroscopic>
127   auto node_macros = node.children("macroscopic");
128   int num_macros = std::distance(node_macros.begin(), node_macros.end());
129 
130   vector<std::string> names;
131   vector<double> densities;
132   if (settings::run_CE && num_macros > 0) {
133     fatal_error("Macroscopic can not be used in continuous-energy mode.");
134   } else if (num_macros > 1) {
135     fatal_error("Only one macroscopic object permitted per material, "
136       + std::to_string(id_));
137   } else if (num_macros == 1) {
138     pugi::xml_node node_nuc = *node_macros.begin();
139 
140     // Check for empty name on nuclide
141     if (!check_for_node(node_nuc, "name")) {
142       fatal_error("No name specified on macroscopic data in material "
143         + std::to_string(id_));
144     }
145 
146     // store nuclide name
147     std::string name = get_node_value(node_nuc, "name", false, true);
148     names.push_back(name);
149 
150     // Set density for macroscopic data
151     if (units == "macro") {
152       densities.push_back(density_);
153     } else {
154       fatal_error("Units can only be macro for macroscopic data " + name);
155     }
156   } else {
157     // Create list of nuclides based on those specified
158     for (auto node_nuc : node.children("nuclide")) {
159       // Check for empty name on nuclide
160       if (!check_for_node(node_nuc, "name")) {
161         fatal_error("No name specified on nuclide in material "
162           + std::to_string(id_));
163       }
164 
165       // store nuclide name
166       std::string name = get_node_value(node_nuc, "name", false, true);
167       names.push_back(name);
168 
169       // Check if no atom/weight percents were specified or if both atom and
170       // weight percents were specified
171       if (units == "macro") {
172         densities.push_back(density_);
173       } else {
174         bool has_ao = check_for_node(node_nuc, "ao");
175         bool has_wo = check_for_node(node_nuc, "wo");
176 
177         if (!has_ao && !has_wo) {
178           fatal_error("No atom or weight percent specified for nuclide: " + name);
179         } else if (has_ao && has_wo) {
180           fatal_error("Cannot specify both atom and weight percents for a "
181             "nuclide: " + name);
182         }
183 
184         // Copy atom/weight percents
185         if (has_ao) {
186           densities.push_back(std::stod(get_node_value(node_nuc, "ao")));
187         } else {
188           densities.push_back(-std::stod(get_node_value(node_nuc, "wo")));
189         }
190       }
191     }
192   }
193 
194   // =======================================================================
195   // READ AND PARSE <isotropic> element
196 
197   vector<std::string> iso_lab;
198   if (check_for_node(node, "isotropic")) {
199     iso_lab = get_node_array<std::string>(node, "isotropic");
200   }
201 
202   // ========================================================================
203   // COPY NUCLIDES TO ARRAYS IN MATERIAL
204 
205   // allocate arrays in Material object
206   auto n = names.size();
207   nuclide_.reserve(n);
208   atom_density_ = xt::empty<double>({n});
209   if (settings::photon_transport) element_.reserve(n);
210 
211   for (int i = 0; i < n; ++i) {
212     const auto& name {names[i]};
213 
214     // Check that this nuclide is listed in the nuclear data library
215     // (cross_sections.xml for CE and the MGXS HDF5 for MG)
216     LibraryKey key {Library::Type::neutron, name};
217     if (data::library_map.find(key) == data::library_map.end()) {
218       fatal_error("Could not find nuclide " + name + " in the "
219                   "nuclear data library.");
220     }
221 
222     // If this nuclide hasn't been encountered yet, we need to add its name
223     // and alias to the nuclide_dict
224     if (data::nuclide_map.find(name) == data::nuclide_map.end()) {
225       int index = data::nuclide_map.size();
226       data::nuclide_map[name] = index;
227       nuclide_.push_back(index);
228     } else {
229       nuclide_.push_back(data::nuclide_map[name]);
230     }
231 
232     // If the corresponding element hasn't been encountered yet and photon
233     // transport will be used, we need to add its symbol to the element_dict
234     if (settings::photon_transport) {
235       std::string element = to_element(name);
236 
237       // Make sure photon cross section data is available
238       LibraryKey key {Library::Type::photon, element};
239       if (data::library_map.find(key) == data::library_map.end()) {
240         fatal_error("Could not find element " + element
241           + " in cross_sections.xml.");
242       }
243 
244       if (data::element_map.find(element) == data::element_map.end()) {
245         int index = data::element_map.size();
246         data::element_map[element] = index;
247         element_.push_back(index);
248       } else {
249         element_.push_back(data::element_map[element]);
250       }
251     }
252 
253     // Copy atom/weight percent
254     atom_density_(i) = densities[i];
255   }
256 
257   if (settings::run_CE) {
258     // By default, isotropic-in-lab is not used
259     if (iso_lab.size() > 0) {
260       p0_.resize(n);
261 
262       // Apply isotropic-in-lab treatment to specified nuclides
263       for (int j = 0; j < n; ++j) {
264         for (const auto& nuc : iso_lab) {
265           if (names[j] == nuc) {
266             p0_[j] = true;
267             break;
268           }
269         }
270       }
271     }
272   }
273 
274   // Check to make sure either all atom percents or all weight percents are
275   // given
276   if (!(xt::all(atom_density_ >= 0.0) || xt::all(atom_density_ <= 0.0))) {
277     fatal_error("Cannot mix atom and weight percents in material "
278       + std::to_string(id_));
279   }
280 
281   // Determine density if it is a sum value
282   if (sum_density) density_ = xt::sum(atom_density_)();
283 
284 
285   if (check_for_node(node, "temperature")) {
286     temperature_ = std::stod(get_node_value(node, "temperature"));
287   }
288 
289   if (check_for_node(node, "volume")) {
290     volume_ = std::stod(get_node_value(node, "volume"));
291   }
292 
293   // =======================================================================
294   // READ AND PARSE <sab> TAG FOR THERMAL SCATTERING DATA
295   if (settings::run_CE) {
296     // Loop over <sab> elements
297 
298     vector<std::string> sab_names;
299     for (auto node_sab : node.children("sab")) {
300       // Determine name of thermal scattering table
301       if (!check_for_node(node_sab, "name")) {
302         fatal_error("Need to specify <name> for thermal scattering table.");
303       }
304       std::string name = get_node_value(node_sab, "name");
305       sab_names.push_back(name);
306 
307       // Read the fraction of nuclei affected by this thermal scattering table
308       double fraction = 1.0;
309       if (check_for_node(node_sab, "fraction")) {
310         fraction = std::stod(get_node_value(node_sab, "fraction"));
311       }
312 
313       // Check that the thermal scattering table is listed in the
314       // cross_sections.xml file
315       LibraryKey key {Library::Type::thermal, name};
316       if (data::library_map.find(key) == data::library_map.end()) {
317         fatal_error("Could not find thermal scattering data " + name +
318           " in cross_sections.xml file.");
319       }
320 
321       // Determine index of thermal scattering data in global
322       // data::thermal_scatt array
323       int index_table;
324       if (data::thermal_scatt_map.find(name) == data::thermal_scatt_map.end()) {
325         index_table = data::thermal_scatt_map.size();
326         data::thermal_scatt_map[name] = index_table;
327       } else {
328         index_table = data::thermal_scatt_map[name];
329       }
330 
331       // Add entry to thermal tables vector. For now, we put the nuclide index
332       // as zero since we don't know which nuclides the table is being applied
333       // to yet (this is assigned in init_thermal)
334       thermal_tables_.push_back({index_table, 0, fraction});
335     }
336   }
337 }
338 
~Material()339 Material::~Material()
340 {
341   model::material_map.erase(id_);
342 }
343 
finalize()344 void Material::finalize()
345 {
346   // Set fissionable if any nuclide is fissionable
347   for (const auto& i_nuc : nuclide_) {
348     if (data::nuclides[i_nuc]->fissionable_) {
349       fissionable_ = true;
350       break;
351     }
352   }
353 
354   // Generate material bremsstrahlung data for electrons and positrons
355   if (settings::photon_transport && settings::electron_treatment == ElectronTreatment::TTB) {
356     this->init_bremsstrahlung();
357   }
358 
359   // Assign thermal scattering tables
360   this->init_thermal();
361 
362   // Normalize density
363   this->normalize_density();
364 }
365 
normalize_density()366 void Material::normalize_density()
367 {
368   bool percent_in_atom = (atom_density_(0) >= 0.0);
369   bool density_in_atom = (density_ >= 0.0);
370 
371   for (int i = 0; i < nuclide_.size(); ++i) {
372     // determine atomic weight ratio
373     int i_nuc = nuclide_[i];
374     double awr = settings::run_CE ?
375       data::nuclides[i_nuc]->awr_ : data::mg.nuclides_[i_nuc].awr;
376 
377     // if given weight percent, convert all values so that they are divided
378     // by awr. thus, when a sum is done over the values, it's actually
379     // sum(w/awr)
380     if (!percent_in_atom) atom_density_(i) = -atom_density_(i) / awr;
381   }
382 
383   // determine normalized atom percents. if given atom percents, this is
384   // straightforward. if given weight percents, the value is w/awr and is
385   // divided by sum(w/awr)
386   atom_density_ /= xt::sum(atom_density_)();
387 
388   // Change density in g/cm^3 to atom/b-cm. Since all values are now in
389   // atom percent, the sum needs to be re-evaluated as 1/sum(x*awr)
390   if (!density_in_atom) {
391     double sum_percent = 0.0;
392     for (int i = 0; i < nuclide_.size(); ++i) {
393       int i_nuc = nuclide_[i];
394       double awr = settings::run_CE ?
395         data::nuclides[i_nuc]->awr_ : data::mg.nuclides_[i_nuc].awr;
396       sum_percent += atom_density_(i)*awr;
397     }
398     sum_percent = 1.0 / sum_percent;
399     density_ = -density_ * N_AVOGADRO / MASS_NEUTRON * sum_percent;
400   }
401 
402   // Calculate nuclide atom densities
403   atom_density_ *= density_;
404 
405   // Calculate density in g/cm^3.
406   density_gpcc_ = 0.0;
407   for (int i = 0; i < nuclide_.size(); ++i) {
408     int i_nuc = nuclide_[i];
409     double awr = settings::run_CE ? data::nuclides[i_nuc]->awr_ : 1.0;
410     density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO;
411   }
412 }
413 
init_thermal()414 void Material::init_thermal()
415 {
416   vector<ThermalTable> tables;
417 
418   std::unordered_set<int> already_checked;
419   for (const auto& table : thermal_tables_) {
420     // Make sure each S(a,b) table only gets checked once
421     if (already_checked.find(table.index_table) != already_checked.end()) {
422       continue;
423     }
424     already_checked.insert(table.index_table);
425 
426     // In order to know which nuclide the S(a,b) table applies to, we need
427     // to search through the list of nuclides for one which has a matching
428     // name
429     bool found = false;
430     for (int j = 0; j < nuclide_.size(); ++j) {
431       const auto& name {data::nuclides[nuclide_[j]]->name_};
432       if (contains(data::thermal_scatt[table.index_table]->nuclides_, name)) {
433         tables.push_back({table.index_table, j, table.fraction});
434         found = true;
435       }
436     }
437 
438     // Check to make sure thermal scattering table matched a nuclide
439     if (!found) {
440       fatal_error("Thermal scattering table " + data::thermal_scatt[
441         table.index_table]->name_  + " did not match any nuclide on material "
442         + std::to_string(id_));
443     }
444   }
445 
446   // Make sure each nuclide only appears in one table.
447   for (int j = 0; j < tables.size(); ++j) {
448     for (int k = j+1; k < tables.size(); ++k) {
449       if (tables[j].index_nuclide == tables[k].index_nuclide) {
450         int index = nuclide_[tables[j].index_nuclide];
451         auto name = data::nuclides[index]->name_;
452         fatal_error(name + " in material " + std::to_string(id_) + " was found "
453           "in multiple thermal scattering tables. Each nuclide can appear in "
454           "only one table per material.");
455       }
456     }
457   }
458 
459   // If there are multiple S(a,b) tables, we need to make sure that the
460   // entries in i_sab_nuclides are sorted or else they won't be applied
461   // correctly in the cross_section module.
462   std::sort(tables.begin(), tables.end(), [](ThermalTable a, ThermalTable b) {
463     return a.index_nuclide < b.index_nuclide;
464   });
465 
466   // Update the list of thermal tables
467   thermal_tables_ = tables;
468 }
469 
collision_stopping_power(double * s_col,bool positron)470 void Material::collision_stopping_power(double* s_col, bool positron)
471 {
472   // Average electron number and average atomic weight
473   double electron_density = 0.0;
474   double mass_density = 0.0;
475 
476   // Log of the mean excitation energy of the material
477   double log_I = 0.0;
478 
479   // Effective number of conduction electrons in the material
480   double n_conduction = 0.0;
481 
482   // Oscillator strength and square of the binding energy for each oscillator
483   // in material
484   vector<double> f;
485   vector<double> e_b_sq;
486 
487   for (int i = 0; i < element_.size(); ++i) {
488     const auto& elm = *data::elements[element_[i]];
489     double awr = data::nuclides[nuclide_[i]]->awr_;
490 
491     // Get atomic density of nuclide given atom/weight percent
492     double atom_density = (atom_density_[0] > 0.0) ?
493       atom_density_[i] : -atom_density_[i] / awr;
494 
495     electron_density += atom_density * elm.Z_;
496     mass_density += atom_density * awr * MASS_NEUTRON;
497     log_I += atom_density * elm.Z_ * std::log(elm.I_);
498 
499     for (int j = 0; j < elm.n_electrons_.size(); ++j) {
500       if (elm.n_electrons_[j] < 0) {
501         n_conduction -= elm.n_electrons_[j] * atom_density;
502         continue;
503       }
504       e_b_sq.push_back(elm.ionization_energy_[j] * elm.ionization_energy_[j]);
505       f.push_back(elm.n_electrons_[j] * atom_density);
506     }
507   }
508   log_I /= electron_density;
509   n_conduction /= electron_density;
510   for (auto& f_i : f) f_i /= electron_density;
511 
512   // Get density in g/cm^3 if it is given in atom/b-cm
513   double density = (density_ < 0.0) ? -density_ : mass_density / N_AVOGADRO;
514 
515   // Calculate the square of the plasma energy
516   double e_p_sq = PLANCK_C * PLANCK_C * PLANCK_C * N_AVOGADRO *
517     electron_density * density / (2.0 * PI * PI * FINE_STRUCTURE *
518     MASS_ELECTRON_EV * mass_density);
519 
520   // Get the Sternheimer adjustment factor
521   double rho = sternheimer_adjustment(f, e_b_sq, e_p_sq, n_conduction, log_I,
522     1.0e-6, 100);
523 
524   // Classical electron radius in cm
525   constexpr double CM_PER_ANGSTROM {1.0e-8};
526   constexpr double r_e = CM_PER_ANGSTROM * PLANCK_C / (2.0 * PI *
527     FINE_STRUCTURE * MASS_ELECTRON_EV);
528 
529   // Constant in expression for collision stopping power
530   constexpr double BARN_PER_CM_SQ {1.0e24};
531   double c = BARN_PER_CM_SQ * 2.0 * PI * r_e * r_e * MASS_ELECTRON_EV *
532     electron_density;
533 
534   // Loop over incident charged particle energies
535   for (int i = 0; i < data::ttb_e_grid.size(); ++i) {
536     double E = data::ttb_e_grid(i);
537 
538     // Get the density effect correction
539     double delta = density_effect(f, e_b_sq, e_p_sq, n_conduction, rho, E,
540       1.0e-6, 100);
541 
542     // Square of the ratio of the speed of light to the velocity of the charged
543     // particle
544     double beta_sq = E * (E + 2.0 * MASS_ELECTRON_EV) / ((E + MASS_ELECTRON_EV)
545       * (E + MASS_ELECTRON_EV));
546 
547     double tau = E / MASS_ELECTRON_EV;
548 
549     double F;
550     if (positron) {
551       double t = tau + 2.0;
552       F = std::log(4.0) - (beta_sq / 12.0) * (23.0 + 14.0 / t + 10.0 / (t * t)
553         + 4.0 / (t * t * t));
554     } else {
555       F = (1.0 - beta_sq) * (1.0 + tau * tau / 8.0 - (2.0 * tau + 1.0) *
556         std::log(2.0));
557     }
558 
559     // Calculate the collision stopping power for this energy
560     s_col[i] = c / beta_sq * (2.0 * (std::log(E) - log_I) + std::log(1.0 + tau
561       / 2.0) + F - delta);
562   }
563 }
564 
init_bremsstrahlung()565 void Material::init_bremsstrahlung()
566 {
567   // Create new object
568   ttb_ = make_unique<Bremsstrahlung>();
569 
570   // Get the size of the energy grids
571   auto n_k = data::ttb_k_grid.size();
572   auto n_e = data::ttb_e_grid.size();
573 
574   // Determine number of elements
575   int n = element_.size();
576 
577   for (int particle = 0; particle < 2; ++particle) {
578     // Loop over logic twice, once for electron, once for positron
579     BremsstrahlungData* ttb = (particle == 0) ? &ttb_->electron : &ttb_->positron;
580     bool positron = (particle == 1);
581 
582     // Allocate arrays for TTB data
583     ttb->pdf = xt::zeros<double>({n_e, n_e});
584     ttb->cdf = xt::zeros<double>({n_e, n_e});
585     ttb->yield = xt::empty<double>({n_e});
586 
587     // Allocate temporary arrays
588     xt::xtensor<double, 1> stopping_power_collision({n_e}, 0.0);
589     xt::xtensor<double, 1> stopping_power_radiative({n_e}, 0.0);
590     xt::xtensor<double, 2> dcs({n_e, n_k}, 0.0);
591 
592     double Z_eq_sq = 0.0;
593     double sum_density = 0.0;
594 
595     // Get the collision stopping power of the material
596     this->collision_stopping_power(stopping_power_collision.data(), positron);
597 
598     // Calculate the molecular DCS and the molecular radiative stopping power using
599     // Bragg's additivity rule.
600     for (int i = 0; i < n; ++i) {
601       // Get pointer to current element
602       const auto& elm = *data::elements[element_[i]];
603       double awr = data::nuclides[nuclide_[i]]->awr_;
604 
605       // Get atomic density and mass density of nuclide given atom/weight percent
606       double atom_density = (atom_density_[0] > 0.0) ?
607         atom_density_[i] : -atom_density_[i] / awr;
608 
609       // Calculate the "equivalent" atomic number Zeq of the material
610       Z_eq_sq += atom_density * elm.Z_ * elm.Z_;
611       sum_density += atom_density;
612 
613       // Accumulate material DCS
614       dcs += (atom_density * elm.Z_ * elm.Z_) * elm.dcs_;
615 
616       // Accumulate material radiative stopping power
617       stopping_power_radiative += atom_density * elm.stopping_power_radiative_;
618     }
619     Z_eq_sq /= sum_density;
620 
621     // Calculate the positron DCS and radiative stopping power. These are
622     // obtained by multiplying the electron DCS and radiative stopping powers by
623     // a factor r, which is a numerical approximation of the ratio of the
624     // radiative stopping powers for positrons and electrons. Source: F. Salvat,
625     // J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for
626     // Monte Carlo Simulation of Electron and Photon Transport," OECD-NEA,
627     // Issy-les-Moulineaux, France (2011).
628     if (positron) {
629       for (int i = 0; i < n_e; ++i) {
630         double t = std::log(1.0 + 1.0e6*data::ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON_EV));
631         double r = 1.0 - std::exp(-1.2359e-1*t + 6.1274e-2*std::pow(t, 2)
632           - 3.1516e-2*std::pow(t, 3) + 7.7446e-3*std::pow(t, 4)
633           - 1.0595e-3*std::pow(t, 5) + 7.0568e-5*std::pow(t, 6)
634           - 1.808e-6*std::pow(t, 7));
635         stopping_power_radiative(i) *= r;
636         auto dcs_i = xt::view(dcs, i, xt::all());
637         dcs_i *= r;
638       }
639     }
640 
641     // Total material stopping power
642     xt::xtensor<double, 1> stopping_power = stopping_power_collision +
643       stopping_power_radiative;
644 
645     // Loop over photon energies
646     xt::xtensor<double, 1> f({n_e}, 0.0);
647     xt::xtensor<double, 1> z({n_e}, 0.0);
648     for (int i = 0; i < n_e - 1; ++i) {
649       double w = data::ttb_e_grid(i);
650 
651       // Loop over incident particle energies
652       for (int j = i; j < n_e; ++j) {
653         double e = data::ttb_e_grid(j);
654 
655         // Reduced photon energy
656         double k = w / e;
657 
658         // Find the lower bounding index of the reduced photon energy
659         int i_k = lower_bound_index(data::ttb_k_grid.cbegin(),
660           data::ttb_k_grid.cend(), k);
661 
662         // Get the interpolation bounds
663         double k_l = data::ttb_k_grid(i_k);
664         double k_r = data::ttb_k_grid(i_k + 1);
665         double x_l = dcs(j, i_k);
666         double x_r = dcs(j, i_k + 1);
667 
668         // Find the value of the DCS using linear interpolation in reduced
669         // photon energy k
670         double x = x_l + (k - k_l)*(x_r - x_l)/(k_r - k_l);
671 
672         // Square of the ratio of the speed of light to the velocity of the
673         // charged particle
674         double beta_sq = e * (e + 2.0 * MASS_ELECTRON_EV) / ((e +
675           MASS_ELECTRON_EV) * (e + MASS_ELECTRON_EV));
676 
677         // Compute the integrand of the PDF
678         f(j) = x / (beta_sq * stopping_power(j) * w);
679       }
680 
681       // Number of points to integrate
682       int n = n_e - i;
683 
684       // Integrate the PDF using cubic spline integration over the incident
685       // particle energy
686       if (n > 2) {
687         spline(n, &data::ttb_e_grid(i), &f(i), &z(i));
688 
689         double c = 0.0;
690         for (int j = i; j < n_e - 1; ++j) {
691           c += spline_integrate(n, &data::ttb_e_grid(i), &f(i), &z(i),
692             data::ttb_e_grid(j), data::ttb_e_grid(j+1));
693 
694           ttb->pdf(j+1,i) = c;
695         }
696 
697       // Integrate the last two points using trapezoidal rule in log-log space
698       } else {
699         double e_l = std::log(data::ttb_e_grid(i));
700         double e_r = std::log(data::ttb_e_grid(i+1));
701         double x_l = std::log(f(i));
702         double x_r = std::log(f(i+1));
703 
704         ttb->pdf(i+1,i) = 0.5*(e_r - e_l)*(std::exp(e_l + x_l) + std::exp(e_r + x_r));
705       }
706     }
707 
708     // Loop over incident particle energies
709     for (int j = 1; j < n_e; ++j) {
710       // Set last element of PDF to small non-zero value to enable log-log
711       // interpolation
712       ttb->pdf(j,j) = std::exp(-500.0);
713 
714       // Loop over photon energies
715       double c = 0.0;
716       for (int i = 0; i < j; ++i) {
717         // Integrate the CDF from the PDF using the trapezoidal rule in log-log
718         // space
719         double w_l = std::log(data::ttb_e_grid(i));
720         double w_r = std::log(data::ttb_e_grid(i+1));
721         double x_l = std::log(ttb->pdf(j,i));
722         double x_r = std::log(ttb->pdf(j,i+1));
723 
724         c += 0.5*(w_r - w_l)*(std::exp(w_l + x_l) + std::exp(w_r + x_r));
725         ttb->cdf(j,i+1) = c;
726       }
727 
728       // Set photon number yield
729       ttb->yield(j) = c;
730     }
731 
732     // Use logarithm of number yield since it is log-log interpolated
733     ttb->yield = xt::where(ttb->yield > 0.0, xt::log(ttb->yield), -500.0);
734   }
735 }
736 
init_nuclide_index()737 void Material::init_nuclide_index()
738 {
739   int n = settings::run_CE ?
740     data::nuclides.size() : data::mg.nuclides_.size();
741   mat_nuclide_index_.resize(n);
742   std::fill(mat_nuclide_index_.begin(), mat_nuclide_index_.end(), C_NONE);
743   for (int i = 0; i < nuclide_.size(); ++i) {
744     mat_nuclide_index_[nuclide_[i]] = i;
745   }
746 }
747 
calculate_xs(Particle & p) const748 void Material::calculate_xs(Particle& p) const
749 {
750   // Set all material macroscopic cross sections to zero
751   p.macro_xs().total = 0.0;
752   p.macro_xs().absorption = 0.0;
753   p.macro_xs().fission = 0.0;
754   p.macro_xs().nu_fission = 0.0;
755 
756   if (p.type() == ParticleType::neutron) {
757     this->calculate_neutron_xs(p);
758   } else if (p.type() == ParticleType::photon) {
759     this->calculate_photon_xs(p);
760   }
761 }
762 
calculate_neutron_xs(Particle & p) const763 void Material::calculate_neutron_xs(Particle& p) const
764 {
765   // Find energy index on energy grid
766   int neutron = static_cast<int>(ParticleType::neutron);
767   int i_grid =
768     std::log(p.E() / data::energy_min[neutron]) / simulation::log_spacing;
769 
770   // Determine if this material has S(a,b) tables
771   bool check_sab = (thermal_tables_.size() > 0);
772 
773   // Initialize position in i_sab_nuclides
774   int j = 0;
775 
776   // Add contribution from each nuclide in material
777   for (int i = 0; i < nuclide_.size(); ++i) {
778     // ======================================================================
779     // CHECK FOR S(A,B) TABLE
780 
781     int i_sab = C_NONE;
782     double sab_frac = 0.0;
783 
784     // Check if this nuclide matches one of the S(a,b) tables specified.
785     // This relies on thermal_tables_ being sorted by .index_nuclide
786     if (check_sab) {
787       const auto& sab {thermal_tables_[j]};
788       if (i == sab.index_nuclide) {
789         // Get index in sab_tables
790         i_sab = sab.index_table;
791         sab_frac = sab.fraction;
792 
793         // If particle energy is greater than the highest energy for the
794         // S(a,b) table, then don't use the S(a,b) table
795         if (p.E() > data::thermal_scatt[i_sab]->energy_max_)
796           i_sab = C_NONE;
797 
798         // Increment position in thermal_tables_
799         ++j;
800 
801         // Don't check for S(a,b) tables if there are no more left
802         if (j == thermal_tables_.size()) check_sab = false;
803       }
804     }
805 
806     // ======================================================================
807     // CALCULATE MICROSCOPIC CROSS SECTION
808 
809     // Determine microscopic cross sections for this nuclide
810     int i_nuclide = nuclide_[i];
811 
812     // Calculate microscopic cross section for this nuclide
813     const auto& micro {p.neutron_xs(i_nuclide)};
814     if (p.E() != micro.last_E || p.sqrtkT() != micro.last_sqrtkT ||
815         i_sab != micro.index_sab || sab_frac != micro.sab_frac) {
816       data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, p);
817     }
818 
819     // ======================================================================
820     // ADD TO MACROSCOPIC CROSS SECTION
821 
822     // Copy atom density of nuclide in material
823     double atom_density = atom_density_(i);
824 
825     // Add contributions to cross sections
826     p.macro_xs().total += atom_density * micro.total;
827     p.macro_xs().absorption += atom_density * micro.absorption;
828     p.macro_xs().fission += atom_density * micro.fission;
829     p.macro_xs().nu_fission += atom_density * micro.nu_fission;
830   }
831 }
832 
calculate_photon_xs(Particle & p) const833 void Material::calculate_photon_xs(Particle& p) const
834 {
835   p.macro_xs().coherent = 0.0;
836   p.macro_xs().incoherent = 0.0;
837   p.macro_xs().photoelectric = 0.0;
838   p.macro_xs().pair_production = 0.0;
839 
840   // Add contribution from each nuclide in material
841   for (int i = 0; i < nuclide_.size(); ++i) {
842     // ========================================================================
843     // CALCULATE MICROSCOPIC CROSS SECTION
844 
845     // Determine microscopic cross sections for this nuclide
846     int i_element = element_[i];
847 
848     // Calculate microscopic cross section for this nuclide
849     const auto& micro {p.photon_xs(i_element)};
850     if (p.E() != micro.last_E) {
851       data::elements[i_element]->calculate_xs(p);
852     }
853 
854     // ========================================================================
855     // ADD TO MACROSCOPIC CROSS SECTION
856 
857     // Copy atom density of nuclide in material
858     double atom_density = atom_density_(i);
859 
860     // Add contributions to material macroscopic cross sections
861     p.macro_xs().total += atom_density * micro.total;
862     p.macro_xs().coherent += atom_density * micro.coherent;
863     p.macro_xs().incoherent += atom_density * micro.incoherent;
864     p.macro_xs().photoelectric += atom_density * micro.photoelectric;
865     p.macro_xs().pair_production += atom_density * micro.pair_production;
866   }
867 }
868 
set_id(int32_t id)869 void Material::set_id(int32_t id)
870 {
871   Expects(id >= 0 || id == C_NONE);
872 
873   // Clear entry in material map if an ID was already assigned before
874   if (id_ != C_NONE) {
875     model::material_map.erase(id_);
876     id_ = C_NONE;
877   }
878 
879   // Make sure no other material has same ID
880   if (model::material_map.find(id) != model::material_map.end()) {
881     throw std::runtime_error{"Two materials have the same ID: " + std::to_string(id)};
882   }
883 
884   // If no ID specified, auto-assign next ID in sequence
885   if (id == C_NONE) {
886     id = 0;
887     for (const auto& m : model::materials) {
888       id = std::max(id, m->id_);
889     }
890     ++id;
891   }
892 
893   // Update ID and entry in material map
894   id_ = id;
895   model::material_map[id] = index_;
896 }
897 
set_density(double density,gsl::cstring_span units)898 void Material::set_density(double density, gsl::cstring_span units)
899 {
900   Expects(density >= 0.0);
901 
902   if (nuclide_.empty()) {
903     throw std::runtime_error{"No nuclides exist in material yet."};
904   }
905 
906   if (units == "atom/b-cm") {
907     // Set total density based on value provided
908     density_ = density;
909 
910     // Determine normalized atom percents
911     double sum_percent = xt::sum(atom_density_)();
912     atom_density_ /= sum_percent;
913 
914     // Recalculate nuclide atom densities based on given density
915     atom_density_ *= density;
916 
917     // Calculate density in g/cm^3.
918     density_gpcc_ = 0.0;
919     for (int i = 0; i < nuclide_.size(); ++i) {
920       int i_nuc = nuclide_[i];
921       double awr = data::nuclides[i_nuc]->awr_;
922       density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO;
923     }
924   } else if (units == "g/cm3" || units == "g/cc") {
925     // Determine factor by which to change densities
926     double previous_density_gpcc = density_gpcc_;
927     double f = density / previous_density_gpcc;
928 
929     // Update densities
930     density_gpcc_ = density;
931     density_ *= f;
932     atom_density_ *= f;
933   } else {
934     throw std::invalid_argument{"Invalid units '" + std::string(units.data())
935       + "' specified."};
936   }
937 }
938 
set_densities(const vector<std::string> & name,const vector<double> & density)939 void Material::set_densities(
940   const vector<std::string>& name, const vector<double>& density)
941 {
942   auto n = name.size();
943   Expects(n > 0);
944   Expects(n == density.size());
945 
946   if (n != nuclide_.size()) {
947     nuclide_.resize(n);
948     atom_density_ = xt::zeros<double>({n});
949     if (settings::photon_transport) element_.resize(n);
950   }
951 
952   double sum_density = 0.0;
953   for (gsl::index i = 0; i < n; ++i) {
954     const auto& nuc {name[i]};
955     if (data::nuclide_map.find(nuc) == data::nuclide_map.end()) {
956       int err = openmc_load_nuclide(nuc.c_str(), nullptr, 0);
957       if (err < 0) throw std::runtime_error{openmc_err_msg};
958     }
959 
960     nuclide_[i] = data::nuclide_map.at(nuc);
961     Expects(density[i] > 0.0);
962     atom_density_(i) = density[i];
963     sum_density += density[i];
964 
965     if (settings::photon_transport) {
966       auto element_name = to_element(nuc);
967       element_[i] = data::element_map.at(element_name);
968     }
969   }
970 
971   // Set total density to the sum of the vector
972   this->set_density(sum_density, "atom/b-cm");
973 
974   // Generate material bremsstrahlung data for electrons and positrons
975   if (settings::photon_transport && settings::electron_treatment == ElectronTreatment::TTB) {
976     this->init_bremsstrahlung();
977   }
978 
979   // Assign S(a,b) tables
980   this->init_thermal();
981 }
982 
volume() const983 double Material::volume() const
984 {
985   if (volume_ < 0.0) {
986     throw std::runtime_error{"Volume for material with ID="
987       + std::to_string(id_) + " not set."};
988   }
989   return volume_;
990 }
991 
temperature() const992 double Material::temperature() const
993 {
994   // If material doesn't have an assigned temperature, use global default
995   return temperature_ >= 0 ? temperature_ : settings::temperature_default;
996 }
997 
to_hdf5(hid_t group) const998 void Material::to_hdf5(hid_t group) const
999 {
1000   hid_t material_group = create_group(group, "material " + std::to_string(id_));
1001 
1002   write_attribute(material_group, "depletable", static_cast<int>(depletable_));
1003   if (volume_ > 0.0) {
1004     write_attribute(material_group, "volume", volume_);
1005   }
1006   if (temperature_ > 0.0) {
1007     write_attribute(material_group, "temperature", temperature_);
1008   }
1009   write_dataset(material_group, "name", name_);
1010   write_dataset(material_group, "atom_density", density_);
1011 
1012   // Copy nuclide/macro name for each nuclide to vector
1013   vector<std::string> nuc_names;
1014   vector<std::string> macro_names;
1015   vector<double> nuc_densities;
1016   if (settings::run_CE) {
1017     for (int i = 0; i < nuclide_.size(); ++i) {
1018       int i_nuc = nuclide_[i];
1019       nuc_names.push_back(data::nuclides[i_nuc]->name_);
1020       nuc_densities.push_back(atom_density_(i));
1021     }
1022   } else {
1023     for (int i = 0; i < nuclide_.size(); ++i) {
1024       int i_nuc = nuclide_[i];
1025       if (data::mg.nuclides_[i_nuc].awr != MACROSCOPIC_AWR) {
1026         nuc_names.push_back(data::mg.nuclides_[i_nuc].name);
1027         nuc_densities.push_back(atom_density_(i));
1028       } else {
1029         macro_names.push_back(data::mg.nuclides_[i_nuc].name);
1030       }
1031     }
1032   }
1033 
1034   // Write vector to 'nuclides'
1035   if (!nuc_names.empty()) {
1036     write_dataset(material_group, "nuclides", nuc_names);
1037     write_dataset(material_group, "nuclide_densities", nuc_densities);
1038   }
1039 
1040   // Write vector to 'macroscopics'
1041   if (!macro_names.empty()) {
1042     write_dataset(material_group, "macroscopics", macro_names);
1043   }
1044 
1045   if (!thermal_tables_.empty()) {
1046     vector<std::string> sab_names;
1047     for (const auto& table : thermal_tables_) {
1048       sab_names.push_back(data::thermal_scatt[table.index_table]->name_);
1049     }
1050     write_dataset(material_group, "sab_names", sab_names);
1051   }
1052 
1053   close_group(material_group);
1054 }
1055 
add_nuclide(const std::string & name,double density)1056 void Material::add_nuclide(const std::string& name, double density)
1057 {
1058   // Check if nuclide is already in material
1059   for (int i = 0; i < nuclide_.size(); ++i) {
1060     int i_nuc = nuclide_[i];
1061     if (data::nuclides[i_nuc]->name_ == name) {
1062       double awr = data::nuclides[i_nuc]->awr_;
1063       density_ += density - atom_density_(i);
1064       density_gpcc_ += (density - atom_density_(i))
1065         * awr * MASS_NEUTRON / N_AVOGADRO;
1066       atom_density_(i) = density;
1067       return;
1068     }
1069   }
1070 
1071   // If nuclide wasn't found, extend nuclide/density arrays
1072   int err = openmc_load_nuclide(name.c_str(), nullptr, 0);
1073   if (err < 0) throw std::runtime_error{openmc_err_msg};
1074 
1075   // Append new nuclide/density
1076   int i_nuc = data::nuclide_map[name];
1077   nuclide_.push_back(i_nuc);
1078 
1079   // Append new element if photon transport is on
1080   if (settings::photon_transport) {
1081     int i_elem = data::element_map[to_element(name)];
1082     element_.push_back(i_elem);
1083   }
1084 
1085   auto n = nuclide_.size();
1086 
1087   // Create copy of atom_density_ array with one extra entry
1088   xt::xtensor<double, 1> atom_density = xt::zeros<double>({n});
1089   xt::view(atom_density, xt::range(0, n-1)) = atom_density_;
1090   atom_density(n-1) = density;
1091   atom_density_ = atom_density;
1092 
1093   density_ += density;
1094   density_gpcc_ += density * data::nuclides[i_nuc]->awr_
1095     * MASS_NEUTRON / N_AVOGADRO;
1096 }
1097 
1098 //==============================================================================
1099 // Non-method functions
1100 //==============================================================================
1101 
sternheimer_adjustment(const vector<double> & f,const vector<double> & e_b_sq,double e_p_sq,double n_conduction,double log_I,double tol,int max_iter)1102 double sternheimer_adjustment(const vector<double>& f,
1103   const vector<double>& e_b_sq, double e_p_sq, double n_conduction,
1104   double log_I, double tol, int max_iter)
1105 {
1106   // Get the total number of oscillators
1107   int n = f.size();
1108 
1109   // Calculate the Sternheimer adjustment factor using Newton's method
1110   double rho = 2.0;
1111   int iter;
1112   for (iter = 0; iter < max_iter; ++iter) {
1113     double rho_0 = rho;
1114 
1115     // Function to find the root of and its derivative
1116     double g = 0.0;
1117     double gp = 0.0;
1118 
1119     for (int i = 0; i < n; ++i) {
1120       // Square of resonance energy of a bound-shell oscillator
1121       double e_r_sq = e_b_sq[i] * rho * rho + 2.0 / 3.0 * f[i] * e_p_sq;
1122       g += f[i] * std::log(e_r_sq);
1123       gp += e_b_sq[i] * f[i] * rho / e_r_sq;
1124     }
1125     // Include conduction electrons
1126     if (n_conduction > 0.0) {
1127       g += n_conduction * std::log(n_conduction * e_p_sq);
1128     }
1129 
1130     // Set the next guess: rho_n+1 = rho_n - g(rho_n)/g'(rho_n)
1131     rho -= (g - 2.0 * log_I) / (2.0 * gp);
1132 
1133     // If the initial guess is too large, rho can be negative
1134     if (rho < 0.0) rho = rho_0 / 2.0;
1135 
1136     // Check for convergence
1137     if (std::abs(rho - rho_0) / rho_0 < tol) break;
1138   }
1139   // Did not converge
1140   if (iter >= max_iter) {
1141     warning("Maximum Newton-Raphson iterations exceeded.");
1142     rho = 1.0e-6;
1143   }
1144   return rho;
1145 }
1146 
density_effect(const vector<double> & f,const vector<double> & e_b_sq,double e_p_sq,double n_conduction,double rho,double E,double tol,int max_iter)1147 double density_effect(const vector<double>& f, const vector<double>& e_b_sq,
1148   double e_p_sq, double n_conduction, double rho, double E, double tol,
1149   int max_iter)
1150 {
1151   // Get the total number of oscillators
1152   int n = f.size();
1153 
1154   // Square of the ratio of the speed of light to the velocity of the charged
1155   // particle
1156   double beta_sq = E * (E + 2.0 * MASS_ELECTRON_EV) / ((E + MASS_ELECTRON_EV) *
1157        (E + MASS_ELECTRON_EV));
1158 
1159   // For nonmetals, delta = 0 for beta < beta_0, where beta_0 is obtained by
1160   // setting the frequency w = 0.
1161   double beta_0_sq = 0.0;
1162   if (n_conduction == 0.0) {
1163     for (int i = 0; i < n; ++i) {
1164       beta_0_sq += f[i] * e_p_sq / (e_b_sq[i] * rho * rho);
1165     }
1166     beta_0_sq = 1.0 / (1.0 + beta_0_sq);
1167   }
1168   double delta = 0.0;
1169   if (beta_sq < beta_0_sq) return delta;
1170 
1171   // Compute the square of the frequency w^2 using Newton's method, with the
1172   // initial guess of w^2 equal to beta^2 * gamma^2
1173   double w_sq = E / MASS_ELECTRON_EV * (E / MASS_ELECTRON_EV + 2);
1174   int iter;
1175   for (iter = 0; iter < max_iter; ++iter) {
1176     double w_sq_0 = w_sq;
1177 
1178     // Function to find the root of and its derivative
1179     double g = 0.0;
1180     double gp = 0.0;
1181 
1182     for (int i = 0; i < n; ++i) {
1183       double c = e_b_sq[i] * rho * rho / e_p_sq + w_sq;
1184       g += f[i] / c;
1185       gp -= f[i] / (c * c);
1186     }
1187     // Include conduction electrons
1188     g += n_conduction / w_sq;
1189     gp -= n_conduction / (w_sq * w_sq);
1190 
1191     // Set the next guess: w_n+1 = w_n - g(w_n)/g'(w_n)
1192     w_sq -= (g + 1.0 - 1.0 / beta_sq) / gp;
1193 
1194     // If the initial guess is too large, w can be negative
1195     if (w_sq < 0.0) w_sq = w_sq_0 / 2.0;
1196 
1197     // Check for convergence
1198     if (std::abs(w_sq - w_sq_0) / w_sq_0 < tol) break;
1199   }
1200   // Did not converge
1201   if (iter >= max_iter) {
1202     warning("Maximum Newton-Raphson iterations exceeded: setting density "
1203         "effect correction to zero.");
1204     return delta;
1205   }
1206 
1207   // Solve for the density effect correction
1208   for (int i = 0; i < n; ++i) {
1209     double l_sq = e_b_sq[i] * rho * rho / e_p_sq + 2.0 / 3.0 * f[i];
1210     delta += f[i] * std::log((l_sq + w_sq)/l_sq);
1211   }
1212   // Include conduction electrons
1213   if (n_conduction > 0.0) {
1214     delta += n_conduction * std::log((n_conduction + w_sq) / n_conduction);
1215   }
1216 
1217   return delta - w_sq * (1.0 - beta_sq);
1218 }
1219 
read_materials_xml()1220 void read_materials_xml()
1221 {
1222   write_message("Reading materials XML file...", 5);
1223 
1224   pugi::xml_document doc;
1225 
1226   bool using_dagmc_mats = false;
1227 #ifdef DAGMC
1228   if (settings::dagmc) {
1229     using_dagmc_mats = read_uwuw_materials(doc);
1230   }
1231 #endif
1232 
1233 
1234   if (!using_dagmc_mats) {
1235     // Check if materials.xml exists
1236     std::string filename = settings::path_input + "materials.xml";
1237     if (!file_exists(filename)) {
1238       fatal_error("Material XML file '" + filename + "' does not exist!");
1239     }
1240 
1241     // Parse materials.xml file and get root element
1242     doc.load_file(filename.c_str());
1243   }
1244 
1245   // Loop over XML material elements and populate the array.
1246   pugi::xml_node root = doc.document_element();
1247   for (pugi::xml_node material_node : root.children("material")) {
1248     model::materials.push_back(make_unique<Material>(material_node));
1249   }
1250   model::materials.shrink_to_fit();
1251 }
1252 
free_memory_material()1253 void free_memory_material()
1254 {
1255   model::materials.clear();
1256   model::material_map.clear();
1257 }
1258 
1259 //==============================================================================
1260 // C API
1261 //==============================================================================
1262 
1263 extern "C" int
openmc_get_material_index(int32_t id,int32_t * index)1264 openmc_get_material_index(int32_t id, int32_t* index)
1265 {
1266   auto it = model::material_map.find(id);
1267   if (it == model::material_map.end()) {
1268     set_errmsg("No material exists with ID=" + std::to_string(id) + ".");
1269     return OPENMC_E_INVALID_ID;
1270   } else {
1271     *index = it->second;
1272     return 0;
1273   }
1274 }
1275 
1276 extern "C" int
openmc_material_add_nuclide(int32_t index,const char * name,double density)1277 openmc_material_add_nuclide(int32_t index, const char* name, double density)
1278 {
1279   int err = 0;
1280   if (index >= 0 && index < model::materials.size()) {
1281     try {
1282       model::materials[index]->add_nuclide(name, density);
1283     } catch (const std::runtime_error& e) {
1284       return OPENMC_E_DATA;
1285     }
1286   } else {
1287     set_errmsg("Index in materials array is out of bounds.");
1288     return OPENMC_E_OUT_OF_BOUNDS;
1289   }
1290   return err;
1291 }
1292 
1293 extern "C" int
openmc_material_get_densities(int32_t index,const int ** nuclides,const double ** densities,int * n)1294 openmc_material_get_densities(int32_t index, const int** nuclides, const double** densities, int* n)
1295 {
1296   if (index >= 0 && index < model::materials.size()) {
1297     auto& mat = model::materials[index];
1298     if (!mat->nuclides().empty()) {
1299       *nuclides = mat->nuclides().data();
1300       *densities = mat->densities().data();
1301       *n = mat->nuclides().size();
1302       return 0;
1303     } else {
1304       set_errmsg("Material atom density array has not been allocated.");
1305       return OPENMC_E_ALLOCATE;
1306     }
1307   } else {
1308     set_errmsg("Index in materials array is out of bounds.");
1309     return OPENMC_E_OUT_OF_BOUNDS;
1310   }
1311 }
1312 
1313 extern "C" int
openmc_material_get_density(int32_t index,double * density)1314 openmc_material_get_density(int32_t index, double* density)
1315 {
1316   if (index >= 0 && index < model::materials.size()) {
1317     auto& mat = model::materials[index];
1318     *density = mat->density_gpcc();
1319     return 0;
1320   } else {
1321     set_errmsg("Index in materials array is out of bounds.");
1322     return OPENMC_E_OUT_OF_BOUNDS;
1323   }
1324 }
1325 
1326 extern "C" int
openmc_material_get_fissionable(int32_t index,bool * fissionable)1327 openmc_material_get_fissionable(int32_t index, bool* fissionable)
1328 {
1329   if (index >= 0 && index < model::materials.size()) {
1330     *fissionable = model::materials[index]->fissionable();
1331     return 0;
1332   } else {
1333     set_errmsg("Index in materials array is out of bounds.");
1334     return OPENMC_E_OUT_OF_BOUNDS;
1335   }
1336 }
1337 
1338 extern "C" int
openmc_material_get_id(int32_t index,int32_t * id)1339 openmc_material_get_id(int32_t index, int32_t* id)
1340 {
1341   if (index >= 0 && index < model::materials.size()) {
1342     *id = model::materials[index]->id();
1343     return 0;
1344   } else {
1345     set_errmsg("Index in materials array is out of bounds.");
1346     return OPENMC_E_OUT_OF_BOUNDS;
1347   }
1348 }
1349 
1350 extern "C" int
openmc_material_get_temperature(int32_t index,double * temperature)1351 openmc_material_get_temperature(int32_t index, double* temperature)
1352 {
1353   if (index < 0 || index >= model::materials.size()) {
1354     set_errmsg("Index in materials array is out of bounds.");
1355     return OPENMC_E_OUT_OF_BOUNDS;
1356   }
1357   *temperature = model::materials[index]->temperature();
1358   return 0;
1359 }
1360 
1361 
1362 extern "C" int
openmc_material_get_volume(int32_t index,double * volume)1363 openmc_material_get_volume(int32_t index, double* volume)
1364 {
1365   if (index >= 0 && index < model::materials.size()) {
1366     try {
1367       *volume = model::materials[index]->volume();
1368     } catch (const std::exception& e) {
1369       set_errmsg(e.what());
1370       return OPENMC_E_UNASSIGNED;
1371     }
1372     return 0;
1373   } else {
1374     set_errmsg("Index in materials array is out of bounds.");
1375     return OPENMC_E_OUT_OF_BOUNDS;
1376   }
1377 }
1378 
1379 extern "C" int
openmc_material_set_density(int32_t index,double density,const char * units)1380 openmc_material_set_density(int32_t index, double density, const char* units)
1381 {
1382   if (index >= 0 && index < model::materials.size()) {
1383     try {
1384       model::materials[index]->set_density(density, units);
1385     } catch (const std::exception& e) {
1386       set_errmsg(e.what());
1387       return OPENMC_E_UNASSIGNED;
1388     }
1389   } else {
1390     set_errmsg("Index in materials array is out of bounds.");
1391     return OPENMC_E_OUT_OF_BOUNDS;
1392   }
1393   return 0;
1394 }
1395 
1396 extern "C" int
openmc_material_set_densities(int32_t index,int n,const char ** name,const double * density)1397 openmc_material_set_densities(int32_t index, int n, const char** name, const double* density)
1398 {
1399   if (index >= 0 && index < model::materials.size()) {
1400     try {
1401       model::materials[index]->set_densities({name, name + n}, {density, density + n});
1402     } catch (const std::exception& e) {
1403       set_errmsg(e.what());
1404       return OPENMC_E_UNASSIGNED;
1405     }
1406   } else {
1407     set_errmsg("Index in materials array is out of bounds.");
1408     return OPENMC_E_OUT_OF_BOUNDS;
1409   }
1410   return 0;
1411 }
1412 
1413 extern "C" int
openmc_material_set_id(int32_t index,int32_t id)1414 openmc_material_set_id(int32_t index, int32_t id)
1415 {
1416   if (index >= 0 && index < model::materials.size()) {
1417     try {
1418       model::materials.at(index)->set_id(id);
1419     } catch (const std::exception& e) {
1420       set_errmsg(e.what());
1421       return OPENMC_E_UNASSIGNED;
1422     }
1423   } else {
1424     set_errmsg("Index in materials array is out of bounds.");
1425     return OPENMC_E_OUT_OF_BOUNDS;
1426   }
1427   return 0;
1428 }
1429 
1430 extern "C" int
openmc_material_get_name(int32_t index,const char ** name)1431 openmc_material_get_name(int32_t index, const char** name) {
1432   if (index < 0 || index >= model::materials.size()) {
1433     set_errmsg("Index in materials array is out of bounds.");
1434     return OPENMC_E_OUT_OF_BOUNDS;
1435   }
1436 
1437   *name = model::materials[index]->name().data();
1438 
1439   return 0;
1440 }
1441 
1442 extern "C" int
openmc_material_set_name(int32_t index,const char * name)1443 openmc_material_set_name(int32_t index, const char* name) {
1444   if (index < 0 || index >= model::materials.size()) {
1445     set_errmsg("Index in materials array is out of bounds.");
1446     return OPENMC_E_OUT_OF_BOUNDS;
1447   }
1448 
1449   model::materials[index]->set_name(name);
1450 
1451   return 0;
1452 }
1453 
1454 extern "C" int
openmc_material_set_volume(int32_t index,double volume)1455 openmc_material_set_volume(int32_t index, double volume)
1456 {
1457   if (index >= 0 && index < model::materials.size()) {
1458     auto& m {model::materials[index]};
1459     if (volume >= 0.0) {
1460       m->volume_ = volume;
1461       return 0;
1462     } else {
1463       set_errmsg("Volume must be non-negative");
1464       return OPENMC_E_INVALID_ARGUMENT;
1465     }
1466   } else {
1467     set_errmsg("Index in materials array is out of bounds.");
1468     return OPENMC_E_OUT_OF_BOUNDS;
1469   }
1470 }
1471 
1472 extern "C" int
openmc_extend_materials(int32_t n,int32_t * index_start,int32_t * index_end)1473 openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end)
1474 {
1475   if (index_start) *index_start = model::materials.size();
1476   if (index_end) *index_end = model::materials.size() + n - 1;
1477   for (int32_t i = 0; i < n; i++) {
1478     model::materials.push_back(make_unique<Material>());
1479   }
1480   return 0;
1481 }
1482 
n_materials()1483 extern "C" size_t n_materials() { return model::materials.size(); }
1484 
1485 } // namespace openmc
1486