Use C++ based count_sites for CMFD. Fix several bugs

CMakeLists.txt

@@ -380,6 +380,7 @@ add_library(libopenmc SHARED
   src/tallies/trigger.F90
   src/tallies/trigger_header.F90
   src/cell.cpp
+  src/cmfd_execute.cpp
   src/distribution.cpp
   src/distribution_angle.cpp
   src/distribution_energy.cpp

src/cmfd_execute.F90

@@ -223,7 +223,7 @@ contains
     integer :: nx       ! maximum number of cells in x direction
     integer :: ny       ! maximum number of cells in y direction
     integer :: nz       ! maximum number of cells in z direction
-    integer :: ng       ! maximum number of energy groups
+    integer(C_INT) :: ng       ! maximum number of energy groups
     integer :: i        ! iteration counter
     integer :: g        ! index for group
     integer :: ijk(3)   ! spatial bin location
@@ -231,12 +231,22 @@ contains
     integer :: mesh_bin ! mesh bin of soruce particle
     integer :: n_groups ! number of energy groups
     real(8) :: norm     ! normalization factor
-    logical :: outside  ! any source sites outside mesh
+    logical(C_BOOL) :: outside  ! any source sites outside mesh
     logical :: in_mesh  ! source site is inside mesh
 #ifdef OPENMC_MPI
     integer :: mpi_err
 #endif
 
+    interface
+      subroutine cmfd_populate_sourcecounts(ng, energies, source_counts, outside) bind(C)
+        import C_INT, C_DOUBLE, C_BOOL
+        integer(C_INT), value :: ng
+        real(C_DOUBLE), intent(in) :: energies
+        real(C_DOUBLE), intent(out) :: source_counts
+        logical(C_BOOL), intent(out) :: outside
+      end subroutine
+    end interface
+
     ! Get maximum of spatial and group indices
     nx = cmfd % indices(1)
     ny = cmfd % indices(2)
@@ -260,8 +270,8 @@ contains
       cmfd%weightfactors = ONE
 
       ! Count bank sites in mesh and reverse due to egrid structure
-      call count_bank_sites(cmfd_mesh, source_bank, cmfd%sourcecounts, &
-           cmfd % egrid, sites_outside=outside, size_bank=work)
+      call cmfd_populate_sourcecounts(ng + 1, cmfd % egrid(1), &
+           cmfd % sourcecounts(1,1), outside)
 
       ! Check for sites outside of the mesh
       if (master .and. outside) then

src/cmfd_execute.cpp

@@ -0,0 +1,34 @@
+#include <algorithm> // for copy
+#include <cstdint>
+#include <iostream>
+
+#include "xtensor/xarray.hpp"
+#include "xtensor/xio.hpp"
+
+#include "openmc/capi.h"
+#include "openmc/mesh.h"
+
+namespace openmc {
+
+extern "C" int index_cmfd_mesh;
+
+extern "C" void
+cmfd_populate_sourcecounts(int n_energy, const double* energies,
+  double* source_counts, bool* outside)
+{
+  // Get pointer to source bank
+  Bank* source_bank;
+  int64_t n;
+  openmc_source_bank(&source_bank, &n);
+
+  // Get source counts in each mesh bin / energy bin
+  auto& m = meshes.at(index_cmfd_mesh);
+  xt::xarray<double> counts = m->count_sites(openmc_work, source_bank, n_energy, energies, outside);
+
+  std::cout << counts << "\n";
+
+  // Copy data from the xarray into the source counts array
+  std::copy(counts.begin(), counts.end(), source_counts);
+}
+
+} // namespace openmc

src/cmfd_header.F90

@@ -1,5 +1,7 @@
 module cmfd_header
 
+  use, intrinsic :: ISO_C_BINDING
+
   use constants,  only: CMFD_NOACCEL, ZERO, ONE
   use mesh_header, only: RegularMesh
   use set_header, only: SetInt
@@ -24,7 +26,7 @@ module cmfd_header
     integer :: mat_dim = CMFD_NOACCEL
 
     ! Energy grid
-    real(8), allocatable :: egrid(:)
+    real(C_DOUBLE), allocatable :: egrid(:)
 
     ! Cross sections
     real(8), allocatable :: totalxs(:,:,:,:)
@@ -53,7 +55,7 @@ module cmfd_header
     real(8), allocatable :: openmc_src(:,:,:,:)
 
     ! Source sites in each mesh box
-    real(8), allocatable :: sourcecounts(:,:)
+    real(C_DOUBLE), allocatable :: sourcecounts(:,:)
 
     ! Weight adjustment factors
     real(8), allocatable :: weightfactors(:,:,:,:)
@@ -95,7 +97,7 @@ module cmfd_header
   ! Main object
   type(cmfd_type), public :: cmfd
 
-  integer, public :: index_cmfd_mesh
+  integer(C_INT), public, bind(C) :: index_cmfd_mesh
   type(RegularMesh), public :: cmfd_mesh
 
   ! Pointers for different tallies

src/eigenvalue.F90

@@ -6,7 +6,6 @@ module eigenvalue
   use constants,   only: ZERO
   use error,       only: fatal_error, warning
   use math,        only: t_percentile
-  use mesh,        only: count_bank_sites
   use message_passing
   use random_lcg,  only: prn, set_particle_seed, advance_prn_seed
   use settings

src/eigenvalue.cpp

@@ -80,7 +80,6 @@ void ufs_count_sites()
     int64_t n;
     openmc_source_bank(&source_bank, &n);
 
-
     // count number of source sites in each ufs mesh cell
     bool sites_outside;
     source_frac = m->count_sites(openmc_work, source_bank, 0, nullptr,

src/mesh.cpp

@@ -659,7 +659,7 @@ xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
   std::size_t m = xt::prod(shape_)();
   std::vector<std::size_t> shape;
   if (n_energy > 0) {
-    shape = {m, static_cast<std::size_t>(n_energy)};
+    shape = {m, static_cast<std::size_t>(n_energy - 1)};
   } else {
     shape = {m};
   }
@@ -670,7 +670,8 @@ xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
 
   for (int64_t i = 0; i < n; ++i) {
     // determine scoring bin for entropy mesh
-    int mesh_bin = get_bin({bank[i].xyz});
+    // TODO: off-by-one
+    int mesh_bin = get_bin({bank[i].xyz}) - 1;
 
     // if outside mesh, skip particle
     if (mesh_bin < 0) {
@@ -680,7 +681,6 @@ xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
 
     if (n_energy > 0) {
       double E = bank[i].E;
-      int e_bin;
       if (E >= energies[0] && E <= energies[n_energy - 1]) {
         // determine energy bin
         int e_bin = lower_bound_index(energies, energies + n_energy, E);

src/message_passing.cpp

@@ -5,7 +5,7 @@ namespace mpi {
 
 int rank {0};
 int n_procs {1};
-bool master {false};
+bool master {true};
 
 #ifdef OPENMC_MPI
 MPI_Comm intracomm;

src/settings.cpp

@@ -520,7 +520,7 @@ void read_settings_xml()
       // If the user did not specify how many mesh cells are to be used in
       // each direction, we automatically determine an appropriate number of
       // cells
-      int n = std::ceil(std::pow(settings::n_particles / 20, 1.0/3.0));
+      int n = std::ceil(std::pow(settings::n_particles / 20.0, 1.0/3.0));
       m.shape_ = {n, n, n};
       m.n_dimension_ = 3;