Convert sample_element to C++ and remove physics.F90

CMakeLists.txt

@@ -332,7 +332,6 @@ add_library(libopenmc SHARED
   src/particle_restart.F90
   src/photon_header.F90
   src/physics_common.F90
-  src/physics.F90
   src/pugixml/pugixml_f.F90
   src/random_lcg.F90
   src/reaction_header.F90

include/openmc/physics.h

@@ -44,7 +44,7 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat);
 void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
   Bank* bank_array, int64_t* bank_size, int64_t bank_capacity);
 
-extern "C" int sample_element(Particle* p);
+int sample_element(Particle* p);
 
 Reaction* sample_fission(int i_nuclide, double E);
 

src/physics.F90

@@ -1,69 +0,0 @@
-module physics
-
-  use constants
-  use error,                  only: fatal_error
-  use material_header,        only: Material, materials
-  use math
-  use message_passing
-  use nuclide_header
-  use particle_header
-  use photon_header
-  use random_lcg,             only: prn
-  use settings
-  use simulation_header
-
-  implicit none
-
-  interface
-    subroutine collision(p) bind(C)
-      import Particle
-      type(Particle), intent(inout) :: p
-    end subroutine
-  end interface
-
-contains
-
-!===============================================================================
-! SAMPLE_ELEMENT
-!===============================================================================
-
-  function sample_element(p) result(i_element) bind(C)
-    type(Particle), intent(in) :: p
-    integer(C_INT) :: i_element
-
-    integer :: i
-    real(8) :: prob
-    real(8) :: cutoff
-    real(8) :: atom_density ! atom density of nuclide in atom/b-cm
-    real(8) :: sigma        ! microscopic total xs for nuclide
-
-    associate (mat => materials(p % material))
-      ! Sample cumulative distribution function
-      cutoff = prn() * material_xs % total
-
-      i = 0
-      prob = ZERO
-      do while (prob < cutoff)
-        i = i + 1
-
-        ! Check to make sure that a nuclide was sampled
-        if (i > mat % n_nuclides) then
-          call particle_write_restart(p)
-          call fatal_error("Did not sample any element during collision.")
-        end if
-
-        ! Find atom density
-        i_element    = mat % element(i)
-        atom_density = mat % atom_density(i)
-
-        ! Determine microscopic cross section
-        sigma = atom_density * micro_photon_xs(i_element) % total
-
-        ! Increment probability to compare to cutoff
-        prob = prob + sigma
-      end do
-    end associate
-
-  end function sample_element
-
-end module physics

src/physics.cpp

@@ -229,10 +229,10 @@ void sample_photon_reaction(Particle* p)
 
   // Sample element within material
   int i_element = sample_element(p);
-  p->event_nuclide = i_element;
   // TODO: off-by-one
-  const auto& micro {simulation::micro_photon_xs[i_element - 1]};
-  const auto& element {data::elements[i_element - 1]};
+  p->event_nuclide = i_element + 1;
+  const auto& micro {simulation::micro_photon_xs[i_element]};
+  const auto& element {data::elements[i_element]};
 
   // Calculate photon energy over electron rest mass equivalent
   double alpha = p->E/MASS_ELECTRON_EV;
@@ -479,32 +479,43 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat)
   }
 }
 
-// int sample_element(Particle* p)
-// {
-//   // Sample cumulative distribution function
-//   double cutoff = prn() * simulation::material_xs.total;
+int sample_element(Particle* p)
+{
+  // Sample cumulative distribution function
+  double cutoff = prn() * simulation::material_xs.total;
 
-//   int i = 0;
-//   double prob = 0.0;
-//   while (prob < cutoff) {
-//     // Check to make sure that a nuclide was sampled
-//     if (i > mat % n_nuclides) {
-//       p->write_restart();
-//       fatal_error("Did not sample any element during collision.");
-//     }
+  // Get pointers to elements, densities
+  int* nuclide;
+  double* density;
+  int n;
+  openmc_material_get_densities(p->material, &nuclide, &density, &n);
+  int* element = material_element(p->material);
 
-//     // Find atom density
-//     int i_element = mat % element(i);
-//     double atom_density = mat % atom_density(i);
+  int i = 0;
+  double prob = 0.0;
+  int i_element;
+  while (prob < cutoff) {
+    // Check to make sure that a nuclide was sampled
+    if (i >= n) {
+      p->write_restart();
+      fatal_error("Did not sample any element during collision.");
+    }
 
-//     // Determine microscopic cross section
-//     double sigma = atom_density * simulation::micro_photon_xs[i_element].total;
+    // Find atom density
+    // TODO: off-by-one
+    i_element = element[i] - 1;
+    double atom_density = density[i];
 
-//     // Increment probability to compare to cutoff
-//     prob += sigma;
-//     ++i;
-//   }
-// }
+    // Determine microscopic cross section
+    double sigma = atom_density * simulation::micro_photon_xs[i_element].total;
+
+    // Increment probability to compare to cutoff
+    prob += sigma;
+    ++i;
+  }
+
+  return i_element;
+}
 
 Reaction* sample_fission(int i_nuclide, double E)
 {

src/tracking.F90

@@ -16,7 +16,6 @@ module tracking
   use mgxs_interface
   use nuclide_header
   use particle_header
-  use physics,            only: collision
   use random_lcg,         only: prn, prn_set_stream
   use settings
   use simulation_header
@@ -33,11 +32,15 @@ module tracking
   implicit none
 
   interface
+    subroutine collision(p) bind(C)
+      import Particle
+      type(Particle), intent(inout) :: p
+    end subroutine
+
     subroutine collision_mg(p) bind(C)
       import Particle, C_DOUBLE
       type(Particle),            intent(inout) :: p
     end subroutine collision_mg
-
   end interface
 
 contains