From ca531024f7b3070586ff1e205f285002f13b364c Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 14 Feb 2012 11:44:36 -0800 Subject: [PATCH] JFNK is now parallelized --- src/Makefile | 2 +- src/cmfd_jacobian_operator.F90 | 94 ++++++++++++++++++++++------------ src/cmfd_snes_solver.F90 | 86 +++++++++++++++++++------------ 3 files changed, 117 insertions(+), 65 deletions(-) diff --git a/src/Makefile b/src/Makefile index afc6b8ab5f..445cc5df02 100644 --- a/src/Makefile +++ b/src/Makefile @@ -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 diff --git a/src/cmfd_jacobian_operator.F90 b/src/cmfd_jacobian_operator.F90 index 745b59fc88..1ada9797ce 100644 --- a/src/cmfd_jacobian_operator.F90 +++ b/src/cmfd_jacobian_operator.F90 @@ -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 diff --git a/src/cmfd_snes_solver.F90 b/src/cmfd_snes_solver.F90 index 9cd8a7d232..d8d1979cf3 100644 --- a/src/cmfd_snes_solver.F90 +++ b/src/cmfd_snes_solver.F90 @@ -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