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