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

openmc-dev / openmc / 20278515727

16 Dec 2025 06:26PM UTC coverage: 81.822% (-0.1%) from 81.92%
20278515727

Pull #3493

github

web-flow
Merge 7f36869a3 into bbfa18d72
Pull Request #3493: Implement vector fitting to replace external `vectfit` package

17013 of 23575 branches covered (72.17%)

Branch coverage included in aggregate %.

188 of 207 new or added lines in 3 files covered. (90.82%)

3101 existing lines in 56 files now uncovered.

54973 of 64404 relevant lines covered (85.36%)

41385524.25 hits per line

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

57.9
/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, uint32_t datalength, 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

UNCOV
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
UNCOV
104
  const char* err = dlerror();
×
UNCOV
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)
82✔
122
  {
82✔
123
    if (!lib_handle)
82!
UNCOV
124
      throw std::runtime_error(
×
UNCOV
125
        "MCPL library handle is null during API binding.");
×
126

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

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

158
    // Try to load mcpl_hdr_add_data (available in MCPL >= 2.1.0)
159
    // Set to nullptr if not available for graceful fallback
160
    try {
161
      hdr_add_data = reinterpret_cast<mcpl_hdr_add_data_fpt>(
82✔
162
        load_symbol_platform("mcpl_hdr_add_data"));
82✔
UNCOV
163
    } catch (const std::runtime_error&) {
×
UNCOV
164
      hdr_add_data = nullptr;
×
UNCOV
165
    }
×
166

167
    // Try to load mcpl_hdr_add_stat_sum (available in MCPL >= 2.1.0)
168
    // Set to nullptr if not available for graceful fallback
169
    try {
170
      hdr_add_stat_sum = reinterpret_cast<mcpl_hdr_add_stat_sum_fpt>(
82✔
171
        load_symbol_platform("mcpl_hdr_add_stat_sum"));
82✔
172
    } catch (const std::runtime_error&) {
×
UNCOV
173
      hdr_add_stat_sum = nullptr;
×
174
    }
×
175
  }
82✔
176
};
177

178
static LibraryHandleType g_mcpl_lib_handle = nullptr;
179
static std::unique_ptr<McplApi> g_mcpl_api;
180
static bool g_mcpl_init_attempted = false;
181
static bool g_mcpl_successfully_loaded = false;
182
static std::string g_mcpl_load_error_msg;
183
static std::once_flag g_mcpl_init_flag;
184

UNCOV
185
void append_error(std::string& existing_msg, const std::string& new_error)
×
186
{
UNCOV
187
  if (!existing_msg.empty()) {
×
UNCOV
188
    existing_msg += "; ";
×
189
  }
UNCOV
190
  existing_msg += new_error;
×
UNCOV
191
}
×
192

193
void initialize_mcpl_interface_impl()
82✔
194
{
195
  g_mcpl_init_attempted = true;
82✔
196
  g_mcpl_load_error_msg.clear();
82✔
197

198
  // Try mcpl-config
199
  if (!g_mcpl_lib_handle) {
82!
200
    FILE* pipe = nullptr;
82✔
201
#ifdef _WIN32
202
    pipe = _popen("mcpl-config --show libpath", "r");
203
#else
204
    pipe = popen("mcpl-config --show libpath 2>/dev/null", "r");
82✔
205
#endif
206
    if (pipe) {
82!
207
      char buffer[512];
208
      if (fgets(buffer, sizeof(buffer), pipe) != nullptr) {
82!
209
        std::string shlibpath = buffer;
82✔
210
        // Remove trailing whitespace
211
        while (!shlibpath.empty() &&
328!
212
               std::isspace(static_cast<unsigned char>(shlibpath.back()))) {
164✔
213
          shlibpath.pop_back();
82✔
214
        }
215

216
        if (!shlibpath.empty()) {
82!
217
#ifdef _WIN32
218
          g_mcpl_lib_handle = LoadLibraryA(shlibpath.c_str());
219
#else
220
          g_mcpl_lib_handle = dlopen(shlibpath.c_str(), RTLD_LAZY);
82✔
221
#endif
222
          if (!g_mcpl_lib_handle) {
82!
UNCOV
223
            append_error(
×
UNCOV
224
              g_mcpl_load_error_msg, fmt::format("From mcpl-config ({}): {}",
×
UNCOV
225
                                       shlibpath, get_last_library_error()));
×
226
          }
227
        }
228
      }
82✔
229
#ifdef _WIN32
230
      _pclose(pipe);
231
#else
232
      pclose(pipe);
82✔
233
#endif
234
    } else { // pipe failed to open
UNCOV
235
      append_error(g_mcpl_load_error_msg,
×
236
        "mcpl-config command not found or failed to execute");
237
    }
238
  }
239

240
  // Try standard library names
241
  if (!g_mcpl_lib_handle) {
82!
242
#ifdef _WIN32
243
    const char* standard_names[] = {"mcpl.dll", "libmcpl.dll"};
244
#else
245
    const char* standard_names[] = {"libmcpl.so", "libmcpl.dylib"};
×
246
#endif
UNCOV
247
    for (const char* name : standard_names) {
×
248
#ifdef _WIN32
249
      g_mcpl_lib_handle = LoadLibraryA(name);
250
#else
251
      g_mcpl_lib_handle = dlopen(name, RTLD_LAZY);
×
252
#endif
UNCOV
253
      if (g_mcpl_lib_handle)
×
254
        break;
×
255
    }
UNCOV
256
    if (!g_mcpl_lib_handle) {
×
UNCOV
257
      append_error(
×
258
        g_mcpl_load_error_msg, fmt::format("Using standard names (e.g. {}): {}",
×
259
                                 standard_names[0], get_last_library_error()));
×
260
    }
261
  }
262

263
  if (!g_mcpl_lib_handle) {
82!
UNCOV
264
    if (mpi::master) {
×
UNCOV
265
      warning(fmt::format("MCPL library could not be loaded. MCPL-dependent "
×
266
                          "features will be unavailable. Load attempts: {}",
267
        g_mcpl_load_error_msg.empty()
×
268
          ? "No specific error during load attempts."
×
269
          : g_mcpl_load_error_msg));
270
    }
271
    g_mcpl_successfully_loaded = false;
×
272
    return;
×
273
  }
274

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

297
void initialize_mcpl_interface_if_needed()
162✔
298
{
299
  std::call_once(g_mcpl_init_flag, initialize_mcpl_interface_impl);
162✔
300
}
162✔
301

302
bool is_mcpl_interface_available()
30✔
303
{
304
  initialize_mcpl_interface_if_needed();
30✔
305
  return g_mcpl_successfully_loaded;
30✔
306
}
307

308
inline void ensure_mcpl_ready_or_fatal()
132✔
309
{
310
  initialize_mcpl_interface_if_needed();
132✔
311
  if (!g_mcpl_successfully_loaded) {
132!
UNCOV
312
    fatal_error("MCPL functionality is required, but the MCPL library is not "
×
313
                "available or failed to initialize. Please ensure MCPL is "
314
                "installed and its library can be found (e.g., via PATH on "
315
                "Windows, LD_LIBRARY_PATH on Linux, or DYLD_LIBRARY_PATH on "
316
                "macOS). You can often install MCPL with 'pip install mcpl' or "
317
                "'conda install mcpl'.");
318
  }
319
}
132✔
320

321
SourceSite mcpl_particle_to_site(const mcpl_particle_repr_t* particle_repr)
28,000✔
322
{
323
  SourceSite site;
28,000✔
324
  switch (particle_repr->pdgcode) {
28,000!
325
  case 2112:
28,000✔
326
    site.particle = ParticleType::neutron;
28,000✔
327
    break;
28,000✔
UNCOV
328
  case 22:
×
UNCOV
329
    site.particle = ParticleType::photon;
×
UNCOV
330
    break;
×
UNCOV
331
  case 11:
×
UNCOV
332
    site.particle = ParticleType::electron;
×
UNCOV
333
    break;
×
UNCOV
334
  case -11:
×
UNCOV
335
    site.particle = ParticleType::positron;
×
UNCOV
336
    break;
×
UNCOV
337
  default:
×
UNCOV
338
    fatal_error(fmt::format(
×
339
      "MCPL: Encountered unexpected PDG code {} when converting to SourceSite.",
UNCOV
340
      particle_repr->pdgcode));
×
341
    break;
342
  }
343

344
  // Copy position and direction
345
  site.r.x = particle_repr->position[0];
28,000✔
346
  site.r.y = particle_repr->position[1];
28,000✔
347
  site.r.z = particle_repr->position[2];
28,000✔
348
  site.u.x = particle_repr->direction[0];
28,000✔
349
  site.u.y = particle_repr->direction[1];
28,000✔
350
  site.u.z = particle_repr->direction[2];
28,000✔
351
  // MCPL stores kinetic energy in [MeV], time in [ms]
352
  site.E = particle_repr->ekin * 1e6;
28,000✔
353
  site.time = particle_repr->time * 1e-3;
28,000✔
354
  site.wgt = particle_repr->weight;
28,000✔
355
  return site;
28,000✔
356
}
357

358
vector<SourceSite> mcpl_source_sites(std::string path)
28✔
359
{
360
  ensure_mcpl_ready_or_fatal();
28✔
361
  vector<SourceSite> sites;
28✔
362

363
  mcpl_file_t* mcpl_file = g_mcpl_api->open_file(path.c_str());
28✔
364
  if (!mcpl_file) {
28!
UNCOV
365
    fatal_error(fmt::format("MCPL: Could not open file '{}'. It might be "
×
366
                            "missing, inaccessible, or not a valid MCPL file.",
367
      path));
368
  }
369

370
  size_t n_particles_in_file = g_mcpl_api->hdr_nparticles(mcpl_file);
28✔
371
  size_t n_skipped = 0;
28✔
372
  if (n_particles_in_file > 0) {
28!
373
    sites.reserve(n_particles_in_file);
28✔
374
  }
375

376
  for (size_t i = 0; i < n_particles_in_file; ++i) {
28,028✔
377
    const mcpl_particle_repr_t* p_repr = g_mcpl_api->read(mcpl_file);
28,000✔
378
    if (!p_repr) {
28,000!
UNCOV
379
      warning(fmt::format("MCPL: Read error or unexpected end of file '{}' "
×
380
                          "after reading {} of {} expected particles.",
UNCOV
381
        path, sites.size(), n_particles_in_file));
×
382
      break;
×
383
    }
384
    if (p_repr->pdgcode == 2112 || p_repr->pdgcode == 22 ||
28,000!
UNCOV
385
        p_repr->pdgcode == 11 || p_repr->pdgcode == -11) {
×
386
      sites.push_back(mcpl_particle_to_site(p_repr));
28,000✔
387
    } else {
UNCOV
388
      n_skipped++;
×
389
    }
390
  }
391

392
  g_mcpl_api->close_file(mcpl_file);
28✔
393

394
  if (n_skipped > 0 && n_particles_in_file > 0) {
28!
UNCOV
395
    double percent_skipped =
×
UNCOV
396
      100.0 * static_cast<double>(n_skipped) / n_particles_in_file;
×
UNCOV
397
    warning(fmt::format(
×
398
      "MCPL: Skipped {} of {} total particles ({:.1f}%) in file '{}' because "
399
      "their type is not supported by OpenMC.",
400
      n_skipped, n_particles_in_file, percent_skipped, path));
401
  }
402

403
  if (sites.empty()) {
28!
UNCOV
404
    if (n_particles_in_file > 0) {
×
UNCOV
405
      fatal_error(fmt::format(
×
406
        "MCPL file '{}' contained {} particles, but none were of the supported "
407
        "types (neutron, photon, electron, positron). OpenMC cannot proceed "
408
        "without source particles.",
409
        path, n_particles_in_file));
410
    } else {
411
      fatal_error(fmt::format(
×
412
        "MCPL file '{}' is empty or contains no particle data.", path));
413
    }
414
  }
415
  return sites;
56✔
UNCOV
416
}
×
417

418
void write_mcpl_source_bank_internal(mcpl_outfile_t* file_id,
94✔
419
  span<SourceSite> local_source_bank,
420
  const vector<int64_t>& bank_index_all_ranks)
421
{
422
  if (mpi::master) {
94✔
423
    if (!file_id) {
90!
UNCOV
424
      fatal_error("MCPL: Internal error - master rank called "
×
425
                  "write_mcpl_source_bank_internal with null file_id.");
426
    }
427
    vector<SourceSite> receive_buffer;
90✔
428

429
    for (int rank_idx = 0; rank_idx < mpi::n_procs; ++rank_idx) {
184✔
430
      size_t num_sites_on_rank = static_cast<size_t>(
431
        bank_index_all_ranks[rank_idx + 1] - bank_index_all_ranks[rank_idx]);
94✔
432
      if (num_sites_on_rank == 0)
94✔
433
        continue;
10✔
434

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

493
void write_mcpl_source_point(const char* filename, span<SourceSite> source_bank,
94✔
494
  const vector<int64_t>& bank_index)
495
{
496
  ensure_mcpl_ready_or_fatal();
94✔
497

498
  std::string filename_(filename);
94✔
499
  const auto extension = get_file_extension(filename_);
94✔
500
  if (extension.empty()) {
94!
UNCOV
501
    filename_.append(".mcpl");
×
502
  } else if (extension != "mcpl") {
94!
UNCOV
503
    warning(fmt::format("Specified filename '{}' has an extension '.{}', but "
×
504
                        "an MCPL file (.mcpl) will be written using this name.",
505
      filename, extension));
506
  }
507

508
  mcpl_outfile_t* file_id = nullptr;
94✔
509

510
  if (mpi::master) {
94✔
511
    file_id = g_mcpl_api->create_outfile(filename_.c_str());
90✔
512
    if (!file_id) {
90!
UNCOV
513
      fatal_error(fmt::format(
×
514
        "MCPL: Failed to create output file '{}'. Check permissions and path.",
515
        filename_));
516
    }
517
    std::string src_line;
90✔
518
    if (VERSION_DEV) {
519
      src_line = fmt::format("OpenMC {}.{}.{}-dev{}", VERSION_MAJOR,
162✔
520
        VERSION_MINOR, VERSION_RELEASE, VERSION_COMMIT_COUNT);
90✔
521
    } else {
522
      src_line = fmt::format(
523
        "OpenMC {}.{}.{}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE);
524
    }
525
    g_mcpl_api->hdr_set_srcname(file_id, src_line.c_str());
90✔
526

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

537
  write_mcpl_source_bank_internal(file_id, source_bank, bank_index);
94✔
538

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

556
      g_mcpl_api->close_outfile(file_id);
90✔
557
    }
558
  }
559
}
94✔
560

561
// Collision track feature with MCPL
562
void write_mcpl_collision_track_internal(mcpl_outfile_t* file_id,
10✔
563
  span<CollisionTrackSite> collision_track_bank,
564
  const vector<int64_t>& bank_index_all_ranks)
565
{
566
  if (mpi::master) {
10!
567
    if (!file_id) {
10!
UNCOV
568
      fatal_error("MCPL: Internal error - master rank called "
×
569
                  "write_mcpl_source_bank_internal with null file_id.");
570
    }
571
    vector<CollisionTrackSite> receive_buffer;
10✔
572
    vector<CollisionTrackSite> all_sites;
10✔
573
    all_sites.reserve(static_cast<size_t>(bank_index_all_ranks.back()));
10✔
574
    vector<std::string> all_blobs;
10✔
575
    all_blobs.reserve(static_cast<size_t>(bank_index_all_ranks.back()));
10✔
576

577
    for (int rank_idx = 0; rank_idx < mpi::n_procs; ++rank_idx) {
20✔
578
      size_t num_sites_on_rank = static_cast<size_t>(
579
        bank_index_all_ranks[rank_idx + 1] - bank_index_all_ranks[rank_idx]);
10✔
580
      if (num_sites_on_rank == 0)
10!
UNCOV
581
        continue;
×
582

583
      span<const CollisionTrackSite> sites_to_process;
10✔
584
#ifdef OPENMC_MPI
585
      if (rank_idx == mpi::rank) {
4!
586
        sites_to_process = openmc::span<const CollisionTrackSite>(
4✔
587
          collision_track_bank.data(), num_sites_on_rank);
4✔
588
      } else {
589
        receive_buffer.resize(num_sites_on_rank);
×
590
        MPI_Recv(receive_buffer.data(), num_sites_on_rank,
×
591
          mpi::collision_track_site, rank_idx, rank_idx, mpi::intracomm,
592
          MPI_STATUS_IGNORE);
593
        sites_to_process = openmc::span<const CollisionTrackSite>(
594
          receive_buffer.data(), num_sites_on_rank);
595
      }
596
#else
597
      sites_to_process = openmc::span<const CollisionTrackSite>(
6✔
598
        collision_track_bank.data(), num_sites_on_rank);
6✔
599
#endif
600

601
      for (const auto& site : sites_to_process) {
610✔
602
        std::ostringstream custom_data_stream;
600✔
603
        custom_data_stream << " dE : " << site.dE
600✔
604
                           << " ; event_mt : " << site.event_mt
600✔
605
                           << " ; delayed_group : " << site.delayed_group
600✔
606
                           << " ; cell_id : " << site.cell_id
600✔
607
                           << " ; nuclide_id : " << site.nuclide_id
600✔
608
                           << " ; material_id : " << site.material_id
600✔
609
                           << " ; universe_id : " << site.universe_id
600✔
610
                           << " ; n_collision : " << site.n_collision
600✔
611
                           << " ; parent_id : " << site.parent_id
600✔
612
                           << " ; progeny_id : " << site.progeny_id;
600✔
613

614
        all_blobs.push_back(custom_data_stream.str());
600✔
615
        all_sites.push_back(site);
600✔
616
      }
600✔
617
    }
618

619
    for (size_t idx = 0; idx < all_blobs.size(); ++idx) {
610✔
620
      const auto& blob = all_blobs[idx];
600✔
621
      std::string key = "blob_" + std::to_string(idx);
600✔
622
      g_mcpl_api->hdr_add_data(file_id, key.c_str(), blob.size(), blob.c_str());
600✔
623
    }
600✔
624

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

664
void write_mcpl_collision_track(const char* filename,
10✔
665
  span<CollisionTrackSite> collision_track_bank,
666
  const vector<int64_t>& bank_index)
667
{
668
  ensure_mcpl_ready_or_fatal();
10✔
669

670
  std::string filename_(filename);
10✔
671
  const auto extension = get_file_extension(filename_);
10✔
672
  if (extension.empty()) {
10!
UNCOV
673
    filename_.append(".mcpl");
×
674
  } else if (extension != "mcpl") {
10!
UNCOV
675
    warning(fmt::format("Specified filename '{}' has an extension '.{}', but "
×
676
                        "an MCPL file (.mcpl) will be written using this name.",
677
      filename, extension));
678
  }
679

680
  mcpl_outfile_t* file_id = nullptr;
10✔
681

682
  if (mpi::master) {
10!
683
    file_id = g_mcpl_api->create_outfile(filename_.c_str());
10✔
684
    if (!file_id) {
10!
UNCOV
685
      fatal_error(fmt::format(
×
686
        "MCPL: Failed to create output file '{}'. Check permissions and path.",
687
        filename_));
688
    }
689
    std::string src_line;
10✔
690
    if (VERSION_DEV) {
691
      src_line = fmt::format("OpenMC {}.{}.{}-dev{}", VERSION_MAJOR,
18✔
692
        VERSION_MINOR, VERSION_RELEASE, VERSION_COMMIT_COUNT);
10✔
693
    } else {
694
      src_line = fmt::format(
695
        "OpenMC {}.{}.{}", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE);
696
    }
697

698
    g_mcpl_api->hdr_set_srcname(file_id, src_line.c_str());
10✔
699
  }
10✔
700
  write_mcpl_collision_track_internal(
10✔
701
    file_id, collision_track_bank, bank_index);
702

703
  if (mpi::master) {
10!
704
    if (file_id) {
10!
705
      g_mcpl_api->close_outfile(file_id);
10✔
706
    }
707
  }
708
}
10✔
709

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