Fix last mesh bugs (all tests pass now)

This commit is contained in:
Paul Romano 2018-09-05 11:23:20 -05:00
parent 7a6d36c279
commit f1fb400f92
6 changed files with 8 additions and 131 deletions

View file

@ -312,7 +312,6 @@ add_library(libopenmc SHARED
src/material_header.F90
src/math.F90
src/matrix_header.F90
src/mesh.F90
src/mesh_header.F90
src/message_passing.F90
src/mgxs_data.F90

View file

@ -41,6 +41,8 @@ _dll.openmc_mesh_set_params.errcheck = _error_handler
_dll.openmc_get_mesh_index.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_mesh_index.restype = c_int
_dll.openmc_get_mesh_index.errcheck = _error_handler
_dll.n_meshes.argtypes = []
_dll.n_meshes.restype = c_int
class Mesh(_FortranObjectWithID):
@ -172,10 +174,10 @@ class _MeshMapping(Mapping):
def __iter__(self):
for i in range(len(self)):
yield Mesh(index=i + 1).id
yield Mesh(index=i).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_meshes').value
return _dll.n_meshes()
def __repr__(self):
return repr(dict(self))

View file

@ -214,7 +214,6 @@ contains
use bank_header, only: source_bank
use constants, only: ZERO, ONE
use error, only: warning, fatal_error
use mesh, only: count_bank_sites
use message_passing
use string, only: to_str

View file

@ -1,106 +0,0 @@
module mesh
use algorithm, only: binary_search
use bank_header, only: bank
use constants
use mesh_header
use message_passing
implicit none
contains
!===============================================================================
! COUNT_BANK_SITES determines the number of fission bank sites in each cell of a
! given mesh as well as an optional energy group structure. This can be used for
! a variety of purposes (Shannon entropy, CMFD, uniform fission source
! weighting)
!===============================================================================
subroutine count_bank_sites(m, bank_array, cnt, energies, size_bank, &
sites_outside)
type(RegularMesh), intent(in) :: m ! mesh to count sites
type(Bank), intent(in) :: bank_array(:) ! fission or source bank
real(8), intent(out) :: cnt(:,:) ! weight of sites in each
! cell and energy group
real(8), intent(in), optional :: energies(:) ! energy grid to search
integer(8), intent(in), optional :: size_bank ! # of bank sites (on each proc)
logical, intent(inout), optional :: sites_outside ! were there sites outside mesh?
real(8), allocatable :: cnt_(:,:)
integer :: i ! loop index for local fission sites
integer :: n_sites ! size of bank array
integer :: n ! number of energy groups / size
integer :: mesh_bin ! mesh bin
integer :: e_bin ! energy bin
#ifdef OPENMC_MPI
integer :: mpi_err ! MPI error code
#endif
logical :: outside ! was any site outside mesh?
! initialize variables
allocate(cnt_(size(cnt,1), size(cnt,2)))
cnt_ = ZERO
outside = .false.
! Set size of bank
if (present(size_bank)) then
n_sites = int(size_bank,4)
else
n_sites = size(bank_array)
end if
! Determine number of energies in group structure
if (present(energies)) then
n = size(energies) - 1
else
n = 1
end if
! loop over fission sites and count how many are in each mesh box
FISSION_SITES: do i = 1, n_sites
! determine scoring bin for entropy mesh
call m % get_bin(bank_array(i) % xyz, mesh_bin)
! if outside mesh, skip particle
if (mesh_bin == NO_BIN_FOUND) then
outside = .true.
cycle
end if
! determine energy bin
if (present(energies)) then
if (bank_array(i) % E < energies(1)) then
e_bin = 1
elseif (bank_array(i) % E > energies(n + 1)) then
e_bin = n
else
e_bin = binary_search(energies, n + 1, bank_array(i) % E)
end if
else
e_bin = 1
end if
! add to appropriate mesh box
cnt_(e_bin, mesh_bin) = cnt_(e_bin, mesh_bin) + bank_array(i) % wgt
end do FISSION_SITES
#ifdef OPENMC_MPI
! collect values from all processors
n = size(cnt_)
call MPI_REDUCE(cnt_, cnt, n, MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err)
! Check if there were sites outside the mesh for any processor
if (present(sites_outside)) then
call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, &
mpi_intracomm, mpi_err)
end if
#else
sites_outside = outside
cnt = cnt_
#endif
end subroutine count_bank_sites
end module mesh

View file

@ -443,7 +443,7 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
}
}
int j = xt::argmin(d)(0);
j = xt::argmin(d)(0);
if (u[j] > 0.0) {
++ijk0(j);
} else {
@ -729,11 +729,11 @@ xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
extern "C" int
openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end)
{
if (!index_start) *index_start = meshes.size();
if (index_start) *index_start = meshes.size();
for (int i = 0; i < n; ++i) {
meshes.emplace_back(new RegularMesh{});
}
if (!index_end) *index_end = meshes.size() - 1;
if (index_end) *index_end = meshes.size() - 1;
return 0;
}
@ -824,6 +824,7 @@ openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width,
*ll = m->lower_left_.data();
*ur = m->upper_right_.data();
*width = m->width_.data();
*n = m->n_dimension_;
return 0;
}

View file

@ -73,25 +73,7 @@ contains
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
integer, parameter :: MAX_SEARCH_ITER = 100 ! Maximum number of times we can
! can loop while trying to find
! the first intersection.
integer :: j ! loop index for direction
integer :: n
integer :: ijk0(3) ! indices of starting coordinates
integer :: ijk1(3) ! indices of ending coordinates
integer :: search_iter ! loop count for intersection search
integer :: bin
real(8) :: uvw(3) ! cosine of angle of particle
real(8) :: xyz0(3) ! starting/intermediate coordinates
real(8) :: xyz1(3) ! ending coordinates of particle
real(8) :: xyz_cross ! coordinates of next boundary
real(8) :: d(3) ! distance to each bounding surface
real(8) :: total_distance ! distance of entire particle track
real(8) :: distance ! distance traveled in mesh cell
logical :: start_in_mesh ! starting coordinates inside mesh?
logical :: end_in_mesh ! ending coordinates inside mesh?
type(RegularMesh) :: m
type(C_PTR) :: ptr_bins, ptr_weights