MGXS Bootstrapping

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

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 <usersguide_fw_cadis>`. 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 
-<usersguide_random_ray_run_modes>`. 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
+<usersguide_random_ray_run_modes>`. 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<usersguide_fw_cadis>` 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<usersguide_fw_cadis>` 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.
 
 ---------------------------------------

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
 

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