FW-CADIS Weight Window Generation with Random Ray (#3273)

Co-authored-by: Olek <[email protected]>
Co-authored-by: Paul Romano <[email protected]>

docs/source/io_formats/settings.rst

@@ -1398,7 +1398,7 @@ mesh-based weight windows.
     *Default*: true
 
   :method:
-    Method used to update weight window values (currently only 'magic' is supported)
+    Method used to update weight window values (one of 'magic' or 'fw_cadis')
 
     *Default*: magic
 

docs/source/methods/index.rst

@@ -20,4 +20,5 @@ Theory and Methodology
     energy_deposition
     parallelization
     cmfd
+    variance_reduction
     random_ray

docs/source/methods/random_ray.rst

@@ -1060,6 +1060,8 @@ random ray and Monte Carlo, however.
   develop the scattering source by way of inactive batches before beginning
   active batches.
 
+.. _adjoint:
+
 ------------------------
 Adjoint Flux Solver Mode
 ------------------------

docs/source/methods/variance_reduction.rst

@@ -0,0 +1,134 @@
+.. _methods_variance_reduction:
+
+==================
+Variance Reduction
+==================
+
+.. _methods_variance_reduction_intro:
+
+------------
+Introduction
+------------
+
+Transport problems can sometimes involve a significant degree of attenuation
+between the source and a detector (tally) region, which can result in a flux
+differential of ten orders of magnitude (or more) throughout the simulation
+domain. As Monte Carlo uncertainties tend to be inversely proportional to the
+physical flux density, it can be extremely difficult to accurately resolve
+tallies in locations that are optically far from the source. This issue is
+particularly common in fixed source simulations, where some tally locations may
+not experience a single scoring event, even after billions of analog histories.
+
+Variance reduction techniques aim to either flatten the global uncertainty
+distribution, such that all regions of phase space have a fairly similar
+uncertainty, or to reduce the uncertainty in specific locations (such as a
+detector). There are two strategies available in OpenMC for variance reduction:
+the Monte Carlo MAGIC method and the FW-CADIS method. Both strategies work by
+developing a weight window mesh that can be utilized by subsequent Monte Carlo
+solves to split particles heading towards areas of lower flux densities while
+terminating particles in higher flux regions---all while maintaining a fair
+game.
+
+------------
+MAGIC Method
+------------
+
+The Method of Automatic Generation of Importances by Calculation, or `MAGIC
+method <https://doi.org/10.1016/j.fusengdes.2011.01.059>`_, is an iterative
+technique that uses spatial flux information :math:`\phi(r)` obtained from a
+normal Monte Carlo solve to produce weight windows :math:`w(r)` that can be
+utilized by a subsequent iteration of Monte Carlo. While the first generation of
+weight windows produced may only help to reduce variance slightly, use of these
+weights to generate another set of weight windows results in a progressively
+improving iterative scheme.
+
+Equation :eq:`magic` defines how the lower bound of weight windows
+:math:`w_{\ell}(r)` are generated with MAGIC using forward flux information.
+Here, we can see that the flux at location :math:`r` is normalized by the
+maximum flux in any group at that location. We can also see that the weights are
+divided by a factor of two, which accounts for the typical :math:`5\times`
+factor separating the lower and upper weight window bounds in OpenMC.
+
+.. math::
+    :label: magic
+
+    w_{\ell}(r) = \frac{\phi(r)}{2\,\text{max}(\phi(r))}
+
+A major advantage of this technique is that it does not require any special
+transport machinery; it simply uses multiple Monte Carlo simulations to
+iteratively improve a set of weight windows (which are typically defined on a
+mesh covering the simulation domain). The downside to this method is that as the
+flux differential increases between areas near and far from the source, it
+requires more outer Monte Carlo iterations, each of which can be expensive in
+itself. Additionally, computation of weight windows based on regular (forward)
+neutron flux tally information does not produce the most numerically effective
+set of weight windows. Nonetheless, MAGIC remains a simple and effective
+technique for generating weight windows.
+
+--------
+FW-CADIS
+--------
+
+As discussed in the previous section, computation of weight windows based on
+regular (forward) neutron flux tally information does not produce the most
+numerically efficient set of weight windows. It is highly preferable to generate
+weight windows based on spatial adjoint flux :math:`\phi^{\dag}(r)`
+information. The adjoint flux is essentially the "reverse" simulation problem,
+where we sample a random point and assume this is where a particle was absorbed,
+and then trace it backwards (upscattering in energy), until we sample the point
+where it was born from.
+
+The Forward-Weighted Consistent Adjoint Driven Importance Sampling method, or
+`FW-CADIS method <https://doi.org/10.13182/NSE12-33>`_, produces weight windows
+for global variance reduction given adjoint flux information throughout the
+entire domain. The weight window lower bound is defined in Equation
+:eq:`fw_cadis`, and also involves a normalization step not shown here.
+
+.. math::
+    :label: fw_cadis
+
+    w_{\ell}(r) = \frac{1}{2\phi^{\dag}(r)}
+
+While the algorithm itself is quite simple, it requires estimates of the global
+adjoint flux distribution, which is difficult to generate directly with Monte
+Carlo transport. Thus, FW-CADIS typically uses an alternative solver (often
+deterministic) that can be more readily adapted for generating adjoint flux
+information, and which is often much cheaper than Monte Carlo given that a rough
+solution is often sufficient for weight window generation.
+
+The FW-CADIS implementation in OpenMC utilizes its own internal random ray
+multigroup transport solver to generate the adjoint source distribution. No
+coupling to any external transport is solver is necessary. The random ray solver
+operates on the same geometry as the Monte Carlo solver, so no redefinition of
+the simulation geometry is required. More details on how the adjoint flux is
+computed are given in the :ref:`adjoint methods section <adjoint>`.
+
+More information on the workflow is available in the :ref:`user guide
+<variance_reduction>`, but generally production of weight windows with FW-CADIS
+involves several stages (some of which are highly automated). These tasks
+include generation of approximate multigroup cross section data for use by the
+random ray solver, running of the random ray solver in normal (forward flux)
+mode to generate a source for the adjoint solver, running of the random ray
+solver in adjoint mode to generate adjoint flux tallies, and finally the
+production of weight windows via the FW-CADIS method. As is discussed in the
+user guide, most of these steps are automated together, making the additional
+burden on the user fairly small.
+
+The major advantage of this technique is that it typically produces much more
+numerically efficient weight windows as compared to those generated with MAGIC,
+sometimes with an order-of-magnitude improvement in the figure of merit
+(Equation :eq:`variance_fom`), which accounts for both the variance and the
+execution time. Another major advantage is that the cost of the random ray
+solver is typically negligible compared to the cost of the subsequent Monte
+Carlo solve itself, making it a very cheap method to deploy. The downside to
+this method is that it introduces a second transport method into the mix (random
+ray), such that there are more free input parameters for the user to know about
+and adjust, potentially making the method more complex to use. However, as many
+of the parameters have natural choices, much of this parameterization can be
+handled automatically behind the scenes without the need for the user to be
+aware of this.
+
+.. math::
+    :label: variance_fom
+
+    \text{FOM} = \frac{1}{\text{Time} \times \sigma^2}

docs/source/usersguide/index.rst

@@ -25,6 +25,7 @@ essential aspects of using OpenMC to perform simulations.
     processing
     parallel
     volume
+    variance_reduction
     random_ray
     troubleshoot
 

docs/source/usersguide/random_ray.rst

@@ -435,10 +435,11 @@ Inputting Multigroup Cross Sections (MGXS)
 Multigroup cross sections for use with OpenMC's random ray solver are input the
 same way as with OpenMC's traditional multigroup Monte Carlo mode. There is more
 information on generating multigroup cross sections via OpenMC in the
-:ref:`multigroup materials <create_mgxs>` user guide. You may also wish to
-use an existing multigroup library. An example of using OpenMC's Python
-interface to generate a correctly formatted ``mgxs.h5`` input file is given
-in the `OpenMC Jupyter notebook collection
+:ref:`multigroup materials <create_mgxs>` user guide. You may also wish to use
+an existing ``mgxs.h5`` MGXS library file, or define your own given a known set
+of cross section data values (e.g., as taken from a benchmark specification). An
+example of using OpenMC's Python interface to generate a correctly formatted
+``mgxs.h5`` input file is given in the `OpenMC Jupyter notebook collection
 <https://nbviewer.org/github/openmc-dev/openmc-notebooks/blob/main/mg-mode-part-i.ipynb>`_.
 
 .. note::
@@ -447,6 +448,184 @@ in the `OpenMC Jupyter notebook collection
     separate materials can be defined each with a separate multigroup dataset
     corresponding to a given temperature.
 
+.. _mgxs_gen:
+
+-------------------------------------------
+Generating Multigroup Cross Sections (MGXS)
+-------------------------------------------
+
+OpenMC is capable of generating multigroup cross sections by way of flux
+collapsing data based on flux solutions obtained from a continuous energy Monte
+Carlo solve. While it is a circular excercise in some respects to use continuous
+energy Monte Carlo to generate cross sections to be used by a reduced-fidelity
+multigroup transport solver, there are many use cases where this is nonetheless
+highly desirable. For instance, generation of a multigroup library may enable
+the same set of approximate multigroup cross section data to be used across a
+variety of problem types (or through a multidimensional parameter sweep of
+design variables) with only modest errors and at greatly reduced cost as
+compared to using only continuous energy Monte Carlo.
+
+We give here a quick summary of how to produce a multigroup cross section data
+file (``mgxs.h5``) from a starting point of a typical continuous energy Monte
+Carlo input file. Notably, continuous energy input files define materials as a
+mixture of nuclides with different densities, whereas multigroup materials are
+simply defined by which name they correspond to in a ``mgxs.h5`` library file.
+
+To generate the cross section data, we begin with a continuous energy Monte
+Carlo input deck and add in the required tallies that will be needed to generate
+our library. In this example, we will specify material-wise cross sections and a
+two group energy decomposition::
+
+  # Define geometry
+  ...
+  ...
+  geometry = openmc.Geometry()
+  ...
+  ...
+
+  # Initialize MGXS library with a finished OpenMC geometry object
+  mgxs_lib = openmc.mgxs.Library(geometry)
+
+  # Pick energy group structure
+  groups = mgxs.EnergyGroups(mgxs.GROUP_STRUCTURES['CASMO-2'])
+  mgxs_lib.energy_groups = groups
+
+  # Disable transport correction
+  mgxs_lib.correction = None
+
+  # Specify needed cross sections for random ray
+  mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission',
+                      'nu-scatter matrix', 'multiplicity matrix', 'chi']
+
+  # Specify a "cell" domain type for the cross section tally filters
+  mgxs_lib.domain_type = "material"
+
+  # Specify the cell domains over which to compute multi-group cross sections
+  mgxs_lib.domains = geom.get_all_materials().values()
+
+  # Do not compute cross sections on a nuclide-by-nuclide basis
+  mgxs_lib.by_nuclide = False
+
+  # Check the library - if no errors are raised, then the library is satisfactory.
+  mgxs_lib.check_library_for_openmc_mgxs()
+
+  # Construct all tallies needed for the multi-group cross section library
+  mgxs_lib.build_library()
+
+  # Create a "tallies.xml" file for the MGXS Library
+  tallies = openmc.Tallies()
+  mgxs_lib.add_to_tallies_file(tallies, merge=True)
+
+  # Export
+  tallies.export_to_xml()
+
+  ...
+
+When selecting an energy decomposition, you can manually define group boundaries
+or pick out a group structure already known to OpenMC (a list of which can be
+found at :class:`openmc.mgxs.GROUP_STRUCTURES`). Once the above input deck has
+been run, the resulting statepoint file will contain the needed flux and
+reaction rate tally data so that a MGXS library file can be generated. Below is
+the postprocessing script needed to generate the ``mgxs.h5`` library file given
+a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g.,
+``summary.h5``) that resulted from running our previous example::
+
+  import openmc
+  import openmc.mgxs as mgxs
+
+  summary = openmc.Summary('summary.h5')
+  geom = summary.geometry
+  mats = summary.materials
+
+  statepoint_filename = 'statepoint.100.h5'
+  sp = openmc.StatePoint(statepoint_filename)
+
+  groups = mgxs.EnergyGroups(mgxs.GROUP_STRUCTURES['CASMO-2'])
+  mgxs_lib = openmc.mgxs.Library(geom)
+  mgxs_lib.energy_groups = groups
+  mgxs_lib.correction = None
+  mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission',
+                          'nu-scatter matrix', 'multiplicity matrix', 'chi']
+
+  # Specify a "cell" domain type for the cross section tally filters
+  mgxs_lib.domain_type = "material"
+
+  # Specify the cell domains over which to compute multi-group cross sections
+  mgxs_lib.domains = geom.get_all_materials().values()
+
+  # Do not compute cross sections on a nuclide-by-nuclide basis
+  mgxs_lib.by_nuclide = False
+
+  # Check the library - if no errors are raised, then the library is satisfactory.
+  mgxs_lib.check_library_for_openmc_mgxs()
+
+  # Construct all tallies needed for the multi-group cross section library
+  mgxs_lib.build_library()
+
+  mgxs_lib.load_from_statepoint(sp)
+
+  names = []
+  for mat in mgxs_lib.domains: names.append(mat.name)
+
+  # Create a MGXS File which can then be written to disk
+  mgxs_file = mgxs_lib.create_mg_library(xs_type='macro', xsdata_names=names)
+
+  # Write the file to disk using the default filename of "mgxs.h5"
+  mgxs_file.export_to_hdf5("mgxs.h5")
+
+Notably, the postprocessing script needs to match the same
+:class:`openmc.mgxs.Library` settings that were used to generate the tallies,
+but otherwise is able to discern the rest of the simulation details from the
+statepoint and summary files. Once the postprocessing script is successfully
+run, the ``mgxs.h5`` file can be loaded by subsequent runs of OpenMC.
+
+If you want to convert continuous energy material objects in an OpenMC input
+deck to multigroup ones from a ``mgxs.h5`` library, you can follow the below
+example. Here we begin with the original continuous energy materials we used to
+generate our MGXS library::
+
+    fuel = openmc.Material(name='UO2 (2.4%)')
+    fuel.set_density('g/cm3', 10.29769)
+    fuel.add_nuclide('U234', 4.4843e-6)
+    fuel.add_nuclide('U235', 5.5815e-4)
+    fuel.add_nuclide('U238', 2.2408e-2)
+    fuel.add_nuclide('O16', 4.5829e-2)
+
+    water = openmc.Material(name='Hot borated water')
+    water.set_density('g/cm3', 0.740582)
+    water.add_nuclide('H1', 4.9457e-2)
+    water.add_nuclide('O16', 2.4672e-2)
+    water.add_nuclide('B10', 8.0042e-6)
+    water.add_nuclide('B11', 3.2218e-5)
+    water.add_s_alpha_beta('c_H_in_H2O')
+
+    materials = openmc.Materials([fuel, water])
+
+Once the ``mgxs.h5`` library file has been generated, we can then manually make
+the necessary edits to the material definitions so that they load from the
+multigroup library instead of defining their isotopic contents, as::
+
+    # Instantiate some Macroscopic Data
+    fuel_data = openmc.Macroscopic('UO2 (2.4%)')
+    water_data = openmc.Macroscopic('Hot borated water')
+
+    # Instantiate some Materials and register the appropriate Macroscopic objects
+    fuel= openmc.Material(name='UO2 (2.4%)')
+    fuel.set_density('macro', 1.0)
+    fuel.add_macroscopic(fuel_data)
+
+    water= openmc.Material(name='Hot borated water')
+    water.set_density('macro', 1.0)
+    water.add_macroscopic(water_data)
+
+    # Instantiate a Materials collection and export to XML
+    materials = openmc.Materials([fuel, water])
+    materials.cross_sections = "mgxs.h5"
+
+In the above example, our ``fuel`` and ``water`` materials will now load MGXS
+data from the ``mgxs.h5`` file instead of loading continuous energy isotopic
+cross section data.
+
 --------------
 Linear Sources
 --------------
@@ -597,9 +776,7 @@ estimator, the following code would be used:
 Adjoint Flux Mode
 -----------------
 
-The adjoint flux random ray solver mode can be enabled as:
-entire
-::
+The adjoint flux random ray solver mode can be enabled as::
 
     settings.random_ray['adjoint'] = True