9789289404

Committed 04 Jul 2024 06:02AM UTC coverage: 84.763% (-0.003%) from 84.766%

Build # 9789289404

Build Type

Pull #3070

github

Committed by

web-flow

Commit Message

Merge d5ce44080 into 1b22dd28d

Pull Request Pull Request #3070: Survival Normalization 0.15.1

Run Details

8 of 13 new or added lines in 4 files covered. (61.54%)

1 existing line in 1 file now uncovered.

49186 of 58028 relevant lines covered (84.76%)

31248475.07 hits per line

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

89.25

/src/physics_mg.cpp

#include "openmc/physics_mg.h"

#include <stdexcept>

#include "xtensor/xarray.hpp"
#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"

namespace openmc {

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

  // Sample the reaction type
  sample_reaction(p);

  // 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 survival biasing is turned on
  if (settings::survival_biasing) {
    // if survival normalization is applicable, use normalized weight cutoff and normalized weight survive
    if (settings::survival_normalization) {
      if (p.wgt() < settings::weight_cutoff*p.wgt0()) {
        russian_roulette(p, settings::weight_survive*p.wgt0());
      } 
    } else if (p.wgt() < settings::weight_cutoff) {
            russian_roulette(p, settings::weight_survive);
    }
  }
}

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.wgt = 1. / weight;
    site.parent_id = p.id();
    site.progeny_id = p.n_progeny()++;

    // 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;

    // 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 {
      p.secondary_bank().push_back(site);
    }

    // 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
    p.nu_bank().emplace_back();
    NuBank* nu_bank_entry = &p.nu_bank().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 "xtensor/xarray.hpp"
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
23	namespace openmc {
24
25	void collision_mg(Particle& p)	1,939,604,940✔
26	{
27	// Add to the collision counter for the particle
28	p.n_collision()++;	1,939,604,940✔
29
30	// Sample the reaction type
31	sample_reaction(p);	1,939,604,940✔
32
33	// Display information about collision
34	if ((settings::verbosity >= 10) \|\| p.trace()) {	1,939,604,940✔
35	write_message(fmt::format(" Energy Group = {}", p.g()), 1);	×
36	}
37	}	1,939,604,940✔
38
39	void sample_reaction(Particle& p)	1,939,604,940✔
40	{
41	// Create fission bank sites. Note that while a fission reaction is sampled,
42	// it never actually "happens", i.e. the weight of the particle does not
43	// change when sampling fission sites. The following block handles all
44	// absorption (including fission)
45
46	if (model::materials[p.material()]->fissionable()) {	1,939,604,940✔
47	if (settings::run_mode == RunMode::EIGENVALUE \|\|	1,937,955,456✔
48	(settings::run_mode == RunMode::FIXED_SOURCE &&	×
49	settings::create_fission_neutrons)) {
50	create_fission_sites(p);	1,937,955,456✔
51	}
52	}
53
54	// If survival biasing is being used, the following subroutine adjusts the
55	// weight of the particle. Otherwise, it checks to see if absorption occurs.
56	if (p.macro_xs().absorption > 0.) {	1,939,604,940✔
57	absorption(p);	1,939,604,940✔
58	}
59	if (!p.alive())	1,939,604,940✔
60	return;	120,885,540✔
61
62	// Sample a scattering event to determine the energy of the exiting neutron
63	scatter(p);	1,818,719,400✔
64
65	// Play Russian roulette if survival biasing is turned on
66	if (settings::survival_biasing) {	1,818,719,400✔
67	// if survival normalization is applicable, use normalized weight cutoff and normalized weight survive
68	if (settings::survival_normalization) {	4,596,048✔
NEW 69	if (p.wgt() < settings::weight_cutoff*p.wgt0()) {	×
NEW 70	russian_roulette(p, settings::weight_survive*p.wgt0());	×
71	}
72	} else if (p.wgt() < settings::weight_cutoff) {	4,596,048✔
73	russian_roulette(p, settings::weight_survive);	155,124✔
74	}
75	}
76	}
77
78	void scatter(Particle& p)	1,818,719,400✔
79	{
80	data::mg.macro_xs_[p.material()].sample_scatter(p.g_last(), p.g(), p.mu(),	2,147,483,647✔
81	p.wgt(), p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a);	1,818,719,400✔
82
83	// Rotate the angle
84	p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed());	1,818,719,400✔
85
86	// Update energy value for downstream compatability (in tallying)
87	p.E() = data::mg.energy_bin_avg_[p.g()];	1,818,719,400✔
88
89	// Set event component
90	p.event() = TallyEvent::SCATTER;	1,818,719,400✔
91	}	1,818,719,400✔
92
93	void create_fission_sites(Particle& p)	1,937,955,456✔
94	{
95	// If uniform fission source weighting is turned on, we increase or decrease
96	// the expected number of fission sites produced
97	double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;	1,937,955,456✔
98
99	// Determine the expected number of neutrons produced
100	double nu_t = p.wgt() / simulation::keff * weight * p.macro_xs().nu_fission /	1,937,955,456✔
101	p.macro_xs().total;	1,937,955,456✔
102
103	// Sample the number of neutrons produced
104	int nu = static_cast<int>(nu_t);	1,937,955,456✔
105	if (prn(p.current_seed()) <= (nu_t - int(nu_t))) {	1,937,955,456✔
106	nu++;	121,148,676✔
107	}
108
109	// If no neutrons were produced then don't continue
110	if (nu == 0)	1,937,955,456✔
111	return;	1,816,805,988✔
112
113	// Initialize the counter of delayed neutrons encountered for each delayed
114	// group.
115	double nu_d[MAX_DELAYED_GROUPS] = {0.};	121,149,468✔
116
117	// Clear out particle's nu fission bank
118	p.nu_bank().clear();	121,149,468✔
119
120	p.fission() = true;	121,149,468✔
121
122	// Determine whether to place fission sites into the shared fission bank
123	// or the secondary particle bank.
124	bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);	121,149,468✔
125
126	// Counter for the number of fission sites successfully stored to the shared
127	// fission bank or the secondary particle bank
128	int n_sites_stored;
129
130	for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) {	242,300,544✔
131	// Initialize fission site object with particle data
132	SourceSite site;	121,151,076✔
133	site.r = p.r();	121,151,076✔
134	site.particle = ParticleType::neutron;	121,151,076✔
135	site.wgt = 1. / weight;	121,151,076✔
136	site.parent_id = p.id();	121,151,076✔
137	site.progeny_id = p.n_progeny()++;	121,151,076✔
138
139	// Sample the cosine of the angle, assuming fission neutrons are emitted
140	// isotropically
141	double mu = 2. * prn(p.current_seed()) - 1.;	121,151,076✔
142
143	// Sample the azimuthal angle uniformly in [0, 2.pi)
144	double phi = 2. * PI * prn(p.current_seed());	121,151,076✔
145	site.u.x = mu;	121,151,076✔
146	site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi);	121,151,076✔
147	site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi);	121,151,076✔
148
149	// Sample secondary energy distribution for the fission reaction
150	int dg;
151	int gout;
152	data::mg.macro_xs_[p.material()].sample_fission_energy(	363,453,228✔
153	p.g(), dg, gout, p.current_seed(), p.mg_xs_cache().t, p.mg_xs_cache().a);	242,302,152✔
154
155	// Store the energy and delayed groups on the fission bank
156	site.E = gout;	121,151,076✔
157
158	// We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest
159	// of the code, 0 is prompt.
160	site.delayed_group = dg + 1;	121,151,076✔
161
162	// Store fission site in bank
163	if (use_fission_bank) {	121,151,076✔
164	int64_t idx = simulation::fission_bank.thread_safe_append(site);	121,151,076✔
165	if (idx == -1) {	121,151,076✔
166	warning(	×
167	"The shared fission bank is full. Additional fission sites created "
168	"in this generation will not be banked. Results may be "
169	"non-deterministic.");
170
171	// Decrement number of particle progeny as storage was unsuccessful.
172	// This step is needed so that the sum of all progeny is equal to the
173	// size of the shared fission bank.
174	p.n_progeny()--;	×
175
176	// Break out of loop as no more sites can be added to fission bank
177	break;	×
178	}
179	} else {
180	p.secondary_bank().push_back(site);	×
181	}
182
183	// Set the delayed group on the particle as well
184	p.delayed_group() = dg + 1;	121,151,076✔
185
186	// Increment the number of neutrons born delayed
187	if (p.delayed_group() > 0) {	121,151,076✔
188	nu_d[dg]++;	1,680✔
189	}
190
191	// Write fission particles to nuBank
192	p.nu_bank().emplace_back();	121,151,076✔
193	NuBank* nu_bank_entry = &p.nu_bank().back();	121,151,076✔
194	nu_bank_entry->wgt = site.wgt;	121,151,076✔
195	nu_bank_entry->E = site.E;	121,151,076✔
196	nu_bank_entry->delayed_group = site.delayed_group;	121,151,076✔
197	}
198
199	// If shared fission bank was full, and no fissions could be added,
200	// set the particle fission flag to false.
201	if (n_sites_stored == 0) {	121,149,468✔
202	p.fission() = false;	×
203	return;	×
204	}
205
206	// Set nu to the number of fission sites successfully stored. If the fission
207	// bank was not found to be full then these values are already equivalent.
208	nu = n_sites_stored;	121,149,468✔
209
210	// Store the total weight banked for analog fission tallies
211	p.n_bank() = nu;	121,149,468✔
212	p.wgt_bank() = nu / weight;	121,149,468✔
213	for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {	1,090,345,212✔
214	p.n_delayed_bank(d) = nu_d[d];	969,195,744✔
215	}
216	}
217
218	void absorption(Particle& p)	1,939,604,940✔
219	{
220	if (settings::survival_biasing) {	1,939,604,940✔
221	// Determine weight absorbed in survival biasing
222	double wgt_absorb = p.wgt() * p.macro_xs().absorption / p.macro_xs().total;	4,596,048✔
223
224	// Adjust weight of particle by the probability of absorption
225	p.wgt() -= wgt_absorb;	4,596,048✔
226
227	// Score implicit absorpion estimate of keff
228	p.keff_tally_absorption() +=	4,596,048✔
229	wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption;	4,596,048✔
230	} else {
231	if (p.macro_xs().absorption > prn(p.current_seed()) * p.macro_xs().total) {	1,935,008,892✔
232	p.keff_tally_absorption() +=	241,771,080✔
233	p.wgt() * p.macro_xs().nu_fission / p.macro_xs().absorption;	120,885,540✔
234	p.wgt() = 0.0;	120,885,540✔
235	p.event() = TallyEvent::ABSORB;	120,885,540✔
236	}
237	}
238	}	1,939,604,940✔
239
240	} // namespace openmc

openmc-dev / openmc / 9789289404

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