Local adjoint source for Random Ray

include/openmc/random_ray/flat_source_domain.h

@@ -40,9 +40,10 @@ class FlatSourceDomain {
   void random_ray_tally();
   virtual void accumulate_iteration_flux();
   void output_to_vtk() const;
-  void convert_external_sources();
+  void convert_external_sources(bool use_adjoint_sources);
   void count_external_source_regions();
-  void set_adjoint_sources();
+  void set_fw_adjoint_sources();
+  void set_local_adjoint_sources();
   void flux_swap();
   virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
   double compute_fixed_source_normalization_factor() const;
@@ -76,6 +77,7 @@ class FlatSourceDomain {
   // Static Data members
   static bool volume_normalized_flux_tallies_;
   static bool adjoint_; // If the user wants outputs based on the adjoint flux
+  static bool cadis_;
   static double
     diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization
                                  // for transport corrected MGXS data

include/openmc/random_ray/random_ray_simulation.h

@@ -21,7 +21,8 @@ class RandomRaySimulation {
   // Methods
   void compute_segment_correction_factors();
   void apply_fixed_sources_and_mesh_domains();
-  void prepare_fixed_sources_adjoint();
+  void prepare_fw_fixed_sources_adjoint();
+  void prepare_local_fixed_sources_adjoint();
   void simulate();
   void output_simulation_results() const;
   void instability_check(

include/openmc/source.h

@@ -34,6 +34,7 @@ class Source;
 namespace model {
 
 extern vector<unique_ptr<Source>> external_sources;
+extern vector<unique_ptr<Source>> adjoint_sources;
 
 // Probability distribution for selecting external sources
 extern DiscreteIndex external_sources_probability;

include/openmc/weight_windows.h

@@ -17,7 +17,7 @@
 
 namespace openmc {
 
-enum class WeightWindowUpdateMethod { MAGIC, FW_CADIS };
+enum class WeightWindowUpdateMethod { MAGIC, FW_CADIS, CADIS };
 
 //==============================================================================
 // Constants

openmc/model/model.py

@@ -265,6 +265,37 @@ def add_kinetics_parameters_tallies(self, num_groups: int | None = None):
             denom_tally = openmc.Tally(name='IFP denominator')
             denom_tally.scores = ['ifp-denominator']
             self.tallies.append(denom_tally)
+
+    # TODO: This should also be incorporated into lower-level calls in 
+    # settings.py, but it requires information about the tallies currently
+    # on the active Model
+    def _assign_cadis_tally_IDs(self):
+        # Verify that all tallies assigned as targets on WeightWindowGenerators 
+        # exist within model.tallies. If this is the case, convert the .targets 
+        # attribute of each WeightWindowGenerator to a sequence of tally IDs.
+        if len(self.settings.weight_window_generators) == 0:
+            return
+
+        for wwg in self.settings.weight_window_generators:
+            # Only proceeds if the "targets" attribute is an openmc.Tallies, 
+            # which means it hasn't been checked against model.tallies.
+            if isinstance(wwg.targets, openmc.Tallies):
+                id_vec = []
+                for tal in wwg.targets:
+                    # check against model tallies
+                    id_next = None
+                    for reference_tal in self.tallies:
+                        if tal == reference_tal:
+                            id_next = reference_tal.id
+                            break
+
+                    if id_next == None:
+                        raise RuntimeError(
+                            f'CADIS target tally {tal.id} not found on model.tallies!')
+                    else:
+                        id_vec.append(id_next)
+
+                wwg.targets = id_vec
 
     @classmethod
     def from_xml(
@@ -576,6 +607,7 @@ def export_to_xml(self, directory: PathLike = '.', remove_surfs: bool = False,
         if not d.is_dir():
             d.mkdir(parents=True, exist_ok=True)
 
+        self._assign_cadis_tally_IDs()
         self.settings.export_to_xml(d)
         self.geometry.export_to_xml(d, remove_surfs=remove_surfs)
 
@@ -634,6 +666,9 @@ def export_to_model_xml(self, path: PathLike = 'model.xml', remove_surfs: bool =
                           "set the Geometry.merge_surfaces attribute instead.")
             self.geometry.merge_surfaces = True
 
+        # Link CADIS WeightWindowGenerator target tallies, if present
+        self._assign_cadis_tally_IDs()
+
         # provide a memo to track which meshes have been written
         mesh_memo = set()
         settings_element = self.settings.to_xml_element(mesh_memo)

openmc/settings.py

@@ -224,6 +224,9 @@ class Settings:
             stabilization, which may be desirable as stronger diagonal stabilization
             also tends to dampen the convergence rate of the solver, thus requiring
             more iterations to converge.
+        :adjoint_source:
+            Source object used to define localized adjoint source/detector response 
+            function.
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -1345,6 +1348,8 @@ def random_ray(self, random_ray: dict):
                 cv.check_type('diagonal stabilization rho', value, Real)
                 cv.check_greater_than('diagonal stabilization rho',
                                       value, 0.0, True)
+            elif key == 'adjoint_source':
+                cv.check_type('adjoint source', value, SourceBase)
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
@@ -1876,8 +1881,9 @@ def _create_random_ray_subelement(self, root, mesh_memo=None):
             element = ET.SubElement(root, "random_ray")
             for key, value in self._random_ray.items():
                 if key == 'ray_source' and isinstance(value, SourceBase):
+                    subelement = ET.SubElement(element, 'ray_source')
                     source_element = value.to_xml_element()
-                    element.append(source_element)
+                    subelement.append(source_element)
                 elif key == 'source_region_meshes':
                     subelement = ET.SubElement(element, 'source_region_meshes')
                     for mesh, domains in value:
@@ -1894,8 +1900,12 @@ def _create_random_ray_subelement(self, root, mesh_memo=None):
                         path = f"./mesh[@id='{mesh.id}']"
                         if root.find(path) is None:
                             root.append(mesh.to_xml_element())
-                            if mesh_memo is not None:
+                            if mesh_memo is not None:    
                                 mesh_memo.add(mesh.id)
+                elif key == 'adjoint_source' and isinstance(value, SourceBase):
+                    subelement = ET.SubElement(element, 'adjoint_source')
+                    adj_source_element = value.to_xml_element()
+                    subelement.append(adj_source_element)
                 elif isinstance(value, bool):
                     subelement = ET.SubElement(element, key)
                     subelement.text = str(value).lower()
@@ -2322,8 +2332,9 @@ def _random_ray_from_xml_element(self, root, meshes=None):
             for child in elem:
                 if child.tag in ('distance_inactive', 'distance_active', 'diagonal_stabilization_rho'):
                     self.random_ray[child.tag] = float(child.text)
-                elif child.tag == 'source':
-                    source = SourceBase.from_xml_element(child)
+                elif child.tag == 'ray_source':
+                    source_element = child.find('source')
+                    source = SourceBase.from_xml_element(source_element)
                     self.random_ray['ray_source'] = source
                 elif child.tag == 'volume_estimator':
                     self.random_ray['volume_estimator'] = child.text
@@ -2337,6 +2348,11 @@ def _random_ray_from_xml_element(self, root, meshes=None):
                     self.random_ray['adjoint'] = (
                         child.text in ('true', '1')
                     )
+                elif child.tag == 'adjoint_source':
+                    for subelem in child.findall('source'):
+                        src = SourceBase.from_xml_element(subelem)
+                        # add newly constructed source object to the list
+                        self.random_ray['adjoint_source'].append(src)
                 elif child.tag == 'sample_method':
                     self.random_ray['sample_method'] = child.text
                 elif child.tag == 'source_region_meshes':

openmc/weight_windows.py

@@ -12,6 +12,7 @@
 import openmc
 from openmc.filter import _PARTICLES
 from openmc.mesh import MeshBase, RectilinearMesh, CylindricalMesh, SphericalMesh, UnstructuredMesh
+from openmc.tally import Tallies
 import openmc.checkvalue as cv
 from openmc.checkvalue import PathLike
 from ._xml import get_elem_list, get_text, clean_indentation
@@ -496,8 +497,10 @@ class WeightWindowGenerator:
         maximum and minimum energy for the data available at runtime.
     particle_type : {'neutron', 'photon'}
         Particle type the weight windows apply to
-    method : {'magic', 'fw_cadis'}
+    method : {'magic', 'fw_cadis', 'cadis'}
         The weight window generation methodology applied during an update.
+    targets : :class:`openmc.Tallies`
+        Target tallies for local variance reduction via CADIS.
     max_realizations : int
         The upper limit for number of tally realizations when generating weight
         windows.
@@ -515,8 +518,10 @@ class WeightWindowGenerator:
         energies in [eV] for a single bin
     particle_type : {'neutron', 'photon'}
         Particle type the weight windows apply to
-    method : {'magic', 'fw_cadis'}
+    method : {'magic', 'fw_cadis', 'cadis'}
         The weight window generation methodology applied during an update.
+    targets : :class:`openmc.Tallies`
+        Target tallies for local variance reduction via CADIS.
     max_realizations : int
         The upper limit for number of tally realizations when generating weight
         windows.
@@ -528,14 +533,15 @@ class WeightWindowGenerator:
         Whether or not to apply weight windows on the fly.
     """
 
-    _MAGIC_PARAMS = {'value': str, 'threshold': float, 'ratio': float}
+    _WWG_PARAMS = {'value': str, 'threshold': float, 'ratio': float}
 
     def __init__(
         self,
         mesh: openmc.MeshBase,
         energy_bounds: Sequence[float] | None = None,
         particle_type: str = 'neutron',
         method: str = 'magic',
+        targets: openmc.Tallies | list[int] | None = None,
         max_realizations: int = 1,
         update_interval: int = 1,
         on_the_fly: bool = True
@@ -548,6 +554,7 @@ def __init__(
             self.energy_bounds = energy_bounds
         self.particle_type = particle_type
         self.method = method
+        self.targets = targets
         self.max_realizations = max_realizations
         self.update_interval = update_interval
         self.on_the_fly = on_the_fly
@@ -601,13 +608,33 @@ def method(self) -> str:
     @method.setter
     def method(self, m: str):
         cv.check_type('generation method', m, str)
-        cv.check_value('generation method', m, ('magic', 'fw_cadis'))
+        cv.check_value('generation method', m, ('magic', 'fw_cadis', 'cadis'))
         self._method = m
         if self._update_parameters is not None:
             try:
                 self._check_update_parameters()
             except (TypeError, KeyError):
                 warnings.warn(f'Update parameters are invalid for the "{m}" method.')
+
+    @property
+    def targets(self) -> openmc.Tallies:
+        return self._targets
+
+    @targets.setter
+    def targets(self, t):
+        if t is None:
+            self._targets = t
+        else:
+            cv.check_type('CADIS target tallies', t, [openmc.Tallies, list])
+            cv.check_greater_than('CADIS target tallies', len(t), 0)
+            if isinstance(t, list):
+                for tally_id in t:
+                    if not isinstance(tally_id, int):
+                        raise TypeError(
+                            "Tally IDs passed to WeightWindowGenerator " \
+                            "must be of type int")
+
+            self._targets = t
 
     @property
     def max_realizations(self) -> int:
@@ -634,14 +661,14 @@ def update_parameters(self) -> dict:
         return self._update_parameters
 
     def _check_update_parameters(self, params: dict):
-        if self.method == 'magic' or self.method == 'fw_cadis':
-            check_params = self._MAGIC_PARAMS
+        if self.method == 'magic' or self.method == 'fw_cadis'  or self.method == 'cadis':
+            check_params = self._WWG_PARAMS
 
         for key, val in params.items():
             if key not in check_params:
                 raise ValueError(f'Invalid param "{key}" for {self.method} '
                                   'weight window generation')
-            cv.check_type(f'weight window generation param: "{key}"', val, self._MAGIC_PARAMS[key])
+            cv.check_type(f'weight window generation param: "{key}"', val, self._WWG_PARAMS[key])
 
     @update_parameters.setter
     def update_parameters(self, params: dict):
@@ -677,8 +704,8 @@ def _sanitize_update_parameters(cls, method: str, update_parameters: dict):
         update_parameters : dict
             The update parameters as-read from the XML node (keys: str, values: str)
         """
-        if method == 'magic' or method == 'fw_cadis':
-            check_params = cls._MAGIC_PARAMS
+        if method == 'magic' or method == 'fw_cadis'  or method == 'cadis':
+            check_params = cls._WWG_PARAMS
 
         for param, param_type in check_params.items():
             if param in update_parameters:
@@ -704,6 +731,20 @@ def to_xml_element(self):
         otf_elem.text = str(self.on_the_fly).lower()
         method_elem = ET.SubElement(element, 'method')
         method_elem.text = self.method
+        if self.targets is not None:
+            if self.method != 'cadis':
+                raise ValueError(
+                    "CADIS update method is required in order to use " \
+                    "target tallies for WeightWindowGenerator.")
+            elif isinstance(self.targets, openmc.Tallies):
+                raise RuntimeError(
+                    "CADIS target tallies must be checked to ensure they are" \
+                    "present on model.tallies. Use model.export_to_xml() or " \
+                    "model.export_to_model_xml() to link CADIS target tallies.")
+            else:
+                targets_elem = ET.SubElement(element, 'targets')
+                targets_elem.text = ' '.join(str(tally_id) for tally_id in self.targets)
+
         if self.update_parameters is not None:
             self._update_parameters_subelement(element)
 
@@ -731,7 +772,7 @@ def from_xml_element(cls, elem: ET.Element, meshes: dict) -> Self:
         mesh_id = int(get_text(elem, 'mesh'))
         mesh = meshes[mesh_id]
 
-        energy_bounds = get_elem_list(elem, "energy_bounds, float")
+        energy_bounds = get_elem_list(elem, "energy_bounds", float)
         particle_type = get_text(elem, 'particle_type')
 
         wwg = cls(mesh, energy_bounds, particle_type)
@@ -740,6 +781,14 @@ def from_xml_element(cls, elem: ET.Element, meshes: dict) -> Self:
         wwg.update_interval = int(get_text(elem, 'update_interval'))
         wwg.on_the_fly = bool(get_text(elem, 'on_the_fly'))
         wwg.method = get_text(elem, 'method')
+        targets_elem = elem.find('targets')
+        if targets_elem is not None:
+            if wwg.method != 'cadis':
+                raise ValueError(
+                    "CADIS update method is required in order to use " \
+                    "target tallies for WeightWindowGenerator.")
+            else:
+                wwg.targets = get_elem_list(elem, "targets")
 
         if elem.find('update_parameters') is not None:
             update_parameters = {}