Random Ray Adjoint Mode (openmc-dev#3191)

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

docs/source/methods/random_ray.rst

@@ -542,7 +542,7 @@ in that cell for the iteration from Equation :eq:`phi_naive` to:
 .. math::
     :label: phi_missed_one
 
-    \phi_{i,g,n}^{missed} = \frac{Q_{i,g,n} }{\Sigma_{t,i,g}} 
+    \phi_{i,g,n}^{missed} = \frac{Q_{i,g,n} }{\Sigma_{t,i,g}}
 
 as the streaming operator has gone to zero. While this is obviously innacurate
 as it ignores transport, for most problems where the region is only occasionally
@@ -1060,6 +1060,49 @@ random ray and Monte Carlo, however.
   develop the scattering source by way of inactive batches before beginning
   active batches.
 
+------------------------
+Adjoint Flux Solver Mode
+------------------------
+
+The random ray solver in OpenMC can also be used to solve for the adjoint flux,
+:math:`\psi^{\dagger}`. In combination with the regular (forward) flux solution,
+the adjoint flux is useful for perturbation methods as well as for computing
+weight windows for subsequent Monte Carlo simulations. The adjoint flux can be
+thought of as the "backwards" flux, representing the flux where a particle is
+born at an absoprtion point (and typical absorption energy), and then undergoes
+transport with a transposed scattering matrix. That is, instead of sampling a
+particle and seeing where it might go as in a standard forward solve, we will
+sample an absorption location and see where the particle that was absorbed there
+might have come from. Notably, for typical neutron absorption at low energy
+levels, this means that adjoint flux particles are typically sampled at a low
+energy and then upscatter (via a transposed scattering matrix) over their
+lifetimes.
+
+In OpenMC, the random ray adjoint solver is implemented simply by transposing
+the scattering matrix, swapping :math:`\nu\Sigma_f` and :math:`\chi`, and then
+running a normal transport solve. When no external fixed source is present, no
+additional changes are needed in the transport process. However, if an external
+fixed forward source is present in the simulation problem, then an additional
+step is taken to compute the accompanying fixed adjoint source. In OpenMC, the
+adjoint flux does *not* represent a response function for a particular detector
+region. Rather, the adjoint flux is the global response, making it appropriate
+for use with weight window generation schemes for global variance reduction.
+Thus, if using a fixed source, the external source for the adjoint mode is
+simply computed as being :math:`1 / \phi`, where :math:`\phi` is the forward
+scalar flux that results from a normal forward solve (which OpenMC will run
+first automatically when in adjoint mode). The adjoint external source will be
+computed for each source region in the simulation mesh, independent of any
+tallies. The adjoint external source is always flat, even when a linear
+scattering and fission source shape is used. When in adjoint mode, all reported
+results (e.g., tallies, eigenvalues, etc.) are derived from the adjoint flux,
+even when the physical meaning is not necessarily obvious. These values are
+still reported, though we emphasize that the primary use case for adjoint mode
+is for producing adjoint flux tallies to support subsequent perturbation studies
+and weight window generation.
+
+Note that the adjoint :math:`k_{eff}` is statistically the same as the forward
+:math:`k_{eff}`, despite the flux distributions taking different shapes.
+
 ---------------------------
 Fundamental Sources of Bias
 ---------------------------

docs/source/usersguide/random_ray.rst

@@ -558,7 +558,7 @@ following methods are currently available in OpenMC:
      - Cons
    * - ``simulation_averaged``
      - Accumulates total active ray lengths in each FSR over all iterations,
-       improving the estimate of the volume in each cell each iteration. 
+       improving the estimate of the volume in each cell each iteration.
      - * Virtually unbiased after several iterations
        * Asymptotically approaches the true analytical volume
        * Typically most efficient in terms of speed vs. accuracy
@@ -593,6 +593,33 @@ estimator, the following code would be used:
 
     settings.random_ray['volume_estimator'] = 'naive'
 
+-----------------
+Adjoint Flux Mode
+-----------------
+
+The adjoint flux random ray solver mode can be enabled as:
+entire
+::
+
+    settings.random_ray['adjoint'] = True
+
+When enabled, OpenMC will first run a forward transport simulation followed by
+an adjoint transport simulation. The purpose of the forward solve is to compute
+the adjoint external source when an external source is present in the
+simulation. Simulation settings (e.g., number of rays, batches, etc.) will be
+identical for both simulations. 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. The initial forward flux
+solution will not be stored or available in the final statepoint file. Those
+wishing to do analysis requiring both the forward and adjoint solutions will
+need to run two separate simulations and load both statepoint files.
+
+.. note::
+    When adjoint mode is selected, OpenMC will always perform a full forward
+    solve and then run a full adjoint solve immediately afterwards. Statepoint
+    and tally results will be derived from the adjoint flux, but will not be
+    labeled any differently.
+
 ---------------------------------------
 Putting it All Together: Example Inputs
 ---------------------------------------

include/openmc/mgxs.h

@@ -86,8 +86,10 @@ class Mgxs {
   bool fissionable; // Is this fissionable
   bool is_isotropic {
     true}; // used to skip search for angle indices if isotropic
+  bool exists_in_model {true}; // Is this present in model
 
   Mgxs() = default;
+  Mgxs(bool exists) : exists_in_model(exists) {}
 
   //! \brief Constructor that loads the Mgxs object from the HDF5 file
   //!

include/openmc/random_ray/flat_source_domain.h

@@ -111,13 +111,17 @@ class FlatSourceDomain {
   virtual void all_reduce_replicated_source_regions();
   void convert_external_sources();
   void count_external_source_regions();
+  void set_adjoint_sources(const vector<double>& forward_flux);
   virtual void flux_swap();
   virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
   double compute_fixed_source_normalization_factor() const;
+  void flatten_xs();
+  void transpose_scattering_matrix();
 
   //----------------------------------------------------------------------------
   // Static Data members
   static bool volume_normalized_flux_tallies_;
+  static bool adjoint_; // If the user wants outputs based on the adjoint flux
 
   //----------------------------------------------------------------------------
   // Static data members
@@ -150,6 +154,19 @@ class FlatSourceDomain {
   vector<float> source_;
   vector<float> external_source_;
   vector<bool> external_source_present_;
+  vector<double> scalar_flux_final_;
+
+  // 2D arrays stored in 1D representing values for all materials x energy
+  // groups
+  int n_materials_;
+  vector<double> sigma_t_;
+  vector<double> nu_sigma_f_;
+  vector<double> sigma_f_;
+  vector<double> chi_;
+
+  // 3D arrays stored in 1D representing values for all materials x energy
+  // groups x energy groups
+  vector<double> sigma_s_;
 
 protected:
   //----------------------------------------------------------------------------
@@ -190,10 +207,6 @@ class FlatSourceDomain {
   vector<int> material_;
   vector<double> volume_naive_;
 
-  // 2D arrays stored in 1D representing values for all source regions x energy
-  // groups
-  vector<float> scalar_flux_final_;
-
   // Volumes for each tally and bin/score combination. This intermediate data
   // structure is used when tallying quantities that must be normalized by
   // volume (i.e., flux). The vector is index by tally index, while the inner 2D

include/openmc/random_ray/random_ray_simulation.h

@@ -20,6 +20,8 @@ class RandomRaySimulation {
   //----------------------------------------------------------------------------
   // Methods
   void compute_segment_correction_factors();
+  void prepare_fixed_sources();
+  void prepare_fixed_sources_adjoint(vector<double>& forward_flux);
   void simulate();
   void reduce_simulation_statistics();
   void output_simulation_results() const;
@@ -30,8 +32,13 @@ class RandomRaySimulation {
     int64_t n_external_source_regions) const;
 
   //----------------------------------------------------------------------------
-  // Data members
+  // Accessors
+  FlatSourceDomain* domain() const { return domain_.get(); }
+
 private:
+  //----------------------------------------------------------------------------
+  // Data members
+
   // Contains all flat source region data
   unique_ptr<FlatSourceDomain> domain_;
 

openmc/deplete/openmc_operator.py

@@ -208,7 +208,7 @@ def _get_burnable_mats(self) -> tuple[list[str], dict[str, float], list[str]]:
                 if nuclide in self.nuclides_with_data or self._decay_nucs:
                     model_nuclides.add(nuclide)
                 else:
-                    msg = (f"Nuclilde {nuclide} in material {mat.id} is not "
+                    msg = (f"Nuclide {nuclide} in material {mat.id} is not "
                            "present in the depletion chain and has no cross "
                            "section data.")
                     warn(msg)

openmc/settings.py

@@ -170,6 +170,9 @@ class Settings:
             cm/cm^3. When disabled, flux tallies will be reported in units
             of cm (i.e., total distance traveled by neutrons in the spatial
             tally region).
+        :adjoint:
+            Whether to run the random ray solver in adjoint mode (bool). The
+            default is 'False'.
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -1113,6 +1116,8 @@ def random_ray(self, random_ray: dict):
                                ('flat', 'linear', 'linear_xy'))
             elif key == 'volume_normalized_flux_tallies':
                 cv.check_type('volume normalized flux tallies', random_ray[key], bool)
+            elif key == 'adjoint':
+                cv.check_type('adjoint', random_ray[key], bool)
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
@@ -1916,6 +1921,10 @@ def _random_ray_from_xml_element(self, root):
                     self.random_ray['volume_normalized_flux_tallies'] = (
                         child.text in ('true', '1')
                     )
+                elif child.tag == 'adjoint':
+                    self.random_ray['adjoint'] = (
+                        child.text in ('true', '1')
+                    )
 
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.

src/mgxs_interface.cpp

@@ -146,7 +146,7 @@ void MgxsInterface::create_macro_xs()
         num_energy_groups_, num_delayed_groups_);
     } else {
       // Preserve the ordering of materials by including a blank entry
-      macro_xs_.emplace_back();
+      macro_xs_.emplace_back(false);
     }
   }
 }