Merge neutral particle reservoirs across MPI processes
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(NeutralParticleReservoir), | intent(inout) | :: | res | Neutral Particle Reservoir |
subroutine parallel_merge_reservoirs(res)
!+ Merge neutral particle reservoirs across MPI processes
type(NeutralParticleReservoir), intent(inout) :: res
!+ Neutral Particle Reservoir
integer :: proc, num_procs, i, j, k
integer(Int32), dimension(:), allocatable :: nums, inds
real(Float64), dimension(:), allocatable :: Ws
real(Float64), dimension(:,:), allocatable :: wghts
real(Float64), dimension(:,:,:), allocatable :: vels
type(NeutralParticleReservoir) :: res1, res2
proc = my_rank() + 1
num_procs = num_ranks()
k = res%k
allocate(nums(num_procs), inds(num_procs), Ws(num_procs))
nums = 0
inds = 0
Ws = 0.d0
nums(proc) = res%n
inds(proc) = res%i
Ws(proc) = res%W
call parallel_sum(nums)
call parallel_sum(inds)
call parallel_sum(Ws)
allocate(vels(3,reservoir_size,num_procs))
allocate(wghts(reservoir_size,num_procs))
vels = 0.d0; wghts = 0.d0
do i=1,min(k,res%n)
vels(:,i,proc) = res%R(i)%v
wghts(i,proc) = res%R(i)%w
enddo
call parallel_sum(vels)
call parallel_sum(wghts)
call init_reservoir(res1)
res1%n = nums(1)
res1%i = inds(1)
res1%W = Ws(1)
do i=1,min(k,res1%n)
res1%R(i) = NeutralParticle(wghts(i,1),vels(:,i,1))
enddo
call init_reservoir(res2)
do i=2,num_procs
res2%n = nums(i)
res2%i = inds(i)
res2%W = Ws(i)
do j=1,min(k,res2%n)
res2%R(j) = NeutralParticle(wghts(j,i),vels(:,j,i))
enddo
call merge_reservoirs(res1,res2)
enddo
deallocate(vels,wghts,nums,inds,Ws)
res = res1
end subroutine parallel_merge_reservoirs