26675018610

Committed 30 May 2026 04:53AM UTC coverage: 81.343% (+0.006%) from 81.337%

Build # 26675018610

Build Type

Pull #3915

github

Committed by

web-flow

Commit Message

Merge b70d83d50 into 46d889613

Pull Request Pull Request #3915: Allow tracklength estimator for neutron heating

Coverage Stats

18010 of 26105 branches covered (68.99%)

Branch coverage included in aggregate %.

9 of 9 new or added lines in 1 file covered. (100.0%)

1782 existing lines in 46 files now uncovered.

59130 of 68728 relevant lines covered (86.03%)

48511807.72 hits per line

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

88.24

/src/physics_mg.cpp

#include "openmc/physics_mg.h"

#include <stdexcept>

#include "openmc/tensor.h"
#include <fmt/core.h>

#include "openmc/bank.h"
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/math_functions.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/particle.h"
#include "openmc/physics_common.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/tally.h"
#include "openmc/weight_windows.h"

namespace openmc {

void collision_mg(Particle& p)
{
  // Add to the collision counter for the particle
  p.n_collision()++;
  p.secondary_bank_index() = p.local_secondary_bank().size();

  // Sample the reaction type
  sample_reaction(p);

  if (settings::weight_windows_on) {
    auto [ww_found, ww] = search_weight_window(p);
    if (!ww_found && p.type() == ParticleType::neutron()) {
      // if the weight window is not valid, apply russian roulette
      // (regardless of weight window collision checkpoint setting)
      apply_russian_roulette(p);
    } else if (settings::weight_window_checkpoint_collision) {
      // if collision checkpointing is on, apply weight window
      apply_weight_window(p, ww);
    }
  }

  // Display information about collision
  if ((settings::verbosity >= 10) || p.trace()) {
    write_message(fmt::format("    Energy Group = {}", p.g()), 1);
  }
}

void sample_reaction(Particle& p)
{
  // Create fission bank sites. Note that while a fission reaction is sampled,
  // it never actually "happens", i.e. the weight of the particle does not
  // change when sampling fission sites. The following block handles all
  // absorption (including fission)

  if (model::materials[p.material()]->fissionable()) {
    if (settings::run_mode == RunMode::EIGENVALUE ||
        (settings::run_mode == RunMode::FIXED_SOURCE &&
          settings::create_fission_neutrons)) {
      create_fission_sites(p);
    }
  }

  // If survival biasing is being used, the following subroutine adjusts the
  // weight of the particle. Otherwise, it checks to see if absorption occurs.
  if (p.macro_xs().absorption > 0.) {
    absorption(p);
  }
  if (!p.alive())
    return;

  // Sample a scattering event to determine the energy of the exiting neutron
  scatter(p);

  // Play russian roulette if there are no weight windows
  if (!settings::weight_windows_on)
    apply_russian_roulette(p);
}

void scatter(Particle& p)
{
  data::mg.macro_xs_[p.material()].sample_scatter(p.g_last(), p.g(), p.mu(),
    p.wgt(), p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a);

  // Rotate the angle
  p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed());

  // Update energy value for downstream compatability (in tallying)
  p.E() = data::mg.energy_bin_avg_[p.g()];

  // Set event component
  p.event() = TallyEvent::SCATTER;
}

void create_fission_sites(Particle& p)
{
  // If uniform fission source weighting is turned on, we increase or decrease
  // the expected number of fission sites produced
  double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;

  // Determine the expected number of neutrons produced
  double nu_t = p.wgt() / simulation::keff * weight * p.macro_xs().nu_fission /
                p.macro_xs().total;

  // Sample the number of neutrons produced
  int nu = static_cast<int>(nu_t);
  if (prn(p.current_seed()) <= (nu_t - int(nu_t))) {
    nu++;
  }

  // If no neutrons were produced then don't continue
  if (nu == 0)
    return;

  // Initialize the counter of delayed neutrons encountered for each delayed
  // group.
  double nu_d[MAX_DELAYED_GROUPS] = {0.};

  // Clear out particle's nu fission bank
  p.nu_bank().clear();

  p.fission() = true;

  // Determine whether to place fission sites into the shared fission bank
  // or the secondary particle bank.
  bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);

  // Counter for the number of fission sites successfully stored to the shared
  // fission bank or the secondary particle bank
  int n_sites_stored;

  for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) {
    // Initialize fission site object with particle data
    SourceSite site;
    site.r = p.r();
    site.particle = ParticleType::neutron();
    site.time = p.time();
    site.wgt = 1. / weight;

    // Sample the cosine of the angle, assuming fission neutrons are emitted
    // isotropically
    double mu = 2. * prn(p.current_seed()) - 1.;

    // Sample the azimuthal angle uniformly in [0, 2.pi)
    double phi = 2. * PI * prn(p.current_seed());
    site.u.x = mu;
    site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi);
    site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi);

    // Sample secondary energy distribution for the fission reaction
    int dg;
    int gout;
    data::mg.macro_xs_[p.material()].sample_fission_energy(
      p.g(), dg, gout, p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a);

    // Store the energy and delayed groups on the fission bank
    site.E = gout;

    // We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest
    // of the code, 0 is prompt.
    site.delayed_group = dg + 1;

    // If delayed product production, sample time of emission
    if (dg != -1) {
      auto& macro_xs = data::mg.macro_xs_[p.material()];
      double decay_rate =
        macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0);
      site.time -= std::log(prn(p.current_seed())) / decay_rate;

      // Reject site if it exceeds time cutoff
      double t_cutoff = settings::time_cutoff[site.particle.transport_index()];
      if (site.time > t_cutoff) {
        continue;
      }
    }

    // Set parent and progeny ID
    site.parent_id = p.current_work();
    site.progeny_id = p.n_progeny()++;

    // Store fission site in bank
    if (use_fission_bank) {
      int64_t idx = simulation::fission_bank.thread_safe_append(site);
      if (idx == -1) {
        warning(
          "The shared fission bank is full. Additional fission sites created "
          "in this generation will not be banked. Results may be "
          "non-deterministic.");

        // Decrement number of particle progeny as storage was unsuccessful.
        // This step is needed so that the sum of all progeny is equal to the
        // size of the shared fission bank.
        p.n_progeny()--;

        // Break out of loop as no more sites can be added to fission bank
        break;
      }
    } else {
      site.wgt_born = p.wgt_born();
      site.wgt_ww_born = p.wgt_ww_born();
      site.n_split = p.n_split();
      p.local_secondary_bank().push_back(site);
      p.n_secondaries()++;
    }

    // Set the delayed group on the particle as well
    p.delayed_group() = dg + 1;

    // Increment the number of neutrons born delayed
    if (p.delayed_group() > 0) {
      nu_d[dg]++;
    }

    // Write fission particles to nuBank
    NuBank& nu_bank_entry = p.nu_bank().emplace_back();
    nu_bank_entry.wgt = site.wgt;
    nu_bank_entry.E = site.E;
    nu_bank_entry.delayed_group = site.delayed_group;
  }

  // If shared fission bank was full, and no fissions could be added,
  // set the particle fission flag to false.
  if (n_sites_stored == 0) {
    p.fission() = false;
    return;
  }

  // Set nu to the number of fission sites successfully stored. If the fission
  // bank was not found to be full then these values are already equivalent.
  nu = n_sites_stored;

  // Store the total weight banked for analog fission tallies
  p.n_bank() = nu;
  p.wgt_bank() = nu / weight;
  for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
    p.n_delayed_bank(d) = nu_d[d];
  }
}

void absorption(Particle& p)
{
  if (settings::survival_biasing) {
    // Determine weight absorbed in survival biasing
    double wgt_absorb = p.wgt() * p.macro_xs().absorption / p.macro_xs().total;

    // Adjust weight of particle by the probability of absorption
    p.wgt() -= wgt_absorb;

    // Score implicit absorpion estimate of keff
    p.keff_tally_absorption() +=
      wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption;
  } else {
    if (p.macro_xs().absorption > prn(p.current_seed()) * p.macro_xs().total) {
      p.keff_tally_absorption() +=
        p.wgt() * p.macro_xs().nu_fission / p.macro_xs().absorption;
      p.wgt() = 0.0;
      p.event() = TallyEvent::ABSORB;
    }
  }
}

} // namespace openmc

1	#include "openmc/physics_mg.h"
2
3	#include <stdexcept>
4
5	#include "openmc/tensor.h"
6	#include <fmt/core.h>
7
8	#include "openmc/bank.h"
9	#include "openmc/constants.h"
10	#include "openmc/eigenvalue.h"
11	#include "openmc/error.h"
12	#include "openmc/material.h"
13	#include "openmc/math_functions.h"
14	#include "openmc/message_passing.h"
15	#include "openmc/mgxs_interface.h"
16	#include "openmc/particle.h"
17	#include "openmc/physics_common.h"
18	#include "openmc/random_lcg.h"
19	#include "openmc/settings.h"
20	#include "openmc/simulation.h"
21	#include "openmc/tallies/tally.h"
22	#include "openmc/weight_windows.h"
23
24	namespace openmc {
25
26	void collision_mg(Particle& p)	1,800,818,030✔
27	{
28	// Add to the collision counter for the particle
29	p.n_collision()++;	1,800,818,030✔
30	p.secondary_bank_index() = p.local_secondary_bank().size();	1,800,818,030✔
31
32	// Sample the reaction type
33	sample_reaction(p);	1,800,818,030✔
34
35	if (settings::weight_windows_on) {	1,800,818,030✔
36	auto [ww_found, ww] = search_weight_window(p);	20,814,640✔
37	if (!ww_found && p.type() == ParticleType::neutron()) {	20,814,640!
38	// if the weight window is not valid, apply russian roulette
39	// (regardless of weight window collision checkpoint setting)
40	apply_russian_roulette(p);	10,780✔
41	} else if (settings::weight_window_checkpoint_collision) {	20,803,860!
42	// if collision checkpointing is on, apply weight window
43	apply_weight_window(p, ww);	20,803,860✔
44	}
45	}
46
47	// Display information about collision
48	if ((settings::verbosity >= 10) \|\| p.trace()) {	1,800,818,030!
49	write_message(fmt::format(" Energy Group = {}", p.g()), 1);	×
50	}
51	}	1,800,818,030✔
52
53	void sample_reaction(Particle& p)	1,800,818,030✔
54	{
55	// Create fission bank sites. Note that while a fission reaction is sampled,
56	// it never actually "happens", i.e. the weight of the particle does not
57	// change when sampling fission sites. The following block handles all
58	// absorption (including fission)
59
60	if (model::materials[p.material()]->fissionable()) {	1,800,818,030✔
61	if (settings::run_mode == RunMode::EIGENVALUE \|\|	1,795,603,689!
62	(settings::run_mode == RunMode::FIXED_SOURCE &&	17,806,866!
63	settings::create_fission_neutrons)) {
64	create_fission_sites(p);	1,795,603,689✔
65	}
66	}
67
68	// If survival biasing is being used, the following subroutine adjusts the
69	// weight of the particle. Otherwise, it checks to see if absorption occurs.
70	if (p.macro_xs().absorption > 0.) {	1,800,818,030!
71	absorption(p);	1,800,818,030✔
72	}
73	if (!p.alive())	1,800,818,030✔
74	return;
75
76	// Sample a scattering event to determine the energy of the exiting neutron
77	scatter(p);	1,686,408,922✔
78
79	// Play russian roulette if there are no weight windows
80	if (!settings::weight_windows_on)	1,686,408,922✔
81	apply_russian_roulette(p);	1,668,643,207✔
82	}
83
84	void scatter(Particle& p)	1,686,408,922✔
85	{
86	data::mg.macro_xs_[p.material()].sample_scatter(p.g_last(), p.g(), p.mu(),	1,686,408,922✔
87	p.wgt(), p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a);	1,686,408,922✔
88
89	// Rotate the angle
90	p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed());	1,686,408,922✔
91
92	// Update energy value for downstream compatability (in tallying)
93	p.E() = data::mg.energy_bin_avg_[p.g()];	1,686,408,922✔
94
95	// Set event component
96	p.event() = TallyEvent::SCATTER;	1,686,408,922✔
97	}	1,686,408,922✔
98
99	void create_fission_sites(Particle& p)	1,795,603,689✔
100	{
101	// If uniform fission source weighting is turned on, we increase or decrease
102	// the expected number of fission sites produced
103	double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;	1,795,603,689!
104
105	// Determine the expected number of neutrons produced
106	double nu_t = p.wgt() / simulation::keff * weight * p.macro_xs().nu_fission /	1,795,603,689✔
107	p.macro_xs().total;	1,795,603,689✔
108
109	// Sample the number of neutrons produced
110	int nu = static_cast<int>(nu_t);	1,795,603,689✔
111	if (prn(p.current_seed()) <= (nu_t - int(nu_t))) {	1,795,603,689✔
112	nu++;	111,887,534✔
113	}
114
115	// If no neutrons were produced then don't continue
116	if (nu == 0)	1,795,603,689✔
117	return;	1,683,715,517✔
118
119	// Initialize the counter of delayed neutrons encountered for each delayed
120	// group.
121	double nu_d[MAX_DELAYED_GROUPS] = {0.};	111,888,172✔
122
123	// Clear out particle's nu fission bank
124	p.nu_bank().clear();	111,888,172✔
125
126	p.fission() = true;	111,888,172✔
127
128	// Determine whether to place fission sites into the shared fission bank
129	// or the secondary particle bank.
130	bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);	111,888,172✔
131
132	// Counter for the number of fission sites successfully stored to the shared
133	// fission bank or the secondary particle bank
134	int n_sites_stored;	111,888,172✔
135
136	for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) {	223,778,071✔
137	// Initialize fission site object with particle data
138	SourceSite site;	111,889,899✔
139	site.r = p.r();	111,889,899✔
140	site.particle = ParticleType::neutron();	111,889,899✔
141	site.time = p.time();	111,889,899✔
142	site.wgt = 1. / weight;	111,889,899✔
143
144	// Sample the cosine of the angle, assuming fission neutrons are emitted
145	// isotropically
146	double mu = 2. * prn(p.current_seed()) - 1.;	111,889,899✔
147
148	// Sample the azimuthal angle uniformly in [0, 2.pi)
149	double phi = 2. * PI * prn(p.current_seed());	111,889,899✔
150	site.u.x = mu;	111,889,899✔
151	site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi);	111,889,899✔
152	site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi);	111,889,899✔
153
154	// Sample secondary energy distribution for the fission reaction
155	int dg;	111,889,899✔
156	int gout;	111,889,899✔
157	data::mg.macro_xs_[p.material()].sample_fission_energy(	111,889,899✔
158	p.g(), dg, gout, p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a);	111,889,899✔
159
160	// Store the energy and delayed groups on the fission bank
161	site.E = gout;	111,889,899✔
162
163	// We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest
164	// of the code, 0 is prompt.
165	site.delayed_group = dg + 1;	111,889,899✔
166
167	// If delayed product production, sample time of emission
168	if (dg != -1) {	111,889,899✔
169	auto& macro_xs = data::mg.macro_xs_[p.material()];	1,507✔
170	double decay_rate =	1,507✔
171	macro_xs.get_xs(MgxsType::DECAY_RATE, 0, nullptr, nullptr, &dg, 0, 0);	1,507✔
172	site.time -= std::log(prn(p.current_seed())) / decay_rate;	1,507✔
173
174	// Reject site if it exceeds time cutoff
175	double t_cutoff = settings::time_cutoff[site.particle.transport_index()];	1,507✔
176	if (site.time > t_cutoff) {	1,507✔
177	continue;	737✔
178	}
179	}
180
181	// Set parent and progeny ID
182	site.parent_id = p.current_work();	111,889,162✔
183	site.progeny_id = p.n_progeny()++;	111,889,162✔
184
185	// Store fission site in bank
186	if (use_fission_bank) {	111,889,162✔
187	int64_t idx = simulation::fission_bank.thread_safe_append(site);	111,089,198✔
188	if (idx == -1) {	111,089,198!
189	warning(	×
190	"The shared fission bank is full. Additional fission sites created "
191	"in this generation will not be banked. Results may be "
192	"non-deterministic.");
193
194	// Decrement number of particle progeny as storage was unsuccessful.
195	// This step is needed so that the sum of all progeny is equal to the
196	// size of the shared fission bank.
197	p.n_progeny()--;	×
198
199	// Break out of loop as no more sites can be added to fission bank
200	break;	×
201	}
202	} else {
203	site.wgt_born = p.wgt_born();	799,964✔
204	site.wgt_ww_born = p.wgt_ww_born();	799,964✔
205	site.n_split = p.n_split();	799,964✔
206	p.local_secondary_bank().push_back(site);	799,964✔
207	p.n_secondaries()++;	799,964✔
208	}
209
210	// Set the delayed group on the particle as well
211	p.delayed_group() = dg + 1;	111,889,162✔
212
213	// Increment the number of neutrons born delayed
214	if (p.delayed_group() > 0) {	111,889,162✔
215	nu_d[dg]++;	770✔
216	}
217
218	// Write fission particles to nuBank
219	NuBank& nu_bank_entry = p.nu_bank().emplace_back();	111,889,162✔
220	nu_bank_entry.wgt = site.wgt;	111,889,162✔
221	nu_bank_entry.E = site.E;	111,889,162✔
222	nu_bank_entry.delayed_group = site.delayed_group;	111,889,162✔
223	}
224
225	// If shared fission bank was full, and no fissions could be added,
226	// set the particle fission flag to false.
227	if (n_sites_stored == 0) {	111,888,172!
UNCOV 228	p.fission() = false;	×
UNCOV 229	return;	×
230	}
231
232	// Set nu to the number of fission sites successfully stored. If the fission
233	// bank was not found to be full then these values are already equivalent.
234	nu = n_sites_stored;	111,888,172✔
235
236	// Store the total weight banked for analog fission tallies
237	p.n_bank() = nu;	111,888,172✔
238	p.wgt_bank() = nu / weight;	111,888,172✔
239	for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {	1,006,993,548✔
240	p.n_delayed_bank(d) = nu_d[d];	895,105,376✔
241	}
242	}
243
244	void absorption(Particle& p)	1,800,818,030✔
245	{
246	if (settings::survival_biasing) {	1,800,818,030✔
247	// Determine weight absorbed in survival biasing
248	double wgt_absorb = p.wgt() * p.macro_xs().absorption / p.macro_xs().total;	4,219,259✔
249
250	// Adjust weight of particle by the probability of absorption
251	p.wgt() -= wgt_absorb;	4,219,259✔
252
253	// Score implicit absorpion estimate of keff
254	p.keff_tally_absorption() +=	4,219,259✔
255	wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption;	4,219,259✔
256	} else {
257	if (p.macro_xs().absorption > prn(p.current_seed()) * p.macro_xs().total) {	1,796,598,771✔
258	p.keff_tally_absorption() +=	114,409,108✔
259	p.wgt() * p.macro_xs().nu_fission / p.macro_xs().absorption;	114,409,108✔
260	p.wgt() = 0.0;	114,409,108✔
261	p.event() = TallyEvent::ABSORB;	114,409,108✔
262	}
263	}
264	}	1,800,818,030✔
265
266	} // namespace openmc

openmc-dev / openmc / 26675018610

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