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

openmc-dev / openmc / 18908039831

29 Oct 2025 12:34PM UTC coverage: 81.875% (+0.06%) from 81.812%
18908039831

Pull #3417

github

web-flow
Merge 0db1f0918 into 4c4176661
Pull Request #3417: Addition of a collision tracking feature

16866 of 23521 branches covered (71.71%)

Branch coverage included in aggregate %.

480 of 522 new or added lines in 13 files covered. (91.95%)

1 existing line in 1 file now uncovered.

54305 of 63405 relevant lines covered (85.65%)

42524539.39 hits per line

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

58.35
/src/mcpl_interface.cpp
1
#include "openmc/mcpl_interface.h"
2

3
#include "openmc/bank.h"
4
#include "openmc/error.h"
5
#include "openmc/file_utils.h"
6
#include "openmc/message_passing.h"
7
#include "openmc/settings.h"
8
#include "openmc/simulation.h"
9
#include "openmc/state_point.h"
10
#include "openmc/vector.h"
11

12
#include <fmt/core.h>
13

14
#include <cctype>
15
#include <cstdint>
16
#include <cstdio>
17
#include <cstdlib>
18
#include <cstring>
19
#include <memory>
20
#include <mutex>
21
#include <sstream>
22
#include <stdexcept>
23
#include <string>
24

25
#ifdef _WIN32
26
#define WIN32_LEAN_AND_MEAN
27
#include <windows.h>
28
#else
29
#include <dlfcn.h>
30
#endif
31

32
// WARNING: These declarations MUST EXACTLY MATCH the structure and function
33
// signatures of the libmcpl being loaded at runtime. Any discrepancy will
34
// likely lead to crashes or incorrect behavior. This is a maintenance risk.
35
// MCPL 2.2.0
36

37
#pragma pack(push, 1)
38
struct mcpl_particle_repr_t {
39
  double ekin;
40
  double polarisation[3];
41
  double position[3];
42
  double direction[3];
43
  double time;
44
  double weight;
45
  int32_t pdgcode;
46
  uint32_t userflags;
47
};
48
#pragma pack(pop)
49

50
// Opaque struct definitions replicating the MCPL C-API to ensure ABI
51
// compatibility without including mcpl.h. These must be kept in sync.
52
struct mcpl_file_t {
53
  void* internal;
54
};
55
struct mcpl_outfile_t {
56
  void* internal;
57
};
58

59
// Function pointer types for the dynamically loaded MCPL library
60
using mcpl_open_file_fpt = mcpl_file_t* (*)(const char* filename);
61
using mcpl_hdr_nparticles_fpt = uint64_t (*)(mcpl_file_t* file_handle);
62
using mcpl_read_fpt = const mcpl_particle_repr_t* (*)(mcpl_file_t* file_handle);
63
using mcpl_close_file_fpt = void (*)(mcpl_file_t* file_handle);
64

65
using mcpl_hdr_add_data_fpt = void (*)(mcpl_outfile_t* file_handle,
66
  const char* key, int32_t ldata, const char* data);
67
using mcpl_create_outfile_fpt = mcpl_outfile_t* (*)(const char* filename);
68
using mcpl_hdr_set_srcname_fpt = void (*)(
69
  mcpl_outfile_t* outfile_handle, const char* srcname);
70
using mcpl_add_particle_fpt = void (*)(
71
  mcpl_outfile_t* outfile_handle, const mcpl_particle_repr_t* particle);
72
using mcpl_close_outfile_fpt = void (*)(mcpl_outfile_t* outfile_handle);
73
using mcpl_hdr_add_stat_sum_fpt = void (*)(
74
  mcpl_outfile_t* outfile_handle, const char* key, double value);
75

76
namespace openmc {
77

78
#ifdef _WIN32
79
using LibraryHandleType = HMODULE;
80
#else
81
using LibraryHandleType = void*;
82
#endif
83

84
std::string get_last_library_error()
×
85
{
86
#ifdef _WIN32
87
  DWORD error_code = GetLastError();
88
  if (error_code == 0)
89
    return "No error reported by system."; // More accurate than "No error."
90
  LPSTR message_buffer = nullptr;
91
  size_t size =
92
    FormatMessageA(FORMAT_MESSAGE_ALLOCATE_BUFFER | FORMAT_MESSAGE_FROM_SYSTEM |
93
                     FORMAT_MESSAGE_IGNORE_INSERTS,
94
      NULL, error_code, MAKELANGID(LANG_NEUTRAL, SUBLANG_DEFAULT),
95
      (LPSTR)&message_buffer, 0, NULL);
96
  std::string message(message_buffer, size);
97
  LocalFree(message_buffer);
98
  while (
99
    !message.empty() && (message.back() == '\n' || message.back() == '\r')) {
100
    message.pop_back();
101
  }
102
  return message;
103
#else
104
  const char* err = dlerror();
×
105
  return err ? std::string(err) : "No error reported by dlerror.";
×
106
#endif
107
}
108

109
struct McplApi {
110
  mcpl_open_file_fpt open_file;
111
  mcpl_hdr_nparticles_fpt hdr_nparticles;
112
  mcpl_read_fpt read;
113
  mcpl_close_file_fpt close_file;
114
  mcpl_create_outfile_fpt create_outfile;
115
  mcpl_hdr_set_srcname_fpt hdr_set_srcname;
116
  mcpl_hdr_add_data_fpt hdr_add_data;
117
  mcpl_add_particle_fpt add_particle;
118
  mcpl_close_outfile_fpt close_outfile;
119
  mcpl_hdr_add_stat_sum_fpt hdr_add_stat_sum;
120

121
  explicit McplApi(LibraryHandleType lib_handle)
92✔
122
  {
92✔
123
    if (!lib_handle)
92!
124
      throw std::runtime_error(
×
125
        "MCPL library handle is null during API binding.");
×
126

127
    auto load_symbol_platform = [lib_handle](const char* name) {
920✔
128
      void* sym = nullptr;
920✔
129
#ifdef _WIN32
130
      sym = (void*)GetProcAddress(lib_handle, name);
131
#else
132
      sym = dlsym(lib_handle, name);
920✔
133
#endif
134
      if (!sym) {
920!
135
        throw std::runtime_error(
×
136
          fmt::format("Failed to load MCPL symbol '{}': {}", name,
×
137
            get_last_library_error()));
×
138
      }
139
      return sym;
920✔
140
    };
92✔
141

142
    open_file = reinterpret_cast<mcpl_open_file_fpt>(
92✔
143
      load_symbol_platform("mcpl_open_file"));
92✔
144
    hdr_nparticles = reinterpret_cast<mcpl_hdr_nparticles_fpt>(
92✔
145
      load_symbol_platform("mcpl_hdr_nparticles"));
92✔
146
    read = reinterpret_cast<mcpl_read_fpt>(load_symbol_platform("mcpl_read"));
92✔
147
    close_file = reinterpret_cast<mcpl_close_file_fpt>(
92✔
148
      load_symbol_platform("mcpl_close_file"));
92✔
149
    create_outfile = reinterpret_cast<mcpl_create_outfile_fpt>(
92✔
150
      load_symbol_platform("mcpl_create_outfile"));
92✔
151
    hdr_set_srcname = reinterpret_cast<mcpl_hdr_set_srcname_fpt>(
92✔
152
      load_symbol_platform("mcpl_hdr_set_srcname"));
92✔
153
    hdr_add_data = reinterpret_cast<mcpl_hdr_add_data_fpt>(
92✔
154
      load_symbol_platform("mcpl_hdr_add_data"));
92✔
155
    add_particle = reinterpret_cast<mcpl_add_particle_fpt>(
92✔
156
      load_symbol_platform("mcpl_add_particle"));
92✔
157
    close_outfile = reinterpret_cast<mcpl_close_outfile_fpt>(
92✔
158
      load_symbol_platform("mcpl_close_outfile"));
92✔
159

160
    // Try to load mcpl_hdr_add_stat_sum (available in MCPL >= 2.1.0)
161
    // Set to nullptr if not available for graceful fallback
162
    try {
163
      hdr_add_stat_sum = reinterpret_cast<mcpl_hdr_add_stat_sum_fpt>(
92✔
164
        load_symbol_platform("mcpl_hdr_add_stat_sum"));
92✔
165
    } catch (const std::runtime_error&) {
×
166
      hdr_add_stat_sum = nullptr;
×
167
    }
×
168
  }
92✔
169
};
170

171
static LibraryHandleType g_mcpl_lib_handle = nullptr;
172
static std::unique_ptr<McplApi> g_mcpl_api;
173
static bool g_mcpl_init_attempted = false;
174
static bool g_mcpl_successfully_loaded = false;
175
static std::string g_mcpl_load_error_msg;
176
static std::once_flag g_mcpl_init_flag;
177

178
void append_error(std::string& existing_msg, const std::string& new_error)
×
179
{
180
  if (!existing_msg.empty()) {
×
181
    existing_msg += "; ";
×
182
  }
183
  existing_msg += new_error;
×
184
}
×
185

186
void initialize_mcpl_interface_impl()
92✔
187
{
188
  g_mcpl_init_attempted = true;
92✔
189
  g_mcpl_load_error_msg.clear();
92✔
190

191
  // Try mcpl-config
192
  if (!g_mcpl_lib_handle) {
92!
193
    FILE* pipe = nullptr;
92✔
194
#ifdef _WIN32
195
    pipe = _popen("mcpl-config --show libpath", "r");
196
#else
197
    pipe = popen("mcpl-config --show libpath 2>/dev/null", "r");
92✔
198
#endif
199
    if (pipe) {
92!
200
      char buffer[512];
201
      if (fgets(buffer, sizeof(buffer), pipe) != nullptr) {
92!
202
        std::string shlibpath = buffer;
92✔
203
        // Remove trailing whitespace
204
        while (!shlibpath.empty() &&
368!
205
               std::isspace(static_cast<unsigned char>(shlibpath.back()))) {
184✔
206
          shlibpath.pop_back();
92✔
207
        }
208

209
        if (!shlibpath.empty()) {
92!
210
#ifdef _WIN32
211
          g_mcpl_lib_handle = LoadLibraryA(shlibpath.c_str());
212
#else
213
          g_mcpl_lib_handle = dlopen(shlibpath.c_str(), RTLD_LAZY);
92✔
214
#endif
215
          if (!g_mcpl_lib_handle) {
92!
216
            append_error(
×
217
              g_mcpl_load_error_msg, fmt::format("From mcpl-config ({}): {}",
×
218
                                       shlibpath, get_last_library_error()));
×
219
          }
220
        }
221
      }
92✔
222
#ifdef _WIN32
223
      _pclose(pipe);
224
#else
225
      pclose(pipe);
92✔
226
#endif
227
    } else { // pipe failed to open
228
      append_error(g_mcpl_load_error_msg,
×
229
        "mcpl-config command not found or failed to execute");
230
    }
231
  }
232

233
  // Try standard library names
234
  if (!g_mcpl_lib_handle) {
92!
235
#ifdef _WIN32
236
    const char* standard_names[] = {"mcpl.dll", "libmcpl.dll"};
237
#else
238
    const char* standard_names[] = {"libmcpl.so", "libmcpl.dylib"};
×
239
#endif
240
    for (const char* name : standard_names) {
×
241
#ifdef _WIN32
242
      g_mcpl_lib_handle = LoadLibraryA(name);
243
#else
244
      g_mcpl_lib_handle = dlopen(name, RTLD_LAZY);
×
245
#endif
246
      if (g_mcpl_lib_handle)
×
247
        break;
×
248
    }
249
    if (!g_mcpl_lib_handle) {
×
250
      append_error(
×
251
        g_mcpl_load_error_msg, fmt::format("Using standard names (e.g. {}): {}",
×
252
                                 standard_names[0], get_last_library_error()));
×
253
    }
254
  }
255

256
  if (!g_mcpl_lib_handle) {
92!
257
    if (mpi::master) {
×
258
      warning(fmt::format("MCPL library could not be loaded. MCPL-dependent "
×
259
                          "features will be unavailable. Load attempts: {}",
260
        g_mcpl_load_error_msg.empty()
×
261
          ? "No specific error during load attempts."
×
262
          : g_mcpl_load_error_msg));
263
    }
264
    g_mcpl_successfully_loaded = false;
×
265
    return;
×
266
  }
267

268
  try {
269
    g_mcpl_api = std::make_unique<McplApi>(g_mcpl_lib_handle);
92✔
270
    g_mcpl_successfully_loaded = true;
92✔
271
    // Do not call dlclose/FreeLibrary at exit. Leaking the handle is safer
272
    // and standard practice for libraries used for the application's lifetime.
273
  } catch (const std::runtime_error& e) {
×
274
    append_error(g_mcpl_load_error_msg,
×
275
      fmt::format(
×
276
        "MCPL library loaded, but failed to bind symbols: {}", e.what()));
×
277
    if (mpi::master) {
×
278
      warning(g_mcpl_load_error_msg);
×
279
    }
280
#ifdef _WIN32
281
    FreeLibrary(g_mcpl_lib_handle);
282
#else
283
    dlclose(g_mcpl_lib_handle);
×
284
#endif
285
    g_mcpl_lib_handle = nullptr;
×
286
    g_mcpl_successfully_loaded = false;
×
287
  }
×
288
}
289

290
void initialize_mcpl_interface_if_needed()
180✔
291
{
292
  std::call_once(g_mcpl_init_flag, initialize_mcpl_interface_impl);
180✔
293
}
180✔
294

295
bool is_mcpl_interface_available()
33✔
296
{
297
  initialize_mcpl_interface_if_needed();
33✔
298
  return g_mcpl_successfully_loaded;
33✔
299
}
300

301
inline void ensure_mcpl_ready_or_fatal()
147✔
302
{
303
  initialize_mcpl_interface_if_needed();
147✔
304
  if (!g_mcpl_successfully_loaded) {
147!
305
    fatal_error("MCPL functionality is required, but the MCPL library is not "
×
306
                "available or failed to initialize. Please ensure MCPL is "
307
                "installed and its library can be found (e.g., via PATH on "
308
                "Windows, LD_LIBRARY_PATH on Linux, or DYLD_LIBRARY_PATH on "
309
                "macOS). You can often install MCPL with 'pip install mcpl' or "
310
                "'conda install mcpl'.");
311
  }
312
}
147✔
313

314
SourceSite mcpl_particle_to_site(const mcpl_particle_repr_t* particle_repr)
32,000✔
315
{
316
  SourceSite site;
32,000✔
317
  switch (particle_repr->pdgcode) {
32,000!
318
  case 2112:
32,000✔
319
    site.particle = ParticleType::neutron;
32,000✔
320
    break;
32,000✔
321
  case 22:
×
322
    site.particle = ParticleType::photon;
×
323
    break;
×
324
  case 11:
×
325
    site.particle = ParticleType::electron;
×
326
    break;
×
327
  case -11:
×
328
    site.particle = ParticleType::positron;
×
329
    break;
×
330
  default:
×
331
    fatal_error(fmt::format(
×
332
      "MCPL: Encountered unexpected PDG code {} when converting to SourceSite.",
333
      particle_repr->pdgcode));
×
334
    break;
335
  }
336

337
  // Copy position and direction
338
  site.r.x = particle_repr->position[0];
32,000✔
339
  site.r.y = particle_repr->position[1];
32,000✔
340
  site.r.z = particle_repr->position[2];
32,000✔
341
  site.u.x = particle_repr->direction[0];
32,000✔
342
  site.u.y = particle_repr->direction[1];
32,000✔
343
  site.u.z = particle_repr->direction[2];
32,000✔
344
  // MCPL stores kinetic energy in [MeV], time in [ms]
345
  site.E = particle_repr->ekin * 1e6;
32,000✔
346
  site.time = particle_repr->time * 1e-3;
32,000✔
347
  site.wgt = particle_repr->weight;
32,000✔
348
  return site;
32,000✔
349
}
350

351
vector<SourceSite> mcpl_source_sites(std::string path)
32✔
352
{
353
  ensure_mcpl_ready_or_fatal();
32✔
354
  vector<SourceSite> sites;
32✔
355

356
  mcpl_file_t* mcpl_file = g_mcpl_api->open_file(path.c_str());
32✔
357
  if (!mcpl_file) {
32!
358
    fatal_error(fmt::format("MCPL: Could not open file '{}'. It might be "
×
359
                            "missing, inaccessible, or not a valid MCPL file.",
360
      path));
361
  }
362

363
  size_t n_particles_in_file = g_mcpl_api->hdr_nparticles(mcpl_file);
32✔
364
  size_t n_skipped = 0;
32✔
365
  if (n_particles_in_file > 0) {
32!
366
    sites.reserve(n_particles_in_file);
32✔
367
  }
368

369
  for (size_t i = 0; i < n_particles_in_file; ++i) {
32,032✔
370
    const mcpl_particle_repr_t* p_repr = g_mcpl_api->read(mcpl_file);
32,000✔
371
    if (!p_repr) {
32,000!
372
      warning(fmt::format("MCPL: Read error or unexpected end of file '{}' "
×
373
                          "after reading {} of {} expected particles.",
374
        path, sites.size(), n_particles_in_file));
×
375
      break;
×
376
    }
377
    if (p_repr->pdgcode == 2112 || p_repr->pdgcode == 22 ||
32,000!
378
        p_repr->pdgcode == 11 || p_repr->pdgcode == -11) {
×
379
      sites.push_back(mcpl_particle_to_site(p_repr));
32,000✔
380
    } else {
381
      n_skipped++;
×
382
    }
383
  }
384

385
  g_mcpl_api->close_file(mcpl_file);
32✔
386

387
  if (n_skipped > 0 && n_particles_in_file > 0) {
32!
388
    double percent_skipped =
×
389
      100.0 * static_cast<double>(n_skipped) / n_particles_in_file;
×
390
    warning(fmt::format(
×
391
      "MCPL: Skipped {} of {} total particles ({:.1f}%) in file '{}' because "
392
      "their type is not supported by OpenMC.",
393
      n_skipped, n_particles_in_file, percent_skipped, path));
394
  }
395

396
  if (sites.empty()) {
32!
397
    if (n_particles_in_file > 0) {
×
398
      fatal_error(fmt::format(
×
399
        "MCPL file '{}' contained {} particles, but none were of the supported "
400
        "types (neutron, photon, electron, positron). OpenMC cannot proceed "
401
        "without source particles.",
402
        path, n_particles_in_file));
403
    } else {
404
      fatal_error(fmt::format(
×
405
        "MCPL file '{}' is empty or contains no particle data.", path));
406
    }
407
  }
408
  return sites;
64✔
409
}
×
410

411
void write_mcpl_source_bank_internal(mcpl_outfile_t* file_id,
104✔
412
  span<SourceSite> local_source_bank,
413
  const vector<int64_t>& bank_index_all_ranks)
414
{
415
  if (mpi::master) {
104✔
416
    if (!file_id) {
99!
417
      fatal_error("MCPL: Internal error - master rank called "
×
418
                  "write_mcpl_source_bank_internal with null file_id.");
419
    }
420
    vector<SourceSite> receive_buffer;
99✔
421

422
    for (int rank_idx = 0; rank_idx < mpi::n_procs; ++rank_idx) {
203✔
423
      size_t num_sites_on_rank = static_cast<size_t>(
424
        bank_index_all_ranks[rank_idx + 1] - bank_index_all_ranks[rank_idx]);
104✔
425
      if (num_sites_on_rank == 0)
104✔
426
        continue;
11✔
427

428
      span<const SourceSite> sites_to_write;
93✔
429
#ifdef OPENMC_MPI
430
      if (rank_idx == mpi::rank) {
45✔
431
        sites_to_write = openmc::span<const SourceSite>(
40✔
432
          local_source_bank.data(), num_sites_on_rank);
40✔
433
      } else {
434
        if (receive_buffer.size() < num_sites_on_rank) {
5!
435
          receive_buffer.resize(num_sites_on_rank);
5✔
436
        }
437
        MPI_Recv(receive_buffer.data(), num_sites_on_rank, mpi::source_site,
5✔
438
          rank_idx, rank_idx, mpi::intracomm, MPI_STATUS_IGNORE);
439
        sites_to_write = openmc::span<const SourceSite>(
5✔
440
          receive_buffer.data(), num_sites_on_rank);
5✔
441
      }
442
#else
443
      sites_to_write = openmc::span<const SourceSite>(
48✔
444
        local_source_bank.data(), num_sites_on_rank);
48✔
445
#endif
446
      for (const auto& site : sites_to_write) {
46,414✔
447
        mcpl_particle_repr_t p_repr {};
46,321✔
448
        p_repr.position[0] = site.r.x;
46,321✔
449
        p_repr.position[1] = site.r.y;
46,321✔
450
        p_repr.position[2] = site.r.z;
46,321✔
451
        p_repr.direction[0] = site.u.x;
46,321✔
452
        p_repr.direction[1] = site.u.y;
46,321✔
453
        p_repr.direction[2] = site.u.z;
46,321✔
454
        p_repr.ekin = site.E * 1e-6;
46,321✔
455
        p_repr.time = site.time * 1e3;
46,321✔
456
        p_repr.weight = site.wgt;
46,321✔
457
        switch (site.particle) {
46,321!
458
        case ParticleType::neutron:
46,321✔
459
          p_repr.pdgcode = 2112;
46,321✔
460
          break;
46,321✔
461
        case ParticleType::photon:
×
462
          p_repr.pdgcode = 22;
×
463
          break;
×
464
        case ParticleType::electron:
×
465
          p_repr.pdgcode = 11;
×
466
          break;
×
467
        case ParticleType::positron:
×
468
          p_repr.pdgcode = -11;
×
469
          break;
×
470
        default:
×
471
          continue;
×
472
        }
473
        g_mcpl_api->add_particle(file_id, &p_repr);
46,321✔
474
      }
475
    }
476
  } else {
99✔
477
#ifdef OPENMC_MPI
478
    if (!local_source_bank.empty()) {
5!
479
      MPI_Send(local_source_bank.data(), local_source_bank.size(),
5✔
480
        mpi::source_site, 0, mpi::rank, mpi::intracomm);
481
    }
482
#endif
483
  }
484
}
104✔
485

486
void write_mcpl_source_point(const char* filename, span<SourceSite> source_bank,
104✔
487
  const vector<int64_t>& bank_index)
488
{
489
  ensure_mcpl_ready_or_fatal();
104✔
490

491
  std::string filename_(filename);
104✔
492
  const auto extension = get_file_extension(filename_);
104✔
493
  if (extension.empty()) {
104!
494
    filename_.append(".mcpl");
×
495
  } else if (extension != "mcpl") {
104!
496
    warning(fmt::format("Specified filename '{}' has an extension '.{}', but "
×
497
                        "an MCPL file (.mcpl) will be written using this name.",
498
      filename, extension));
499
  }
500

501
  mcpl_outfile_t* file_id = nullptr;
104✔
502

503
  if (mpi::master) {
104✔
504
    file_id = g_mcpl_api->create_outfile(filename_.c_str());
99✔
505
    if (!file_id) {
99!
506
      fatal_error(fmt::format(
×
507
        "MCPL: Failed to create output file '{}'. Check permissions and path.",
508
        filename_));
509
    }
510
    std::string src_line;
99✔
511
    if (VERSION_DEV) {
512
      src_line = fmt::format("OpenMC {}.{}.{}-dev{}", VERSION_MAJOR,
180✔
513
        VERSION_MINOR, VERSION_RELEASE, VERSION_COMMIT_COUNT);
99✔
514
    } else {
515
      src_line = fmt::format(
516
        "OpenMC {}.{}.{}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE);
517
    }
518
    g_mcpl_api->hdr_set_srcname(file_id, src_line.c_str());
99✔
519

520
    // Initialize stat:sum with -1 to indicate incomplete file (issue #3514)
521
    // This follows MCPL >= 2.1.0 convention for tracking simulation statistics
522
    // The -1 value indicates "not available" if file creation is interrupted
523
    if (g_mcpl_api->hdr_add_stat_sum) {
99!
524
      // Using key "openmc_np1" following tkittel's recommendation
525
      // Initial value of -1 prevents misleading values in case of crashes
526
      g_mcpl_api->hdr_add_stat_sum(file_id, "openmc_np1", -1.0);
99✔
527
    }
528
  }
99✔
529

530
  write_mcpl_source_bank_internal(file_id, source_bank, bank_index);
104✔
531

532
  if (mpi::master) {
104✔
533
    if (file_id) {
99!
534
      // Update stat:sum with actual particle count before closing (issue #3514)
535
      // This represents the original number of source particles in the
536
      // simulation (not the number of particles in the file)
537
      if (g_mcpl_api->hdr_add_stat_sum) {
99!
538
        // Calculate total source particles from active batches
539
        // Per issue #3514: this should be the original number of source
540
        // particles, not the number written to the file
541
        int64_t total_source_particles =
99✔
542
          static_cast<int64_t>(settings::n_batches - settings::n_inactive) *
99✔
543
          settings::gen_per_batch * settings::n_particles;
99✔
544
        // Update with actual count - this overwrites the initial -1 value
545
        g_mcpl_api->hdr_add_stat_sum(
99✔
546
          file_id, "openmc_np1", static_cast<double>(total_source_particles));
547
      }
548

549
      g_mcpl_api->close_outfile(file_id);
99✔
550
    }
551
  }
552
}
104✔
553

554
// Collision track feature with MCPL
555
void write_mcpl_collision_track_internal(mcpl_outfile_t* file_id,
11✔
556
  span<CollisionTrackSite> collision_track_bank,
557
  const vector<int64_t>& bank_index_all_ranks)
558
{
559
  if (mpi::master) {
11!
560
    if (!file_id) {
11!
NEW
561
      fatal_error("MCPL: Internal error - master rank called "
×
562
                  "write_mcpl_source_bank_internal with null file_id.");
563
    }
564
    vector<CollisionTrackSite> receive_buffer;
11✔
565
    vector<CollisionTrackSite> all_sites;
11✔
566
    all_sites.reserve(static_cast<size_t>(bank_index_all_ranks.back()));
11✔
567
    vector<std::string> all_blobs;
11✔
568
    all_blobs.reserve(static_cast<size_t>(bank_index_all_ranks.back()));
11✔
569

570
    for (int rank_idx = 0; rank_idx < mpi::n_procs; ++rank_idx) {
22✔
571
      size_t num_sites_on_rank = static_cast<size_t>(
572
        bank_index_all_ranks[rank_idx + 1] - bank_index_all_ranks[rank_idx]);
11✔
573
      if (num_sites_on_rank == 0)
11!
NEW
574
        continue;
×
575

576
      span<const CollisionTrackSite> sites_to_process;
11✔
577
#ifdef OPENMC_MPI
578
      if (rank_idx == mpi::rank) {
5!
579
        sites_to_process = openmc::span<const CollisionTrackSite>(
5✔
580
          collision_track_bank.data(), num_sites_on_rank);
5✔
581
      } else {
582
        receive_buffer.resize(num_sites_on_rank);
×
583
        MPI_Recv(receive_buffer.data(), num_sites_on_rank,
×
584
          mpi::collision_track_site, rank_idx, rank_idx, mpi::intracomm,
585
          MPI_STATUS_IGNORE);
586
        sites_to_process = openmc::span<const CollisionTrackSite>(
587
          receive_buffer.data(), num_sites_on_rank);
588
      }
589
#else
590
      sites_to_process = openmc::span<const CollisionTrackSite>(
6✔
591
        collision_track_bank.data(), num_sites_on_rank);
6✔
592
#endif
593

594
      for (const auto& site : sites_to_process) {
671✔
595
        std::ostringstream custom_data_stream;
660✔
596
        custom_data_stream << " dE : " << site.dE
660✔
597
                           << " ; event_mt : " << site.event_mt
660✔
598
                           << " ; delayed_group : " << site.delayed_group
660✔
599
                           << " ; cell_id : " << site.cell_id
660✔
600
                           << " ; nuclide_id : " << site.nuclide_id
660✔
601
                           << " ; material_id : " << site.material_id
660✔
602
                           << " ; universe_id : " << site.universe_id
660✔
603
                           << " ; n_collision : " << site.n_collision
660✔
604
                           << " ; parent_id : " << site.parent_id
660✔
605
                           << " ; progeny_id : " << site.progeny_id;
660✔
606

607
        all_blobs.push_back(custom_data_stream.str());
660✔
608
        all_sites.push_back(site);
660✔
609
      }
660✔
610
    }
611

612
    for (size_t idx = 0; idx < all_blobs.size(); ++idx) {
671✔
613
      const auto& blob = all_blobs[idx];
660✔
614
      std::string key = "blob_" + std::to_string(idx);
660✔
615
      g_mcpl_api->hdr_add_data(file_id, key.c_str(), blob.size(), blob.c_str());
660✔
616
    }
660✔
617

618
    for (const auto& site : all_sites) {
671✔
619
      mcpl_particle_repr_t p_repr {};
660✔
620
      p_repr.position[0] = site.r.x;
660✔
621
      p_repr.position[1] = site.r.y;
660✔
622
      p_repr.position[2] = site.r.z;
660✔
623
      p_repr.direction[0] = site.u.x;
660✔
624
      p_repr.direction[1] = site.u.y;
660✔
625
      p_repr.direction[2] = site.u.z;
660✔
626
      p_repr.ekin = site.E * 1e-6;
660✔
627
      p_repr.time = site.time * 1e3;
660✔
628
      p_repr.weight = site.wgt;
660✔
629
      switch (site.particle) {
660!
630
      case ParticleType::neutron:
660✔
631
        p_repr.pdgcode = 2112;
660✔
632
        break;
660✔
NEW
633
      case ParticleType::photon:
×
NEW
634
        p_repr.pdgcode = 22;
×
NEW
635
        break;
×
NEW
636
      case ParticleType::electron:
×
NEW
637
        p_repr.pdgcode = 11;
×
NEW
638
        break;
×
NEW
639
      case ParticleType::positron:
×
NEW
640
        p_repr.pdgcode = -11;
×
NEW
641
        break;
×
NEW
642
      default:
×
NEW
643
        continue;
×
644
      }
645
      g_mcpl_api->add_particle(file_id, &p_repr);
660✔
646
    }
647
  } else {
11✔
648
#ifdef OPENMC_MPI
649
    if (!collision_track_bank.empty()) {
×
650
      MPI_Send(collision_track_bank.data(), collision_track_bank.size(),
651
        mpi::collision_track_site, 0, mpi::rank, mpi::intracomm);
652
    }
653
#endif
654
  }
655
}
11✔
656

657
void write_mcpl_collision_track(const char* filename,
11✔
658
  span<CollisionTrackSite> collision_track_bank,
659
  const vector<int64_t>& bank_index)
660
{
661
  ensure_mcpl_ready_or_fatal();
11✔
662

663
  std::string filename_(filename);
11✔
664
  const auto extension = get_file_extension(filename_);
11✔
665
  if (extension.empty()) {
11!
NEW
666
    filename_.append(".mcpl");
×
667
  } else if (extension != "mcpl") {
11!
NEW
668
    warning(fmt::format("Specified filename '{}' has an extension '.{}', but "
×
669
                        "an MCPL file (.mcpl) will be written using this name.",
670
      filename, extension));
671
  }
672

673
  mcpl_outfile_t* file_id = nullptr;
11✔
674

675
  if (mpi::master) {
11!
676
    file_id = g_mcpl_api->create_outfile(filename_.c_str());
11✔
677
    if (!file_id) {
11!
NEW
678
      fatal_error(fmt::format(
×
679
        "MCPL: Failed to create output file '{}'. Check permissions and path.",
680
        filename_));
681
    }
682
    std::string src_line;
11✔
683
    if (VERSION_DEV) {
684
      src_line = fmt::format("OpenMC {}.{}.{}-dev{}", VERSION_MAJOR,
20✔
685
        VERSION_MINOR, VERSION_RELEASE, VERSION_COMMIT_COUNT);
11✔
686
    } else {
687
      src_line = fmt::format(
688
        "OpenMC {}.{}.{}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE);
689
    }
690

691
    g_mcpl_api->hdr_set_srcname(file_id, src_line.c_str());
11✔
692
  }
11✔
693
  write_mcpl_collision_track_internal(
11✔
694
    file_id, collision_track_bank, bank_index);
695

696
  if (mpi::master) {
11!
697
    if (file_id) {
11!
698
      g_mcpl_api->close_outfile(file_id);
11✔
699
    }
700
  }
701
}
11✔
702

703
} // 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