diff --git a/docs/source/pythonapi/examples.rst b/docs/source/pythonapi/examples.rst index e7ef523c0e3..9efbd29c710 100644 --- a/docs/source/pythonapi/examples.rst +++ b/docs/source/pythonapi/examples.rst @@ -11,6 +11,7 @@ Simple Models :template: myfunction.rst openmc.examples.slab_mg + openmc.examples.sphere_with_shielded_pocket Reactor Models -------------- diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index 06d4e50ed7f..adf5227f5f4 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -672,7 +672,9 @@ of these methods is given below: both spatial and resonance self shielding effects - * Potentially slower as the full geometry must be run * If a material is only present far from the source and doesn't get tallied - to in the CE simulation, the MGXS will be zero for that material. + to in the CE simulation, the MGXS will be zero for that material. This + can be mitigated by supplying weight windows via the + ``weight_windows_file`` argument (see :ref:`mgxs_bootstrap`). * - ``stochastic_slab`` - * Medium Fidelity * Runs a CE simulation with a greatly simplified geometry, where materials @@ -757,6 +759,53 @@ simulation, and if more fidelity is needed the user may wish to follow the instructions below or experiment with transport correction techniques to improve the fidelity of the generated MGXS data. +.. _mgxs_bootstrap: + +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +Bootstrapping Material-Wise MGXS with Weight Windows +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The ``"material_wise"`` method runs a continuous energy simulation of the +original geometry, so it produces the highest fidelity cross sections of the +three methods. However, it has a notable weakness: if a material only appears +far from the source (for example, a detector or structural material located +outside a thick shield), an analog continuous energy simulation may be unable to +transport any particles to that material. No tallies are scored there, and the +resulting cross sections for that material are zero. This situation is common in +shielding problems. + +This limitation can be overcome by "bootstrapping" the cross section generation +with weight windows. The idea is to first cheaply produce a set of weight +windows that cover the entire problem and then reuse them to push particles into +the far regions during the higher fidelity ``"material_wise"`` solve. The weight +windows are generated using the ``"stochastic_slab"`` method (which produces +cross sections for *all* materials regardless of their location) together with +the random ray solver and a :class:`~openmc.WeightWindowGenerator`, exactly as +described in the :ref:`FW-CADIS user guide `. The resulting +``weight_windows.h5`` file is then passed to a second, higher fidelity +``"material_wise"`` cross section generation via the ``weight_windows_file`` +argument:: + + # First, generate weight windows with the stochastic slab method and random + # ray (see the FW-CADIS user guide), producing a weight_windows.h5 file. + ... + + # Then, bootstrap a higher fidelity material-wise library, applying those + # weight windows during the continuous energy solve so that particles can + # reach materials far from the source. + model.convert_to_multigroup( + method="material_wise", + weight_windows_file="weight_windows.h5", + overwrite_mgxs_library=True, + ) + +The ``weight_windows_file`` argument is only used with the ``"material_wise"`` +method, as the ``"stochastic_slab"`` and ``"infinite_medium"`` methods use +simplified surrogate geometries that are incompatible with a weight window mesh +defined over the original geometry (and do not need weight windows, since they +already tally all materials). A warning is issued and the argument is ignored if +it is supplied to another method. + ~~~~~~~~~~~~ The Hard Way ~~~~~~~~~~~~ @@ -1075,9 +1124,9 @@ The adjoint flux random ray solver mode can be enabled as:: settings.random_ray['adjoint'] = True -When enabled, OpenMC will first run a forward transport simulation if there are -no user-specified adjoint sources present, followed by an adjoint transport -simulation. Fixed adjoint sources can be specified on the +When enabled, OpenMC will first run a forward transport simulation if there are +no user-specified adjoint sources present, followed by an adjoint transport +simulation. Fixed adjoint sources can be specified on the :attr:`openmc.Settings.random_ray` dictionary as follows:: # Geometry definition @@ -1090,21 +1139,21 @@ simulation. Fixed adjoint sources can be specified on the energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths) adj_source = openmc.IndependentSource( - energy=energy_distribution, + energy=energy_distribution, constraints={'domains': [detector_cell]} ) # Add to random_ray dict settings.random_ray['adjoint_source'] = adj_source -The same constraints apply to the user-defined adjoint source as to the forward -source, described in the :ref:`Fixed Source and Eigenvalue section -`. If this source is not provided, a forward -solve must take place to compute the adjoint external source when a forward -external source is present in the problem. Simulation settings (e.g., number of -rays, batches, etc.) will be identical for both calculations. At the -conclusion of the run, all results (e.g., tallies, plots, etc.) will be -derived from the adjoint flux rather than the forward flux but are not labeled +The same constraints apply to the user-defined adjoint source as to the forward +source, described in the :ref:`Fixed Source and Eigenvalue section +`. If this source is not provided, a forward +solve must take place to compute the adjoint external source when a forward +external source is present in the problem. Simulation settings (e.g., number of +rays, batches, etc.) will be identical for both calculations. At the +conclusion of the run, all results (e.g., tallies, plots, etc.) will be +derived from the adjoint flux rather than the forward flux but are not labeled any differently. When an initial forward solve is performed (i.e., when no user-specified adjoint source is present), its output files are also written to disk with a ``forward`` infix, so they are not overwritten by the subsequent @@ -1116,10 +1165,10 @@ generating FW-CADIS weight windows, no weight window file is written for the forward solve, as only the final adjoint-derived weight windows are meaningful. .. note:: - Use of the automated - :ref:`FW-CADIS weight window generator` is not - currently compatible with user-defined adjoint sources. Instead, the - initial forward calculation is used to assign "forward-weighted" adjoint + Use of the automated + :ref:`FW-CADIS weight window generator` is not + currently compatible with user-defined adjoint sources. Instead, the + initial forward calculation is used to assign "forward-weighted" adjoint sources to the tally regions of interest. --------------------------------------- diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 26b4a224c0e..a1d94e5819e 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -75,6 +75,19 @@ constexpr double ZERO_FLUX_CUTOFF {1e-22}; // value will be converted to pure void. constexpr double MINIMUM_MACRO_XS {1e-6}; +// Relative dead band applied to weight window comparisons: particles split +// only above upper * (1 + tol) and roulette only below lower * (1 - tol). +// Weight window arithmetic can land a particle's weight exactly back on a +// bound value (e.g., a roulette survivor is assigned survival_ratio * lower +// and a later split divides that back down), in which case the branch taken +// would be decided by the last ulp of the bound. Since window data carries +// ulp-level noise from non-associative parallel reductions in the solver that +// generated it, transport results would otherwise be chaotically sensitive to +// bit-level differences in the weight window file. Treating weights within +// the band as inside the window is statistically negligible, and weight +// window games are unbiased regardless of where the thresholds sit. +constexpr double WEIGHT_WINDOW_REL_TOL {1e-9}; + // ============================================================================ // MATH AND PHYSICAL CONSTANTS diff --git a/openmc/examples.py b/openmc/examples.py index 6895bd94b53..6c18b81e41e 100644 --- a/openmc/examples.py +++ b/openmc/examples.py @@ -1314,17 +1314,17 @@ def fill_cube(N, n_1, n_2, fill_1, fill_2, fill_3): def random_ray_three_region_cube_with_detectors() -> openmc.Model: """Create a three region cube model with two external tally regions. - This is an adaptation of the simple monoenergetic problem of a cube with - three concentric cubic regions. The innermost region is near void (with - Sigma_t around 10^-5) and contains an external isotropic source term, the - middle region is a mild scatterer (with Sigma_t around 10^-3), and the - outer region of the cube is a scatterer and absorber (with Sigma_t around + This is an adaptation of the simple monoenergetic problem of a cube with + three concentric cubic regions. The innermost region is near void (with + Sigma_t around 10^-5) and contains an external isotropic source term, the + middle region is a mild scatterer (with Sigma_t around 10^-3), and the + outer region of the cube is a scatterer and absorber (with Sigma_t around 1). - Two cubic "detector" regions are found outside this geometry, one along the - y-axis near z=0, and the other in the upper right corner of the system. - The size of each detector is scaled to be equal to that of the source - region. The model returned by this function contains cell tallies on each + Two cubic "detector" regions are found outside this geometry, one along the + y-axis near z=0, and the other in the upper right corner of the system. + The size of each detector is scaled to be equal to that of the source + region. The model returned by this function contains cell tallies on each detector. Returns @@ -1498,29 +1498,29 @@ def fill_cube(N, n_1, n_2, fill_1, fill_2, fill_3): fill=absorber_mat, region=detector2_region ) - + external_x = ( - +x_high & +y_low & +z_low & -x_outer & + +x_high & +y_low & +z_low & -x_outer & ((-y_outer & -z_high) | (-y_high & +z_high & -z_outer)) ) external_y = ( - +y_high & -y_outer & + +y_high & -y_outer & ( - (+detector1_right & -x_high & +z_low & -z_outer) | - (-detector1_right & +x_low & +detector1_top & -z_outer) | + (+detector1_right & -x_high & +z_low & -z_outer) | + (-detector1_right & +x_low & +detector1_top & -z_outer) | (+x_high & -x_outer & +z_low & -z_high) ) ) external_z = ( - +x_low & +y_low & +z_high & -z_outer & + +x_low & +y_low & +z_high & -z_outer & ((-y_outer & -x_high) | (-y_high & +x_high & -x_outer)) ) - external_cell = openmc.Cell(fill=cavity_mat, - region=(external_x | external_y | external_z), + external_cell = openmc.Cell(fill=cavity_mat, + region=(external_x | external_y | external_z), name='outside cube') root = openmc.Universe( - name='root universe', + name='root universe', cells=[cube_domain, detector1, detector2, external_cell] ) @@ -1604,8 +1604,8 @@ def fill_cube(N, n_1, n_2, fill_1, fill_2, fill_3): source_tally.estimator = estimator # Instantiate a Tallies collection and export to XML - tallies = openmc.Tallies([detector1_tally, - detector2_tally, + tallies = openmc.Tallies([detector1_tally, + detector2_tally, absorber_tally, cavity_tally, source_tally]) @@ -1619,3 +1619,116 @@ def fill_cube(N, n_1, n_2, fill_1, fill_2, fill_3): model.tallies = tallies return model + + +def sphere_with_shielded_pocket() -> openmc.Model: + """Create a continuous energy deep-shielding model with a far detector pocket. + + A concrete sphere is centered at the origin. A 2 MeV isotropic neutron + source sits in a small air cavity just inside the sphere surface on the -x + side, and a small steel pocket is embedded flush with the surface on the + +x axis, so roughly a meter of concrete separates the source from the + pocket while only a few centimeters of concrete back the cavity. The + sphere is enclosed in a vacuum-bounded box, with a void gap in between, so + that solvers which sample uniformly over a rectangular domain (e.g., + random ray) can be applied directly. The geometry is designed for testing + weight window and variance reduction workflows: + + - The probability that an analog source neutron reaches the steel pocket is + ~4e-5 (the product of the concrete attenuation and the pocket's small + solid angle), so an analog simulation with a few hundred histories + essentially never tallies the steel, while even crude global weight + windows allow particles to reach it reliably. + - Because the cavity sits near the surface, deep shielding (and thus a wide + weight window dynamic range) exists only within the small solid angle + subtended by the pocket, which keeps weight window splitting cheap and + convergent and the whole model fast enough for regression testing. + + Returns + ------- + model : openmc.Model + A deep-shielding model with a steel pocket behind a thick concrete + shield + + """ + model = openmc.Model() + + ########################################################################### + # Materials (few nuclides, to keep data loading cheap in multi-solve tests) + + air = openmc.Material(name='Air') + air.set_density('g/cm3', 0.001225) + air.add_nuclide('N14', 0.79, 'ao') + air.add_nuclide('O16', 0.21, 'ao') + + concrete = openmc.Material(name='Concrete') + concrete.set_density('g/cm3', 2.3) + concrete.add_nuclide('H1', 0.168759, 'ao') + concrete.add_nuclide('O16', 0.562489, 'ao') + concrete.add_nuclide('Si28', 0.203031, 'ao') + concrete.add_nuclide('Ca40', 0.044849, 'ao') + concrete.add_nuclide('Al27', 0.020872, 'ao') + + steel = openmc.Material(name='Steel') + steel.set_density('g/cm3', 7.87) + steel.add_nuclide('Fe56', 1.0) + + ########################################################################### + # Geometry + + # ~92 cm of concrete separates the cavity from the pocket face, while only + # ~6 cm of concrete backs the cavity on the -x side, so deep shielding is + # confined to the solid angle subtended by the pocket. + r_sphere = 66.0 + box_half_width = 70.0 + cavity_center_x = -54.0 + cavity_half_width = 6.0 + pocket_inner_face = 44.0 + pocket_half_width = 4.0 + + sphere = openmc.Sphere(r=r_sphere) + outer_box = openmc.model.RectangularParallelepiped( + -box_half_width, box_half_width, + -box_half_width, box_half_width, + -box_half_width, box_half_width, boundary_type='vacuum') + cavity_box = openmc.model.RectangularParallelepiped( + cavity_center_x - cavity_half_width, cavity_center_x + cavity_half_width, + -cavity_half_width, cavity_half_width, + -cavity_half_width, cavity_half_width) + # The pocket box extends past the sphere surface and is clipped by it, so + # the pocket sits flush with (and just inside) the outer surface. + pocket_box = openmc.model.RectangularParallelepiped( + pocket_inner_face, r_sphere + 1.0, + -pocket_half_width, pocket_half_width, + -pocket_half_width, pocket_half_width) + + cavity_cell = openmc.Cell(name='cavity', fill=air, region=-cavity_box) + pocket_cell = openmc.Cell(name='pocket', fill=steel, + region=-pocket_box & -sphere) + concrete_cell = openmc.Cell( + name='concrete', fill=concrete, + region=-sphere & +cavity_box & ~(-pocket_box & -sphere)) + void_cell = openmc.Cell(name='void', region=+sphere & -outer_box) + + model.geometry = openmc.Geometry( + [cavity_cell, pocket_cell, concrete_cell, void_cell]) + + ########################################################################### + # Source and settings + + source = openmc.IndependentSource() + source.space = openmc.stats.Box( + [cavity_center_x - cavity_half_width, + -cavity_half_width, -cavity_half_width], + [cavity_center_x + cavity_half_width, + cavity_half_width, cavity_half_width]) + source.constraints = {'domains': [cavity_cell]} + source.angle = openmc.stats.Isotropic() + source.energy = openmc.stats.delta_function(2.0e6) + + model.settings.run_mode = 'fixed source' + model.settings.source = source + model.settings.particles = 1000 + model.settings.batches = 10 + + return model diff --git a/openmc/model/model.py b/openmc/model/model.py index 3284942885e..82b0ae4dad8 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -2526,6 +2526,7 @@ def _isothermal_materialwise_mgxs( directory: PathLike, temperature_settings: dict, temperature: float | None = None, + weight_windows_file: PathLike | None = None, ) -> dict[str, openmc.XSdata]: """Generate a material-wise MGXS library for the model by running the original continuous energy OpenMC simulation. If a temperature is @@ -2556,6 +2557,10 @@ def _isothermal_materialwise_mgxs( The isothermal temperature value to apply to the materials in the input model. If not specified, defaults to the temperatures in the materials. + weight_windows_file : PathLike, optional + Path to a weight windows file to load and apply during the + continuous energy solve. See :meth:`Model.convert_to_multigroup` + for details. Returns ------- @@ -2576,6 +2581,14 @@ def _isothermal_materialwise_mgxs( model.settings.output = {'summary': True, 'tallies': False} model.settings.temperature = temperature_settings + # If a weight window file was provided, load and apply the weight + # windows during the continuous energy solve. This allows materials far + # from the source -- which an analog simulation may struggle to reach -- + # to still be tallied, and thus obtain nonzero cross sections. + if weight_windows_file is not None: + model.settings.weight_windows_file = weight_windows_file + model.settings.weight_windows_on = True + # Generate MGXS mgxs_lib = Model._auto_generate_mgxs_lib( model, groups, correction, directory) @@ -2602,6 +2615,7 @@ def _generate_material_wise_mgxs( directory: PathLike, temperatures: Sequence[float] | None = None, temperature_settings: dict | None = None, + weight_windows_file: PathLike | None = None, ) -> None: """Generate a material-wise MGXS library for the model by running the original continuous energy OpenMC simulation of the full material @@ -2634,6 +2648,10 @@ def _generate_material_wise_mgxs( A dictionary of temperature settings to use when generating MGXS. Valid entries for temperature_settings are the same as the valid entries in openmc.Settings.temperature_settings. + weight_windows_file : PathLike, optional + Path to a weight windows file to load and apply during the + continuous energy solve. See :meth:`Model.convert_to_multigroup` + for details. """ temp_settings = {} if temperature_settings is None: @@ -2648,7 +2666,8 @@ def _generate_material_wise_mgxs( nparticles, correction, directory, - temp_settings + temp_settings, + weight_windows_file=weight_windows_file ).values() # Write the file to disk. @@ -2667,7 +2686,8 @@ def _generate_material_wise_mgxs( correction, directory, temp_settings, - temperature + temperature, + weight_windows_file=weight_windows_file ) # Unpack the isothermal XSData objects and build a single XSData object per material. @@ -2695,6 +2715,7 @@ def convert_to_multigroup( source_energy: openmc.stats.Univariate | None = None, temperatures: Sequence[float] | None = None, temperature_settings: dict | None = None, + weight_windows_file: PathLike | None = None, ): """Convert all materials from continuous energy to multigroup. @@ -2747,10 +2768,50 @@ def convert_to_multigroup( A dictionary of temperature settings to use when generating MGXS. Valid entries for temperature_settings are the same as the valid entries in openmc.Settings.temperature_settings. + weight_windows_file : PathLike, optional + Path to a weight windows file (e.g., ``"weight_windows.h5"``) to + load and apply during the continuous energy MGXS generation + simulation. Applying weight windows allows the simulation to obtain + tallies -- and thus nonzero cross sections -- for materials located + far from the source, for example behind a thick shield, which an + analog simulation may be unable to reach. A typical use case is to + first generate weight windows with the ``"stochastic_slab"`` method + and the random ray solver, then "bootstrap" a higher-fidelity + ``"material_wise"`` library by passing those weight windows here. + This argument is only used with the ``"material_wise"`` method; a + warning is issued and the argument is ignored for the + ``"stochastic_slab"`` and ``"infinite_medium"`` methods. + + .. versionadded:: 0.15.4 """ if not isinstance(groups, openmc.mgxs.EnergyGroups): groups = openmc.mgxs.EnergyGroups(groups) + # Weight windows are only applicable to the "material_wise" method, + # which runs a continuous energy simulation of the original geometry + # (and is therefore compatible with a weight window mesh defined over + # that geometry). The "stochastic_slab" and "infinite_medium" methods + # use simplified surrogate geometries for which the provided weight + # windows are neither applicable nor needed. + if weight_windows_file is not None: + if method == "material_wise": + # Resolve to an absolute path now: MGXS generation runs in a + # temporary working directory, so a path relative to the current + # working directory would not be found at run time. + weight_windows_file = Path(weight_windows_file).resolve() + if not weight_windows_file.is_file(): + raise FileNotFoundError( + f'Weight windows file "{weight_windows_file}" could ' + 'not be found.' + ) + else: + warnings.warn( + 'The "weight_windows_file" argument is only applicable to ' + 'the "material_wise" MGXS generation method and will be ' + f'ignored for the "{method}" method.' + ) + weight_windows_file = None + # Do all work (including MGXS generation) in a temporary directory # to avoid polluting the working directory with residual XML files with TemporaryDirectory() as tmpdir: @@ -2791,7 +2852,7 @@ def convert_to_multigroup( elif method == "material_wise": self._generate_material_wise_mgxs( groups, nparticles, mgxs_path, correction, tmpdir, - temperatures, temperature_settings) + temperatures, temperature_settings, weight_windows_file) elif method == "stochastic_slab": self._generate_stochastic_slab_mgxs( groups, nparticles, mgxs_path, correction, tmpdir, source_energy, diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 00fff99a7c9..6d6f2c743de 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -392,6 +392,11 @@ void RandomRaySimulation::simulate() // Begin main simulation timer simulation::time_total.start(); + // Reset per-solve accumulators, as simulate() may run more than once on the + // same object (e.g. forward then adjoint when generating weight windows) + avg_miss_rate_ = 0.0; + total_geometric_intersections_ = 0; + // Random ray power iteration loop while (simulation::current_batch < settings::n_batches) { // Initialize the current batch diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 39e46026cf0..416d1421c3e 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -945,6 +945,11 @@ void apply_weight_windows(Particle& p) if (!settings::weight_windows_on) return; + // Random ray rays are not Monte Carlo particles and must not be biased by + // weight windows; the solver generates weight windows but never applies them + if (settings::solver_type == SolverType::RANDOM_RAY) + return; + // WW on photon and neutron only if (!p.type().is_neutron() && !p.type().is_photon()) return; @@ -1003,14 +1008,25 @@ void apply_weight_window(Particle& p, WeightWindow weight_window) 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) { + // If the particle's weight is above the weight window, split it until the + // resulting particles are within the window. The comparisons use a relative + // dead band so that the branch taken is insensitive to bit-level differences + // in the window bounds (see WEIGHT_WINDOW_REL_TOL). + if (weight > weight_window.upper_weight * (1.0 + WEIGHT_WINDOW_REL_TOL)) { // 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); + // The (1 - WEIGHT_WINDOW_REL_TOL) factor keeps the number of splits + // stable when the weight-to-bound ratio sits within rounding of an exact + // integer, which the weight window arithmetic itself can produce (e.g., a + // roulette survivor assigned weight * max_split, later split against an + // upper bound that is an exact multiple of the same lower bound). Ratios + // within the dead band of an integer consistently round down; the lower + // clamp of 2 preserves the guarantee that the split branch always splits. + double n_split = + std::max(2.0, std::ceil((1.0 - WEIGHT_WINDOW_REL_TOL) * weight / + weight_window.upper_weight)); double max_split = weight_window.max_split; n_split = std::min(n_split, max_split); @@ -1024,7 +1040,8 @@ void apply_weight_window(Particle& p, WeightWindow weight_window) // remaining weight is applied to current particle p.wgt() = weight / n_split; - } else if (weight <= weight_window.lower_weight) { + } else if (weight < + weight_window.lower_weight * (1.0 - WEIGHT_WINDOW_REL_TOL)) { // 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); diff --git a/tests/regression_tests/random_ray_auto_convert_bootstrap/__init__.py b/tests/regression_tests/random_ray_auto_convert_bootstrap/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/random_ray_auto_convert_bootstrap/inputs_true.dat b/tests/regression_tests/random_ray_auto_convert_bootstrap/inputs_true.dat new file mode 100644 index 00000000000..87043f116d6 --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_bootstrap/inputs_true.dat @@ -0,0 +1,81 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 20 + 10 + + + -60.0 -6.0 -6.0 -48.0 6.0 6.0 + + + + 2000000.0 1.0 + + + cell + 1 + + + true + weight_windows.h5 + + true + true + + 100000 + + + + 1 2 3 4 + + + 193 + flux + + + diff --git a/tests/regression_tests/random_ray_auto_convert_bootstrap/results_true.dat b/tests/regression_tests/random_ray_auto_convert_bootstrap/results_true.dat new file mode 100644 index 00000000000..f5d1b050692 --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_bootstrap/results_true.dat @@ -0,0 +1,9 @@ +tally 1: +9.500337E+01 +9.709273E+02 +2.417257E-04 +1.153221E-08 +9.046900E+02 +8.904339E+04 +1.034722E+02 +1.134807E+03 diff --git a/tests/regression_tests/random_ray_auto_convert_bootstrap/test.py b/tests/regression_tests/random_ray_auto_convert_bootstrap/test.py new file mode 100644 index 00000000000..52b15c0c7b3 --- /dev/null +++ b/tests/regression_tests/random_ray_auto_convert_bootstrap/test.py @@ -0,0 +1,101 @@ +import copy +import os + +import numpy as np +import openmc +import openmc.mgxs +from openmc.examples import sphere_with_shielded_pocket + +from tests.testing_harness import TolerantPyAPITestHarness + + +GROUPS = 'CASMO-4' + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + for f in ('mgxs.h5', 'weight_windows.h5'): + if os.path.exists(f): + os.remove(f) + + +def generate_weight_windows(model, mesh): + """Convert a multigroup model to random ray and run a FW-CADIS weight + window generation solve, producing weight_windows.h5.""" + model.convert_to_random_ray() + rr = model.settings.random_ray + rr['source_region_meshes'] = [(mesh, [model.geometry.root_universe])] + rr['distance_inactive'] = 100.0 + rr['distance_active'] = 400.0 + rr['source_shape'] = 'flat' + rr['sample_method'] = 'halton' + rr['volume_estimator'] = 'hybrid' + model.settings.particles = 1000 + model.settings.batches = 20 + model.settings.inactive = 15 + model.settings.weight_window_generators = openmc.WeightWindowGenerator( + method='fw_cadis', mesh=mesh, max_realizations=20, + energy_bounds=openmc.mgxs.EnergyGroups(GROUPS).group_edges) + model.run() + + +def test_random_ray_auto_convert_bootstrap(): + # Tests the weight window bootstrapped MGXS generation workflow, which + # consists of five OpenMC runs: a stochastic_slab MGXS generation, a random + # ray FW-CADIS weight window generation, a material_wise MGXS generation + # using those weight windows (which reach the steel pocket that an analog + # solve essentially never tallies), a second FW-CADIS generation using the + # bootstrapped library, and a final continuous energy Monte Carlo solve + # using the improved weight windows, whose tallies are compared against + # the reference results. + openmc.reset_auto_ids() + model = sphere_with_shielded_pocket() + + # Used by the continuous energy solves below; these have no effect on the + # random ray solves (random ray generates weight windows but never applies + # them). + model.settings.weight_window_checkpoints = { + 'collision': True, 'surface': True} + model.settings.max_history_splits = 100_000 + + # Overlay a ~10 cm source region / weight window mesh on the geometry + bbox = model.geometry.bounding_box + mesh = openmc.RegularMesh() + mesh.dimension = np.round( + (np.asarray(bbox.upper_right) - np.asarray(bbox.lower_left)) / 10.0 + ).astype(int) + mesh.lower_left = bbox.lower_left + mesh.upper_right = bbox.upper_right + + # Generate weight windows covering the whole problem with the + # stochastic_slab method and a random ray FW-CADIS solve + slab_model = copy.deepcopy(model) + slab_model.convert_to_multigroup( + method='stochastic_slab', groups=GROUPS, nparticles=50, + overwrite_mgxs_library=True, mgxs_path='mgxs.h5') + generate_weight_windows(slab_model, mesh) + + # Bootstrap the material-wise MGXS generation with those weight windows, + # then regenerate the weight windows from the higher-fidelity library + boot_model = copy.deepcopy(model) + boot_model.convert_to_multigroup( + method='material_wise', groups=GROUPS, nparticles=1, + overwrite_mgxs_library=True, mgxs_path='mgxs.h5', + weight_windows_file='weight_windows.h5') + generate_weight_windows(boot_model, mesh) + + # Run the continuous energy model with the improved weight windows, + # tallying the flux in every region + model.settings.weight_windows_file = 'weight_windows.h5' + model.settings.weight_windows_on = True + model.settings.particles = 20 + model.settings.batches = 10 + tally = openmc.Tally(name='flux') + tally.filters = [ + openmc.CellFilter(list(model.geometry.get_all_cells().values()))] + tally.scores = ['flux'] + model.tallies = openmc.Tallies([tally]) + + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/weightwindows/local/results_true.dat b/tests/regression_tests/weightwindows/local/results_true.dat index d3dfa119d84..ebff4edacbd 100644 --- a/tests/regression_tests/weightwindows/local/results_true.dat +++ b/tests/regression_tests/weightwindows/local/results_true.dat @@ -1 +1 @@ -a78972fe3c0dadfc256cfb139c5ca8c7b634738e1895ad2d6286ed5e2567c6dc518e499363908e9058432684148a706a179a39110f703d5322491184a1d0c3e4 \ No newline at end of file +41d8ba482783b724c61bf058a67ff701109076ac460f0be69b3a33c7a9e3e2aa5b5773442d33561b71c69500a9535f1ffb0f5e81f6ce364288c5094bbbc0a292 \ No newline at end of file diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index efa203b510a..c9b24e3059f 100644 --- a/tests/unit_tests/weightwindows/test.py +++ b/tests/unit_tests/weightwindows/test.py @@ -249,7 +249,7 @@ def test_photon_heating(run_in_tmpdir, shared_secondary): model.settings.run_mode = 'fixed source' model.settings.batches = 5 - model.settings.particles = 101 + model.settings.particles = 100 model.settings.shared_secondary_bank = shared_secondary tally = openmc.Tally()