30942440215

Committed 04 Aug 2026 07:15PM UTC coverage: 81.481% (+0.06%) from 81.425%

Build # 30942440215

Build Type

Pull #4044

github

Committed by

web-flow

Commit Message

Merge 5f06a886a into 8202ef6fb

Pull Request Pull Request #4044: Hybrid delta tracking

Coverage Stats

18963 of 27411 branches covered (69.18%)

Branch coverage included in aggregate %.

721 of 782 new or added lines in 20 files covered. (92.2%)

3 existing lines in 1 file now uncovered.

60943 of 70656 relevant lines covered (86.25%)

51279318.47 hits per line

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

84.42

/src/event.cpp

#include "openmc/event.h"

#include "openmc/bank.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/timer.h"

namespace openmc {

//==============================================================================
// Global variables
//==============================================================================

namespace simulation {

SharedArray<EventQueueItem> calculate_fuel_xs_queue;
SharedArray<EventQueueItem> calculate_nonfuel_xs_queue;
SharedArray<EventQueueItem> calculate_majorant_xs_queue;
SharedArray<EventQueueItem> advance_particle_queue;
SharedArray<EventQueueItem> surface_crossing_queue;
SharedArray<EventQueueItem> collision_queue;

vector<Particle> particles;

} // namespace simulation

//==============================================================================
// Non-member functions
//==============================================================================

void init_event_queues(int64_t n_particles)
{
  simulation::calculate_fuel_xs_queue.reserve(n_particles);
  simulation::calculate_nonfuel_xs_queue.reserve(n_particles);
  if (settings::delta_tracking) {
    simulation::calculate_majorant_xs_queue.reserve(n_particles);
  }
  simulation::advance_particle_queue.reserve(n_particles);
  simulation::surface_crossing_queue.reserve(n_particles);
  simulation::collision_queue.reserve(n_particles);

  simulation::particles.resize(n_particles);
}

void free_event_queues(void)
{
  simulation::calculate_fuel_xs_queue.clear();
  simulation::calculate_nonfuel_xs_queue.clear();
  if (settings::delta_tracking) {
    simulation::calculate_majorant_xs_queue.clear();
  }
  simulation::advance_particle_queue.clear();
  simulation::surface_crossing_queue.clear();
  simulation::collision_queue.clear();

  simulation::particles.clear();
}

void dispatch_xs_event(int64_t buffer_idx)
{
  Particle& p = simulation::particles[buffer_idx];
  if (p.material() == MATERIAL_VOID ||
      !model::materials[p.material()]->fissionable()) {
    simulation::calculate_nonfuel_xs_queue.thread_safe_append({p, buffer_idx});
  } else {
    simulation::calculate_fuel_xs_queue.thread_safe_append({p, buffer_idx});
  }
}

void process_death_events(int64_t n_particles)
{
  simulation::time_event_death.start();
#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < n_particles; i++) {
    Particle& p = simulation::particles[i];
    p.event_death();
  }
  simulation::time_event_death.stop();
}

//==============================================================================
// Functions for surface tracking
//==============================================================================

void process_init_events(int64_t n_particles, int64_t source_offset)
{
  simulation::time_event_init.start();
#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < n_particles; i++) {
    initialize_particle_track(
      simulation::particles[i], source_offset + i + 1, false);
    dispatch_xs_event(i);
  }
  simulation::time_event_init.stop();
}

void process_calculate_xs_events(SharedArray<EventQueueItem>& queue)
{
  simulation::time_event_calculate_xs.start();

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

  int64_t offset = simulation::advance_particle_queue.size();

#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < queue.size(); i++) {
    Particle* p = &simulation::particles[queue[i].idx];
    p->event_calculate_xs();

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

  simulation::advance_particle_queue.resize(offset + queue.size());

  queue.resize(0);

  simulation::time_event_calculate_xs.stop();
}

void process_advance_particle_events()
{
  simulation::time_event_advance_particle.start();

#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {
    int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
    Particle& p = simulation::particles[buffer_idx];
    p.event_advance();
    if (!p.alive())
      continue;
    if (p.collision_distance() > p.boundary().distance()) {
      simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx});
    } else {
      simulation::collision_queue.thread_safe_append({p, buffer_idx});
    }
  }

  simulation::advance_particle_queue.resize(0);

  simulation::time_event_advance_particle.stop();
}

void process_surface_crossing_events()
{
  simulation::time_event_surface_crossing.start();

#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < simulation::surface_crossing_queue.size(); i++) {
    int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;
    Particle& p = simulation::particles[buffer_idx];
    p.event_cross_surface();
    p.event_check_limit_and_revive();
    if (p.alive())
      dispatch_xs_event(buffer_idx);
  }

  simulation::surface_crossing_queue.resize(0);

  simulation::time_event_surface_crossing.stop();
}

void process_collision_events()
{
  simulation::time_event_collision.start();

#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < simulation::collision_queue.size(); i++) {
    int64_t buffer_idx = simulation::collision_queue[i].idx;
    Particle& p = simulation::particles[buffer_idx];
    p.event_collide();
    p.event_check_limit_and_revive();
    if (p.alive())
      dispatch_xs_event(buffer_idx);
  }

  simulation::collision_queue.resize(0);

  simulation::time_event_collision.stop();
}

void process_transport_events()
{
  while (true) {
    int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
      simulation::calculate_nonfuel_xs_queue.size(),
      simulation::advance_particle_queue.size(),
      simulation::surface_crossing_queue.size(),
      simulation::collision_queue.size()});

    if (max == 0) {
      break;
    } else if (max == simulation::calculate_fuel_xs_queue.size()) {
      process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
    } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
      process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
    } else if (max == simulation::advance_particle_queue.size()) {
      process_advance_particle_events();
    } else if (max == simulation::surface_crossing_queue.size()) {
      process_surface_crossing_events();
    } else if (max == simulation::collision_queue.size()) {
      process_collision_events();
    }
  }
}

void process_init_secondary_events(int64_t n_particles, int64_t offset,
  const SharedArray<SourceSite>& shared_secondary_bank)
{
  simulation::time_event_init.start();
#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < n_particles; i++) {
    initialize_particle_track(simulation::particles[i], offset + i + 1, true);
    const SourceSite& site = shared_secondary_bank[offset + i];
    simulation::particles[i].event_revive_from_secondary(site);
    if (simulation::particles[i].alive()) {
      dispatch_xs_event(i);
    }
  }
  simulation::time_event_init.stop();
}

//==============================================================================
// Functions for hybrid delta tracking
//==============================================================================

void process_hybrid_init_events(int64_t n_particles, int64_t source_offset)
{
  simulation::time_event_init.start();
#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < n_particles; i++) {
    initialize_particle_track(
      simulation::particles[i], source_offset + i + 1, false);

    simulation::calculate_majorant_xs_queue[i] = {simulation::particles[i], i};
  }
  simulation::calculate_majorant_xs_queue.resize(n_particles);
  simulation::time_event_init.stop();
}

void process_hybrid_calculate_xs_events(SharedArray<EventQueueItem>& queue)
{
  simulation::time_event_calculate_xs.start();

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

#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < queue.size(); i++) {
    int64_t buffer_idx = queue[i].idx;
    Particle& p = simulation::particles[buffer_idx];
    p.event_calculate_xs();
    if (!p.alive()) {
      continue;
    }

    // Particles always require a delta collision event after cross sections
    // are computed if running delta tracking. Otherwise, an advance event is
    // required for surface tracking.
    if (p.delta_tracking()) {
      simulation::collision_queue.thread_safe_append({p, buffer_idx});
    } else {
      simulation::advance_particle_queue.thread_safe_append({p, buffer_idx});
    }
  }

  queue.resize(0);

  simulation::time_event_calculate_xs.stop();
}

void process_hybrid_calculate_majorant_events()
{
  simulation::time_event_calculate_majorant_xs.start();

#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < simulation::calculate_majorant_xs_queue.size(); i++) {
    int64_t buffer_idx = simulation::calculate_majorant_xs_queue[i].idx;
    Particle& p = simulation::particles[buffer_idx];

    p.event_update_majorant();

    // Need to compute cross sections before advancing when running surface
    // tracking. Otherwise, we can go right to advancing with delta tracking.
    if (!p.delta_tracking()) {
      dispatch_xs_event(buffer_idx);
    } else {
      simulation::advance_particle_queue.thread_safe_append({p, buffer_idx});
    }
  }

  simulation::calculate_majorant_xs_queue.resize(0);

  simulation::time_event_calculate_majorant_xs.stop();
}

void process_hybrid_advance_particle_events()
{
  simulation::time_event_advance_particle.start();

#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {
    int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
    Particle& p = simulation::particles[buffer_idx];
    p.event_delta_advance();
    if (!p.alive()) {
      continue;
    }

    if (p.collision_distance() < p.boundary().distance()) {
      // We need to compute cross sections prior to processing a collision when
      // running delta tracking. Otherwise, we can go right to a collision when
      // running surface tracking.
      if (p.delta_tracking()) {
        dispatch_xs_event(buffer_idx);
      } else {
        simulation::collision_queue.thread_safe_append({p, buffer_idx});
      }
    } else {
      simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx});
    }
  }

  simulation::advance_particle_queue.resize(0);

  simulation::time_event_advance_particle.stop();
}

void process_hybrid_surface_crossing_events()
{
  simulation::time_event_surface_crossing.start();

#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < simulation::surface_crossing_queue.size(); i++) {
    int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;
    Particle& p = simulation::particles[buffer_idx];
    p.event_cross_surface();
    p.event_check_limit_and_revive();
    if (p.alive())
      simulation::calculate_majorant_xs_queue.thread_safe_append(
        {p, buffer_idx});
  }

  simulation::surface_crossing_queue.resize(0);

  simulation::time_event_surface_crossing.stop();
}

void process_hybrid_collision_events()
{
  simulation::time_event_collision.start();

#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < simulation::collision_queue.size(); i++) {
    int64_t buffer_idx = simulation::collision_queue[i].idx;
    Particle& p = simulation::particles[buffer_idx];

    if (p.delta_tracking()) {
      // Check to ensure the majorant is valid.
      if (p.kill_invalid_maj()) {
        continue;
      }

      // Perform rejection sampling. If this is true, a real collision is
      // processed prior to going back to majorant calculations.
      if (prn(p.current_seed()) < (p.macro_xs().total / p.majorant())) {
        p.event_collide();
      }
    } else {
      p.event_collide();
    }

    // Update the tracking type based on the chosen hybrid scheme.
    p.update_tracking_type();

    p.event_check_limit_and_revive();
    if (p.alive()) {
      simulation::calculate_majorant_xs_queue.thread_safe_append(
        {p, buffer_idx});
    }
  }

  simulation::collision_queue.resize(0);

  simulation::time_event_collision.stop();
}

void process_hybrid_transport_events()
{
  while (true) {
    int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
      simulation::calculate_nonfuel_xs_queue.size(),
      simulation::calculate_majorant_xs_queue.size(),
      simulation::advance_particle_queue.size(),
      simulation::surface_crossing_queue.size(),
      simulation::collision_queue.size()});

    if (max == 0) {
      break;
    } else if (max == simulation::calculate_fuel_xs_queue.size()) {
      process_hybrid_calculate_xs_events(simulation::calculate_fuel_xs_queue);
    } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
      process_hybrid_calculate_xs_events(
        simulation::calculate_nonfuel_xs_queue);
    } else if (max == simulation::calculate_majorant_xs_queue.size()) {
      process_hybrid_calculate_majorant_events();
    } else if (max == simulation::advance_particle_queue.size()) {
      process_hybrid_advance_particle_events();
    } else if (max == simulation::surface_crossing_queue.size()) {
      process_hybrid_surface_crossing_events();
    } else if (max == simulation::collision_queue.size()) {
      process_hybrid_collision_events();
    }
  }
}

void process_hybrid_init_secondary_events(int64_t n_particles, int64_t offset,
  const SharedArray<SourceSite>& shared_secondary_bank)
{
  simulation::time_event_init.start();
#pragma omp parallel for schedule(runtime)
  for (int64_t i = 0; i < n_particles; i++) {
    Particle& p = simulation::particles[i];
    initialize_particle_track(p, offset + i + 1, true);
    const SourceSite& site = shared_secondary_bank[offset + i];
    p.event_revive_from_secondary(site);

    if (p.alive()) {
      simulation::calculate_majorant_xs_queue.thread_safe_append({p, i});
    }
  }
  simulation::time_event_init.stop();
}

} // namespace openmc

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> calculate_majorant_xs_queue;
21	SharedArray<EventQueueItem> advance_particle_queue;
22	SharedArray<EventQueueItem> surface_crossing_queue;
23	SharedArray<EventQueueItem> collision_queue;
24
25	vector<Particle> particles;
26
27	} // namespace simulation
28
29	//==============================================================================
30	// Non-member functions
31	//==============================================================================
32
33	void init_event_queues(int64_t n_particles)	375✔
34	{
35	simulation::calculate_fuel_xs_queue.reserve(n_particles);	375✔
36	simulation::calculate_nonfuel_xs_queue.reserve(n_particles);	375✔
37	if (settings::delta_tracking) {	375✔
38	simulation::calculate_majorant_xs_queue.reserve(n_particles);	98✔
39	}
40	simulation::advance_particle_queue.reserve(n_particles);	375✔
41	simulation::surface_crossing_queue.reserve(n_particles);	375✔
42	simulation::collision_queue.reserve(n_particles);	375✔
43
44	simulation::particles.resize(n_particles);	375✔
45	}	375✔
46
47	void free_event_queues(void)	×
48	{
49	simulation::calculate_fuel_xs_queue.clear();	×
50	simulation::calculate_nonfuel_xs_queue.clear();	×
NEW 51	if (settings::delta_tracking) {	×
NEW 52	simulation::calculate_majorant_xs_queue.clear();	×
53	}
54	simulation::advance_particle_queue.clear();	×
55	simulation::surface_crossing_queue.clear();	×
56	simulation::collision_queue.clear();	×
57
58	simulation::particles.clear();	×
59	}	×
60
61	void dispatch_xs_event(int64_t buffer_idx)	306,170,368✔
62	{
63	Particle& p = simulation::particles[buffer_idx];	306,170,368✔
64	if (p.material() == MATERIAL_VOID \|\|	306,170,368✔
65	!model::materials[p.material()]->fissionable()) {	287,221,276✔
66	simulation::calculate_nonfuel_xs_queue.thread_safe_append({p, buffer_idx});	103,968,268✔
67	} else {
68	simulation::calculate_fuel_xs_queue.thread_safe_append({p, buffer_idx});	202,202,100✔
69	}
70	}	306,170,368✔
71
72	void process_death_events(int64_t n_particles)	4,906✔
73	{
74	simulation::time_event_death.start();	4,906✔
75	#pragma omp parallel for schedule(runtime)	4,546✔
76	for (int64_t i = 0; i < n_particles; i++) {	300,360✔
77	Particle& p = simulation::particles[i];	300,000✔
78	p.event_death();	300,000✔
79	}
80	simulation::time_event_death.stop();	4,906✔
81	}	4,906✔
82
83	//==============================================================================
84	// Functions for surface tracking
85	//==============================================================================
86
87	void process_init_events(int64_t n_particles, int64_t source_offset)	3,242✔
88	{
89	simulation::time_event_init.start();	3,242✔
90	#pragma omp parallel for schedule(runtime)	3,242✔
91	for (int64_t i = 0; i < n_particles; i++) {	×
92	initialize_particle_track(
93	simulation::particles[i], source_offset + i + 1, false);
94	dispatch_xs_event(i);
95	}
96	simulation::time_event_init.stop();	3,242✔
97	}	3,242✔
98
99	void process_calculate_xs_events(SharedArray<EventQueueItem>& queue)	1,399,162✔
100	{
101	simulation::time_event_calculate_xs.start();	1,399,162✔
102
103	// TODO: If using C++17, we could perform a parallel sort of the queue by
104	// particle type, material type, and then energy, in order to improve cache
105	// locality and reduce thread divergence on GPU. However, the parallel
106	// algorithms typically require linking against an additional library (Intel
107	// TBB). Prior to C++17, std::sort is a serial only operation, which in this
108	// case makes it too slow to be practical for most test problems.
109	//
110	// std::sort(std::execution::par_unseq, queue.data(), queue.data() +
111	// queue.size());
112
113	int64_t offset = simulation::advance_particle_queue.size();	1,399,162✔
114
115	#pragma omp parallel for schedule(runtime)	1,399,162✔
116	for (int64_t i = 0; i < queue.size(); i++) {	×
117	Particle* p = &simulation::particles[queue[i].idx];
118	p->event_calculate_xs();
119
120	// After executing a calculate_xs event, particles will
121	// always require an advance event. Therefore, we don't need to use
122	// the protected enqueuing function.
123	simulation::advance_particle_queue[offset + i] = queue[i];
124	}
125
126	simulation::advance_particle_queue.resize(offset + queue.size());	1,399,162✔
127
128	queue.resize(0);	1,399,162✔
129
130	simulation::time_event_calculate_xs.stop();	1,399,162✔
131	}	1,399,162✔
132
133	void process_advance_particle_events()	1,394,191✔
134	{
135	simulation::time_event_advance_particle.start();	1,394,191✔
136
137	#pragma omp parallel for schedule(runtime)	1,394,191✔
138	for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {	×
139	int64_t buffer_idx = simulation::advance_particle_queue[i].idx;
140	Particle& p = simulation::particles[buffer_idx];
141	p.event_advance();
142	if (!p.alive())	×
143	continue;
144	if (p.collision_distance() > p.boundary().distance()) {	×
145	simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx});
146	} else {
147	simulation::collision_queue.thread_safe_append({p, buffer_idx});
148	}
149	}
150
151	simulation::advance_particle_queue.resize(0);	1,394,191✔
152
153	simulation::time_event_advance_particle.stop();	1,394,191✔
154	}	1,394,191✔
155
156	void process_surface_crossing_events()	267,591✔
157	{
158	simulation::time_event_surface_crossing.start();	267,591✔
159
160	#pragma omp parallel for schedule(runtime)	267,591✔
161	for (int64_t i = 0; i < simulation::surface_crossing_queue.size(); i++) {	×
162	int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;
163	Particle& p = simulation::particles[buffer_idx];
164	p.event_cross_surface();
165	p.event_check_limit_and_revive();
166	if (p.alive())	×
167	dispatch_xs_event(buffer_idx);
168	}
169
170	simulation::surface_crossing_queue.resize(0);	267,591✔
171
172	simulation::time_event_surface_crossing.stop();	267,591✔
173	}	267,591✔
174
175	void process_collision_events()	1,146,934✔
176	{
177	simulation::time_event_collision.start();	1,146,934✔
178
179	#pragma omp parallel for schedule(runtime)	1,146,934✔
180	for (int64_t i = 0; i < simulation::collision_queue.size(); i++) {	×
181	int64_t buffer_idx = simulation::collision_queue[i].idx;
182	Particle& p = simulation::particles[buffer_idx];
183	p.event_collide();
184	p.event_check_limit_and_revive();
185	if (p.alive())	×
186	dispatch_xs_event(buffer_idx);
187	}
188
189	simulation::collision_queue.resize(0);	1,146,934✔
190
191	simulation::time_event_collision.stop();	1,146,934✔
192	}	1,146,934✔
193
194	void process_transport_events()	3,926✔
195	{
196	while (true) {	4,211,804✔
197	int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),	4,211,804✔
198	simulation::calculate_nonfuel_xs_queue.size(),	4,211,804✔
199	simulation::advance_particle_queue.size(),	4,211,804✔
200	simulation::surface_crossing_queue.size(),	4,211,804✔
201	simulation::collision_queue.size()});	4,211,804✔
202
203	if (max == 0) {	4,211,804✔
204	break;
205	} else if (max == simulation::calculate_fuel_xs_queue.size()) {	4,207,878✔
206	process_calculate_xs_events(simulation::calculate_fuel_xs_queue);	465,856✔
207	} else if (max == simulation::calculate_nonfuel_xs_queue.size()) {	3,742,022✔
208	process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);	933,306✔
209	} else if (max == simulation::advance_particle_queue.size()) {	2,808,716✔
210	process_advance_particle_events();	1,394,191✔
211	} else if (max == simulation::surface_crossing_queue.size()) {	1,414,525✔
212	process_surface_crossing_events();	267,591✔
213	} else if (max == simulation::collision_queue.size()) {	1,146,934!
214	process_collision_events();	1,146,934✔
215	}
216	}
217	}	3,926✔
218
219	void process_init_secondary_events(int64_t n_particles, int64_t offset,	684✔
220	const SharedArray<SourceSite>& shared_secondary_bank)
221	{
222	simulation::time_event_init.start();	684✔
223	#pragma omp parallel for schedule(runtime)	684✔
224	for (int64_t i = 0; i < n_particles; i++) {	×
225	initialize_particle_track(simulation::particles[i], offset + i + 1, true);
226	const SourceSite& site = shared_secondary_bank[offset + i];
227	simulation::particles[i].event_revive_from_secondary(site);
228	if (simulation::particles[i].alive()) {	×
229	dispatch_xs_event(i);
230	}
231	}
232	simulation::time_event_init.stop();	684✔
233	}	684✔
234
235	//==============================================================================
236	// Functions for hybrid delta tracking
237	//==============================================================================
238
239	void process_hybrid_init_events(int64_t n_particles, int64_t source_offset)	980✔
240	{
241	simulation::time_event_init.start();	980✔
242	#pragma omp parallel for schedule(runtime)	620✔
243	for (int64_t i = 0; i < n_particles; i++) {	300,360✔
244	initialize_particle_track(	300,000✔
245	simulation::particles[i], source_offset + i + 1, false);	300,000✔
246
247	simulation::calculate_majorant_xs_queue[i] = {simulation::particles[i], i};	300,000✔
248	}
249	simulation::calculate_majorant_xs_queue.resize(n_particles);	980✔
250	simulation::time_event_init.stop();	980✔
251	}	980✔
252
253	void process_hybrid_calculate_xs_events(SharedArray<EventQueueItem>& queue)	1,173,997✔
254	{
255	simulation::time_event_calculate_xs.start();	1,173,997✔
256
257	// TODO: If using C++17, we could perform a parallel sort of the queue by
258	// particle type, material type, and then energy, in order to improve cache
259	// locality and reduce thread divergence on GPU. However, the parallel
260	// algorithms typically require linking against an additional library (Intel
261	// TBB). Prior to C++17, std::sort is a serial only operation, which in this
262	// case makes it too slow to be practical for most test problems.
263	//
264	// std::sort(std::execution::par_unseq, queue.data(), queue.data() +
265	// queue.size());
266
267	#pragma omp parallel for schedule(runtime)	739,838✔
268	for (int64_t i = 0; i < queue.size(); i++) {	30,848,454✔
269	int64_t buffer_idx = queue[i].idx;	30,414,295✔
270	Particle& p = simulation::particles[buffer_idx];	30,414,295✔
271	p.event_calculate_xs();	30,414,295✔
272	if (!p.alive()) {	30,414,295!
273	continue;
274	}
275
276	// Particles always require a delta collision event after cross sections
277	// are computed if running delta tracking. Otherwise, an advance event is
278	// required for surface tracking.
279	if (p.delta_tracking()) {	30,414,295✔
280	simulation::collision_queue.thread_safe_append({p, buffer_idx});	19,577,530✔
281	} else {
282	simulation::advance_particle_queue.thread_safe_append({p, buffer_idx});	10,836,765✔
283	}
284	}
285
286	queue.resize(0);	1,173,997✔
287
288	simulation::time_event_calculate_xs.stop();	1,173,997✔
289	}	1,173,997✔
290
291	void process_hybrid_calculate_majorant_events()	1,184,481✔
292	{
293	simulation::time_event_calculate_majorant_xs.start();	1,184,481✔
294
295	#pragma omp parallel for schedule(runtime)	745,269✔
296	for (int64_t i = 0; i < simulation::calculate_majorant_xs_queue.size(); i++) {	34,844,362✔
297	int64_t buffer_idx = simulation::calculate_majorant_xs_queue[i].idx;	34,405,150✔
298	Particle& p = simulation::particles[buffer_idx];	34,405,150✔
299
300	p.event_update_majorant();	34,405,150✔
301
302	// Need to compute cross sections before advancing when running surface
303	// tracking. Otherwise, we can go right to advancing with delta tracking.
304	if (!p.delta_tracking()) {	34,405,150✔
305	dispatch_xs_event(buffer_idx);	10,836,765✔
306	} else {
307	simulation::advance_particle_queue.thread_safe_append({p, buffer_idx});	23,568,385✔
308	}
309	}
310
311	simulation::calculate_majorant_xs_queue.resize(0);	1,184,481✔
312
313	simulation::time_event_calculate_majorant_xs.stop();	1,184,481✔
314	}	1,184,481✔
315
316	void process_hybrid_advance_particle_events()	1,133,933✔
317	{
318	simulation::time_event_advance_particle.start();	1,133,933✔
319
320	#pragma omp parallel for schedule(runtime)	715,323✔
321	for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) {	34,823,760✔
322	int64_t buffer_idx = simulation::advance_particle_queue[i].idx;	34,405,150✔
323	Particle& p = simulation::particles[buffer_idx];	34,405,150✔
324	p.event_delta_advance();	34,405,150✔
325	if (!p.alive()) {	34,405,150!
326	continue;
327	}
328
329	if (p.collision_distance() < p.boundary().distance()) {	34,405,150✔
330	// We need to compute cross sections prior to processing a collision when
331	// running delta tracking. Otherwise, we can go right to a collision when
332	// running surface tracking.
333	if (p.delta_tracking()) {	23,946,145✔
334	dispatch_xs_event(buffer_idx);	19,577,530✔
335	} else {
336	simulation::collision_queue.thread_safe_append({p, buffer_idx});	4,368,615✔
337	}
338	} else {
339	simulation::surface_crossing_queue.thread_safe_append({p, buffer_idx});	10,459,005✔
340	}
341	}
342
343	simulation::advance_particle_queue.resize(0);	1,133,933✔
344
345	simulation::time_event_advance_particle.stop();	1,133,933✔
346	}	1,133,933✔
347
348	void process_hybrid_surface_crossing_events()	253,647✔
349	{
350	simulation::time_event_surface_crossing.start();	253,647✔
351
352	#pragma omp parallel for schedule(runtime)	153,221✔
353	for (int64_t i = 0; i < simulation::surface_crossing_queue.size(); i++) {	10,559,431✔
354	int64_t buffer_idx = simulation::surface_crossing_queue[i].idx;	10,459,005✔
355	Particle& p = simulation::particles[buffer_idx];	10,459,005✔
356	p.event_cross_surface();	10,459,005✔
357	p.event_check_limit_and_revive();	10,459,005✔
358	if (p.alive())	10,459,005!
359	simulation::calculate_majorant_xs_queue.thread_safe_append(	10,459,005✔
360	{p, buffer_idx});
361	}
362
363	simulation::surface_crossing_queue.resize(0);	253,647✔
364
365	simulation::time_event_surface_crossing.stop();	253,647✔
366	}	253,647✔
367
368	void process_hybrid_collision_events()	943,743✔
369	{
370	simulation::time_event_collision.start();	943,743✔
371
372	#pragma omp parallel for schedule(runtime)	600,186✔
373	for (int64_t i = 0; i < simulation::collision_queue.size(); i++) {	24,289,702✔
374	int64_t buffer_idx = simulation::collision_queue[i].idx;	23,946,145✔
375	Particle& p = simulation::particles[buffer_idx];	23,946,145✔
376
377	if (p.delta_tracking()) {	23,946,145✔
378	// Check to ensure the majorant is valid.
379	if (p.kill_invalid_maj()) {	19,577,530!
380	continue;
381	}
382
383	// Perform rejection sampling. If this is true, a real collision is
384	// processed prior to going back to majorant calculations.
385	if (prn(p.current_seed()) < (p.macro_xs().total / p.majorant())) {	19,577,530✔
386	p.event_collide();	9,299,615✔
387	}
388	} else {
389	p.event_collide();	4,368,615✔
390	}
391
392	// Update the tracking type based on the chosen hybrid scheme.
393	p.update_tracking_type();	23,946,145✔
394
395	p.event_check_limit_and_revive();	23,946,145✔
396	if (p.alive()) {	23,946,145✔
397	simulation::calculate_majorant_xs_queue.thread_safe_append(	23,646,145✔
398	{p, buffer_idx});
399	}
400	}
401
402	simulation::collision_queue.resize(0);	943,743✔
403
404	simulation::time_event_collision.stop();	943,743✔
405	}	943,743✔
406
407	void process_hybrid_transport_events()	980✔
408	{
409	while (true) {	4,690,781✔
410	int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),	4,690,781✔
411	simulation::calculate_nonfuel_xs_queue.size(),	4,690,781✔
412	simulation::calculate_majorant_xs_queue.size(),	4,690,781✔
413	simulation::advance_particle_queue.size(),	4,690,781✔
414	simulation::surface_crossing_queue.size(),	4,690,781✔
415	simulation::collision_queue.size()});	4,690,781✔
416
417	if (max == 0) {	4,690,781✔
418	break;
419	} else if (max == simulation::calculate_fuel_xs_queue.size()) {	4,689,801✔
420	process_hybrid_calculate_xs_events(simulation::calculate_fuel_xs_queue);	280,549✔
421	} else if (max == simulation::calculate_nonfuel_xs_queue.size()) {	4,409,252✔
422	process_hybrid_calculate_xs_events(	893,448✔
423	simulation::calculate_nonfuel_xs_queue);
424	} else if (max == simulation::calculate_majorant_xs_queue.size()) {	3,515,804✔
425	process_hybrid_calculate_majorant_events();	1,184,481✔
426	} else if (max == simulation::advance_particle_queue.size()) {	2,331,323✔
427	process_hybrid_advance_particle_events();	1,133,933✔
428	} else if (max == simulation::surface_crossing_queue.size()) {	1,197,390✔
429	process_hybrid_surface_crossing_events();	253,647✔
430	} else if (max == simulation::collision_queue.size()) {	943,743!
431	process_hybrid_collision_events();	943,743✔
432	}
433	}
434	}	980✔
435
NEW 436	void process_hybrid_init_secondary_events(int64_t n_particles, int64_t offset,	×
437	const SharedArray<SourceSite>& shared_secondary_bank)
438	{
NEW 439	simulation::time_event_init.start();	×
440	#pragma omp parallel for schedule(runtime)
441	for (int64_t i = 0; i < n_particles; i++) {	×
442	Particle& p = simulation::particles[i];
443	initialize_particle_track(p, offset + i + 1, true);
444	const SourceSite& site = shared_secondary_bank[offset + i];
445	p.event_revive_from_secondary(site);
446
447	if (p.alive()) {	×
448	simulation::calculate_majorant_xs_queue.thread_safe_append({p, i});
449	}
450	}
NEW 451	simulation::time_event_init.stop();	×
NEW 452	}	×
453
454	} // namespace openmc

openmc-dev / openmc / 30942440215

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

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