crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

JFNK is now parallelized

Browse source
This commit is contained in:
Bryan Herman 2012-02-14 11:44:36 -08:00
parent 80777a7e36
commit ca531024f7
3 changed files with 117 additions and 65 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/Makefile
View file

@ -19,7 +19,7 @@ F90X = /home/bherman/petsc/release/bin/mpif90
F90XFLAGS = -Wall -Wno-unused-variable
COMPILER = petsc
DEBUG = no
PROFILE = yes
PROFILE = no
OPTIMIZE = no
USE_MPI = no
USE_HDF5 = yes

94
src/cmfd_jacobian_operator.F90
View file

@ -87,7 +87,7 @@ contains
end if
! add 1 for eigenvalue row
this%n = this%n + 1
! this%n = this%n + 1
end subroutine get_J_indices
@ -97,13 +97,11 @@ contains
subroutine preallocate_jacobian_matrix(this,ctx)
use global, only: cmfd,cmfd_coremap
use global, only: cmfd,cmfd_coremap,rank,n_procs
type(jacobian_operator) :: this
type(operators) :: ctx
integer :: rank ! rank of processor
integer :: sizen ! number of procs
integer :: i ! iteration counter for x
integer :: j ! iteration counter for y
integer :: k ! iteration counter for z
@ -123,25 +121,24 @@ contains
integer :: neig_mat_idx ! matrix index of neighbor cell
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)
! get local problem size
n = this%n
! determine local size, divide evenly between all other procs
this%localn = n/(sizen)
this%localn = n/(n_procs)
! add 1 more if less proc id is less than mod
if (rank < mod(n,sizen)) this%localn = this%localn + 1
if (rank < mod(n,n_procs)) this%localn = this%localn + 1
! add another 1 on last proc
if (rank == n_procs - 1) this%localn = this%localn + 1
! determine local starting row
row_start = 0
if (rank < mod(n,sizen)) then
row_start = rank*(n/sizen+1)
if (rank < mod(n,n_procs)) then
row_start = rank*(n/n_procs+1)
else
row_start = min(mod(n,sizen)*(n/sizen+1)+(rank - mod(n,sizen))*(n/sizen),n)
row_start = min(mod(n,n_procs)*(n/n_procs+1)+(rank - mod(n,n_procs))*(n/n_procs),n)
end if
! determine local final row
@ -154,18 +151,23 @@ contains
this % o_nnz = 0
! start with pattern from loss matrix
this % d_nnz = ctx%loss%d_nnz
this % o_nnz = ctx%loss%o_nnz
if (rank == n_procs - 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
this % d_nnz = ctx%loss%d_nnz
this % o_nnz = ctx%loss%o_nnz
end if
! append -F*phi term for last processor will take care of 1 for lambda
if (row_end == n - 1) then
if (rank == n_procs - 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 (row_end == n - 1) then
if (rank == n_procs - 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))
@ -179,6 +181,8 @@ contains
subroutine build_jacobian_matrix(snes,x,jac,jac_prec,flag,ctx,ierr)
use global, only: rank,n_procs
! formal variables
SNES :: snes ! the snes context
Vec :: x ! the solution vector
@ -197,12 +201,16 @@ contains
integer :: irow ! row counter
integer :: row_start! starting local row on process
integer :: row_end ! ending local row on process
integer, allocatable :: dims(:) ! vec of starting and ending rows
integer, allocatable :: dims1(:) ! vec of sizes on each proc
integer, allocatable :: nnzv(:) ! vector of number of nonzeros for jac
integer, allocatable :: cols(:) ! vector of column numbers
real(8) :: lambda ! eigenvalue
real(8), pointer :: xptr(:) ! pointer to solution vector
real(8), pointer :: sptr(:) ! pointer to source vector
real(8), allocatable :: vals(:) ! vector of row values
real(8), allocatable :: phi_tmp(:) ! temp buffer for flux
print *,'In Jacobian'
! create operators
call build_loss_matrix(ctx%loss)
@ -224,21 +232,18 @@ contains
call VecGetArrayF90(x,xptr,ierr)
! create petsc vector for flux
call VecCreate(PETSC_COMM_SELF,phi,ierr)
call VecSetSizes(phi,PETSC_DECIDE,n,ierr)
call VecSetFromOptions(phi,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,ctx%loss%localn,PETSC_DECIDE,phi,ierr)
! extract flux and eigenvalue
! extract flux and eigenvalue
call VecPlaceArray(phi,xptr,ierr)
lambda = xptr(n+1)
if (rank == n_procs - 1) lambda = xptr(size(xptr))
call MPI_BCAST(lambda,1,MPI_REAL8,n_procs-1,MPI_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)
! create tmp petsc vector for source
call VecCreate(PETSC_COMM_SELF,source,ierr)
call VecSetSizes(source,PETSC_DECIDE,n,ierr)
call VecSetFromOptions(source,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,ctx%loss%localn,PETSC_DECIDE,source,ierr)
! perform math (-F*phi --> source)
call MatMult(ctx%prod%F,phi,source,ierr)
@ -263,13 +268,33 @@ contains
call MatRestoreRow(ctx%loss%M,irow,ncols,cols,vals,ierr)
! insert source value
call MatSetValue(jac_prec,irow,n,sptr(irow+1),INSERT_VALUES,ierr)
call MatSetValue(jac_prec,irow,n,sptr(irow-row_start+1),INSERT_VALUES,ierr)
end do
! allocate space for flux vector buffer
if (rank == n_procs - 1) then
! temporary receive buffer
if (.not. allocated(phi_tmp)) allocate(phi_tmp(0:n-1))
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))
call VecGetOwnershipRanges(phi,dims,ierr)
do k = 0,n_procs-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)
! set values in last row of matrix
if (row_end - 1 == n) then
call MatSetValues(jac_prec,1,n,n,(/(k,k=0,n-1)/),xptr(1:n),INSERT_VALUES,ierr)
if (rank == n_procs - 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)
end if
@ -293,6 +318,12 @@ contains
! deallocate all temporary space
if (allocated(cols)) deallocate(cols)
if (allocated(vals)) deallocate(vals)
if (allocated(phi_tmp)) deallocate(phi_tmp)
if (allocated(dims)) deallocate(dims)
if (allocated(dims1)) deallocate(dims1)
! print jacobian out
call print_J_operator(jac_prec)
end subroutine build_jacobian_matrix
@ -300,18 +331,17 @@ contains
! PRINT_J_OPERATOR
!===============================================================================
subroutine print_J_operator(this)
subroutine print_J_operator(jac)
use global, only: path_input
Mat :: jac
PetscViewer :: viewer
type(jacobian_operator) :: this
! write out matrix in binary file (debugging)
call PetscViewerBinaryOpen(PETSC_COMM_WORLD,trim(path_input)//'jacobian.bin'&
& ,FILE_MODE_WRITE,viewer,ierr)
call MatView(this%J,viewer,ierr)
call MatView(jac,viewer,ierr)
call PetscViewerDestroy(viewer,ierr)
end subroutine print_J_operator

86
src/cmfd_snes_solver.F90
View file

@ -64,7 +64,7 @@ print *,'solving system'
subroutine init_data()
use global, only: cmfd
use global, only: cmfd,rank,n_procs
integer :: k ! implied do counter
integer :: n ! problem size
@ -79,25 +79,26 @@ print *,'solving system'
call init_J_operator(jac_prec,ctx)
! get problem size
n = jac_prec%n - 1
n = jac_prec%n
print *,cmfd%phi
print *,cmfd%keff
! create PETSc vectors
call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n+1,resvec,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n+1,xvec,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,jac_prec%localn,PETSC_DECIDE,resvec,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,jac_prec%localn,PETSC_DECIDE,xvec,ierr)
! get the local dimensions for each process
call VecGetOwnershipRange(xvec,row_start,row_end,ierr)
if (row_end == n + 1) row_end = n
if (rank == n_procs - 1) row_end = n
! set flux in guess
call VecSetValues(xvec,row_end-row_start-1,(/(k,k=row_start,row_end-1)/), &
call VecSetValues(xvec,row_end-row_start,(/(k,k=row_start,row_end-1)/), &
& cmfd%phi(row_start+1:row_end),INSERT_VALUES,ierr)
call VecAssemblyBegin(xvec,ierr)
call VecAssemblyEnd(xvec,ierr)
! set keff in guess
if (row_end == n) then
if (rank == n_procs - 1) then
call VecGetArrayF90(xvec,xptr,ierr)
xptr(size(xptr)) = 1.0_8/cmfd%keff
call VecRestoreArrayF90(xvec,xptr,ierr)
@ -172,6 +173,8 @@ print *,'solving system'
subroutine compute_nonlinear_residual(snes,x,res,ierr)
use global, only: rank,n_procs,path_input
! arguments
SNES :: snes ! nonlinear solver context
Vec :: x ! independent vector
@ -183,50 +186,39 @@ print *,'solving system'
Vec :: rphi ! flux part of residual
Vec :: phiM ! M part of residual flux calc
integer :: n ! problem size
integer :: row_start ! local row start
integer :: row_end ! local row end
integer :: n_procs ! number of processors
real(8) :: lambda ! eigenvalue
real(8) :: reslamb ! residual for lambda
real(8), pointer :: xptr(:) ! pointer to solution vector
real(8), pointer :: rptr(:) ! pointer to residual vector
PetscViewer :: viewer
print *,'In residual'
! get number of processors
call MPI_COMM_SIZE(MPI_COMM_WORLD,n_procs,ierr)
! create operators
call build_loss_matrix(ctx%loss)
call build_prod_matrix(ctx%prod)
! get problem size
call VecGetSize(x,n,ierr)
n = n - 1 ! subtract off last row
! get the local dimensions for each process
call VecGetOwnershipRange(x,row_start,row_end,ierr)
n = ctx%loss%n
! get pointers to vectors
call VecGetArrayF90(x,xptr,ierr)
call VecGetArrayF90(res,rptr,ierr)
! create petsc vector for flux
call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n,phi,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n,rphi,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,ctx%loss%localn,PETSC_DECIDE,phi,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,ctx%loss%localn,PETSC_DECIDE,rphi,ierr)
! extract flux and place in petsc vector
call VecPlaceArray(phi,xptr,ierr)
call VecPlaceArray(rphi,rptr,ierr)
! extract eigenvalue and broadcast (going to want to make this more general in future)
if (row_end == n+1) then
lambda = xptr(size(xptr))
end if
if (rank == n_procs - 1) lambda = xptr(size(xptr))
call MPI_BCAST(lambda,1,MPI_REAL8,n_procs-1,MPI_COMM_WORLD,ierr)
! create operators
call build_loss_matrix(ctx%loss)
call build_prod_matrix(ctx%prod)
! create new petsc vectors to perform math
call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,n,phiM,ierr)
call VecCreateMPI(PETSC_COMM_WORLD,ctx%loss%localn,PETSC_DECIDE,phiM,ierr)
! calculate flux part of residual vector
call MatMult(ctx%loss%M,phi,phiM,ierr)
@ -237,7 +229,7 @@ print *,'In residual'
call VecDot(phi,phi,reslamb,ierr)
! map to ptr
if (row_end == n+1) rptr(size(rptr)) = 0.5_8 - 0.5_8*reslamb
if (rank == n_procs) rptr(size(rptr)) = 0.5_8 - 0.5_8*reslamb
! reset arrays that are not used
call VecResetArray(phi,ierr)
@ -252,6 +244,17 @@ print *,'In residual'
call VecDestroy(phiM,ierr)
call VecDestroy(rphi,ierr)
! write out matrix in binary file (debugging)
call PetscViewerBinaryOpen(PETSC_COMM_WORLD,trim(path_input)//'residual.bin' &
& ,FILE_MODE_WRITE,viewer,ierr)
call VecView(res,viewer,ierr)
call PetscViewerDestroy(viewer,ierr)
call PetscViewerBinaryOpen(PETSC_COMM_WORLD,trim(path_input)//'x.bin' &
& ,FILE_MODE_WRITE,viewer,ierr)
call VecView(x,viewer,ierr)
call PetscViewerDestroy(viewer,ierr)
end subroutine compute_nonlinear_residual
!===============================================================================
@ -260,9 +263,12 @@ print *,'In residual'
subroutine extract_results()
use global, only: cmfd
use global, only: cmfd,rank,n_procs
integer :: n ! problem size
integer :: row_start ! local row start
integer :: row_end ! local row end
real(8),allocatable :: mybuf(:) ! temp buffer
PetscScalar, pointer :: xptr(:) ! pointer to eigenvector info
! get problem size
@ -270,13 +276,29 @@ print *,'In residual'
! also allocate in cmfd object
if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n))
if (.not. allocated(mybuf)) allocate(mybuf(n))
! get ownership range
call VecGetOwnershipRange(xvec,row_start,row_end,ierr)
! resize the last proc
if (rank == n_procs - 1) row_end = row_end - 1
! convert petsc phi_object to cmfd_obj
call VecGetArrayF90(xvec,xptr,ierr)
cmfd%phi = xptr(1:n)
cmfd%phi(row_start+1:row_end) = xptr(1:row_end - row_start)
! reduce result to all
mybuf = 0.0_8
call MPI_ALLREDUCE(cmfd%phi,mybuf,n,MPI_REAL8,MPI_SUM,MPI_COMM_WORLD,ierr)
! move buffer to object and deallocate
cmfd%phi = mybuf
if(allocated(mybuf)) deallocate(mybuf)
! save eigenvalue
cmfd%keff = 1.0_8 / xptr(n+1)
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)
call VecRestoreArrayF90(xvec,xptr,ierr)
end subroutine extract_results