diff --git a/include/openmc/physics_common.h b/include/openmc/physics_common.h index e38a3c7f8..b0e395c1e 100644 --- a/include/openmc/physics_common.h +++ b/include/openmc/physics_common.h @@ -13,5 +13,9 @@ namespace openmc { //! \param[in] weight_survive Weight assigned to particles that survive void russian_roulette(Particle& p, double weight_survive); +//! \brief Performs the global russian roulette operation on a particle +//! \param[in,out] p Particle object +void apply_russian_roulette(Particle& p); + } // namespace openmc #endif // OPENMC_PHYSICS_COMMON_H diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 97abcfe82..42846f9d1 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -25,17 +25,6 @@ enum class WeightWindowUpdateMethod { MAGIC, FW_CADIS }; constexpr double DEFAULT_WEIGHT_CUTOFF {1.0e-38}; // default low weight cutoff -//============================================================================== -// Non-member functions -//============================================================================== - -//! Apply weight windows to a particle -//! \param[in] p Particle to apply weight windows to -void apply_weight_windows(Particle& p); - -//! Free memory associated with weight windows -void free_memory_weight_windows(); - //============================================================================== // Global variables //============================================================================== @@ -137,7 +126,7 @@ public: //! Retrieve the weight window for a particle //! \param[in] p Particle to get weight window for - WeightWindow get_weight_window(const Particle& p) const; + std::pair get_weight_window(const Particle& p) const; std::array bounds_size() const; @@ -236,6 +225,26 @@ public: double ratio_ {5.0}; // search_weight_window(const Particle& p); + //! Finalize variance reduction objects after all inputs have been read void finalize_variance_reduction(); diff --git a/src/physics.cpp b/src/physics.cpp index c9737205c..72b4a6611 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -64,8 +64,17 @@ void collision(Particle& p) fatal_error("Unsupported particle PDG for collision sampling."); } - if (settings::weight_window_checkpoint_collision) - apply_weight_windows(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 for neutrons + // (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); + } + } // Kill particle if energy falls below cutoff int type = p.type().transport_index(); @@ -156,18 +165,9 @@ void sample_neutron_reaction(Particle& p) advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE)); } - // Play russian roulette if survival biasing is turned on - if (settings::survival_biasing) { - // 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); - } - } + // Play russian roulette if there are no weight windows + if (!settings::weight_windows_on) + apply_russian_roulette(p); } void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) diff --git a/src/physics_common.cpp b/src/physics_common.cpp index e25ae6b97..10760ce2d 100644 --- a/src/physics_common.cpp +++ b/src/physics_common.cpp @@ -18,4 +18,20 @@ void russian_roulette(Particle& p, double weight_survive) } } +void apply_russian_roulette(Particle& p) +{ + // Exit if survival biasing is turned off + if (!settings::survival_biasing) + return; + + // 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); + } +} } // namespace openmc diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 32ae10a60..6991c7365 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -32,8 +32,17 @@ void collision_mg(Particle& p) // Sample the reaction type sample_reaction(p); - if (settings::weight_window_checkpoint_collision) - apply_weight_windows(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()) { @@ -67,18 +76,9 @@ void sample_reaction(Particle& p) // 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.wgt_born()) { - russian_roulette(p, settings::weight_survive * p.wgt_born()); - } - } else if (p.wgt() < settings::weight_cutoff) { - russian_roulette(p, settings::weight_survive); - } - } + // Play russian roulette if there are no weight windows + if (!settings::weight_windows_on) + apply_russian_roulette(p); } void scatter(Particle& p) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index b881ced71..1a6e7c0be 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -843,7 +843,7 @@ void FlatSourceDomain::output_to_vtk() const voxel_positions[z * Ny * Nx + y * Nx + x] = sample; if (variance_reduction::weight_windows.size() == 1) { - WeightWindow ww = + auto [ww_found, ww] = variance_reduction::weight_windows[0]->get_weight_window(p); float weight = ww.lower_weight; weight_windows[z * Ny * Nx + y * Nx + x] = weight; diff --git a/src/settings.cpp b/src/settings.cpp index b6a9ab953..ee8edd4a5 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1266,6 +1266,13 @@ void read_settings_xml(pugi::xml_node root) } } + if (weight_windows_on) { + if (!weight_window_checkpoint_surface && + !weight_window_checkpoint_collision) + fatal_error( + "Weight Windows are enabled but there are no valid checkpoints."); + } + if (check_for_node(root, "use_decay_photons")) { settings::use_decay_photons = get_node_value_bool(root, "use_decay_photons"); diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index c05169132..c00872e56 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -44,108 +44,6 @@ openmc::vector> weight_windows_generators; } // namespace variance_reduction -//============================================================================== -// Non-member functions -//============================================================================== - -void apply_weight_windows(Particle& p) -{ - if (!settings::weight_windows_on) - return; - - // WW on photon and neutron only - if (!p.type().is_neutron() && !p.type().is_photon()) - return; - - // skip dead or no energy - if (p.E() <= 0 || !p.alive()) - return; - - bool in_domain = false; - // TODO: this is a linear search - should do something more clever - WeightWindow weight_window; - for (const auto& ww : variance_reduction::weight_windows) { - weight_window = ww->get_weight_window(p); - if (weight_window.is_valid()) - break; - } - - // If particle has not yet had its birth weight window value set, set it to - // the current weight window (or 1.0 if not born in a weight window). - if (p.wgt_ww_born() == -1.0) { - if (weight_window.is_valid()) { - p.wgt_ww_born() = - (weight_window.lower_weight + weight_window.upper_weight) / 2; - } else { - p.wgt_ww_born() = 1.0; - } - } - - // particle is not in any of the ww domains, do nothing - if (!weight_window.is_valid()) - return; - - // Normalize weight windows based on particle's starting weight - // and the value of the weight window the particle was born in. - weight_window.scale(p.wgt_born() / p.wgt_ww_born()); - - // get the paramters - double weight = p.wgt(); - - // first check to see if particle should be killed for weight cutoff - if (p.wgt() < weight_window.weight_cutoff) { - p.wgt() = 0.0; - return; - } - - // check if particle is far above current weight window - // only do this if the factor is not already set on the particle and a - // maximum lower bound ratio is specified - if (p.ww_factor() == 0.0 && weight_window.max_lb_ratio > 1.0 && - p.wgt() > weight_window.lower_weight * weight_window.max_lb_ratio) { - p.ww_factor() = - p.wgt() / (weight_window.lower_weight * weight_window.max_lb_ratio); - } - - // move weight window closer to the particle weight if needed - if (p.ww_factor() > 1.0) - weight_window.scale(p.ww_factor()); - - // if particle's weight is above the weight window split until they are within - // the window - if (weight > weight_window.upper_weight) { - // do not further split the particle if above the limit - if (p.n_split() >= settings::max_history_splits) - return; - - double n_split = std::ceil(weight / weight_window.upper_weight); - double max_split = weight_window.max_split; - n_split = std::min(n_split, max_split); - - p.n_split() += n_split; - - // Create secondaries and divide weight among all particles - int i_split = std::round(n_split); - for (int l = 0; l < i_split - 1; l++) { - p.split(weight / n_split); - } - // remaining weight is applied to current particle - p.wgt() = weight / n_split; - - } else if (weight <= weight_window.lower_weight) { - // if the particle weight is below the window, play Russian roulette - double weight_survive = - std::min(weight * weight_window.max_split, weight_window.survival_weight); - russian_roulette(p, weight_survive); - } // else particle is in the window, continue as normal -} - -void free_memory_weight_windows() -{ - variance_reduction::ww_map.clear(); - variance_reduction::weight_windows.clear(); -} - //============================================================================== // WeightWindowSettings implementation //============================================================================== @@ -375,27 +273,28 @@ void WeightWindows::set_mesh(const Mesh* mesh) set_mesh(model::mesh_map[mesh->id_]); } -WeightWindow WeightWindows::get_weight_window(const Particle& p) const +std::pair WeightWindows::get_weight_window( + const Particle& p) const { // check for particle type if (particle_type_ != p.type()) { - return {}; + return {false, {}}; } + // particle energy + double E = p.E(); + + // check to make sure energy is in range, expects sorted energy values + if (E < energy_bounds_.front() || E > energy_bounds_.back()) + return {false, {}}; + // Get mesh index for particle's position const auto& mesh = this->mesh(); int mesh_bin = mesh->get_bin(p.r()); // particle is outside the weight window mesh if (mesh_bin < 0) - return {}; - - // particle energy - double E = p.E(); - - // check to make sure energy is in range, expects sorted energy values - if (E < energy_bounds_.front() || E > energy_bounds_.back()) - return {}; + return {false, {}}; // get the mesh bin in energy group int energy_bin = @@ -410,7 +309,7 @@ WeightWindow WeightWindows::get_weight_window(const Particle& p) const ww.max_lb_ratio = max_lb_ratio_; ww.max_split = max_split_; ww.weight_cutoff = weight_cutoff_; - return ww; + return {true, ww}; } std::array WeightWindows::bounds_size() const @@ -1023,6 +922,115 @@ void WeightWindowsGenerator::update() const // Non-member functions //============================================================================== +std::pair search_weight_window(const Particle& p) +{ + // TODO: this is a linear search - should do something more clever + for (const auto& ww : variance_reduction::weight_windows) { + auto [ww_found, weight_window] = ww->get_weight_window(p); + if (ww_found) + return {true, weight_window}; + } + return {false, {}}; +} + +void apply_weight_windows(Particle& p) +{ + if (!settings::weight_windows_on) + return; + + // WW on photon and neutron only + if (!p.type().is_neutron() && !p.type().is_photon()) + return; + + // skip dead or no energy + if (p.E() <= 0 || !p.alive()) + return; + + auto [ww_found, ww] = search_weight_window(p); + if (ww_found && ww.is_valid()) { + apply_weight_window(p, ww); + } else { + if (p.wgt_ww_born() == -1.0) + p.wgt_ww_born() = 1.0; + } +} + +void apply_weight_window(Particle& p, WeightWindow weight_window) +{ + if (!weight_window.is_valid()) + return; + + // skip dead or no energy + if (p.E() <= 0 || !p.alive()) + return; + + // If particle has not yet had its birth weight window value set, set it to + // the current weight window. + if (p.wgt_ww_born() == -1.0) + p.wgt_ww_born() = + (weight_window.lower_weight + weight_window.upper_weight) / 2; + + // Normalize weight windows based on particle's starting weight + // and the value of the weight window the particle was born in. + weight_window.scale(p.wgt_born() / p.wgt_ww_born()); + + // get the paramters + double weight = p.wgt(); + + // first check to see if particle should be killed for weight cutoff + if (p.wgt() < weight_window.weight_cutoff) { + p.wgt() = 0.0; + return; + } + + // check if particle is far above current weight window + // only do this if the factor is not already set on the particle and a + // maximum lower bound ratio is specified + if (p.ww_factor() == 0.0 && weight_window.max_lb_ratio > 1.0 && + p.wgt() > weight_window.lower_weight * weight_window.max_lb_ratio) { + p.ww_factor() = + p.wgt() / (weight_window.lower_weight * weight_window.max_lb_ratio); + } + + // move weight window closer to the particle weight if needed + if (p.ww_factor() > 1.0) + weight_window.scale(p.ww_factor()); + + // if particle's weight is above the weight window split until they are within + // the window + if (weight > weight_window.upper_weight) { + // do not further split the particle if above the limit + if (p.n_split() >= settings::max_history_splits) + return; + + double n_split = std::ceil(weight / weight_window.upper_weight); + double max_split = weight_window.max_split; + n_split = std::min(n_split, max_split); + + p.n_split() += n_split; + + // Create secondaries and divide weight among all particles + int i_split = std::round(n_split); + for (int l = 0; l < i_split - 1; l++) { + p.split(weight / n_split); + } + // remaining weight is applied to current particle + p.wgt() = weight / n_split; + + } else if (weight <= weight_window.lower_weight) { + // if the particle weight is below the window, play Russian roulette + double weight_survive = + std::min(weight * weight_window.max_split, weight_window.survival_weight); + russian_roulette(p, weight_survive); + } // else particle is in the window, continue as normal +} + +void free_memory_weight_windows() +{ + variance_reduction::ww_map.clear(); + variance_reduction::weight_windows.clear(); +} + void finalize_variance_reduction() { for (const auto& wwg : variance_reduction::weight_windows_generators) { diff --git a/tests/regression_tests/weightwindows/survival_biasing/__init__.py b/tests/regression_tests/weightwindows/survival_biasing/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/weightwindows/survival_biasing/inputs_true.dat b/tests/regression_tests/weightwindows/survival_biasing/inputs_true.dat new file mode 100644 index 000000000..d5aa135d4 --- /dev/null +++ b/tests/regression_tests/weightwindows/survival_biasing/inputs_true.dat @@ -0,0 +1,68 @@ + + + + + + + + + + + + + + + + + + + + + + + fixed source + 50 + 5 + + + + 0.01 1.0 + + + + + + + 14100000.0 1.0 + + + true + + 1 + neutron + 0.46135961568957096 0.2874485390202603 0.13256013950494605 0.052765209609807295 0.020958440832685638 0.006317250035749876 0.002627037175682506 0.00030032343592437284 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4611178493390162 0.28624862560842024 0.13999870622621136 0.05508479408959756 0.018281241958254993 0.0055577764872361555 0.0015265432226293397 0.00024298772352636966 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.46294957913063317 0.2857734367546051 0.13365623190336945 0.05034742166227103 0.017525851812913266 0.005096725451851077 0.0026297640951531403 0.00033732424392026144 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.45772598750472554 0.281494921471495 0.13604397962762246 0.054792881472295364 0.019802376898755525 0.0073351745579128495 0.002067392946066426 9.529474085926143e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4739762983490318 0.28419061537046764 0.12757509901507888 0.0516920293019733 0.01919167874787235 0.007593035964707848 0.0017488176314107277 9.678267909681988e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.47560199098558276 0.27728389146956994 0.12760802572535623 0.05251965811713552 0.022498960644951632 0.01063372459325151 0.004242071054914633 0.00045603112202665183 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.47659199471934693 0.28886666944252215 0.1258637573398161 0.05818862375955657 0.020476313720744762 0.007143193244018263 0.0022364083022515273 0.0003953123781408474 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4693487369606823 0.27180862139490125 0.12833455570423483 0.055146743559938344 0.020667403529470333 0.007891508723137146 0.001643551705407358 0.0006501416781353278 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4538724931023337 0.2824491421084296 0.1275341706351566 0.055027501141893705 0.02305114640508087 0.008599303158607198 0.0025082611194161804 0.000518530311231696 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.45258444753220123 0.26839102466484255 0.12991633918797083 0.05442263709947767 0.019695895223471486 0.005353204961887518 0.0015074698293840157 0.0001598680898180058 8.872047880811405e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4482467985982804 0.2813276470183359 0.13449860790768647 0.056103278190533436 0.026274320334196962 0.006884071908429303 0.0033099390738735458 0.0005337454397674922 0.0001386380399465575 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4875748950139487 0.2908586705316248 0.14212060658392278 0.06068108009637272 0.019467970468789477 0.00637102427946399 0.0028510425266191687 0.0012184421778115488 0.00029404541567146187 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.5 0.2852101551297824 0.13622307172321643 0.058302519884339946 0.023124734915152625 0.007401527839638796 0.002527635652743992 0.0008509612315162482 0.00018213334574554768 7.815977984795461e-06 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.47706397688984176 0.28313436859495966 0.13170624744221976 0.0525453512766549 0.021111753158319105 0.0067451713955609645 0.003204270539718037 0.0010308107242263192 0.0002759210144794204 0.00013727805365387318 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.48641988659284513 0.2816173231471242 0.131541850061199 0.054793200248327956 0.016517701691156576 0.006757591781856257 0.0025425350750908644 0.0006383278085839443 2.279421641064226e-05 0.0003597432472224371 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4307755159245372 0.2652713784080849 0.12753857957535653 0.05470396852881577 0.02049817309420563 0.00566559741247352 0.0007213846765465147 6.848024591311009e-05 -1.0 8.609308814924014e-05 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.43631324024956025 0.2639077825530567 0.13368737962742414 0.05213531603720763 0.020628860469743833 0.008268384844678654 0.0028599221013934375 0.00013059757129982714 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.45926287238219 0.2635824148883888 0.13104280996799478 0.055379812186041905 0.019255042561941074 0.007252095690246872 0.0018026593488140996 0.0003215836458542421 9.14840176000331e-06 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4574378998400604 0.28714971179927723 0.1363507911051337 0.05010315954654347 0.017963278989571074 0.0058998176297971605 0.0021393135666168 0.00015677575033083482 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 0.4821442782978812 0.27704132233083506 0.14149074035662307 0.05546042815834048 0.016800649382269998 0.004096804603054233 0.0019161108398578026 0.0005289600253155307 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 + 2.306798078447855 1.4372426951013015 0.6628006975247303 0.2638260480490365 0.10479220416342819 0.03158625017874938 0.013135185878412529 0.0015016171796218643 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.3055892466950807 1.4312431280421012 0.6999935311310568 0.2754239704479878 0.09140620979127496 0.027788882436180776 0.007632716113146698 0.0012149386176318483 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.314747895653166 1.4288671837730256 0.6682811595168472 0.25173710831135515 0.08762925906456634 0.025483627259255386 0.013148820475765701 0.0016866212196013073 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.288629937523628 1.407474607357475 0.6802198981381123 0.2739644073614768 0.09901188449377762 0.036675872789564246 0.01033696473033213 0.00047647370429630714 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.369881491745159 1.4209530768523382 0.6378754950753944 0.2584601465098665 0.09595839373936176 0.03796517982353924 0.008744088157053638 0.00048391339548409945 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.3780099549279137 1.3864194573478497 0.6380401286267812 0.26259829058567763 0.11249480322475816 0.053168622966257545 0.021210355274573163 0.0022801556101332593 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.382959973596735 1.4443333472126108 0.6293187866990806 0.29094311879778284 0.10238156860372381 0.035715966220091315 0.011182041511257637 0.001976561890704237 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.3467436848034113 1.3590431069745064 0.6416727785211741 0.2757337177996917 0.10333701764735166 0.03945754361568573 0.00821775852703679 0.0032507083906766388 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.2693624655116684 1.412245710542148 0.637670853175783 0.27513750570946854 0.11525573202540434 0.04299651579303599 0.012541305597080901 0.00259265155615848 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.262922237661006 1.3419551233242126 0.6495816959398542 0.2721131854973884 0.09847947611735743 0.02676602480943759 0.007537349146920078 0.000799340449090029 0.0004436023940405703 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.2412339929914022 1.4066382350916795 0.6724930395384323 0.28051639095266717 0.1313716016709848 0.03442035954214652 0.016549695369367727 0.0026687271988374613 0.0006931901997327875 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.4378744750697434 1.4542933526581239 0.7106030329196139 0.3034054004818636 0.09733985234394739 0.03185512139731995 0.014255212633095843 0.006092210889057744 0.0014702270783573093 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.5 1.4260507756489118 0.6811153586160822 0.2915125994216997 0.11562367457576313 0.03700763919819398 0.012638178263719959 0.004254806157581241 0.0009106667287277384 3.9079889923977307e-05 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.385319884449209 1.4156718429747983 0.6585312372110987 0.2627267563832745 0.10555876579159552 0.03372585697780482 0.016021352698590185 0.005154053621131596 0.001379605072397102 0.000686390268269366 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.432099432964226 1.4080866157356209 0.657709250305995 0.2739660012416398 0.08258850845578287 0.03378795890928128 0.012712675375454322 0.0031916390429197216 0.0001139710820532113 0.0017987162361121855 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.1538775796226863 1.3263568920404245 0.6376928978767826 0.27351984264407886 0.10249086547102815 0.028327987062367603 0.0036069233827325737 0.0003424012295655504 -5.0 0.0004304654407462007 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.181566201247801 1.3195389127652835 0.6684368981371207 0.2606765801860381 0.10314430234871916 0.041341924223393264 0.014299610506967188 0.0006529878564991358 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.29631436191095 1.317912074441944 0.6552140498399739 0.27689906093020955 0.09627521280970537 0.03626047845123436 0.009013296744070498 0.0016079182292712106 4.574200880001655e-05 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.287189499200302 1.435748558996386 0.6817539555256685 0.25051579773271737 0.08981639494785537 0.029499088148985803 0.010696567833083998 0.0007838787516541741 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 2.410721391489406 1.3852066116541752 0.7074537017831153 0.2773021407917024 0.08400324691134999 0.020484023015271167 0.009580554199289014 0.0026448001265776534 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 -5.0 + 3.0 + 10 + 1e-38 + + + 20 20 1 + 0.0 0.0 0.0 + 160.0 160.0 160.0 + + + true + true + + + + + 1 + + + 1 + flux + + + diff --git a/tests/regression_tests/weightwindows/survival_biasing/results_true.dat b/tests/regression_tests/weightwindows/survival_biasing/results_true.dat new file mode 100644 index 000000000..b457a245c --- /dev/null +++ b/tests/regression_tests/weightwindows/survival_biasing/results_true.dat @@ -0,0 +1 @@ +386e507008ed3c72c6e1a101aafc01cfaaff2c0b10555558d1eab6332441f7b17264b55002d721151bac2f3e7c1a8aac4c5ed243f5270533d171850f985206ed \ No newline at end of file diff --git a/tests/regression_tests/weightwindows/survival_biasing/test.py b/tests/regression_tests/weightwindows/survival_biasing/test.py new file mode 100644 index 000000000..92604e338 --- /dev/null +++ b/tests/regression_tests/weightwindows/survival_biasing/test.py @@ -0,0 +1,81 @@ +import pytest +import numpy as np + +import openmc +from openmc.stats import Discrete, Point + +from tests.testing_harness import HashedPyAPITestHarness + + +@pytest.fixture +def model(): + # Material + w = openmc.Material(name='Tungsten') + w.add_element('W', 1.0) + w.set_density('g/cm3', 19.25) + + materials = openmc.Materials([w]) + + # Geometry surfaces + x0 = openmc.XPlane(x0=0.0, boundary_type='reflective') + x1 = openmc.XPlane(x0=160.0, boundary_type='vacuum') + y0 = openmc.YPlane(y0=0.0, boundary_type='reflective') + y1 = openmc.YPlane(y0=160.0, boundary_type='reflective') + z0 = openmc.ZPlane(z0=0.0, boundary_type='reflective') + z1 = openmc.ZPlane(z0=160.0, boundary_type='reflective') + + region = +x0 & -x1 & +y0 & -y1 & +z0 & -z1 + cell = openmc.Cell(region=region, fill=w) + root = openmc.Universe(cells=[cell]) + geometry = openmc.Geometry(root) + + # Source: planar on x=0, mono-directional along +x, 14.1 MeV neutrons + space = openmc.stats.CartesianIndependent( + openmc.stats.Discrete([0.01], [1.0]), + openmc.stats.Uniform(0.0, 160.0), + openmc.stats.Uniform(0.0, 160.0), + ) + angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0)) + energy = openmc.stats.Discrete([14.1e6], [1.0]) + + source = openmc.Source(space=space, angle=angle, energy=energy) + + settings = openmc.Settings() + settings.run_mode = 'fixed source' + settings.batches = 5 + settings.particles = 50 + settings.source = source + + model = openmc.Model(geometry=geometry, materials=materials, settings=settings) + + # Mesh tally: 1 cm voxels, flux only + mesh = openmc.RegularMesh() + mesh.dimension = (20, 20, 1) + mesh.lower_left = (0.0, 0.0, 0.0) + mesh.upper_right = (160.0, 160.0, 160.0) + + mesh_filter = openmc.MeshFilter(mesh) + flux_tally = openmc.Tally(name='flux') + flux_tally.filters = [mesh_filter] + flux_tally.scores = ['flux'] + tallies = openmc.Tallies([flux_tally]) + model.tallies = tallies + + lower_ww_bounds = np.loadtxt('ww_n.txt') + + weight_windows = openmc.WeightWindows(mesh, + lower_ww_bounds, + upper_bound_ratio=5.0, + particle_type='neutron') + + model.settings.weight_windows = weight_windows + model.settings.weight_window_checkpoints = {'surface': True, + 'collision': True} + model.settings.survival_biasing = True + + return model + + +def test_weight_windows_with_survival_biasing(model): + harness = HashedPyAPITestHarness('statepoint.5.h5', model) + harness.main() diff --git a/tests/regression_tests/weightwindows/survival_biasing/ww_n.txt b/tests/regression_tests/weightwindows/survival_biasing/ww_n.txt new file mode 100644 index 000000000..aea1ff132 --- /dev/null +++ b/tests/regression_tests/weightwindows/survival_biasing/ww_n.txt @@ -0,0 +1,400 @@ +4.613596156895709566e-01 +4.611178493390161726e-01 +4.629495791306331709e-01 +4.577259875047255400e-01 +4.739762983490318216e-01 +4.756019909855827565e-01 +4.765919947193469342e-01 +4.693487369606823001e-01 +4.538724931023336850e-01 +4.525844475322012284e-01 +4.482467985982804271e-01 +4.875748950139486837e-01 +5.000000000000000000e-01 +4.770639768898417565e-01 +4.864198865928451299e-01 +4.307755159245372223e-01 +4.363132402495602524e-01 +4.592628723821899905e-01 +4.574378998400603913e-01 +4.821442782978812014e-01 +2.874485390202602964e-01 +2.862486256084202374e-01 +2.857734367546050924e-01 +2.814949214714950187e-01 +2.841906153704676363e-01 +2.772838914695699430e-01 +2.888666694425221504e-01 +2.718086213949012508e-01 +2.824491421084295850e-01 +2.683910246648425479e-01 +2.813276470183359024e-01 +2.908586705316247856e-01 +2.852101551297823723e-01 +2.831343685949596622e-01 +2.816173231471241767e-01 +2.652713784080849013e-01 +2.639077825530566912e-01 +2.635824148883887941e-01 +2.871497117992772297e-01 +2.770413223308350603e-01 +1.325601395049460507e-01 +1.399987062262113557e-01 +1.336562319033694490e-01 +1.360439796276224633e-01 +1.275750990150788799e-01 +1.276080257253562333e-01 +1.258637573398161125e-01 +1.283345557042348262e-01 +1.275341706351565962e-01 +1.299163391879708251e-01 +1.344986079076864738e-01 +1.421206065839227817e-01 +1.362230717232164323e-01 +1.317062474422197593e-01 +1.315418500611990060e-01 +1.275385795753565255e-01 +1.336873796274241355e-01 +1.310428099679947778e-01 +1.363507911051337063e-01 +1.414907403566230681e-01 +5.276520960980729535e-02 +5.508479408959755796e-02 +5.034742166227103299e-02 +5.479288147229536415e-02 +5.169202930197330098e-02 +5.251965811713552035e-02 +5.818862375955657223e-02 +5.514674355993834376e-02 +5.502750114189370462e-02 +5.442263709947767203e-02 +5.610327819053343573e-02 +6.068108009637272066e-02 +5.830251988433994559e-02 +5.254535127665489747e-02 +5.479320024832795566e-02 +5.470396852881576760e-02 +5.213531603720762686e-02 +5.537981218604190459e-02 +5.010315954654347148e-02 +5.546042815834047873e-02 +2.095844083268563751e-02 +1.828124195825499287e-02 +1.752585181291326649e-02 +1.980237689875552487e-02 +1.919167874787235106e-02 +2.249896064495163217e-02 +2.047631372074476194e-02 +2.066740352947033302e-02 +2.305114640508086968e-02 +1.969589522347148583e-02 +2.627432033419696208e-02 +1.946797046878947710e-02 +2.312473491515262478e-02 +2.111175315831910482e-02 +1.651770169115657563e-02 +2.049817309420563088e-02 +2.062886046974383297e-02 +1.925504256194107353e-02 +1.796327898957107358e-02 +1.680064938226999774e-02 +6.317250035749876098e-03 +5.557776487236155451e-03 +5.096725451851077254e-03 +7.335174557912849461e-03 +7.593035964707848043e-03 +1.063372459325150933e-02 +7.143193244018263000e-03 +7.891508723137145853e-03 +8.599303158607197670e-03 +5.353204961887517849e-03 +6.884071908429303423e-03 +6.371024279463989581e-03 +7.401527839638795750e-03 +6.745171395560964518e-03 +6.757591781856257113e-03 +5.665597412473520264e-03 +8.268384844678653561e-03 +7.252095690246871881e-03 +5.899817629797160512e-03 +4.096804603054233183e-03 +2.627037175682505905e-03 +1.526543222629339657e-03 +2.629764095153140288e-03 +2.067392946066426082e-03 +1.748817631410727689e-03 +4.242071054914632773e-03 +2.236408302251527251e-03 +1.643551705407358000e-03 +2.508261119416180414e-03 +1.507469829384015672e-03 +3.309939073873545759e-03 +2.851042526619168658e-03 +2.527635652743991969e-03 +3.204270539718037138e-03 +2.542535075090864363e-03 +7.213846765465147101e-04 +2.859922101393437468e-03 +1.802659348814099642e-03 +2.139313566616799812e-03 +1.916110839857802610e-03 +3.003234359243728432e-04 +2.429877235263696591e-04 +3.373242439202614419e-04 +9.529474085926143009e-05 +9.678267909681988429e-05 +4.560311220266518293e-04 +3.953123781408473961e-04 +6.501416781353277531e-04 +5.185303112316960250e-04 +1.598680898180058108e-04 +5.337454397674922446e-04 +1.218442177811548824e-03 +8.509612315162482432e-04 +1.030810724226319166e-03 +6.383278085839443235e-04 +6.848024591311008757e-05 +1.305975712998271444e-04 +3.215836458542421176e-04 +1.567757503308348204e-04 +5.289600253155306818e-04 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +8.872047880811405188e-05 +1.386380399465574921e-04 +2.940454156714618719e-04 +1.821333457455476831e-04 +2.759210144794204192e-04 +2.279421641064226044e-05 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +9.148401760003310179e-06 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +7.815977984795461318e-06 +1.372780536538731832e-04 +3.597432472224371168e-04 +8.609308814924013556e-05 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00 +-1.000000000000000000e+00