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

openmc-dev / openmc / 30830380864

03 Aug 2026 04:03PM UTC coverage: 81.48% (+0.05%) from 81.429%
30830380864

Pull #3971

github

web-flow
Merge 13d0297e2 into 5982acdf8
Pull Request #3971: Delta tracking

18827 of 27241 branches covered (69.11%)

Branch coverage included in aggregate %.

597 of 645 new or added lines in 20 files covered. (92.56%)

60852 of 70549 relevant lines covered (86.25%)

50385682.17 hits per line

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

83.46
/src/event.cpp
1
#include "openmc/event.h"
2

3
#include "openmc/bank.h"
4
#include "openmc/error.h"
5
#include "openmc/material.h"
6
#include "openmc/settings.h"
7
#include "openmc/simulation.h"
8
#include "openmc/timer.h"
9

10
namespace openmc {
11

12
//==============================================================================
13
// Global variables
14
//==============================================================================
15

16
namespace simulation {
17

18
SharedArray<EventQueueItem> calculate_fuel_xs_queue;
19
SharedArray<EventQueueItem> calculate_nonfuel_xs_queue;
20
SharedArray<EventQueueItem> advance_particle_queue;
21
SharedArray<EventQueueItem> surface_crossing_queue;
22
SharedArray<EventQueueItem> collision_queue;
23

24
vector<Particle> particles;
25

26
} // namespace simulation
27

28
//==============================================================================
29
// Non-member functions
30
//==============================================================================
31

32
void init_event_queues(int64_t n_particles)
315✔
33
{
34
  simulation::calculate_fuel_xs_queue.reserve(n_particles);
315✔
35
  simulation::calculate_nonfuel_xs_queue.reserve(n_particles);
315✔
36
  simulation::advance_particle_queue.reserve(n_particles);
315✔
37
  simulation::surface_crossing_queue.reserve(n_particles);
315✔
38
  simulation::collision_queue.reserve(n_particles);
315✔
39

40
  simulation::particles.resize(n_particles);
315✔
41
}
315✔
42

43
void free_event_queues(void)
×
44
{
45
  simulation::calculate_fuel_xs_queue.clear();
×
46
  simulation::calculate_nonfuel_xs_queue.clear();
×
47
  simulation::advance_particle_queue.clear();
×
48
  simulation::surface_crossing_queue.clear();
×
49
  simulation::collision_queue.clear();
×
50

51
  simulation::particles.clear();
×
52
}
×
53

54
void dispatch_xs_event(int64_t buffer_idx)
267,072,265✔
55
{
56
  Particle& p = simulation::particles[buffer_idx];
267,072,265✔
57
  if (p.material() == MATERIAL_VOID ||
267,072,265✔
58
      !model::materials[p.material()]->fissionable()) {
252,806,434✔
59
    simulation::calculate_nonfuel_xs_queue.thread_safe_append({p, buffer_idx});
78,240,137✔
60
  } else {
61
    simulation::calculate_fuel_xs_queue.thread_safe_append({p, buffer_idx});
188,832,128✔
62
  }
63
}
267,072,265✔
64

65
void process_death_events(int64_t n_particles)
4,306✔
66
{
67
  simulation::time_event_death.start();
4,306✔
68
#pragma omp parallel for schedule(runtime)
4,186✔
69
  for (int64_t i = 0; i < n_particles; i++) {
100,120✔
70
    Particle& p = simulation::particles[i];
100,000✔
71
    p.event_death();
100,000✔
72
  }
73
  simulation::time_event_death.stop();
4,306✔
74
}
4,306✔
75

76
//==============================================================================
77
// Functions for surface tracking
78
//==============================================================================
79

80
void process_init_events(int64_t n_particles, int64_t source_offset)
3,242✔
81
{
82
  simulation::time_event_init.start();
3,242✔
83
#pragma omp parallel for schedule(runtime)
3,242✔
84
  for (int64_t i = 0; i < n_particles; i++) {
×
85
    initialize_particle_track(
86
      simulation::particles[i], source_offset + i + 1, false);
87
    dispatch_xs_event(i);
88
  }
89
  simulation::time_event_init.stop();
3,242✔
90
}
3,242✔
91

92
void process_calculate_xs_events(SharedArray<EventQueueItem>& queue)
1,399,162✔
93
{
94
  simulation::time_event_calculate_xs.start();
1,399,162✔
95

96
  // TODO: If using C++17, we could perform a parallel sort of the queue by
97
  // particle type, material type, and then energy, in order to improve cache
98
  // locality and reduce thread divergence on GPU. However, the parallel
99
  // algorithms typically require linking against an additional library (Intel
100
  // TBB). Prior to C++17, std::sort is a serial only operation, which in this
101
  // case makes it too slow to be practical for most test problems.
102
  //
103
  // std::sort(std::execution::par_unseq, queue.data(), queue.data() +
104
  // queue.size());
105

106
  int64_t offset = simulation::advance_particle_queue.size();
1,399,162✔
107

108
#pragma omp parallel for schedule(runtime)
1,399,162✔
109
  for (int64_t i = 0; i < queue.size(); i++) {
×
110
    Particle* p = &simulation::particles[queue[i].idx];
111
    p->event_calculate_xs();
112

113
    // After executing a calculate_xs event, particles will
114
    // always require an advance event. Therefore, we don't need to use
115
    // the protected enqueuing function.
116
    simulation::advance_particle_queue[offset + i] = queue[i];
117
  }
118

119
  simulation::advance_particle_queue.resize(offset + queue.size());
1,399,162✔
120

121
  queue.resize(0);
1,399,162✔
122

123
  simulation::time_event_calculate_xs.stop();
1,399,162✔
124
}
1,399,162✔
125

126
void process_advance_particle_events()
1,394,191✔
127
{
128
  simulation::time_event_advance_particle.start();
1,394,191✔
129

130
#pragma omp parallel for schedule(runtime)
1,394,191✔
131
  for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {
×
132
    int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
133
    Particle& p = simulation::particles[buffer_idx];
134
    p.event_advance();
135
    if (!p.alive())
×
136
      continue;
137
    if (p.collision_distance() > p.boundary().distance()) {
×
138
      simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx});
139
    } else {
140
      simulation::collision_queue.thread_safe_append({p, buffer_idx});
141
    }
142
  }
143

144
  simulation::advance_particle_queue.resize(0);
1,394,191✔
145

146
  simulation::time_event_advance_particle.stop();
1,394,191✔
147
}
1,394,191✔
148

149
void process_surface_crossing_events()
267,591✔
150
{
151
  simulation::time_event_surface_crossing.start();
267,591✔
152

153
#pragma omp parallel for schedule(runtime)
267,591✔
154
  for (int64_t i = 0; i < simulation::surface_crossing_queue.size(); i++) {
×
155
    int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;
156
    Particle& p = simulation::particles[buffer_idx];
157
    p.event_cross_surface();
158
    p.event_check_limit_and_revive();
159
    if (p.alive())
×
160
      dispatch_xs_event(buffer_idx);
161
  }
162

163
  simulation::surface_crossing_queue.resize(0);
267,591✔
164

165
  simulation::time_event_surface_crossing.stop();
267,591✔
166
}
267,591✔
167

168
void process_collision_events()
1,146,934✔
169
{
170
  simulation::time_event_collision.start();
1,146,934✔
171

172
#pragma omp parallel for schedule(runtime)
1,146,934✔
173
  for (int64_t i = 0; i < simulation::collision_queue.size(); i++) {
×
174
    int64_t buffer_idx = simulation::collision_queue[i].idx;
175
    Particle& p = simulation::particles[buffer_idx];
176
    p.event_collide();
177
    p.event_check_limit_and_revive();
178
    if (p.alive())
×
179
      dispatch_xs_event(buffer_idx);
180
  }
181

182
  simulation::collision_queue.resize(0);
1,146,934✔
183

184
  simulation::time_event_collision.stop();
1,146,934✔
185
}
1,146,934✔
186

187
void process_transport_events()
3,926✔
188
{
189
  while (true) {
4,211,804✔
190
    int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
4,211,804✔
191
      simulation::calculate_nonfuel_xs_queue.size(),
4,211,804✔
192
      simulation::advance_particle_queue.size(),
4,211,804✔
193
      simulation::surface_crossing_queue.size(),
4,211,804✔
194
      simulation::collision_queue.size()});
4,211,804✔
195

196
    if (max == 0) {
4,211,804✔
197
      break;
198
    } else if (max == simulation::calculate_fuel_xs_queue.size()) {
4,207,878✔
199
      process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
465,856✔
200
    } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
3,742,022✔
201
      process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
933,306✔
202
    } else if (max == simulation::advance_particle_queue.size()) {
2,808,716✔
203
      process_advance_particle_events();
1,394,191✔
204
    } else if (max == simulation::surface_crossing_queue.size()) {
1,414,525✔
205
      process_surface_crossing_events();
267,591✔
206
    } else if (max == simulation::collision_queue.size()) {
1,146,934!
207
      process_collision_events();
1,146,934✔
208
    }
209
  }
210
}
3,926✔
211

212
void process_init_secondary_events(int64_t n_particles, int64_t offset,
684✔
213
  const SharedArray<SourceSite>& shared_secondary_bank)
214
{
215
  simulation::time_event_init.start();
684✔
216
#pragma omp parallel for schedule(runtime)
684✔
217
  for (int64_t i = 0; i < n_particles; i++) {
×
218
    initialize_particle_track(simulation::particles[i], offset + i + 1, true);
219
    const SourceSite& site = shared_secondary_bank[offset + i];
220
    simulation::particles[i].event_revive_from_secondary(site);
221
    if (simulation::particles[i].alive()) {
×
222
      dispatch_xs_event(i);
223
    }
224
  }
225
  simulation::time_event_init.stop();
684✔
226
}
684✔
227

228
//==============================================================================
229
// Functions for delta tracking
230
//==============================================================================
231

232
void process_delta_init_events(int64_t n_particles, int64_t source_offset)
380✔
233
{
234
  simulation::time_event_init.start();
380✔
235
#pragma omp parallel for schedule(runtime)
260✔
236
  for (int64_t i = 0; i < n_particles; i++) {
100,120✔
237
    simulation::particles[i].delta_tracking() = true;
100,000✔
238
    initialize_particle_track(
100,000✔
239
      simulation::particles[i], source_offset + i + 1, false);
100,000✔
240

241
    simulation::advance_particle_queue[i] = {simulation::particles[i], i};
100,000✔
242
  }
243
  simulation::advance_particle_queue.resize(n_particles);
380✔
244
  simulation::time_event_init.stop();
380✔
245
}
380✔
246

247
void process_delta_calculate_xs_events(SharedArray<EventQueueItem>& queue)
739,299✔
248
{
249
  simulation::time_event_calculate_xs.start();
739,299✔
250

251
  // TODO: If using C++17, we could perform a parallel sort of the queue by
252
  // particle type, material type, and then energy, in order to improve cache
253
  // locality and reduce thread divergence on GPU. However, the parallel
254
  // algorithms typically require linking against an additional library (Intel
255
  // TBB). Prior to C++17, std::sort is a serial only operation, which in this
256
  // case makes it too slow to be practical for most test problems.
257
  //
258
  // std::sort(std::execution::par_unseq, queue.data(), queue.data() +
259
  // queue.size());
260

261
  int64_t offset = simulation::collision_queue.size();
739,299✔
262

263
  simulation::collision_queue.resize(offset + queue.size());
739,299✔
264

265
#pragma omp parallel for schedule(runtime)
482,075✔
266
  for (int64_t i = 0; i < queue.size(); i++) {
12,899,654✔
267
    Particle* p = &simulation::particles[queue[i].idx];
12,642,430✔
268
    p->event_calculate_xs();
12,642,430✔
269

270
    // After executing a calculate_xs event in delta tracking, particles will
271
    // always require a collision event. Therefore, we don't need to use
272
    // the protected enqueuing function.
273
    simulation::collision_queue[offset + i] = queue[i];
12,642,430✔
274
  }
275

276
  queue.resize(0);
739,299✔
277

278
  simulation::time_event_calculate_xs.stop();
739,299✔
279
}
739,299✔
280

281
void process_delta_advance_particle_events()
737,706✔
282
{
283
  simulation::time_event_advance_particle.start();
737,706✔
284

285
#pragma omp parallel for schedule(runtime)
480,534✔
286
  for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {
14,804,312✔
287
    int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
14,547,140✔
288
    Particle& p = simulation::particles[buffer_idx];
14,547,140✔
289
    p.event_delta_advance();
14,547,140✔
290
    if (!p.alive())
14,547,140!
291
      continue;
292

293
    if (p.collision_distance() < p.boundary().distance()) {
14,547,140✔
294
      // We need to compute cross sections prior to processing a collision.
295
      dispatch_xs_event(buffer_idx);
12,642,430✔
296
    } else {
297
      simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx});
1,904,710✔
298
    }
299
  }
300

301
  simulation::advance_particle_queue.resize(0);
737,706✔
302

303
  simulation::time_event_advance_particle.stop();
737,706✔
304
}
737,706✔
305

306
void process_delta_surface_crossing_events()
51,472✔
307
{
308
  simulation::time_event_surface_crossing.start();
51,472✔
309

310
#pragma omp parallel for schedule(runtime)
33,470✔
311
  for (int64_t i = 0; i < simulation::surface_crossing_queue.size(); i++) {
1,922,712✔
312
    int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;
1,904,710✔
313
    Particle& p = simulation::particles[buffer_idx];
1,904,710✔
314
    p.event_cross_surface();
1,904,710✔
315
    p.event_check_limit_and_revive();
1,904,710✔
316
    if (p.alive())
1,904,710!
317
      simulation::advance_particle_queue.thread_safe_append({p, buffer_idx});
1,904,710✔
318
  }
319

320
  simulation::surface_crossing_queue.resize(0);
51,472✔
321

322
  simulation::time_event_surface_crossing.stop();
51,472✔
323
}
51,472✔
324

325
void process_delta_collision_events()
694,025✔
326
{
327
  simulation::time_event_collision.start();
694,025✔
328

329
#pragma omp parallel for schedule(runtime)
452,223✔
330
  for (int64_t i = 0; i < simulation::collision_queue.size(); i++) {
12,884,232✔
331
    int64_t buffer_idx = simulation::collision_queue[i].idx;
12,642,430✔
332
    Particle& p = simulation::particles[buffer_idx];
12,642,430✔
333

334
    if (p.kill_invalid_maj()) {
12,642,430!
335
      continue;
336
    }
337

338
    if (prn(p.current_seed()) < (p.macro_xs().total / p.majorant())) {
12,642,430✔
339
      // Real collision, need to process the collision prior to enqueuing an
340
      // advance event.
341
      p.event_collide();
4,556,135✔
342
    }
343

344
    p.event_check_limit_and_revive();
12,642,430✔
345
    if (p.alive()) {
12,642,430✔
346
      simulation::advance_particle_queue.thread_safe_append({p, buffer_idx});
12,542,430✔
347
    }
348
  }
349

350
  simulation::collision_queue.resize(0);
694,025✔
351

352
  simulation::time_event_collision.stop();
694,025✔
353
}
694,025✔
354

355
void process_delta_transport_events()
380✔
356
{
357
  while (true) {
2,222,882✔
358
    int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
2,222,882✔
359
      simulation::calculate_nonfuel_xs_queue.size(),
2,222,882✔
360
      simulation::advance_particle_queue.size(),
2,222,882✔
361
      simulation::surface_crossing_queue.size(),
2,222,882✔
362
      simulation::collision_queue.size()});
2,222,882✔
363

364
    if (max == 0) {
2,222,882✔
365
      break;
366
    } else if (max == simulation::calculate_fuel_xs_queue.size()) {
2,222,502✔
367
      process_delta_calculate_xs_events(simulation::calculate_fuel_xs_queue);
122,646✔
368
    } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
2,099,856✔
369
      process_delta_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
616,653✔
370
    } else if (max == simulation::advance_particle_queue.size()) {
1,483,203✔
371
      process_delta_advance_particle_events();
737,706✔
372
    } else if (max == simulation::surface_crossing_queue.size()) {
745,497✔
373
      process_delta_surface_crossing_events();
51,472✔
374
    } else if (max == simulation::collision_queue.size()) {
694,025!
375
      process_delta_collision_events();
694,025✔
376
    }
377
  }
378
}
380✔
379

NEW
380
void process_delta_init_secondary_events(int64_t n_particles, int64_t offset,
×
381
  const SharedArray<SourceSite>& shared_secondary_bank)
382
{
NEW
383
  simulation::time_event_init.start();
×
384
#pragma omp parallel for schedule(runtime)
385
  for (int64_t i = 0; i < n_particles; i++) {
×
386
    simulation::particles[i].delta_tracking() = true;
387
    initialize_particle_track(simulation::particles[i], offset + i + 1, true);
388
    const SourceSite& site = shared_secondary_bank[offset + i];
389
    simulation::particles[i].event_revive_from_secondary(site);
390

391
    if (simulation::particles[i].alive()) {
×
392
      simulation::particles[i].event_calculate_xs();
393
      simulation::advance_particle_queue.thread_safe_append(
394
        {simulation::particles[i], i});
395
    }
396
  }
NEW
397
  simulation::time_event_init.stop();
×
NEW
398
}
×
399

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

© 2026 Coveralls, Inc