• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 13134187828

04 Feb 2025 11:13AM UTC coverage: 84.945% (+0.08%) from 84.869%
13134187828

Pull #3252

github

web-flow
Merge 6f397fdd5 into 59c398be8
Pull Request #3252: Adding vtkhdf option to write vtk data

80 of 92 new or added lines in 1 file covered. (86.96%)

1267 existing lines in 48 files now uncovered.

50221 of 59122 relevant lines covered (84.94%)

35221818.06 hits per line

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

86.5
/src/cross_sections.cpp
1
#include "openmc/cross_sections.h"
2

3
#include "openmc/capi.h"
4
#include "openmc/constants.h"
5
#include "openmc/container_util.h"
6
#include "openmc/error.h"
7
#include "openmc/file_utils.h"
8
#include "openmc/geometry_aux.h"
9
#include "openmc/hdf5_interface.h"
10
#include "openmc/material.h"
11
#include "openmc/message_passing.h"
12
#include "openmc/mgxs_interface.h"
13
#include "openmc/nuclide.h"
14
#include "openmc/photon.h"
15
#include "openmc/settings.h"
16
#include "openmc/simulation.h"
17
#include "openmc/string_utils.h"
18
#include "openmc/thermal.h"
19
#include "openmc/timer.h"
20
#include "openmc/wmp.h"
21
#include "openmc/xml_interface.h"
22

23
#include "pugixml.hpp"
24

25
#include <cstdlib> // for getenv
26
#include <filesystem>
27
#include <unordered_set>
28

29
namespace openmc {
30

31
//==============================================================================
32
// Global variable declarations
33
//==============================================================================
34

35
namespace data {
36

37
std::map<LibraryKey, std::size_t> library_map;
38
vector<Library> libraries;
39
} // namespace data
40

41
//==============================================================================
42
// Library methods
43
//==============================================================================
44

45
Library::Library(pugi::xml_node node, const std::string& directory)
5,275,698✔
46
{
47
  // Get type of library
48
  if (check_for_node(node, "type")) {
5,275,698✔
49
    auto type = get_node_value(node, "type");
5,275,698✔
50
    if (type == "neutron") {
5,275,698✔
51
      type_ = Type::neutron;
2,310,183✔
52
    } else if (type == "thermal") {
2,965,515✔
53
      type_ = Type::thermal;
109,412✔
54
    } else if (type == "photon") {
2,856,103✔
55
      type_ = Type::photon;
546,100✔
56
    } else if (type == "wmp") {
2,310,003✔
57
      type_ = Type::wmp;
2,310,003✔
58
    } else {
UNCOV
59
      fatal_error("Unrecognized library type: " + type);
×
60
    }
61
  } else {
5,275,698✔
UNCOV
62
    fatal_error("Missing library type");
×
63
  }
64

65
  // Get list of materials
66
  if (check_for_node(node, "materials")) {
5,275,698✔
67
    materials_ = get_node_array<std::string>(node, "materials");
5,275,698✔
68
  }
69

70
  // determine path of cross section table
71
  if (!check_for_node(node, "path")) {
5,275,698✔
UNCOV
72
    fatal_error("Missing library path");
×
73
  }
74
  std::string path = get_node_value(node, "path");
5,275,698✔
75

76
  if (starts_with(path, "/")) {
5,275,698✔
UNCOV
77
    path_ = path;
×
78
  } else if (ends_with(directory, "/")) {
5,275,698✔
UNCOV
79
    path_ = directory + path;
×
80
  } else if (!directory.empty()) {
5,275,698✔
81
    path_ = directory + "/" + path;
5,263,734✔
82
  } else {
83
    path_ = path;
11,964✔
84
  }
85

86
  if (!file_exists(path_)) {
5,275,698✔
UNCOV
87
    warning("Cross section library " + path_ + " does not exist.");
×
88
  }
89
}
5,275,698✔
90

91
//==============================================================================
92
// Non-member functions
93
//==============================================================================
94

95
void read_cross_sections_xml()
1,501✔
96
{
97
  pugi::xml_document doc;
1,501✔
98
  std::string filename = settings::path_input + "materials.xml";
1,501✔
99
  // Check if materials.xml exists
100
  if (!file_exists(filename)) {
1,501✔
UNCOV
101
    fatal_error("Material XML file '" + filename + "' does not exist.");
×
102
  }
103
  // Parse materials.xml file
104
  doc.load_file(filename.c_str());
1,501✔
105

106
  auto root = doc.document_element();
1,501✔
107

108
  read_cross_sections_xml(root);
1,501✔
109
}
1,501✔
110

111
void read_cross_sections_xml(pugi::xml_node root)
6,583✔
112
{
113
  // Find cross_sections.xml file -- the first place to look is the
114
  // materials.xml file. If no file is found there, then we check the
115
  // OPENMC_CROSS_SECTIONS environment variable
116
  if (!check_for_node(root, "cross_sections")) {
6,583✔
117
    // No cross_sections.xml file specified in settings.xml, check
118
    // environment variable
119
    if (settings::run_CE) {
5,449✔
120
      char* envvar = std::getenv("OPENMC_CROSS_SECTIONS");
5,449✔
121
      if (!envvar) {
5,449✔
UNCOV
122
        fatal_error(
×
123
          "No cross_sections.xml file was specified in "
124
          "materials.xml or in the OPENMC_CROSS_SECTIONS"
125
          " environment variable. OpenMC needs such a file to identify "
126
          "where to find data libraries. Please consult the"
127
          " user's guide at https://docs.openmc.org/ for "
128
          "information on how to set up data libraries.");
129
      }
130
      settings::path_cross_sections = envvar;
5,449✔
131
    } else {
132
      char* envvar = std::getenv("OPENMC_MG_CROSS_SECTIONS");
×
133
      if (!envvar) {
×
UNCOV
134
        fatal_error(
×
135
          "No mgxs.h5 file was specified in "
136
          "materials.xml or in the OPENMC_MG_CROSS_SECTIONS environment "
137
          "variable. OpenMC needs such a file to identify where to "
138
          "find MG cross section libraries. Please consult the user's "
139
          "guide at https://docs.openmc.org for information on "
140
          "how to set up MG cross section libraries.");
141
      }
UNCOV
142
      settings::path_cross_sections = envvar;
×
143
    }
144
  } else {
145
    settings::path_cross_sections = get_node_value(root, "cross_sections");
1,134✔
146

147
    // If no '/' found, the file is probably in the input directory
148
    auto pos = settings::path_cross_sections.rfind("/");
1,134✔
149
    if (pos == std::string::npos && !settings::path_input.empty()) {
1,134✔
150
      settings::path_cross_sections =
UNCOV
151
        settings::path_input + "/" + settings::path_cross_sections;
×
152
    }
153
  }
154

155
  // Now that the cross_sections.xml or mgxs.h5 has been located, read it in
156
  if (settings::run_CE) {
6,583✔
157
    read_ce_cross_sections_xml();
5,641✔
158
  } else {
159
    data::mg.read_header(settings::path_cross_sections);
942✔
160
    put_mgxs_header_data_to_globals();
942✔
161
  }
162

163
  // Establish mapping between (type, material) and index in libraries
164
  int i = 0;
6,583✔
165
  for (const auto& lib : data::libraries) {
5,284,716✔
166
    for (const auto& name : lib.materials_) {
10,556,266✔
167
      LibraryKey key {lib.type_, name};
5,278,133✔
168
      data::library_map.insert({key, i});
5,278,133✔
169
    }
5,278,133✔
170
    ++i;
5,278,133✔
171
  }
172

173
  // Check that 0K nuclides are listed in the cross_sections.xml file
174
  for (const auto& name : settings::res_scat_nuclides) {
6,634✔
175
    LibraryKey key {Library::Type::neutron, name};
51✔
176
    if (data::library_map.find(key) == data::library_map.end()) {
51✔
UNCOV
177
      fatal_error("Could not find resonant scatterer " + name +
×
178
                  " in cross_sections.xml file!");
179
    }
180
  }
51✔
181
}
6,583✔
182

183
void read_ce_cross_sections(const vector<vector<double>>& nuc_temps,
5,641✔
184
  const vector<vector<double>>& thermal_temps)
185
{
186
  std::unordered_set<std::string> already_read;
5,641✔
187

188
  // Construct a vector of nuclide names because we haven't loaded nuclide data
189
  // yet, but we need to know the name of the i-th nuclide
190
  vector<std::string> nuclide_names(data::nuclide_map.size());
5,641✔
191
  vector<std::string> thermal_names(data::thermal_scatt_map.size());
5,641✔
192
  for (const auto& kv : data::nuclide_map) {
33,137✔
193
    nuclide_names[kv.second] = kv.first;
27,496✔
194
  }
195
  for (const auto& kv : data::thermal_scatt_map) {
6,769✔
196
    thermal_names[kv.second] = kv.first;
1,128✔
197
  }
198

199
  // Read cross sections
200
  for (const auto& mat : model::materials) {
18,232✔
201
    for (int i_nuc : mat->nuclide_) {
72,105✔
202
      // Find name of corresponding nuclide. Because we haven't actually loaded
203
      // data, we don't have the name available, so instead we search through
204
      // all key/value pairs in nuclide_map
205
      std::string& name = nuclide_names[i_nuc];
59,514✔
206

207
      // If we've already read this nuclide, skip it
208
      if (already_read.find(name) != already_read.end())
59,514✔
209
        continue;
32,018✔
210

211
      const auto& temps = nuc_temps[i_nuc];
27,496✔
212
      int err = openmc_load_nuclide(name.c_str(), temps.data(), temps.size());
27,496✔
213
      if (err < 0)
27,496✔
UNCOV
214
        throw std::runtime_error {openmc_err_msg};
×
215

216
      already_read.insert(name);
27,496✔
217
    }
218
  }
219

220
  // Perform final tasks -- reading S(a,b) tables, normalizing densities
221
  for (auto& mat : model::materials) {
18,232✔
222
    for (const auto& table : mat->thermal_tables_) {
14,489✔
223
      // Get name of S(a,b) table
224
      int i_table = table.index_table;
1,898✔
225
      std::string& name = thermal_names[i_table];
1,898✔
226

227
      if (already_read.find(name) == already_read.end()) {
1,898✔
228
        LibraryKey key {Library::Type::thermal, name};
1,128✔
229
        int idx = data::library_map[key];
1,128✔
230
        std::string& filename = data::libraries[idx].path_;
1,128✔
231

232
        write_message(6, "Reading {} from {}", name, filename);
1,128✔
233

234
        // Open file and make sure version matches
235
        hid_t file_id = file_open(filename, 'r');
1,128✔
236
        check_data_version(file_id);
1,128✔
237

238
        // Read thermal scattering data from HDF5
239
        hid_t group = open_group(file_id, name.c_str());
1,128✔
240
        data::thermal_scatt.push_back(
1,128✔
241
          make_unique<ThermalScattering>(group, thermal_temps[i_table]));
2,256✔
242
        close_group(group);
1,128✔
243
        file_close(file_id);
1,128✔
244

245
        // Add name to dictionary
246
        already_read.insert(name);
1,128✔
247
      }
1,128✔
248
    } // thermal_tables_
249

250
    // Finish setting up materials (normalizing densities, etc.)
251
    mat->finalize();
12,591✔
252
  } // materials
253

254
  if (settings::photon_transport &&
5,641✔
255
      settings::electron_treatment == ElectronTreatment::TTB) {
266✔
256
    // Take logarithm of energies since they are log-log interpolated
257
    data::ttb_e_grid = xt::log(data::ttb_e_grid);
254✔
258
  }
259

260
  // Show minimum/maximum temperature
261
  write_message(
5,641✔
262
    4, "Minimum neutron data temperature: {} K", data::temperature_min);
263
  write_message(
5,641✔
264
    4, "Maximum neutron data temperature: {} K", data::temperature_max);
265

266
  // If the user wants multipole, make sure we found a multipole library.
267
  if (settings::temperature_multipole) {
5,641✔
268
    bool mp_found = false;
34✔
269
    for (const auto& nuc : data::nuclides) {
34✔
270
      if (nuc->multipole_) {
34✔
271
        mp_found = true;
34✔
272
        break;
34✔
273
      }
274
    }
275
    if (mpi::master && !mp_found) {
34✔
UNCOV
276
      warning("Windowed multipole functionality is turned on, but no multipole "
×
277
              "libraries were found. Make sure that windowed multipole data is "
278
              "present in your cross_sections.xml file.");
279
    }
280
  }
281
}
5,641✔
282

283
void read_ce_cross_sections_xml()
5,641✔
284
{
285
  // Check if cross_sections.xml exists
286
  std::filesystem::path filename(settings::path_cross_sections);
5,641✔
287
  if (!std::filesystem::exists(filename)) {
5,641✔
UNCOV
288
    fatal_error(
×
UNCOV
289
      "Cross sections XML file '" + filename.string() + "' does not exist.");
×
290
  }
291

292
  if (std::filesystem::is_directory(filename)) {
5,641✔
UNCOV
293
    fatal_error("OPENMC_CROSS_SECTIONS is set to a directory. "
×
294
                "It should be set to an XML file.");
295
  }
296

297
  write_message("Reading cross sections XML file...", 5);
5,641✔
298

299
  // Parse cross_sections.xml file
300
  pugi::xml_document doc;
5,641✔
301
  auto result = doc.load_file(filename.c_str());
5,641✔
302
  if (!result) {
5,641✔
UNCOV
303
    fatal_error("Error processing cross_sections.xml file.");
×
304
  }
305
  auto root = doc.document_element();
5,641✔
306

307
  std::string directory;
5,641✔
308
  if (check_for_node(root, "directory")) {
5,641✔
309
    // Copy directory information if present
UNCOV
310
    directory = get_node_value(root, "directory");
×
311
  } else {
312
    // If no directory is listed in cross_sections.xml, by default select the
313
    // directory in which the cross_sections.xml file resides
314
    if (filename.has_parent_path()) {
5,641✔
315
      directory = filename.parent_path().string();
5,449✔
316
    } else {
317
      directory = settings::path_input;
192✔
318
    }
319
  }
320

321
  for (const auto& node_library : root.children("library")) {
5,281,339✔
322
    data::libraries.emplace_back(node_library, directory);
5,275,698✔
323
  }
324

325
  // Make sure file was not empty
326
  if (data::libraries.empty()) {
5,641✔
UNCOV
327
    fatal_error(
×
328
      "No cross section libraries present in cross_sections.xml file.");
329
  }
330
}
5,641✔
331

332
void finalize_cross_sections()
6,869✔
333
{
334
  if (settings::run_mode != RunMode::PLOTTING) {
6,869✔
335
    simulation::time_read_xs.start();
6,583✔
336
    if (settings::run_CE) {
6,583✔
337
      // Determine desired temperatures for each nuclide and S(a,b) table
338
      double_2dvec nuc_temps(data::nuclide_map.size());
5,641✔
339
      double_2dvec thermal_temps(data::thermal_scatt_map.size());
5,641✔
340
      get_temperatures(nuc_temps, thermal_temps);
5,641✔
341

342
      // Read continuous-energy cross sections from HDF5
343
      read_ce_cross_sections(nuc_temps, thermal_temps);
5,641✔
344
    } else {
5,641✔
345
      // Create material macroscopic data for MGXS
346
      set_mg_interface_nuclides_and_temps();
942✔
347
      data::mg.init();
942✔
348
      mark_fissionable_mgxs_materials();
942✔
349
    }
350
    simulation::time_read_xs.stop();
6,583✔
351
  }
352
}
6,869✔
353

354
void library_clear()
6,988✔
355
{
356
  data::libraries.clear();
6,988✔
357
  data::library_map.clear();
6,988✔
358
}
6,988✔
359

360
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc