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
@@ -84,6 +86,8 @@ class FlatSourceDomain {
   static std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
     mesh_domain_map_;
 
+  static std::vector<size_t> cadis_targets_;
+
   //----------------------------------------------------------------------------
   // Static data members
   static RandomRayVolumeEstimator volume_estimator_;

include/openmc/random_ray/random_ray_simulation.h

@@ -19,9 +19,10 @@ class RandomRaySimulation {
 
   //----------------------------------------------------------------------------
   // Methods
+  void compute_segment_correction_factors();
   void apply_fixed_sources_and_mesh_domains();
-  void prepare_fixed_sources_adjoint();
-  void prepare_adjoint_simulation();
+  void prepare_fw_fixed_sources_adjoint();
+  void prepare_local_fixed_sources_adjoint();
   void simulate();
   void output_simulation_results() const;
   void instability_check(
@@ -34,15 +35,9 @@ class RandomRaySimulation {
   // Accessors
   FlatSourceDomain* domain() const { return domain_.get(); }
 
-  //----------------------------------------------------------------------------
-  // Public data members
-
-  // Flag for adjoint simulation;
-  bool adjoint_needed_;
-
 private:
   //----------------------------------------------------------------------------
-  // Private data members
+  // Data members
 
   // Contains all flat source region data
   unique_ptr<FlatSourceDomain> domain_;
@@ -57,9 +52,6 @@ class RandomRaySimulation {
   // Number of energy groups
   int negroups_;
 
-  // Toggle for first simulation
-  bool is_first_simulation_;
-
 }; // class RandomRaySimulation
 
 //============================================================================
@@ -69,7 +61,6 @@ class RandomRaySimulation {
 void openmc_run_random_ray();
 void validate_random_ray_inputs();
 void openmc_reset_random_ray();
-void print_adjoint_header();
 
 } // namespace openmc
 

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
@@ -235,6 +235,9 @@ class WeightWindowsGenerator {
   double threshold_ {1.0}; //<! Relative error threshold for values used to
                            // update weight windows
   double ratio_ {5.0};     //<! ratio of lower to upper weight window bounds
+
+  // CADIS target tallies
+  std::vector<size_t> targets_;
 };
 
 //! Finalize variance reduction objects after all inputs have been read

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

@@ -225,6 +225,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
@@ -1356,6 +1359,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')
@@ -1887,11 +1892,12 @@ 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()
                     if source_element.find('bias') is not None:
                         raise RuntimeError(
                             "Ray source distributions should not be biased.")
-                    element.append(source_element)
+                    subelement.append(source_element)
 
                 elif key == 'source_region_meshes':
                     subelement = ET.SubElement(element, 'source_region_meshes')
@@ -1909,8 +1915,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()
@@ -2337,8 +2347,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)
                     if child.find('bias') is not None:
                         raise RuntimeError(
                             "Ray source distributions should not be biased.")
@@ -2355,6 +2366,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':