1 #include "openmc/settings.h"
2
3 #include <cmath> // for ceil, pow
4 #include <limits> // for numeric_limits
5 #include <string>
6
7 #include <fmt/core.h>
8 #ifdef _OPENMP
9 #include <omp.h>
10 #endif
11
12 #include "openmc/capi.h"
13 #include "openmc/constants.h"
14 #include "openmc/container_util.h"
15 #include "openmc/distribution.h"
16 #include "openmc/distribution_multi.h"
17 #include "openmc/distribution_spatial.h"
18 #include "openmc/eigenvalue.h"
19 #include "openmc/error.h"
20 #include "openmc/file_utils.h"
21 #include "openmc/mesh.h"
22 #include "openmc/message_passing.h"
23 #include "openmc/output.h"
24 #include "openmc/random_lcg.h"
25 #include "openmc/simulation.h"
26 #include "openmc/source.h"
27 #include "openmc/string_utils.h"
28 #include "openmc/tallies/trigger.h"
29 #include "openmc/volume_calc.h"
30 #include "openmc/xml_interface.h"
31
32 namespace openmc {
33
34 //==============================================================================
35 // Global variables
36 //==============================================================================
37
38 namespace settings {
39
40 // Default values for boolean flags
41 bool assume_separate {false};
42 bool check_overlaps {false};
43 bool cmfd_run {false};
44 bool confidence_intervals {false};
45 bool create_fission_neutrons {true};
46 bool dagmc {false};
47 bool delayed_photon_scaling {true};
48 bool entropy_on {false};
49 bool event_based {false};
50 bool legendre_to_tabular {true};
51 bool material_cell_offsets {true};
52 bool output_summary {true};
53 bool output_tallies {true};
54 bool particle_restart_run {false};
55 bool photon_transport {false};
56 bool reduce_tallies {true};
57 bool res_scat_on {false};
58 bool restart_run {false};
59 bool run_CE {true};
60 bool source_latest {false};
61 bool source_separate {false};
62 bool source_write {true};
63 bool surf_source_write {false};
64 bool surf_source_read {false};
65 bool survival_biasing {false};
66 bool temperature_multipole {false};
67 bool trigger_on {false};
68 bool trigger_predict {false};
69 bool ufs_on {false};
70 bool urr_ptables_on {true};
71 bool write_all_tracks {false};
72 bool write_initial_source {false};
73
74 std::string path_cross_sections;
75 std::string path_input;
76 std::string path_output;
77 std::string path_particle_restart;
78 std::string path_sourcepoint;
79 std::string path_statepoint;
80
81 int32_t n_inactive {0};
82 int32_t max_lost_particles {10};
83 double rel_max_lost_particles {1.0e-6};
84 int32_t gen_per_batch {1};
85 int64_t n_particles {-1};
86
87 int64_t max_particles_in_flight {100000};
88
89 ElectronTreatment electron_treatment {ElectronTreatment::TTB};
90 array<double, 4> energy_cutoff {0.0, 1000.0, 0.0, 0.0};
91 int legendre_to_tabular_points {C_NONE};
92 int max_order {0};
93 int n_log_bins {8000};
94 int n_batches;
95 int n_max_batches;
96 ResScatMethod res_scat_method {ResScatMethod::rvs};
97 double res_scat_energy_min {0.01};
98 double res_scat_energy_max {1000.0};
99 vector<std::string> res_scat_nuclides;
100 RunMode run_mode {RunMode::UNSET};
101 std::unordered_set<int> sourcepoint_batch;
102 std::unordered_set<int> statepoint_batch;
103 std::unordered_set<int> source_write_surf_id;
104 int64_t max_surface_particles;
105 TemperatureMethod temperature_method {TemperatureMethod::NEAREST};
106 double temperature_tolerance {10.0};
107 double temperature_default {293.6};
108 array<double, 2> temperature_range {0.0, 0.0};
109 int trace_batch;
110 int trace_gen;
111 int64_t trace_particle;
112 vector<array<int, 3>> track_identifiers;
113 int trigger_batch_interval {1};
114 int verbosity {7};
115 double weight_cutoff {0.25};
116 double weight_survive {1.0};
117
118 } // namespace settings
119
120 //==============================================================================
121 // Functions
122 //==============================================================================
123
get_run_parameters(pugi::xml_node node_base)124 void get_run_parameters(pugi::xml_node node_base)
125 {
126 using namespace settings;
127 using namespace pugi;
128
129 // Check number of particles
130 if (!check_for_node(node_base, "particles")) {
131 fatal_error("Need to specify number of particles.");
132 }
133
134 // Get number of particles if it wasn't specified as a command-line argument
135 if (n_particles == -1) {
136 n_particles = std::stoll(get_node_value(node_base, "particles"));
137 }
138
139 // Get maximum number of in flight particles for event-based mode
140 if (check_for_node(node_base, "max_particles_in_flight")) {
141 max_particles_in_flight = std::stoll(get_node_value(node_base,
142 "max_particles_in_flight"));
143 }
144
145 // Get number of basic batches
146 if (check_for_node(node_base, "batches")) {
147 n_batches = std::stoi(get_node_value(node_base, "batches"));
148 }
149 if (!trigger_on) n_max_batches = n_batches;
150
151 // Get max number of lost particles
152 if (check_for_node(node_base, "max_lost_particles")) {
153 max_lost_particles = std::stoi(get_node_value(node_base, "max_lost_particles"));
154 }
155
156 // Get relative number of lost particles
157 if (check_for_node(node_base, "rel_max_lost_particles")) {
158 rel_max_lost_particles = std::stod(get_node_value(node_base, "rel_max_lost_particles"));
159 }
160
161 // Get number of inactive batches
162 if (run_mode == RunMode::EIGENVALUE) {
163 if (check_for_node(node_base, "inactive")) {
164 n_inactive = std::stoi(get_node_value(node_base, "inactive"));
165 }
166 if (check_for_node(node_base, "generations_per_batch")) {
167 gen_per_batch = std::stoi(get_node_value(node_base, "generations_per_batch"));
168 }
169
170 // Preallocate space for keff and entropy by generation
171 int m = settings::n_max_batches * settings::gen_per_batch;
172 simulation::k_generation.reserve(m);
173 simulation::entropy.reserve(m);
174
175 // Get the trigger information for keff
176 if (check_for_node(node_base, "keff_trigger")) {
177 xml_node node_keff_trigger = node_base.child("keff_trigger");
178
179 if (check_for_node(node_keff_trigger, "type")) {
180 auto temp = get_node_value(node_keff_trigger, "type", true, true);
181 if (temp == "std_dev") {
182 keff_trigger.metric = TriggerMetric::standard_deviation;
183 } else if (temp == "variance") {
184 keff_trigger.metric = TriggerMetric::variance;
185 } else if (temp == "rel_err") {
186 keff_trigger.metric = TriggerMetric::relative_error;
187 } else {
188 fatal_error("Unrecognized keff trigger type " + temp);
189 }
190 } else {
191 fatal_error("Specify keff trigger type in settings XML");
192 }
193
194 if (check_for_node(node_keff_trigger, "threshold")) {
195 keff_trigger.threshold = std::stod(get_node_value(
196 node_keff_trigger, "threshold"));
197 if (keff_trigger.threshold <= 0) {
198 fatal_error("keff trigger threshold must be positive");
199 }
200 } else {
201 fatal_error("Specify keff trigger threshold in settings XML");
202 }
203 }
204 }
205 }
206
read_settings_xml()207 void read_settings_xml()
208 {
209 using namespace settings;
210 using namespace pugi;
211
212 // Check if settings.xml exists
213 std::string filename = path_input + "settings.xml";
214 if (!file_exists(filename)) {
215 if (run_mode != RunMode::PLOTTING) {
216 fatal_error(fmt::format(
217 "Settings XML file '{}' does not exist! In order "
218 "to run OpenMC, you first need a set of input files; at a minimum, this "
219 "includes settings.xml, geometry.xml, and materials.xml. Please consult "
220 "the user's guide at https://docs.openmc.org for further "
221 "information.", filename
222 ));
223 } else {
224 // The settings.xml file is optional if we just want to make a plot.
225 return;
226 }
227 }
228
229 // Parse settings.xml file
230 xml_document doc;
231 auto result = doc.load_file(filename.c_str());
232 if (!result) {
233 fatal_error("Error processing settings.xml file.");
234 }
235
236 // Get root element
237 xml_node root = doc.document_element();
238
239 // Verbosity
240 if (check_for_node(root, "verbosity")) {
241 verbosity = std::stoi(get_node_value(root, "verbosity"));
242 }
243
244 // DAGMC geometry check
245 if (check_for_node(root, "dagmc")) {
246 dagmc = get_node_value_bool(root, "dagmc");
247 }
248
249 #ifndef DAGMC
250 if (dagmc) {
251 fatal_error("DAGMC mode unsupported for this build of OpenMC");
252 }
253 #endif
254
255 // To this point, we haven't displayed any output since we didn't know what
256 // the verbosity is. Now that we checked for it, show the title if necessary
257 if (mpi::master) {
258 if (verbosity >= 2) title();
259 }
260 write_message("Reading settings XML file...", 5);
261
262 // Find if a multi-group or continuous-energy simulation is desired
263 if (check_for_node(root, "energy_mode")) {
264 std::string temp_str = get_node_value(root, "energy_mode", true, true);
265 if (temp_str == "mg" || temp_str == "multi-group") {
266 run_CE = false;
267 } else if (temp_str == "ce" || temp_str == "continuous-energy") {
268 run_CE = true;
269 }
270 }
271
272 // Look for deprecated cross_sections.xml file in settings.xml
273 if (check_for_node(root, "cross_sections")) {
274 warning("Setting cross_sections in settings.xml has been deprecated."
275 " The cross_sections are now set in materials.xml and the "
276 "cross_sections input to materials.xml and the OPENMC_CROSS_SECTIONS"
277 " environment variable will take precendent over setting "
278 "cross_sections in settings.xml.");
279 path_cross_sections = get_node_value(root, "cross_sections");
280 }
281
282 if (!run_CE) {
283 // Scattering Treatments
284 if (check_for_node(root, "max_order")) {
285 max_order = std::stoi(get_node_value(root, "max_order"));
286 } else {
287 // Set to default of largest int - 1, which means to use whatever is
288 // contained in library. This is largest int - 1 because for legendre
289 // scattering, a value of 1 is added to the order; adding 1 to the largest
290 // int gets you the largest negative integer, which is not what we want.
291 max_order = std::numeric_limits<int>::max() - 1;
292 }
293 }
294
295 // Check for a trigger node and get trigger information
296 if (check_for_node(root, "trigger")) {
297 xml_node node_trigger = root.child("trigger");
298
299 // Check if trigger(s) are to be turned on
300 trigger_on = get_node_value_bool(node_trigger, "active");
301
302 if (trigger_on) {
303 if (check_for_node(node_trigger, "max_batches") ){
304 n_max_batches = std::stoi(get_node_value(node_trigger, "max_batches"));
305 } else {
306 fatal_error("<max_batches> must be specified with triggers");
307 }
308
309 // Get the batch interval to check triggers
310 if (!check_for_node(node_trigger, "batch_interval")){
311 trigger_predict = true;
312 } else {
313 trigger_batch_interval = std::stoi(get_node_value(node_trigger, "batch_interval"));
314 if (trigger_batch_interval <= 0) {
315 fatal_error("Trigger batch interval must be greater than zero");
316 }
317 }
318 }
319 }
320
321 // Check run mode if it hasn't been set from the command line
322 xml_node node_mode;
323 if (run_mode == RunMode::UNSET) {
324 if (check_for_node(root, "run_mode")) {
325 std::string temp_str = get_node_value(root, "run_mode", true, true);
326 if (temp_str == "eigenvalue") {
327 run_mode = RunMode::EIGENVALUE;
328 } else if (temp_str == "fixed source") {
329 run_mode = RunMode::FIXED_SOURCE;
330 } else if (temp_str == "plot") {
331 run_mode = RunMode::PLOTTING;
332 } else if (temp_str == "particle restart") {
333 run_mode = RunMode::PARTICLE;
334 } else if (temp_str == "volume") {
335 run_mode = RunMode::VOLUME;
336 } else {
337 fatal_error("Unrecognized run mode: " + temp_str);
338 }
339
340 // Assume XML specifies <particles>, <batches>, etc. directly
341 node_mode = root;
342 } else {
343 warning("<run_mode> should be specified.");
344
345 // Make sure that either eigenvalue or fixed source was specified
346 node_mode = root.child("eigenvalue");
347 if (node_mode) {
348 run_mode = RunMode::EIGENVALUE;
349 } else {
350 node_mode = root.child("fixed_source");
351 if (node_mode) {
352 run_mode = RunMode::FIXED_SOURCE;
353 } else {
354 fatal_error("<eigenvalue> or <fixed_source> not specified.");
355 }
356 }
357 }
358 }
359
360 if (run_mode == RunMode::EIGENVALUE || run_mode == RunMode::FIXED_SOURCE) {
361 // Read run parameters
362 get_run_parameters(node_mode);
363
364 // Check number of active batches, inactive batches, max lost particles and particles
365 if (n_batches <= n_inactive) {
366 fatal_error("Number of active batches must be greater than zero.");
367 } else if (n_inactive < 0) {
368 fatal_error("Number of inactive batches must be non-negative.");
369 } else if (n_particles <= 0) {
370 fatal_error("Number of particles must be greater than zero.");
371 } else if (max_lost_particles <= 0) {
372 fatal_error("Number of max lost particles must be greater than zero.");
373 } else if (rel_max_lost_particles <= 0.0 || rel_max_lost_particles >= 1.0) {
374 fatal_error("Relative max lost particles must be between zero and one.");
375 }
376 }
377
378 // Copy random number seed if specified
379 if (check_for_node(root, "seed")) {
380 auto seed = std::stoll(get_node_value(root, "seed"));
381 openmc_set_seed(seed);
382 }
383
384 // Check for electron treatment
385 if (check_for_node(root, "electron_treatment")) {
386 auto temp_str = get_node_value(root, "electron_treatment", true, true);
387 if (temp_str == "led") {
388 electron_treatment = ElectronTreatment::LED;
389 } else if (temp_str == "ttb") {
390 electron_treatment = ElectronTreatment::TTB;
391 } else {
392 fatal_error("Unrecognized electron treatment: " + temp_str + ".");
393 }
394 }
395
396 // Check for photon transport
397 if (check_for_node(root, "photon_transport")) {
398 photon_transport = get_node_value_bool(root, "photon_transport");
399
400 if (!run_CE && photon_transport) {
401 fatal_error("Photon transport is not currently supported in "
402 "multigroup mode");
403 }
404 }
405
406 // Number of bins for logarithmic grid
407 if (check_for_node(root, "log_grid_bins")) {
408 n_log_bins = std::stoi(get_node_value(root, "log_grid_bins"));
409 if (n_log_bins < 1) {
410 fatal_error("Number of bins for logarithmic grid must be greater "
411 "than zero.");
412 }
413 }
414
415 // Number of OpenMP threads
416 if (check_for_node(root, "threads")) {
417 if (mpi::master) warning("The <threads> element has been deprecated. Use "
418 "the OMP_NUM_THREADS environment variable to set the number of threads.");
419 }
420
421 // ==========================================================================
422 // EXTERNAL SOURCE
423
424 // Get point to list of <source> elements and make sure there is at least one
425 for (pugi::xml_node node : root.children("source")) {
426 if (check_for_node(node, "file")) {
427 auto path = get_node_value(node, "file", false, true);
428 model::external_sources.push_back(make_unique<FileSource>(path));
429 } else if (check_for_node(node, "library")) {
430 // Get shared library path and parameters
431 auto path = get_node_value(node, "library", false, true);
432 std::string parameters;
433 if (check_for_node(node, "parameters")) {
434 parameters = get_node_value(node, "parameters", false, true);
435 }
436
437 // Create custom source
438 model::external_sources.push_back(
439 make_unique<CustomSourceWrapper>(path, parameters));
440 } else {
441 model::external_sources.push_back(make_unique<IndependentSource>(node));
442 }
443 }
444
445 // Check if the user has specified to read surface source
446 if (check_for_node(root, "surf_source_read")) {
447 surf_source_read = true;
448 // Get surface source read node
449 xml_node node_ssr = root.child("surf_source_read");
450
451 std::string path = "surface_source.h5";
452 // Check if the user has specified different file for surface source reading
453 if (check_for_node(node_ssr, "path")) {
454 path = get_node_value(node_ssr, "path", false, true);
455 }
456 model::external_sources.push_back(make_unique<FileSource>(path));
457 }
458
459 // If no source specified, default to isotropic point source at origin with Watt spectrum
460 if (model::external_sources.empty()) {
461 model::external_sources.push_back(make_unique<IndependentSource>(
462 UPtrSpace {new SpatialPoint({0.0, 0.0, 0.0})},
463 UPtrAngle {new Isotropic()}, UPtrDist {new Watt(0.988e6, 2.249e-6)}));
464 }
465
466 // Check if we want to write out source
467 if (check_for_node(root, "write_initial_source")) {
468 write_initial_source = get_node_value_bool(root, "write_initial_source");
469 }
470
471 // Survival biasing
472 if (check_for_node(root, "survival_biasing")) {
473 survival_biasing = get_node_value_bool(root, "survival_biasing");
474 }
475
476 // Probability tables
477 if (check_for_node(root, "ptables")) {
478 urr_ptables_on = get_node_value_bool(root, "ptables");
479 }
480
481 // Cutoffs
482 if (check_for_node(root, "cutoff")) {
483 xml_node node_cutoff = root.child("cutoff");
484 if (check_for_node(node_cutoff, "weight")) {
485 weight_cutoff = std::stod(get_node_value(node_cutoff, "weight"));
486 }
487 if (check_for_node(node_cutoff, "weight_avg")) {
488 weight_survive = std::stod(get_node_value(node_cutoff, "weight_avg"));
489 }
490 if (check_for_node(node_cutoff, "energy_neutron")) {
491 energy_cutoff[0] = std::stod(get_node_value(node_cutoff, "energy_neutron"));
492 } else if (check_for_node(node_cutoff, "energy")) {
493 warning("The use of an <energy> cutoff is deprecated and should "
494 "be replaced by <energy_neutron>.");
495 energy_cutoff[0] = std::stod(get_node_value(node_cutoff, "energy"));
496 }
497 if (check_for_node(node_cutoff, "energy_photon")) {
498 energy_cutoff[1] = std::stod(get_node_value(node_cutoff, "energy_photon"));
499 }
500 if (check_for_node(node_cutoff, "energy_electron")) {
501 energy_cutoff[2] = std::stof(get_node_value(node_cutoff, "energy_electron"));
502 }
503 if (check_for_node(node_cutoff, "energy_positron")) {
504 energy_cutoff[3] = std::stod(get_node_value(node_cutoff, "energy_positron"));
505 }
506 }
507
508 // Particle trace
509 if (check_for_node(root, "trace")) {
510 auto temp = get_node_array<int64_t>(root, "trace");
511 if (temp.size() != 3) {
512 fatal_error("Must provide 3 integers for <trace> that specify the "
513 "batch, generation, and particle number.");
514 }
515 trace_batch = temp.at(0);
516 trace_gen = temp.at(1);
517 trace_particle = temp.at(2);
518 }
519
520 // Particle tracks
521 if (check_for_node(root, "track")) {
522 // Get values and make sure there are three per particle
523 auto temp = get_node_array<int>(root, "track");
524 if (temp.size() % 3 != 0) {
525 fatal_error("Number of integers specified in 'track' is not "
526 "divisible by 3. Please provide 3 integers per particle to be "
527 "tracked.");
528 }
529
530 // Reshape into track_identifiers
531 int n_tracks = temp.size() / 3;
532 for (int i = 0; i < n_tracks; ++i) {
533 track_identifiers.push_back({temp[3*i], temp[3*i + 1],
534 temp[3*i + 2]});
535 }
536 }
537
538 // Read meshes
539 read_meshes(root);
540
541 // Shannon Entropy mesh
542 if (check_for_node(root, "entropy_mesh")) {
543 int temp = std::stoi(get_node_value(root, "entropy_mesh"));
544 if (model::mesh_map.find(temp) == model::mesh_map.end()) {
545 fatal_error(fmt::format(
546 "Mesh {} specified for Shannon entropy does not exist.", temp));
547 }
548
549 auto* m = dynamic_cast<RegularMesh*>(
550 model::meshes[model::mesh_map.at(temp)].get());
551 if (!m) fatal_error("Only regular meshes can be used as an entropy mesh");
552 simulation::entropy_mesh = m;
553
554 // Turn on Shannon entropy calculation
555 entropy_on = true;
556
557 } else if (check_for_node(root, "entropy")) {
558 fatal_error("Specifying a Shannon entropy mesh via the <entropy> element "
559 "is deprecated. Please create a mesh using <mesh> and then reference "
560 "it by specifying its ID in an <entropy_mesh> element.");
561 }
562
563 // Uniform fission source weighting mesh
564 if (check_for_node(root, "ufs_mesh")) {
565 auto temp = std::stoi(get_node_value(root, "ufs_mesh"));
566 if (model::mesh_map.find(temp) == model::mesh_map.end()) {
567 fatal_error(fmt::format("Mesh {} specified for uniform fission site "
568 "method does not exist.", temp));
569 }
570
571 auto* m = dynamic_cast<RegularMesh*>(model::meshes[model::mesh_map.at(temp)].get());
572 if (!m) fatal_error("Only regular meshes can be used as a UFS mesh");
573 simulation::ufs_mesh = m;
574
575 // Turn on uniform fission source weighting
576 ufs_on = true;
577
578 } else if (check_for_node(root, "uniform_fs")) {
579 fatal_error("Specifying a UFS mesh via the <uniform_fs> element "
580 "is deprecated. Please create a mesh using <mesh> and then reference "
581 "it by specifying its ID in a <ufs_mesh> element.");
582 }
583
584 // Check if the user has specified to write state points
585 if (check_for_node(root, "state_point")) {
586
587 // Get pointer to state_point node
588 auto node_sp = root.child("state_point");
589
590 // Determine number of batches at which to store state points
591 if (check_for_node(node_sp, "batches")) {
592 // User gave specific batches to write state points
593 auto temp = get_node_array<int>(node_sp, "batches");
594 for (const auto& b : temp) {
595 statepoint_batch.insert(b);
596 }
597 } else {
598 // If neither were specified, write state point at last batch
599 statepoint_batch.insert(n_batches);
600 }
601 } else {
602 // If no <state_point> tag was present, by default write state point at
603 // last batch only
604 statepoint_batch.insert(n_batches);
605 }
606
607 // Check if the user has specified to write source points
608 if (check_for_node(root, "source_point")) {
609 // Get source_point node
610 xml_node node_sp = root.child("source_point");
611
612 // Determine batches at which to store source points
613 if (check_for_node(node_sp, "batches")) {
614 // User gave specific batches to write source points
615 auto temp = get_node_array<int>(node_sp, "batches");
616 for (const auto& b : temp) {
617 sourcepoint_batch.insert(b);
618 }
619 } else {
620 // If neither were specified, write source points with state points
621 sourcepoint_batch = statepoint_batch;
622 }
623
624 // Check if the user has specified to write binary source file
625 if (check_for_node(node_sp, "separate")) {
626 source_separate = get_node_value_bool(node_sp, "separate");
627 }
628 if (check_for_node(node_sp, "write")) {
629 source_write = get_node_value_bool(node_sp, "write");
630 }
631 if (check_for_node(node_sp, "overwrite_latest")) {
632 source_latest = get_node_value_bool(node_sp, "overwrite_latest");
633 source_separate = source_latest;
634 }
635 } else {
636 // If no <source_point> tag was present, by default we keep source bank in
637 // statepoint file and write it out at statepoints intervals
638 source_separate = false;
639 sourcepoint_batch = statepoint_batch;
640 }
641
642 // Check if the user has specified to write surface source
643 if (check_for_node(root, "surf_source_write")) {
644 surf_source_write = true;
645 // Get surface source write node
646 xml_node node_ssw = root.child("surf_source_write");
647
648 // Determine surface ids at which crossing particles are to be banked
649 if (check_for_node(node_ssw, "surface_ids")) {
650 auto temp = get_node_array<int>(node_ssw, "surface_ids");
651 for (const auto& b : temp) {
652 source_write_surf_id.insert(b);
653 }
654 }
655
656 // Get maximum number of particles to be banked per surface
657 if (check_for_node(node_ssw, "max_particles")) {
658 max_surface_particles = std::stoll(get_node_value(node_ssw, "max_particles"));
659 }
660 }
661
662 // If source is not seperate and is to be written out in the statepoint file,
663 // make sure that the sourcepoint batch numbers are contained in the
664 // statepoint list
665 if (!source_separate) {
666 for (const auto& b : sourcepoint_batch) {
667 if (!contains(statepoint_batch, b)) {
668 fatal_error("Sourcepoint batches are not a subset of statepoint batches.");
669 }
670 }
671 }
672
673 // Check if the user has specified to not reduce tallies at the end of every
674 // batch
675 if (check_for_node(root, "no_reduce")) {
676 reduce_tallies = !get_node_value_bool(root, "no_reduce");
677 }
678
679 // Check if the user has specified to use confidence intervals for
680 // uncertainties rather than standard deviations
681 if (check_for_node(root, "confidence_intervals")) {
682 confidence_intervals = get_node_value_bool(root, "confidence_intervals");
683 }
684
685 // Check for output options
686 if (check_for_node(root, "output")) {
687 // Get pointer to output node
688 pugi::xml_node node_output = root.child("output");
689
690 // Check for summary option
691 if (check_for_node(node_output, "summary")) {
692 output_summary = get_node_value_bool(node_output, "summary");
693 }
694
695 // Check for ASCII tallies output option
696 if (check_for_node(node_output, "tallies")) {
697 output_tallies = get_node_value_bool(node_output, "tallies");
698 }
699
700 // Set output directory if a path has been specified
701 if (check_for_node(node_output, "path")) {
702 path_output = get_node_value(node_output, "path");
703 if (!ends_with(path_output, "/")) {
704 path_output += "/";
705 }
706 }
707 }
708
709 // Resonance scattering parameters
710 if (check_for_node(root, "resonance_scattering")) {
711 xml_node node_res_scat = root.child("resonance_scattering");
712
713 // See if resonance scattering is enabled
714 if (check_for_node(node_res_scat, "enable")) {
715 res_scat_on = get_node_value_bool(node_res_scat, "enable");
716 } else {
717 res_scat_on = true;
718 }
719
720 // Determine what method is used
721 if (check_for_node(node_res_scat, "method")) {
722 auto temp = get_node_value(node_res_scat, "method", true, true);
723 if (temp == "rvs") {
724 res_scat_method = ResScatMethod::rvs;
725 } else if (temp == "dbrc") {
726 res_scat_method = ResScatMethod::dbrc;
727 } else {
728 fatal_error("Unrecognized resonance elastic scattering method: "
729 + temp + ".");
730 }
731 }
732
733 // Minimum energy for resonance scattering
734 if (check_for_node(node_res_scat, "energy_min")) {
735 res_scat_energy_min = std::stod(get_node_value(node_res_scat, "energy_min"));
736 }
737 if (res_scat_energy_min < 0.0) {
738 fatal_error("Lower resonance scattering energy bound is negative");
739 }
740
741 // Maximum energy for resonance scattering
742 if (check_for_node(node_res_scat, "energy_max")) {
743 res_scat_energy_max = std::stod(get_node_value(node_res_scat, "energy_max"));
744 }
745 if (res_scat_energy_max < res_scat_energy_min) {
746 fatal_error("Upper resonance scattering energy bound is below the "
747 "lower resonance scattering energy bound.");
748 }
749
750 // Get resonance scattering nuclides
751 if (check_for_node(node_res_scat, "nuclides")) {
752 res_scat_nuclides = get_node_array<std::string>(node_res_scat, "nuclides");
753 }
754 }
755
756 // Get volume calculations
757 for (pugi::xml_node node_vol : root.children("volume_calc")) {
758 model::volume_calcs.emplace_back(node_vol);
759 }
760
761 // Get temperature settings
762 if (check_for_node(root, "temperature_default")) {
763 temperature_default = std::stod(get_node_value(root, "temperature_default"));
764 }
765 if (check_for_node(root, "temperature_method")) {
766 auto temp = get_node_value(root, "temperature_method", true, true);
767 if (temp == "nearest") {
768 temperature_method = TemperatureMethod::NEAREST;
769 } else if (temp == "interpolation") {
770 temperature_method = TemperatureMethod::INTERPOLATION;
771 } else {
772 fatal_error("Unknown temperature method: " + temp);
773 }
774 }
775 if (check_for_node(root, "temperature_tolerance")) {
776 temperature_tolerance = std::stod(get_node_value(root, "temperature_tolerance"));
777 }
778 if (check_for_node(root, "temperature_multipole")) {
779 temperature_multipole = get_node_value_bool(root, "temperature_multipole");
780 }
781 if (check_for_node(root, "temperature_range")) {
782 auto range = get_node_array<double>(root, "temperature_range");
783 temperature_range[0] = range.at(0);
784 temperature_range[1] = range.at(1);
785 }
786
787 // Check for tabular_legendre options
788 if (check_for_node(root, "tabular_legendre")) {
789 // Get pointer to tabular_legendre node
790 xml_node node_tab_leg = root.child("tabular_legendre");
791
792 // Check for enable option
793 if (check_for_node(node_tab_leg, "enable")) {
794 legendre_to_tabular = get_node_value_bool(node_tab_leg, "enable");
795 }
796
797 // Check for the number of points
798 if (check_for_node(node_tab_leg, "num_points")) {
799 legendre_to_tabular_points = std::stoi(get_node_value(
800 node_tab_leg, "num_points"));
801 if (legendre_to_tabular_points <= 1 && !run_CE) {
802 fatal_error("The 'num_points' subelement/attribute of the "
803 "<tabular_legendre> element must contain a value greater than 1");
804 }
805 }
806 }
807
808 // Check whether create fission sites
809 if (run_mode == RunMode::FIXED_SOURCE) {
810 if (check_for_node(root, "create_fission_neutrons")) {
811 create_fission_neutrons = get_node_value_bool(root, "create_fission_neutrons");
812 }
813 }
814
815 // Check whether to scale fission photon yields
816 if (check_for_node(root, "delayed_photon_scaling")) {
817 delayed_photon_scaling = get_node_value_bool(root, "delayed_photon_scaling");
818 }
819
820 // Check whether to use event-based parallelism
821 if (check_for_node(root, "event_based")) {
822 event_based = get_node_value_bool(root, "event_based");
823 }
824
825 // Check whether material cell offsets should be generated
826 if (check_for_node(root, "material_cell_offsets")) {
827 material_cell_offsets = get_node_value_bool(root, "material_cell_offsets");
828 }
829 }
830
free_memory_settings()831 void free_memory_settings() {
832 settings::statepoint_batch.clear();
833 settings::sourcepoint_batch.clear();
834 settings::source_write_surf_id.clear();
835 settings::res_scat_nuclides.clear();
836 }
837
838 //==============================================================================
839 // C API functions
840 //==============================================================================
841
842 extern "C" int
openmc_set_n_batches(int32_t n_batches,bool set_max_batches,bool add_statepoint_batch)843 openmc_set_n_batches(int32_t n_batches, bool set_max_batches,
844 bool add_statepoint_batch)
845 {
846 if (settings::n_inactive >= n_batches) {
847 set_errmsg("Number of active batches must be greater than zero.");
848 return OPENMC_E_INVALID_ARGUMENT;
849 }
850
851 if (simulation::current_batch >= n_batches) {
852 set_errmsg("Number of batches must be greater than current batch.");
853 return OPENMC_E_INVALID_ARGUMENT;
854 }
855
856 if (!settings::trigger_on) {
857 // Set n_batches and n_max_batches to same value
858 settings::n_batches = n_batches;
859 settings::n_max_batches = n_batches;
860 } else {
861 // Set n_batches and n_max_batches based on value of set_max_batches
862 if (set_max_batches) {
863 settings::n_max_batches = n_batches;
864 } else {
865 settings::n_batches = n_batches;
866 }
867 }
868
869 // Update size of k_generation and entropy
870 int m = settings::n_max_batches * settings::gen_per_batch;
871 simulation::k_generation.reserve(m);
872 simulation::entropy.reserve(m);
873
874 // Add value of n_batches to statepoint_batch
875 if (add_statepoint_batch &&
876 !(contains(settings::statepoint_batch, n_batches)))
877 settings::statepoint_batch.insert(n_batches);
878
879 return 0;
880 }
881
882 extern "C" int
openmc_get_n_batches(int * n_batches,bool get_max_batches)883 openmc_get_n_batches(int* n_batches, bool get_max_batches)
884 {
885 *n_batches = get_max_batches ? settings::n_max_batches : settings::n_batches;
886
887 return 0;
888 }
889
890 } // namespace openmc
891