diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 5174cfd8b..0d5367727 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -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. diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 6b318385e..15ae57893 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -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 diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 12835fdef..f3aff08b0 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -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? diff --git a/openmc/settings.py b/openmc/settings.py index 882a17b68..2f054ebd1 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -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') diff --git a/src/physics.cpp b/src/physics.cpp index e04f54c0b..72c04a5ca 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -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); } } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index e381b27c9..c97a3d6f3 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -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); } } diff --git a/src/settings.cpp b/src/settings.cpp index cd925be70..67fdfea66 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -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")); diff --git a/src/simulation.cpp b/src/simulation.cpp index a3c6495ae..74a7dbe63 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -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(); diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index 02a476251..397d70434 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -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,