First Collision Source Method (FCS)

docs/source/methods/random_ray.rst

@@ -994,6 +994,114 @@ random ray and Monte Carlo, however.
   develop the scattering source by way of inactive batches before beginning
   active batches.
 
+
+.. _usersguide_fixed_source_first_collision_source_methods:
+
+----------------------------------
+First Collision Source Method
+----------------------------------
+
+In cases where the fixed source is not a well-defined volumetric source (e.g. point source),
+there is the need to preprocess the source into volumetric sources to be computed by the
+random ray solver. One possible way is through the First Collision Source Method (FCS), which works as a 
+preconditioner of the source, distributing the source contribution throughout the domain. This is
+not a new method (`Alcouffe <Alcouffe-1989_>`_), and has been recently implemented in other 
+neutron transport codes (`Ragusa <Ragusa-2016_>`_, `Falabino <Falabino-2022_>`_).
+
+The FCS works by generating uncollided neutron angular fluxes :math:`\psi^{\text{un}}_{g}` at the source that travels through
+the domain interacting with the media, being naturally attenuated in the process. Neutrons that experience
+any collision are treated as first collided neutrons and will be used to estimate the volumetric neutron 
+fixed source at that cell. Once the FCS preprocess stage is complete, the fixed volumetric source will be
+added to each iteration of the random ray solver. The remaining uncollided flux that has not interacted 
+at any point of the preprocessing stage is added to the final solution at the end of the neutron transport code.
+
+The FCS has a very similar mathematical formulation to the regular Random Ray method. The formulation for
+the method with flat sources will be provided. The neutron transport equation for the uncollided rays can be
+described as in Equation :eq:`moc_final`, however without any pre-existing source terms:
+
+.. math::
+    :label: fcs_moc_final
+
+    \psi_g^{\text{un}}(s, \mathbf{\Omega}) = \psi_g^{\text{un}}(\mathbf{r_0}, \mathbf{\Omega}) e^{-\int_0^s ds^\prime \Sigma_{t_g}(s^\prime)}.
+
+The analytical solution for this ODE is a simple attenuation problem:
+
+.. math::
+    :label: fcs_moc_sol
+
+    \psi_g^{\text{un}}(s) = \psi_g(0) e^{-\Sigma_{t,i,g} s} .
+
+For convenience, we can also write this equation in terms of the incoming and
+outgoing angular flux (:math:`\psi_g^{in}` and :math:`\psi_g^{out}`) 
+
+.. math::
+    :label: fcs_fsr_attenuation_in_out
+
+    \psi_g^{out} = \psi_g^{in} e^{-\Sigma_{t,i,g} \ell_r}.
+
+
+Equation :eq:`fcs_fsr_attenuation_in_out` can be rearranged in terms of :math:`\Delta \psi_{r,g}`:
+
+.. math::
+    :label: fcs_difference
+
+    \Delta\psi_{r,g}^{\text{un}} =\psi_{r,g}^{in} - \psi_{r,g}^{out} = \psi_{r,g}^{in}\big(1-e^{-\Sigma_{t,i,g}l_r}\big) .
+
+The average angular flux can be computed substitute :eq:`average` into :eq:`fcs_moc_final`
+and obtain:
+
+.. math::
+    :label: fcs_average_solved
+
+    \overline{\psi}_{r,i,g}^{\text{un}} = - \frac{\psi_{r,g}^{out} - \psi_{r,g}^{in}}{\ell_r \Sigma_{t,i,g}}.
+
+Which can also be described in terms of :math:`\Delta \psi_{r,g}`:
+
+.. math::
+    :label: fcs_average_solved_difference
+
+    \overline{\psi}_{r,i,g}^{\text{un}} = \frac{\Delta \psi_{r,g}^{\text{un}}}{\ell_r \Sigma_{t,i,g}}.
+
+Similarly to Equation :eq:`discretized`, the scalar flux in cell :math:`i` can be defined as the summation
+of the contributions of the angular fluxes traveling through it. However, in this method, there is the need
+to differentiate how the volume is estimated. In the Random Ray method, the rays are generated uniformly
+across the domain and the volume estimation is unbiased. Nonetheless, angular fluxes have a specific and non-well-distributed
+birth location.
+
+.. math::
+    :label: fcs_discretized
+
+    \phi^{\text{un}}(i,g) = \frac{\int_{V_i} \int_{4\pi} \psi(r, \Omega) d\Omega d\mathbf{r}}{\int_{V_i} d\mathbf{r}} = \overline{\overline{\psi}}_{i,g} \approx \frac{\sum\limits_{r=1}^{N_i} \ell_r \overline{\psi}_{r,i,g}^{un}}{V_i} .
+
+To avoid bias in the volume estimation, it is added an initial volume estimation stage
+to compute the volume of each cell before running FCS, using the same method presented in 
+the Random Ray method. 
+
+In the case of flat sources, the volume estimation can be improved with the random ray solver, 
+allowing a fast and rough initial estimate. However, linear sources require a more careful procedure.
+Linear source requires calculating angular flux moments, that depend on the estimated
+cell's spatial moments and centroids. If a poor estimation is made, inaccurate gradients will be
+carried throughout the calculation and impact the accuracy of the method. 
+
+Substituting Equation :eq:`fcs_average_solved_difference` 
+into Equation :eq:`fcs_discretized`, we obtain the equation for the uncollided angular flux:
+
+.. math::
+    :label: fcs_uncollided_flux
+
+    \phi^{\text{un}}(i,g) =  \frac{\sum\limits_{r=1}^{N_i} \Delta \psi_{r,i,g}^{\text{un}}}{\Sigma_{t,i,g} V_i} .
+
+The fixed-source term can be calculated as the first collided flux that undergoes scattering events:
+
+.. math::
+    :label: fcs_first_collided_flux
+
+        Q_\text{fixed}(i,g) = \phi^{\text{FCS}}(i,g) = \sum_{g'}^{G}\Sigma_s(i,g,g') \phi^{\text{un}} (i,g') .
+
+The fixed source :math:`Q_\text{fixed}` will be treated as a fixed volumetric source for all remaining Random Ray
+iterations, as presented in Eq :eq:`fixed_source_update`.
+
+
 ---------------------------
 Fundamental Sources of Bias
 ---------------------------
@@ -1048,6 +1156,9 @@ in random ray particle transport are:
 .. _Cosgrove-2023: https://doi.org/10.1080/00295639.2023.2270618
 .. _Ferrer-2016: https://doi.org/10.13182/NSE15-6
 .. _Gunow-2018: https://dspace.mit.edu/handle/1721.1/119030
+.. _Alcouffe-1989: https://doi.org/10.13182/NSE90-A23749
+.. _Ragusa-2016: https://www.osti.gov/biblio/1364492
+.. _Falabino-2022: https://doi.org/10.1016/j.jcp.2022.111156
 
 .. only:: html
 

docs/source/usersguide/random_ray.rst

@@ -475,6 +475,41 @@ as::
 which will greatly improve the quality of the linear source term in 2D
 simulations.
 
+-----------------------------
+First Collision Source Method
+-----------------------------
+
+The First Collision Source method (FCS) is supported within the fixed source random ray
+solver. The preprocessing stage can be toggle by setting the ``first_collision_source``
+field in the :attr:`openmc.Settings.random_ray` dictionary to ``True`` as::
+
+    settings.random_ray['first_collision_source'] = True
+
+The method has two input variables: the number of rays for initial volume estimation
+and number of uncollided rays used in FCS.
+
+The user can define the amount of uncollided rays used in the FCS mode by setting the
+``first_collision_rays`` field in the :attr:`openmc.Settings.random_ray` dictionary 
+to ``value`` as::
+
+        settings.random_ray['first_collision_rays'] = 1000
+
+If not provided, the solver will run the default procedure to iteratively generate more
+rays until:
+
+* All cells are hit.
+* No extra cells were hit with additional rays.
+* Maximum value reached (default :math:`= 1e7` rays).
+
+The user can also define the amount of rays for initial volume estimation by setting the
+``first_collision_volume_rays`` field in the :attr:`openmc.Settings.random_ray` dictionary 
+to ``value`` as::
+
+        settings.random_ray['first_collision_volume_rays'] = 2000
+
+If not provided, the solver will use the same amount of rays used for a regular batch in
+the random ray solver (:attr:`settings::n_particles`).
+
 ---------------------------------
 Fixed Source and Eigenvalue Modes
 ---------------------------------

include/openmc/distribution.h

@@ -1,4 +1,4 @@
 //! \file distribution.h
 //! Univariate probability distributions
 
 #ifndef OPENMC_DISTRIBUTION_H
@@ -56,13 +56,15 @@
   // Properties
   const vector<double>& prob() const { return prob_; }
   const vector<size_t>& alias() const { return alias_; }
+  const vector<double>& prob_actual() const { return prob_actual_; }
   double integral() const { return integral_; }
 
 private:
   vector<double> prob_; //!< Probability of accepting the uniformly sampled bin,
                         //!< mapped to alias method table
   vector<size_t> alias_; //!< Alias table
   double integral_;      //!< Integral of distribution
+  vector<double> prob_actual_; //!< actual probability before the Vose algorithm
 
   //! Normalize distribution so that probabilities sum to unity
   void normalize();
@@ -91,6 +93,7 @@
   const vector<double>& x() const { return x_; }
   const vector<double>& prob() const { return di_.prob(); }
   const vector<size_t>& alias() const { return di_.alias(); }
+  const vector<double>& prob_actual() const { return di_.prob_actual(); }
 
 private:
   vector<double> x_; //!< Possible outcomes

include/openmc/particle_data.h

@@ -51,6 +51,7 @@ struct SourceSite {
   ParticleType particle;
   int64_t parent_id;
   int64_t progeny_id;
+  int source_id;
 };
 
 //! State of a particle used for particle track files

include/openmc/random_ray/flat_source_domain.h

@@ -113,6 +113,19 @@ class FlatSourceDomain {
   virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const;
   double compute_fixed_source_normalization_factor() const;
 
+  virtual void compute_first_collided_flux();
+  virtual void normalize_uncollided_scalar_flux(double number_of_particles);
+  virtual void update_volume_uncollided_flux();
+  virtual void compute_first_collided_external_source();
+  virtual void compute_uncollided_scalar_flux();
+  int64_t check_fsr_hits();
+  virtual void uncollided_moments();
+  virtual void batch_reset_fc();
+  virtual double evaluate_uncollided_flux_at_point(
+    Position r, int64_t sr, int g) const;
+  virtual double evaluate_external_source_at_point(
+    Position r, int64_t sr, int g) const;
+
   //----------------------------------------------------------------------------
   // Static Data members
   static bool volume_normalized_flux_tallies_;
@@ -144,6 +157,8 @@ class FlatSourceDomain {
   vector<float> scalar_flux_new_;
   vector<float> source_;
   vector<float> external_source_;
+  vector<float> scalar_uncollided_flux_;
+  vector<float> scalar_first_collided_flux_;
 
 protected:
   //----------------------------------------------------------------------------

include/openmc/random_ray/linear_source_domain.h

@@ -41,13 +41,30 @@ class LinearSourceDomain : public FlatSourceDomain {
   void flux_swap() override;
   double evaluate_flux_at_point(Position r, int64_t sr, int g) const override;
 
+  void compute_first_collided_external_source() override;
+  void compute_uncollided_scalar_flux() override;
+  void compute_first_collided_flux() override;
+  void normalize_uncollided_scalar_flux(double number_of_particles) override;
+  void update_volume_uncollided_flux() override;
+  void uncollided_moments() override;
+  void batch_reset_fc() override;
+  double evaluate_uncollided_flux_at_point(
+    Position r, int64_t sr, int g) const override;
+  double evaluate_external_source_at_point(
+    Position r, int64_t sr, int g) const override;
+
   //----------------------------------------------------------------------------
   // Public Data members
 
   vector<MomentArray> source_gradients_;
   vector<MomentArray> flux_moments_old_;
   vector<MomentArray> flux_moments_new_;
   vector<MomentArray> flux_moments_t_;
+  // First Collided Method
+  vector<MomentArray> flux_moments_uncollided_;
+  vector<MomentArray> flux_moments_first_collided_;
+  vector<MomentArray> external_source_gradients_;
+
   vector<Position> centroid_;
   vector<Position> centroid_iteration_;
   vector<Position> centroid_t_;

include/openmc/random_ray/random_ray.h

@@ -1,4 +1,4 @@
 #ifndef OPENMC_RANDOM_RAY_H
 #define OPENMC_RANDOM_RAY_H
 
 #include "openmc/memory.h"
@@ -20,28 +20,35 @@
   //----------------------------------------------------------------------------
   // Constructors
   RandomRay();
-  RandomRay(uint64_t ray_id, FlatSourceDomain* domain);
+  RandomRay(uint64_t ray_id, FlatSourceDomain* domain, bool uncollided_ray);
 
   //----------------------------------------------------------------------------
   // Methods
   void event_advance_ray();
   void attenuate_flux(double distance, bool is_active);
   void attenuate_flux_flat_source(double distance, bool is_active);
   void attenuate_flux_linear_source(double distance, bool is_active);
-
-  void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain);
+  void event_advance_ray_first_collided();
+  void initialize_ray(uint64_t ray_id, FlatSourceDomain* domain,bool uncollided_ray);
   uint64_t transport_history_based_single_ray();
-
+  
   //----------------------------------------------------------------------------
   // Static data members
   static double distance_inactive_;          // Inactive (dead zone) ray length
   static double distance_active_;            // Active ray length
   static unique_ptr<Source> ray_source_;     // Starting source for ray sampling
   static RandomRaySourceShape source_shape_; // Flag for linear source
 
+  static bool first_collision_source_;
+  static int first_collision_rays_;
+  static int first_collision_volume_rays_;
+
+  static bool no_volume_calc; // Flag for FCS flux calculation 
+
   //----------------------------------------------------------------------------
   // Public data members
   vector<float> angular_flux_;
+  vector<float> angular_flux_initial_;
 
 private:
   //----------------------------------------------------------------------------
@@ -50,6 +57,7 @@
   vector<MomentArray> delta_moments_;
 
   int negroups_;
+
   FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source
                                        // data needed for ray transport
   double distance_travelled_ {0};

include/openmc/random_ray/random_ray_simulation.h

@@ -1,4 +1,4 @@
 #ifndef OPENMC_RANDOM_RAY_SIMULATION_H
 #define OPENMC_RANDOM_RAY_SIMULATION_H
 
 #include "openmc/random_ray/flat_source_domain.h"
@@ -27,9 +27,13 @@
   void print_results_random_ray(uint64_t total_geometric_intersections,
     double avg_miss_rate, int negroups, int64_t n_source_regions,
     int64_t n_external_source_regions) const;
+  void first_collision_source_simulation();
 
   //----------------------------------------------------------------------------
   // Data members
+  // Initial volume estimation timer - First Collided Source Method  
+  double time_volume_fc;
+
 private:
   // Contains all flat source region data
   unique_ptr<FlatSourceDomain> domain_;

include/openmc/timer.h

@@ -32,6 +32,7 @@ extern Timer time_event_surface_crossing;
 extern Timer time_event_collision;
 extern Timer time_event_death;
 extern Timer time_update_src;
+extern Timer time_first_collided;
 
 } // namespace simulation
 

openmc/__init__.py

@@ -40,5 +40,7 @@
 
 from . import examples
 
+#__version__ = '0.14.1-dev'
+
 
 __version__ = importlib.metadata.version("openmc")

openmc/lib/core.py

@@ -26,7 +26,8 @@ class _SourceSite(Structure):
                 ('surf_id', c_int),
                 ('particle', c_int),
                 ('parent_id', c_int64),
-                ('progeny_id', c_int64)]
+                ('progeny_id', c_int64),
+                ('source_id', c_int)]
 
 
 # Define input type for numpy arrays that will be passed into C++ functions

openmc/settings.py

@@ -164,6 +164,19 @@ 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).
+        :first_collision_source:
+            Boolean that indicates if first collided source method (FSC) is 
+            used to initialize the external source input. Default is 'False'.
+        :first_collision_rays:
+            Number of rays (int) used to generate the first collided source.
+            If not provided, the method will automatically run enough rays to
+            adequately converge the first collided source, via an iterative 
+            doubling process. If provided, the method will run the exact
+            prescribed amount.
+        :first_collision_volume_rays:
+            Number of rays (int) used to estimate the initial volume for the
+            first collied source calculation. If not provided, it will use the
+            same number of rays as the main random ray solver stage.
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -1096,6 +1109,14 @@ 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 == 'first_collision_source':
+                cv.check_type('first_collision_source', random_ray[key], bool)
+            elif key == 'first_collision_rays':
+                cv.check_type('first_collision_rays', random_ray[key], int)
+                cv.check_greater_than('first_collision_rays',random_ray[key], 0)
+            elif key == 'first_collision_volume_rays':
+                cv.check_type('first_collision_volume_rays', random_ray[key], int)
+                cv.check_greater_than('first_collision_volume_rays',random_ray[key], 0)                                
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
@@ -1895,6 +1916,12 @@ def _random_ray_from_xml_element(self, root):
                     self.random_ray['volume_normalized_flux_tallies'] = (
                         child.text in ('true', '1')
                     )
+                elif child.tag == 'first_collision_source':
+                    self.random_ray['first_collision_source'] = (
+                        child.text in ('true', '1')
+                    )
+                elif child.tag in ('first_collision_rays', 'first_collision_volume_rays'):
+                    self.random_ray[child.tag] = int(child.text) 
 
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.