Simplify interface for sample_nuclide, fix iso_in_lab test

include/openmc/physics.h

@@ -37,7 +37,12 @@ void sample_electron_reaction(Particle* p);
 //! MeV) are created and travel in opposite directions.
 void sample_positron_reaction(Particle* p);
 
-void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat);
+//! Sample a nuclide based on their total cross sections and densities within
+//! the current material
+//!
+//! \param[in] p Particle
+//! \return Index in the data::nuclides vector
+int sample_nuclide(const Particle* p);
 
 //! Determine the average total, prompt, and delayed neutrons produced from
 //! fission and creates appropriate bank sites.
@@ -52,7 +57,7 @@ void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
 
 void absorption(Particle* p, int i_nuclide);
 
-void scatter(Particle*, int i_nuclide, int i_nuc_mat);
+void scatter(Particle*, int i_nuclide);
 
 //! Treats the elastic scattering of a neutron with a target.
 void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,

src/material.cpp

@@ -203,7 +203,7 @@ Material::Material(pugi::xml_node node)
   atom_density_ = xt::empty<double>({n});
   if (settings::photon_transport) element_.reserve(n);
 
-  for (int i = 0; i < names.size(); ++i) {
+  for (int i = 0; i < n; ++i) {
     const auto& name {names[i]};
 
     // Check that this nuclide is listed in the cross_sections.xml file
@@ -244,8 +244,7 @@ Material::Material(pugi::xml_node node)
       }
     }
 
-    // Copy name and atom/weight percent
-    //mat % names(j) = name
+    // Copy atom/weight percent
     atom_density_(i) = densities[i];
   }
 
@@ -255,9 +254,12 @@ Material::Material(pugi::xml_node node)
       p0_.resize(n);
 
       // Apply isotropic-in-lab treatment to specified nuclides
-      for (const auto& nuc : iso_lab) {
-        for (int j = 0; j < names.size(); ++j) {
-          if (names[j] == nuc) p0_[j] = true;
+      for (int j = 0; j < n; ++j) {
+        for (const auto& nuc : iso_lab) {
+          if (names[j] == nuc) {
+            p0_[j] = true;
+            break;
+          }
         }
       }
     }

src/physics.cpp

@@ -74,9 +74,8 @@ void collision(Particle* p)
 
 void sample_neutron_reaction(Particle* p)
 {
-  int i_nuclide;
-  int i_nuc_mat;
-  sample_nuclide(p, SCORE_TOTAL, &i_nuclide, &i_nuc_mat);
+  // Sample a nuclide within the material
+  int i_nuclide = sample_nuclide(p);
 
   // Save which nuclide particle had collision with
   // TODO: off-by-one
@@ -120,7 +119,7 @@ void sample_neutron_reaction(Particle* p)
 
   // Sample a scattering reaction and determine the secondary energy of the
   // exiting neutron
-  scatter(p, i_nuclide, i_nuc_mat);
+  scatter(p, i_nuclide);
 
   // Advance URR seed stream 'N' times after energy changes
   if (p->E != p->last_E) {
@@ -420,61 +419,30 @@ void sample_positron_reaction(Particle* p)
   p->alive = false;
 }
 
-void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat)
+int sample_nuclide(const Particle* p)
 {
   // Sample cumulative distribution function
-  double cutoff;
-  switch (mt) {
-  case SCORE_TOTAL:
-    cutoff = prn() * simulation::material_xs.total;
-    break;
-  case SCORE_SCATTER:
-    cutoff = prn() * (simulation::material_xs.total -
-      simulation::material_xs.absorption);
-    break;
-  case SCORE_FISSION:
-    cutoff = prn() * simulation::material_xs.fission;
-    break;
-  }
+  double cutoff = prn() * simulation::material_xs.total;
 
   // Get pointers to nuclide/density arrays
   // TODO: off-by-one
   const auto& mat {model::materials[p->material - 1]};
   int n = mat->nuclide_.size();
 
-  *i_nuc_mat = 0;
   double prob = 0.0;
-  while (prob < cutoff) {
-    // Check to make sure that a nuclide was sampled
-    if (*i_nuc_mat >= n) {
-      p->write_restart();
-      fatal_error("Did not sample any nuclide during collision.");
-    }
-
+  for (int i = 0; i < n; ++i) {
     // Get atom density
-    *i_nuclide = mat->nuclide_[*i_nuc_mat];
-    double atom_density = mat->atom_density_[*i_nuc_mat];
-
-    // Determine microscopic cross section
-    double sigma;
-    switch (mt) {
-    case SCORE_TOTAL:
-      sigma = atom_density * simulation::micro_xs[*i_nuclide].total;
-      break;
-    case SCORE_SCATTER:
-      sigma = atom_density * (simulation::micro_xs[*i_nuclide].total -
-        simulation::micro_xs[*i_nuclide].absorption);
-        break;
-    case SCORE_FISSION:
-      sigma = atom_density * simulation::micro_xs[*i_nuclide].fission;
-      break;
-    }
+    int i_nuclide = mat->nuclide_[i];
+    double atom_density = mat->atom_density_[i];
 
     // Increment probability to compare to cutoff
-    prob += sigma;
-
-    ++(*i_nuc_mat);
+    prob += atom_density * simulation::micro_xs[i_nuclide].total;
+    if (prob >= cutoff) return i_nuclide;
   }
+
+  // If we reach here, no nuclide was sampled
+  p->write_restart();
+  throw std::runtime_error{"Did not sample any nuclide during collision."};
 }
 
 int sample_element(Particle* p)
@@ -621,7 +589,7 @@ void absorption(Particle* p, int i_nuclide)
   }
 }
 
-void scatter(Particle* p, int i_nuclide, int i_nuc_mat)
+void scatter(Particle* p, int i_nuclide)
 {
   // copy incoming direction
   Direction u_old {p->coord[0].uvw};
@@ -709,6 +677,7 @@ void scatter(Particle* p, int i_nuclide, int i_nuc_mat)
   // TODO: off-by-one
   const auto& mat {model::materials[p->material - 1]};
   if (!mat->p0_.empty()) {
+    int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
     if (mat->p0_[i_nuc_mat]) {
       // Sample isotropic-in-lab outgoing direction
       double mu = 2.0*prn() - 1.0;