diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 008d299708..59217c2753 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -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 # include # include -#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 diff --git a/src/cmfd_jacobian_operator.F90 b/src/cmfd_jacobian_operator.F90 index e47468bd0c..f58450b7db 100644 --- a/src/cmfd_jacobian_operator.F90 +++ b/src/cmfd_jacobian_operator.F90 @@ -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) diff --git a/src/cmfd_loss_operator.F90 b/src/cmfd_loss_operator.F90 index 38a6852602..ded0a613e0 100644 --- a/src/cmfd_loss_operator.F90 +++ b/src/cmfd_loss_operator.F90 @@ -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 diff --git a/src/cmfd_prod_operator.F90 b/src/cmfd_prod_operator.F90 index 93e1188399..85f355ab09 100644 --- a/src/cmfd_prod_operator.F90 +++ b/src/cmfd_prod_operator.F90 @@ -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 diff --git a/src/cmfd_slepc_solver.F90 b/src/cmfd_slepc_solver.F90 index d76afce478..9ce096f119 100644 --- a/src/cmfd_slepc_solver.F90 +++ b/src/cmfd_slepc_solver.F90 @@ -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 diff --git a/src/cmfd_snes_solver.F90 b/src/cmfd_snes_solver.F90 index 8eeba47535..439edaf5e8 100644 --- a/src/cmfd_snes_solver.F90 +++ b/src/cmfd_snes_solver.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index bb2c872af8..49c303ce29 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -236,6 +236,9 @@ module global ! Flag for coremap accelerator logical :: cmfd_coremap = .FALSE. + ! number of processors for cmfd + integer :: n_procs_cmfd + contains !===============================================================================