From 829b4c4f913e7aeb069ec35de16c046687d4f83c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Sep 2018 06:35:38 -0500 Subject: [PATCH] Use C++ based count_sites for CMFD. Fix several bugs --- CMakeLists.txt | 1 + src/cmfd_execute.F90 | 18 ++++++++++++++---- src/cmfd_execute.cpp | 34 ++++++++++++++++++++++++++++++++++ src/cmfd_header.F90 | 8 +++++--- src/eigenvalue.F90 | 1 - src/eigenvalue.cpp | 1 - src/mesh.cpp | 6 +++--- src/message_passing.cpp | 2 +- src/settings.cpp | 2 +- 9 files changed, 59 insertions(+), 14 deletions(-) create mode 100644 src/cmfd_execute.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 713bf7c077..f4ff9ae56d 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 74264a036e..6709f782ac 100644 --- a/src/cmfd_execute.F90 +++ b/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 diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp new file mode 100644 index 0000000000..7774cfe1d2 --- /dev/null +++ b/src/cmfd_execute.cpp @@ -0,0 +1,34 @@ +#include // for copy +#include +#include + +#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 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 diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index eafba7b8ad..e5e5cc423e 100644 --- a/src/cmfd_header.F90 +++ b/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 diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 14fe59110e..bd2a20b413 100644 --- a/src/eigenvalue.F90 +++ b/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 diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index ef4159e782..fa995e6a03 100644 --- a/src/eigenvalue.cpp +++ b/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, diff --git a/src/mesh.cpp b/src/mesh.cpp index 87d5f1e289..58f166f625 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -659,7 +659,7 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, std::size_t m = xt::prod(shape_)(); std::vector shape; if (n_energy > 0) { - shape = {m, static_cast(n_energy)}; + shape = {m, static_cast(n_energy - 1)}; } else { shape = {m}; } @@ -670,7 +670,8 @@ xt::xarray 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 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); diff --git a/src/message_passing.cpp b/src/message_passing.cpp index ec2fc2d688..dc287e3b3a 100644 --- a/src/message_passing.cpp +++ b/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; diff --git a/src/settings.cpp b/src/settings.cpp index db6aee7f54..7da8b31669 100644 --- a/src/settings.cpp +++ b/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;