Merge pull request #1974 from paulromano/survival-biasing-fix

Remove wgt_absorb data member from Particle

include/openmc/particle_data.h

@@ -244,7 +244,6 @@ class ParticleData {
   Position r_last_;         //!< previous coordinates
   Direction u_last_;        //!< previous direction coordinates
   double wgt_last_ {1.0};   //!< pre-collision particle weight
-  double wgt_absorb_ {0.0}; //!< weight absorbed for survival biasing
 
   // What event took place
   bool fission_ {false};  //!< did particle cause implicit fission
@@ -375,8 +374,6 @@ class ParticleData {
   const Position& u_last() const { return u_last_; }
   double& wgt_last() { return wgt_last_; }
   const double& wgt_last() const { return wgt_last_; }
-  double& wgt_absorb() { return wgt_absorb_; }
-  const double& wgt_absorb() const { return wgt_absorb_; }
 
   bool& fission() { return fission_; }
   TallyEvent& event() { return event_; }

include/openmc/physics_common.h

@@ -9,7 +9,9 @@
 namespace openmc {
 
 //! \brief Performs the russian roulette operation for a particle
-void russian_roulette(Particle& p);
+//! \param[in,out] p  Particle object
+//! \param[in] weight_survive Weight assigned to particles that survive
+void russian_roulette(Particle& p, double weight_survive);
 
 } // namespace openmc
 #endif // OPENMC_PHYSICS_COMMON_H

openmc/surface.py

@@ -2633,7 +2633,7 @@ def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False,
             memo = {}
 
         # If rotated surface not in memo, add it
-        key = (self.surface, tuple(rotation), tuple(pivot), order, inplace)
+        key = (self.surface, tuple(np.ravel(rotation)), tuple(pivot), order, inplace)
         if key not in memo:
             memo[key] = self.surface.rotate(rotation, pivot=pivot, order=order,
                                             inplace=inplace)

src/physics.cpp

@@ -134,8 +134,6 @@ void sample_neutron_reaction(Particle& p)
 
   if (p.neutron_xs(i_nuclide).absorption > 0.0) {
     absorption(p, i_nuclide);
-  } else {
-    p.wgt_absorb() = 0.0;
   }
   if (!p.alive())
     return;
@@ -153,9 +151,9 @@ void sample_neutron_reaction(Particle& p)
 
   // Play russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    russian_roulette(p);
-    if (!p.alive())
-      return;
+    if (p.wgt() < settings::weight_cutoff) {
+      russian_roulette(p, settings::weight_survive);
+    }
   }
 }
 
@@ -626,16 +624,15 @@ void absorption(Particle& p, int i_nuclide)
 {
   if (settings::survival_biasing) {
     // Determine weight absorbed in survival biasing
-    p.wgt_absorb() = p.wgt() * p.neutron_xs(i_nuclide).absorption /
-                     p.neutron_xs(i_nuclide).total;
+    const double wgt_absorb = p.wgt() * p.neutron_xs(i_nuclide).absorption /
+                              p.neutron_xs(i_nuclide).total;
 
     // Adjust weight of particle by probability of absorption
-    p.wgt() -= p.wgt_absorb();
-    p.wgt_last() = p.wgt();
+    p.wgt() -= wgt_absorb;
 
     // Score implicit absorption estimate of keff
     if (settings::run_mode == RunMode::EIGENVALUE) {
-      p.keff_tally_absorption() += p.wgt_absorb() *
+      p.keff_tally_absorption() += wgt_absorb *
                                    p.neutron_xs(i_nuclide).nu_fission /
                                    p.neutron_xs(i_nuclide).absorption;
     }
@@ -952,19 +949,17 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
 
       // cdf value at upper bound attainable energy
       double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) /
-          (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]);
-      double cdf_up =
-        nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]);
+                 (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]);
+      double cdf_up = nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]);
 
       while (true) {
         // directly sample Maxwellian
         double E_t = -kT * std::log(prn(seed));
 
         // sample a relative energy using the xs cdf
         double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low);
-        int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low, 
-                                        nuc.xs_cdf_.begin() + i_E_up+2, 
-                                        cdf_rel);
+        int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low,
+          nuc.xs_cdf_.begin() + i_E_up + 2, cdf_rel);
         double E_rel = nuc.energy_0K_[i_E_low + i_E_rel];
         double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] -
                      nuc.xs_cdf_[i_E_low + i_E_rel]) /

src/physics_common.cpp

@@ -9,17 +9,13 @@ namespace openmc {
 // RUSSIAN_ROULETTE
 //==============================================================================
 
-void russian_roulette(Particle& p)
+void russian_roulette(Particle& p, double weight_survive)
 {
-  if (p.wgt() < settings::weight_cutoff) {
-    if (prn(p.current_seed()) < p.wgt() / settings::weight_survive) {
-      p.wgt() = settings::weight_survive;
-      p.wgt_last() = p.wgt();
-    } else {
-      p.wgt() = 0.;
-      p.wgt_last() = 0.;
-      p.alive() = false;
-    }
+  if (weight_survive * prn(p.current_seed()) < p.wgt()) {
+    p.wgt() = weight_survive;
+  } else {
+    p.wgt() = 0.;
+    p.alive() = false;
   }
 }
 

src/physics_mg.cpp

@@ -55,8 +55,6 @@ void sample_reaction(Particle& p)
   // weight of the particle. Otherwise, it checks to see if absorption occurs.
   if (p.macro_xs().absorption > 0.) {
     absorption(p);
-  } else {
-    p.wgt_absorb() = 0.;
   }
   if (!p.alive())
     return;
@@ -66,9 +64,9 @@ void sample_reaction(Particle& p)
 
   // Play Russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
-    russian_roulette(p);
-    if (!p.alive())
-      return;
+    if (p.wgt() < settings::weight_cutoff) {
+      russian_roulette(p, settings::weight_survive);
+    }
   }
 }
 
@@ -216,15 +214,14 @@ void absorption(Particle& p)
 {
   if (settings::survival_biasing) {
     // Determine weight absorbed in survival biasing
-    p.wgt_absorb() = p.wgt() * p.macro_xs().absorption / p.macro_xs().total;
+    double wgt_absorb = p.wgt() * p.macro_xs().absorption / p.macro_xs().total;
 
     // Adjust weight of particle by the probability of absorption
-    p.wgt() -= p.wgt_absorb();
-    p.wgt_last() = p.wgt();
+    p.wgt() -= wgt_absorb;
 
     // Score implicit absorpion estimate of keff
     p.keff_tally_absorption() +=
-      p.wgt_absorb() * p.macro_xs().nu_fission / p.macro_xs().absorption;
+      wgt_absorb * p.macro_xs().nu_fission / p.macro_xs().absorption;
   } else {
     if (p.macro_xs().absorption > prn(p.current_seed()) * p.macro_xs().total) {
       p.keff_tally_absorption() +=