crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

added a new MPI Communicator for CMFD run which is a subset of MPI_COMM_WORLD,it is set to PETSC_COMM_WORLD, appropriate changes were made to routines to have PETSC_COMM_WORLD, and number of processor for CMFD was added as a global var

Browse source
This commit is contained in:
Bryan Herman 2012-02-21 17:25:03 -08:00
parent 6371a19341
commit 8897963f45
7 changed files with 109 additions and 70 deletions
Download patch file Download diff file Expand all files Collapse all files

93
src/cmfd_execute.F90
View file

@ -3,7 +3,7 @@ module cmfd_execute
use cmfd_data, only: set_up_cmfd
use cmfd_output, only: write_cmfd_vtk
use global, only: cmfd,cmfd_only,time_cmfd,master,rank,mpi_err, &
& current_cycle,n_inactive,n_cycles
& current_cycle,n_inactive,n_cycles,n_procs,n_procs_cmfd
use timing, only: timer_start,timer_stop,timer_reset
@ -19,7 +19,6 @@ module cmfd_execute
# include <finclude/petsc.h90>
# include <finclude/slepcsys.h>
# include <finclude/slepceps.h>
#endif
contains
@ -30,12 +29,19 @@ contains
subroutine execute_cmfd()
integer :: ierr ! petsc error code
integer :: myrank
! initialize mpi communicator
call petsc_init_mpi()
if (rank < 6) then
#ifdef PETSC
! initialize slepc/petsc (communicates to world)
if(current_cycle == n_inactive + 1) call SlepcInitialize &
& (PETSC_NULL_CHARACTER,ierr)
#endif
! set global variable for number of procs in cmfd calc
call MPI_COMM_SIZE(PETSC_COMM_WORLD,n_procs_cmfd,ierr)
! only run if master process
if (master) then
@ -51,10 +57,8 @@ contains
! broadcast cmfd object to all procs
call cmfd_bcast()
#ifdef PETSC
! execute snes solver
call cmfd_snes_execute()
#endif
! stop timer
call timer_stop(time_cmfd)
@ -62,13 +66,13 @@ contains
! write vtk file
!if(.not. cmfd_only) call write_cmfd_vtk()
#ifdef PETSC
! finalize slepc
! finalize slepc
if (current_cycle == n_cycles) call SlepcFinalize(ierr)
end if
! sync up procs
call MPI_Barrier(MPI_COMM_WORLD,mpi_err)
#endif
end subroutine execute_cmfd
@ -97,32 +101,69 @@ contains
call allocate_cmfd(cmfd)
! sync up procs
call MPI_Barrier(MPI_COMM_WORLD,mpi_err)
call MPI_Barrier(PETSC_COMM_WORLD,mpi_err)
! broadcast all data
call MPI_BCAST(cmfd%flux,ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%totalxs,ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%p1scattxs,ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%scattxs,ng*ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%nfissxs,ng*ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%diffcof,ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%dtilde,6*ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%dhat,6*ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%hxyz,3*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%current,12*ng*nx*ny*nz,MPI_REAL8,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%flux,ng*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%totalxs,ng*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%p1scattxs,ng*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%scattxs,ng*ng*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%nfissxs,ng*ng*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%diffcof,ng*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%dtilde,6*ng*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%dhat,6*ng*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%hxyz,3*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%current,12*ng*nx*ny*nz,MPI_REAL8,0,PETSC_COMM_WORLD,mpi_err)
! broadcast coremap info
if (cmfd_coremap) then
call MPI_BCAST(cmfd%coremap,nx*ny*nz,MPI_INT,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%mat_dim,1,MPI_INT,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%coremap,nx*ny*nz,MPI_INT,0,PETSC_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%mat_dim,1,MPI_INT,0,PETSC_COMM_WORLD,mpi_err)
if (.not. allocated(cmfd % indexmap)) allocate &
& (cmfd % indexmap(cmfd % mat_dim,3))
call MPI_BCAST(cmfd%indexmap,cmfd%mat_dim*3,MPI_INT,0,MPI_COMM_WORLD,mpi_err)
call MPI_BCAST(cmfd%indexmap,cmfd%mat_dim*3,MPI_INT,0,PETSC_COMM_WORLD,mpi_err)
end if
! sync up procs
call MPI_Barrier(MPI_COMM_WORLD,mpi_err)
end subroutine cmfd_bcast
!===============================================================================
! PETSC_INIT_MPI
!===============================================================================
subroutine petsc_init_mpi()
integer :: new_comm ! new communicator
integer :: orig_group ! original MPI group for MPI_COMM_WORLD
integer :: new_group ! new MPI group subset of orig_group
integer,allocatable :: ranks(:) ! ranks to include for petsc
integer :: k ! iteration counter
integer :: myrank=9999
! set ranks 0-6 or min
if (n_procs >= 6) then
if (.not. allocated(ranks)) allocate(ranks(0:5))
ranks = (/0,1,2,3,4,5/)
else
if (.not. allocated(ranks)) allocate(ranks(0:n_procs-1))
ranks = (/(k,k=0,n_procs-1)/)
end if
! get the original mpi group
call MPI_COMM_GROUP(MPI_COMM_WORLD,orig_group,mpi_err)
! new group init
call MPI_GROUP_INCL(orig_group,size(ranks),ranks,new_group,mpi_err)
! create new communicator
call MPI_COMM_CREATE(MPI_COMM_WORLD,new_group,new_comm,mpi_err)
! deallocate ranks
if (allocated(ranks)) deallocate(ranks)
! set PETSC_COMM_WORLD to this subset
PETSC_COMM_WORLD = new_comm
end subroutine petsc_init_mpi
#endif
end module cmfd_execute

40
src/cmfd_jacobian_operator.F90
View file

@ -97,7 +97,7 @@ contains
subroutine preallocate_jacobian_matrix(this,ctx)
use global, only: cmfd,cmfd_coremap,rank,n_procs
use global, only: cmfd,cmfd_coremap,rank,n_procs_cmfd
type(jacobian_operator) :: this
type(operators) :: ctx
@ -125,20 +125,20 @@ contains
n = this%n
! determine local size, divide evenly between all other procs
this%localn = n/(n_procs)
this%localn = n/(n_procs_cmfd)
! add 1 more if less proc id is less than mod
if (rank < mod(n,n_procs)) this%localn = this%localn + 1
if (rank < mod(n,n_procs_cmfd)) this%localn = this%localn + 1
! add another 1 on last proc
if (rank == n_procs - 1) this%localn = this%localn + 1
if (rank == n_procs_cmfd - 1) this%localn = this%localn + 1
! determine local starting row
row_start = 0
if (rank < mod(n,n_procs)) then
row_start = rank*(n/n_procs+1)
if (rank < mod(n,n_procs_cmfd)) then
row_start = rank*(n/n_procs_cmfd+1)
else
row_start = min(mod(n,n_procs)*(n/n_procs+1)+(rank - mod(n,n_procs))*(n/n_procs),n)
row_start = min(mod(n,n_procs_cmfd)*(n/n_procs_cmfd+1)+(rank - mod(n,n_procs_cmfd))*(n/n_procs_cmfd),n)
end if
! determine local final row
@ -151,7 +151,7 @@ contains
this % o_nnz = 0
! start with pattern from loss matrix
if (rank == n_procs - 1) then
if (rank == n_procs_cmfd - 1) then
this % d_nnz(row_start:row_end-1) = ctx%loss%d_nnz
this % o_nnz(row_start:row_end-1) = ctx%loss%o_nnz
else
@ -160,14 +160,14 @@ contains
end if
! append -F*phi term for last processor will take care of 1 for lambda
if (rank == n_procs - 1) then
if (rank == n_procs_cmfd - 1) then
this%d_nnz = this%d_nnz + 1
else
this%o_nnz = this%o_nnz + 1
end if
! do last row which has all filled (already did lower left corner above)
if (rank == n_procs - 1) then
if (rank == n_procs_cmfd - 1) then
this % d_nnz(row_end) = this % d_nnz(row_end) + (row_end - row_start)
this % o_nnz(row_end) = this % o_nnz(row_end) + ((this%n-1) - (row_end - &
& row_start))
@ -181,7 +181,7 @@ contains
subroutine build_jacobian_matrix(snes,x,jac,jac_prec,flag,ctx,ierr)
use global, only: rank,n_procs
use global, only: rank,n_procs_cmfd
! formal variables
SNES :: snes ! the snes context
@ -235,8 +235,8 @@ contains
! extract flux and eigenvalue
call VecPlaceArray(phi,xptr,ierr)
if (rank == n_procs - 1) lambda = xptr(size(xptr))
call MPI_BCAST(lambda,1,MPI_REAL8,n_procs-1,MPI_COMM_WORLD,ierr)
if (rank == n_procs_cmfd - 1) lambda = xptr(size(xptr))
call MPI_BCAST(lambda,1,MPI_REAL8,n_procs_cmfd-1,PETSC_COMM_WORLD,ierr)
! compute math (M-lambda*F) M is overwritten here
call MatAXPY(ctx%loss%M,-1.0_8*lambda,ctx%prod%F,SUBSET_NONZERO_PATTERN,ierr)
@ -272,26 +272,26 @@ contains
end do
! allocate space for flux vector buffer
if (rank == n_procs - 1) then
! if (rank == n_procs_cmfd - 1) then
! temporary receive buffer
if (.not. allocated(phi_tmp)) allocate(phi_tmp(0:n-1))
end if
! end if
! get size on each proc
if (.not. allocated(dims)) allocate(dims(0:n_procs))
if (.not. allocated(dims1)) allocate(dims1(0:n_procs-1))
if (.not. allocated(dims)) allocate(dims(0:n_procs_cmfd))
if (.not. allocated(dims1)) allocate(dims1(0:n_procs_cmfd-1))
call VecGetOwnershipRanges(phi,dims,ierr)
do k = 0,n_procs-1
do k = 0,n_procs_cmfd-1
dims1(k) = dims(k+1) - dims(k)
end do
! gather data on all procs (will truncate xptr if needed for last proc)
call MPI_GATHERV(xptr,dims1(rank),MPI_REAL8,phi_tmp,dims1,dims(0:n_procs-1),MPI_REAL8,n_procs-1,MPI_COMM_WORLD,ierr)
call MPI_GATHERV(xptr,dims1(rank),MPI_REAL8,phi_tmp,dims1,dims(0:n_procs_cmfd-1),MPI_REAL8,n_procs_cmfd-1,PETSC_COMM_WORLD,ierr)
! set values in last row of matrix
if (rank == n_procs - 1) then
if (rank == n_procs_cmfd - 1) then
phi_tmp = -1.0_8*phi_tmp ! negate the transpose
call MatSetValues(jac_prec,1,n,n,(/(k,k=0,n-1)/),phi_tmp,INSERT_VALUES,ierr)
call MatSetValue(jac_prec,n,n,1.0_8,INSERT_VALUES,ierr)

4
src/cmfd_loss_operator.F90
View file

@ -110,8 +110,8 @@ contains
integer :: scatt_mat_idx ! matrix index for h-->g scattering terms
! get rank and max rank of procs
call MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD,sizen,ierr)
call MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr)
call MPI_COMM_SIZE(PETSC_COMM_WORLD,sizen,ierr)
! get local problem size
n = this%n

4
src/cmfd_prod_operator.F90
View file

@ -101,8 +101,8 @@ contains
integer :: hmat_idx ! index in matrix for energy group h
! get rank and max rank of procs
call MPI_COMM_RANK(MPI_COMM_WORLD,rank,ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD,sizen,ierr)
call MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr)
call MPI_COMM_SIZE(PETSC_COMM_WORLD,sizen,ierr)
! get local problem size
n = this%n

2
src/cmfd_slepc_solver.F90
View file

@ -177,7 +177,7 @@ call timer_start(time_cmfd)
! reduce result to master
mybuf = 0.0_8
call MPI_ALLREDUCE(cmfd%phi,mybuf,n,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,ierr)
call MPI_ALLREDUCE(cmfd%phi,mybuf,n,MPI_REAL8,MPI_SUM,PETSC_COMM_WORLD,ierr)
! move buffer to object and deallocate
cmfd%phi = mybuf

33
src/cmfd_snes_solver.F90
View file

@ -34,9 +34,6 @@ contains
subroutine cmfd_snes_execute()
use global, only: time_snes,master
use timing
! call slepc solver
call cmfd_slepc_execute()
@ -45,12 +42,10 @@ contains
! initialize solver
call init_solver()
call timer_reset(time_snes)
call timer_start(time_snes)
! solve the system
call SNESSolve(snes,PETSC_NULL,xvec,ierr)
call timer_stop(time_snes)
if(master) print *,'SNES Solution Time:',time_snes%elapsed
! extracts results to cmfd object
call extract_results()
@ -65,7 +60,7 @@ if(master) print *,'SNES Solution Time:',time_snes%elapsed
subroutine init_data()
use global, only: cmfd,rank,n_procs
use global, only: cmfd,rank,n_procs_cmfd
integer :: k ! implied do counter
integer :: n ! problem size
@ -89,7 +84,7 @@ if(master) print *,'SNES Solution Time:',time_snes%elapsed
! get the local dimensions for each process
call VecGetOwnershipRange(xvec,row_start,row_end,ierr)
if (rank == n_procs - 1) row_end = n
if (rank == n_procs_cmfd - 1) row_end = n
! set flux in guess
call VecSetValues(xvec,row_end-row_start,(/(k,k=row_start,row_end-1)/), &
@ -98,7 +93,7 @@ if(master) print *,'SNES Solution Time:',time_snes%elapsed
call VecAssemblyEnd(xvec,ierr)
! set keff in guess
if (rank == n_procs - 1) then
if (rank == n_procs_cmfd - 1) then
call VecGetArrayF90(xvec,xptr,ierr)
xptr(size(xptr)) = 1.0_8/cmfd%keff
call VecRestoreArrayF90(xvec,xptr,ierr)
@ -172,7 +167,7 @@ if(master) print *,'SNES Solution Time:',time_snes%elapsed
subroutine compute_nonlinear_residual(snes,x,res,ierr)
use global, only: rank,n_procs,path_input
use global, only: rank,n_procs_cmfd,path_input
! arguments
SNES :: snes ! nonlinear solver context
@ -212,8 +207,8 @@ PetscViewer :: viewer
call VecPlaceArray(rphi,rptr,ierr)
! extract eigenvalue and broadcast (going to want to make this more general in future)
if (rank == n_procs - 1) lambda = xptr(size(xptr))
call MPI_BCAST(lambda,1,MPI_REAL8,n_procs-1,MPI_COMM_WORLD,ierr)
if (rank == n_procs_cmfd - 1) lambda = xptr(size(xptr))
call MPI_BCAST(lambda,1,MPI_REAL8,n_procs_cmfd-1,PETSC_COMM_WORLD,ierr)
! create new petsc vectors to perform math
call VecCreateMPI(PETSC_COMM_WORLD,ctx%loss%localn,PETSC_DECIDE,phiM,ierr)
@ -227,7 +222,7 @@ PetscViewer :: viewer
call VecDot(phi,phi,reslamb,ierr)
! map to ptr
if (rank == n_procs) rptr(size(rptr)) = 0.5_8 - 0.5_8*reslamb
if (rank == n_procs_cmfd) rptr(size(rptr)) = 0.5_8 - 0.5_8*reslamb
! reset arrays that are not used
call VecResetArray(phi,ierr)
@ -261,7 +256,7 @@ PetscViewer :: viewer
subroutine extract_results()
use global, only: cmfd,rank,n_procs
use global, only: cmfd,rank,n_procs_cmfd
integer :: n ! problem size
integer :: row_start ! local row start
@ -280,7 +275,7 @@ PetscViewer :: viewer
call VecGetOwnershipRange(xvec,row_start,row_end,ierr)
! resize the last proc
if (rank == n_procs - 1) row_end = row_end - 1
if (rank == n_procs_cmfd - 1) row_end = row_end - 1
! convert petsc phi_object to cmfd_obj
call VecGetArrayF90(xvec,xptr,ierr)
@ -288,15 +283,15 @@ PetscViewer :: viewer
! reduce result to all
mybuf = 0.0_8
call MPI_ALLREDUCE(cmfd%phi,mybuf,n,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,ierr)
call MPI_ALLREDUCE(cmfd%phi,mybuf,n,MPI_REAL8,MPI_SUM,PETSC_COMM_WORLD,ierr)
! move buffer to object and deallocate
cmfd%phi = mybuf
if(allocated(mybuf)) deallocate(mybuf)
! save eigenvalue
if(rank == n_procs - 1) cmfd%keff = 1.0_8 / xptr(size(xptr))
call MPI_BCAST(cmfd%keff,1,MPI_REAL8,n_procs-1,MPI_COMM_WORLD,ierr)
if(rank == n_procs_cmfd - 1) cmfd%keff = 1.0_8 / xptr(size(xptr))
call MPI_BCAST(cmfd%keff,1,MPI_REAL8,n_procs_cmfd-1,PETSC_COMM_WORLD,ierr)
call VecRestoreArrayF90(xvec,xptr,ierr)
end subroutine extract_results

3
src/global.F90
View file

@ -236,6 +236,9 @@ module global
! Flag for coremap accelerator
logical :: cmfd_coremap = .FALSE.
! number of processors for cmfd
integer :: n_procs_cmfd
contains
!===============================================================================