Add option for survival biasing source normalization (#3070)

Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Perry 2025-03-06 16:00:29 -08:00 committed by GitHub
parent 9bfce4ee1c
commit e8c9134ff6
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
9 changed files with 76 additions and 26 deletions

View file

@ -81,6 +81,13 @@ time.
*Default*: 1.0
:survival_normalization:
If this element is set to "true", this will enable the use of survival
biasing source normalization, whereby the weight parameters, weight and
weight_avg, are multiplied per history by the start weight of said history.
*Default*: false
:energy_neutron:
The energy under which neutrons will be killed.

View file

@ -434,6 +434,7 @@ private:
int g_last_;
double wgt_ {1.0};
double wgt_born_ {1.0};
double mu_;
double time_ {0.0};
double time_last_ {0.0};
@ -532,12 +533,20 @@ public:
int& g_last() { return g_last_; }
const int& g_last() const { return g_last_; }
// Statistic weight of particle. Setting to zero
// indicates that the particle is dead.
// Statistic weight of particle. Setting to zero indicates that the particle
// is dead.
double& wgt() { return wgt_; }
double wgt() const { return wgt_; }
// Statistic weight of particle at birth
double& wgt_born() { return wgt_born_; }
double wgt_born() const { return wgt_born_; }
// Statistic weight of particle at last collision
double& wgt_last() { return wgt_last_; }
const double& wgt_last() const { return wgt_last_; }
// Whether particle is alive
bool alive() const { return wgt_ != 0.0; }
// Polar scattering angle after a collision

View file

@ -61,6 +61,7 @@ extern bool surf_source_write; //!< write surface source file?
extern bool surf_mcpl_write; //!< write surface mcpl file?
extern bool surf_source_read; //!< read surface source file?
extern bool survival_biasing; //!< use survival biasing?
extern bool survival_normalization; //!< use survival normalization?
extern bool temperature_multipole; //!< use multipole data?
extern "C" bool trigger_on; //!< tally triggers enabled?
extern bool trigger_predict; //!< predict batches for triggers?

View file

@ -50,15 +50,19 @@ class Settings:
Indicate whether fission neutrons should be created or not.
cutoff : dict
Dictionary defining weight cutoff, energy cutoff and time cutoff. The
dictionary may have ten keys, 'weight', 'weight_avg', 'energy_neutron',
'energy_photon', 'energy_electron', 'energy_positron', 'time_neutron',
'time_photon', 'time_electron', and 'time_positron'. Value for 'weight'
should be a float indicating weight cutoff below which particle undergo
Russian roulette. Value for 'weight_avg' should be a float indicating
weight assigned to particles that are not killed after Russian roulette.
Value of energy should be a float indicating energy in eV below which
particle type will be killed. Value of time should be a float in
seconds. Particles will be killed exactly at the specified time.
dictionary may have the following keys, 'weight', 'weight_avg',
'survival_normalization', 'energy_neutron', 'energy_photon',
'energy_electron', 'energy_positron', 'time_neutron', 'time_photon',
'time_electron', and 'time_positron'. Value for 'weight' should be a
float indicating weight cutoff below which particle undergo Russian
roulette. Value for 'weight_avg' should be a float indicating weight
assigned to particles that are not killed after Russian roulette. Value
of energy should be a float indicating energy in eV below which particle
type will be killed. Value of time should be a float in seconds.
Particles will be killed exactly at the specified time. Value for
'survival_normalization' is a bool indicating whether or not the weight
cutoff parameters will be applied relative to the particle's starting
weight or to its current weight.
delayed_photon_scaling : bool
Indicate whether to scale the fission photon yield by (EGP + EGD)/EGP
where EGP is the energy release of prompt photons and EGD is the energy
@ -162,20 +166,20 @@ class Settings:
'naive', 'simulation_averaged', or 'hybrid'.
The default is 'hybrid'.
:source_shape:
Assumed shape of the source distribution within each source
region. Options are 'flat' (default), 'linear', or 'linear_xy'.
Assumed shape of the source distribution within each source region.
Options are 'flat' (default), 'linear', or 'linear_xy'.
:volume_normalized_flux_tallies:
Whether to normalize flux tallies by volume (bool). The default
is 'False'. When enabled, flux tallies will be reported in units of
cm/cm^3. When disabled, flux tallies will be reported in units
of cm (i.e., total distance traveled by neutrons in the spatial
tally region).
Whether to normalize flux tallies by volume (bool). The default is
'False'. When enabled, flux tallies will be reported in units of
cm/cm^3. When disabled, flux tallies will be reported in units of cm
(i.e., total distance traveled by neutrons in the spatial tally
region).
:adjoint:
Whether to run the random ray solver in adjoint mode (bool). The
default is 'False'.
:sample_method:
Sampling method for the ray starting location and direction of travel.
Options are `prng` (default) or 'halton`.
Sampling method for the ray starting location and direction of
travel. Options are `prng` (default) or 'halton`.
.. versionadded:: 0.15.0
resonance_scattering : dict
@ -899,6 +903,8 @@ class Settings:
cv.check_type('average survival weight', cutoff[key], Real)
cv.check_greater_than('average survival weight',
cutoff[key], 0.0)
elif key == 'survival_normalization':
cv.check_type('survival normalization', cutoff[key], bool)
elif key in ['energy_neutron', 'energy_photon', 'energy_electron',
'energy_positron']:
cv.check_type('energy cutoff', cutoff[key], Real)
@ -1364,7 +1370,8 @@ class Settings:
element = ET.SubElement(root, "cutoff")
for key, value in self._cutoff.items():
subelement = ET.SubElement(element, key)
subelement.text = str(value)
subelement.text = str(value) if key != 'survival_normalization' \
else str(value).lower()
def _create_entropy_mesh_subelement(self, root, mesh_memo=None):
if self.entropy_mesh is None:
@ -1797,10 +1804,14 @@ class Settings:
self.cutoff = {}
for key in ('energy_neutron', 'energy_photon', 'energy_electron',
'energy_positron', 'weight', 'weight_avg', 'time_neutron',
'time_photon', 'time_electron', 'time_positron'):
'time_photon', 'time_electron', 'time_positron',
'survival_normalization'):
value = get_text(elem, key)
if value is not None:
self.cutoff[key] = float(value)
if key == 'survival_normalization':
self.cutoff[key] = value in ('true', '1')
else:
self.cutoff[key] = float(value)
def _entropy_mesh_from_xml_element(self, root, meshes):
text = get_text(root, 'entropy_mesh')

View file

@ -153,7 +153,13 @@ void sample_neutron_reaction(Particle& p)
// Play russian roulette if survival biasing is turned on
if (settings::survival_biasing) {
if (p.wgt() < settings::weight_cutoff) {
// if survival normalization is on, use normalized weight cutoff and
// normalized weight survive
if (settings::survival_normalization) {
if (p.wgt() < settings::weight_cutoff * p.wgt_born()) {
russian_roulette(p, settings::weight_survive * p.wgt_born());
}
} else if (p.wgt() < settings::weight_cutoff) {
russian_roulette(p, settings::weight_survive);
}
}

View file

@ -68,7 +68,13 @@ void sample_reaction(Particle& p)
// Play Russian roulette if survival biasing is turned on
if (settings::survival_biasing) {
if (p.wgt() < settings::weight_cutoff) {
// if survival normalization is applicable, use normalized weight cutoff and
// normalized weight survive
if (settings::survival_normalization) {
if (p.wgt() < settings::weight_cutoff * p.wgt_born()) {
russian_roulette(p, settings::weight_survive * p.wgt_born());
}
} else if (p.wgt() < settings::weight_cutoff) {
russian_roulette(p, settings::weight_survive);
}
}

View file

@ -69,6 +69,7 @@ bool surf_source_write {false};
bool surf_mcpl_write {false};
bool surf_source_read {false};
bool survival_biasing {false};
bool survival_normalization {false};
bool temperature_multipole {false};
bool trigger_on {false};
bool trigger_predict {false};
@ -624,6 +625,10 @@ void read_settings_xml(pugi::xml_node root)
if (check_for_node(node_cutoff, "weight_avg")) {
weight_survive = std::stod(get_node_value(node_cutoff, "weight_avg"));
}
if (check_for_node(node_cutoff, "survival_normalization")) {
survival_normalization =
get_node_value_bool(node_cutoff, "survival_normalization");
}
if (check_for_node(node_cutoff, "energy_neutron")) {
energy_cutoff[0] =
std::stod(get_node_value(node_cutoff, "energy_neutron"));

View file

@ -586,6 +586,9 @@ void initialize_history(Particle& p, int64_t index_source)
// Reset weight window ratio
p.ww_factor() = 0.0;
// set particle history start weight
p.wgt_born() = p.wgt();
// Reset pulse_height_storage
std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0);
@ -610,7 +613,7 @@ void initialize_history(Particle& p, int64_t index_source)
write_message("Simulating Particle {}", p.id());
}
// Add paricle's starting weight to count for normalizing tallies later
// Add particle's starting weight to count for normalizing tallies later
#pragma omp atomic
simulation::total_weight += p.wgt();

View file

@ -26,6 +26,7 @@ def test_export_to_xml(run_in_tmpdir):
s.plot_seed = 100
s.survival_biasing = True
s.cutoff = {'weight': 0.25, 'weight_avg': 0.5, 'energy_neutron': 1.0e-5,
'survival_normalization': True,
'energy_photon': 1000.0, 'energy_electron': 1.0e-5,
'energy_positron': 1.0e-5, 'time_neutron': 1.0e-5,
'time_photon': 1.0e-5, 'time_electron': 1.0e-5,
@ -99,6 +100,7 @@ def test_export_to_xml(run_in_tmpdir):
assert s.seed == 17
assert s.survival_biasing
assert s.cutoff == {'weight': 0.25, 'weight_avg': 0.5,
'survival_normalization': True,
'energy_neutron': 1.0e-5, 'energy_photon': 1000.0,
'energy_electron': 1.0e-5, 'energy_positron': 1.0e-5,
'time_neutron': 1.0e-5, 'time_photon': 1.0e-5,