diff --git a/src/Makefile b/src/Makefile index c0027f849f..3259848ede 100644 --- a/src/Makefile +++ b/src/Makefile @@ -18,7 +18,7 @@ include OBJECTS F90X = /home/bherman/petsc/release/bin/mpif90 F90XFLAGS = -Wall -Wno-unused-variable COMPILER = petsc -DEBUG = no +DEBUG = yes PROFILE = no OPTIMIZE = no USE_MPI = no diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 2fc2e8907b..053a14b0d0 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -38,7 +38,7 @@ contains #endif ! only run if master process -! if (master) then + if (master) then ! begin timer call timer_start(time_cmfd) @@ -46,6 +46,11 @@ contains ! set up cmfd if(.not. cmfd_only) call set_up_cmfd() + end if + + ! broadcast cmfd object to all procs + call cmfd_bcast() + #ifdef PETSC ! execute snes solver call cmfd_snes_execute() @@ -57,8 +62,6 @@ contains ! write vtk file !if(.not. cmfd_only) call write_cmfd_vtk() -! end if - #ifdef PETSC ! finalize slepc if (current_cycle == n_cycles) call SlepcFinalize(ierr) @@ -69,4 +72,44 @@ contains end subroutine execute_cmfd +!============================================================================== +! EXECUTE_CMFD +!============================================================================== + + subroutine cmfd_bcast() + + use cmfd_header, only: allocate_cmfd + + integer :: nx ! number of mesh cells in x direction + integer :: ny ! number of mesh cells in y direction + integer :: nz ! number of mesh cells in z direction + integer :: ng ! number of energy groups + integer rank + + ! extract spatial and energy indices from object + nx = cmfd % indices(1) + ny = cmfd % indices(2) + nz = cmfd % indices(3) + ng = cmfd % indices(4) + + ! initialize data + call allocate_cmfd(cmfd) + + ! 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) + + ! sync up procs + call MPI_Barrier(MPI_COMM_WORLD,mpi_err) + + end subroutine cmfd_bcast + end module cmfd_execute diff --git a/src/cmfd_loss_operator.F90 b/src/cmfd_loss_operator.F90 index 664cc97da2..dc5a3be0d7 100644 --- a/src/cmfd_loss_operator.F90 +++ b/src/cmfd_loss_operator.F90 @@ -36,11 +36,11 @@ contains call get_M_indices(this) ! set up M operator - call MatCreateSeqAIJ(PETSC_COMM_SELF,this%n,this%n,this%nnz, & - & PETSC_NULL_INTEGER,this%M,ierr) + call MatCreateMPIAIJ(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,this%n,this%n,this%nnz, & + & PETSC_NULL_INTEGER,this%nnz,PETSC_NULL_INTEGER,this%M,ierr) call MatSetOption(this%M,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE,ierr) call MatSetOption(this%M,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr) - call MatSetOption(this%M,MAT_USE_HASH_TABLE,PETSC_TRUE,ierr) +! call MatSetOption(this%M,MAT_USE_HASH_TABLE,PETSC_TRUE,ierr) end subroutine init_M_operator @@ -78,7 +78,7 @@ contains subroutine build_loss_matrix(this) - use global, only: cmfd,cmfd_coremap + use global, only: cmfd,cmfd_coremap,mpi_err type(loss_operator) :: this @@ -89,7 +89,6 @@ contains integer :: g ! iteration counter for groups integer :: l ! iteration counter for leakages integer :: h ! energy group when doing scattering - integer :: cell_mat_idx ! matrix index of current cell integer :: neig_mat_idx ! matrix index of neighbor cell integer :: scatt_mat_idx ! matrix index for h-->g scattering terms integer :: bound(6) ! vector for comparing when looking for bound @@ -98,6 +97,9 @@ contains integer :: neig_idx(3) ! spatial indices of neighbour integer :: shift_idx ! parameter to shift index by +1 or -1 integer :: kount ! integer for counting values in vector + integer :: row_start ! the first local row on the processor + integer :: row_finish ! the last local row on the processor + integer :: irow ! iteration counter over row real(8) :: totxs ! total macro cross section real(8) :: scattxsgg ! scattering macro cross section g-->g real(8) :: scattxshg ! scattering macro cross section h-->g @@ -109,158 +111,135 @@ contains real(8) :: jnet ! net leakage from jo real(8) :: val ! temporary variable before saving to matrix - ! initialize matrix for building - call MatAssemblyBegin(this%M,MAT_FLUSH_ASSEMBLY,ierr) - ! create single vector of these indices for boundary calculation nxyz(1,:) = (/1,nx/) nxyz(2,:) = (/1,ny/) nxyz(3,:) = (/1,nz/) + ! get row bounds for this processor + call MatGetOwnershipRange(this%M,row_start,row_finish,ierr) + ! begin iteration loops - ZLOOP: do k = 1,nz + ROWS: do irow = row_start,row_finish-1 - YLOOP: do j = 1,ny + ! get indices for that row + call matrix_to_indices(irow,g,i,j,k) - XLOOP: do i = 1,nx + ! retrieve cell data + totxs = cmfd%totalxs(g,i,j,k) + scattxsgg = cmfd%scattxs(g,g,i,j,k) + dtilde = cmfd%dtilde(:,g,i,j,k) + hxyz = cmfd%hxyz(:,i,j,k) - ! check if not including reflector + ! check and get dhat + if (allocated(cmfd%dhat)) then + dhat = cmfd%dhat(:,g,i,j,k) + else + dhat = 0.0_8 + end if + + ! create boundary vector + bound = (/i,i,j,j,k,k/) + + ! begin loop over leakages + ! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z + LEAK: do l = 1,6 + + ! define (x,y,z) and (-,+) indices + xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3 + dir_idx = 2 - mod(l,2) ! -=1, +=2 + + ! calculate spatial indices of neighbor + neig_idx = (/i,j,k/) ! begin with i,j,k + shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1 + neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx) + + ! check for global boundary + if (bound(l) /= nxyz(xyz_idx,dir_idx)) then + + ! check for core map if (cmfd_coremap) then - ! check if at a reflector - if (cmfd % coremap(i,j,k) == 99999) then - cycle + ! check that neighbor is not reflector + if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) /= & + & 99999) then + + ! compute leakage coefficient for neighbor + jn = -dtilde(l) + shift_idx*dhat(l) + + ! get neighbor matrix index + call indices_to_matrix(g,neig_idx(1),neig_idx(2),neig_idx(3), & + & neig_mat_idx) + + ! compute value and record to bank + val = jn/hxyz(xyz_idx) + + ! record value in matrix + call MatSetValue(this%M,irow,neig_mat_idx-1,val, & + & INSERT_VALUES,ierr) + end if + else + + ! compute leakage coefficient for neighbor + jn = -dtilde(l) + shift_idx*dhat(l) + + ! get neighbor matrix index + call indices_to_matrix(g,neig_idx(1),neig_idx(2),neig_idx(3), & + & neig_mat_idx) + + ! compute value and record to bank + val = jn/hxyz(xyz_idx) + + ! record value in matrix + call MatSetValue(this%M,irow,neig_mat_idx-1,val, & + & INSERT_VALUES,ierr) + end if - GROUP: do g = 1,ng + end if - ! get matrix index of cell - cell_mat_idx = get_matrix_idx(g,i,j,k) + ! compute leakage coefficient for target + jo(l) = shift_idx*dtilde(l) + dhat(l) - ! retrieve cell data - totxs = cmfd%totalxs(g,i,j,k) - scattxsgg = cmfd%scattxs(g,g,i,j,k) - dtilde = cmfd%dtilde(:,g,i,j,k) - hxyz = cmfd%hxyz(:,i,j,k) + end do LEAK - ! check and get dhat - if (allocated(cmfd%dhat)) then - dhat = cmfd%dhat(:,g,i,j,k) - else - dhat = 0.0 - end if + ! calate net leakage coefficient for target + jnet = (jo(2) - jo(1))/hxyz(1) + (jo(4) - jo(3))/hxyz(2) + & + & (jo(6) - jo(5))/hxyz(3) - ! create boundary vector - bound = (/i,i,j,j,k,k/) + ! calculate loss of neutrons + val = jnet + totxs - scattxsgg - ! begin loop over leakages - ! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z - LEAK: do l = 1,6 + ! record diagonal term + call MatSetValue(this%M,irow,irow,val,INSERT_VALUES,ierr) - ! define (x,y,z) and (-,+) indices - xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3 - dir_idx = 2 - mod(l,2) ! -=1, +=2 + ! begin loop over off diagonal in-scattering + SCATTR: do h = 1,ng - ! calculate spatial indices of neighbor - neig_idx = (/i,j,k/) ! begin with i,j,k - shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1 - neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx) + ! cycle though if h=g + if (h == g) then + cycle + end if - ! check for global boundary - if (bound(l) /= nxyz(xyz_idx,dir_idx)) then + ! get neighbor matrix index + call indices_to_matrix(h,i,j,k,scatt_mat_idx) - ! check for core map - if (cmfd_coremap) then + ! get scattering macro xs + scattxshg = cmfd%scattxs(h,g,i,j,k) - ! check that neighbor is not reflector - if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) /= & - & 99999) then + ! record value in matrix (negate it) + val = -scattxshg - ! compute leakage coefficient for neighbor - jn = -dtilde(l) + shift_idx*dhat(l) + call MatSetValue(this%M,irow,scatt_mat_idx-1,val, INSERT_VALUES,ierr) - ! get neighbor matrix index - neig_mat_idx = get_matrix_idx(g,neig_idx(1),neig_idx(2), & - & neig_idx(3)) + end do SCATTR - ! compute value and record to bank - val = jn/hxyz(xyz_idx) + end do ROWS - ! record value in matrix - call MatSetValue(this%M,cell_mat_idx-1,neig_mat_idx-1,val, & - & INSERT_VALUES,ierr) - - end if - - else - - ! compute leakage coefficient for neighbor - jn = -dtilde(l) + shift_idx*dhat(l) - - ! get neighbor matrix index - neig_mat_idx = get_matrix_idx(g,neig_idx(1),neig_idx(2), & - & neig_idx(3)) - - ! compute value and record to bank - val = jn/hxyz(xyz_idx) - - ! record value in matrix - call MatSetValue(this%M,cell_mat_idx-1,neig_mat_idx-1,val, & - & INSERT_VALUES,ierr) - - end if - - end if - - ! compute leakage coefficient for target - jo(l) = shift_idx*dtilde(l) + dhat(l) - - end do LEAK - - ! calate net leakage coefficient for target - jnet = (jo(2) - jo(1))/hxyz(1) + (jo(4) - jo(3))/hxyz(2) + & - & (jo(6) - jo(5))/hxyz(3) - - ! calculate loss of neutrons - val = jnet + totxs - scattxsgg - - ! record diagonal term - call MatSetValue(this%M,cell_mat_idx-1,cell_mat_idx-1,val, & - & INSERT_VALUES,ierr) - - ! begin loop over off diagonal in-scattering - SCATTR: do h = 1,ng - - ! cycle though if h=g - if (h == g) then - cycle - end if - - ! get matrix index of in-scatter - scatt_mat_idx = get_matrix_idx(h,i,j,k) - - ! get scattering macro xs - scattxshg = cmfd%scattxs(h,g,i,j,k) - - ! record value in matrix (negate it) - val = -scattxshg - - call MatSetValue(this%M,cell_mat_idx-1,scatt_mat_idx-1,val, & - & INSERT_VALUES,ierr) - - end do SCATTR - - end do GROUP - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - ! finalize matrix assembly + ! assemble matrix + call MatAssemblyBegin(this%M,MAT_FLUSH_ASSEMBLY,ierr) call MatAssemblyEnd(this%M,MAT_FINAL_ASSEMBLY,ierr) ! print out operator to file @@ -269,40 +248,53 @@ contains end subroutine build_loss_matrix !=============================================================================== -! GET_MATRIX_IDX takes (x,y,z,g) indices and computes location in matrix +! INDICES_TO_MATRIX takes (x,y,z,g) indices and computes location in matrix !=============================================================================== - function get_matrix_idx(g,i,j,k) + subroutine indices_to_matrix(g,i,j,k,matidx) use global, only: cmfd,cmfd_coremap - ! arguments - integer :: get_matrix_idx ! the index location in matrix + integer :: matidx ! the index location in matrix integer :: i ! current x index integer :: j ! current y index integer :: k ! current z index integer :: g ! current group index - ! local variables - integer :: nidx ! index in matrix - ! check if coremap is used if (cmfd_coremap) then ! get idx from core map - nidx = ng*(cmfd % coremap(i,j,k)) - (ng - g) + matidx = ng*(cmfd % coremap(i,j,k)) - (ng - g) else ! compute index - nidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1) + matidx = g + ng*(i - 1) + ng*nx*(j - 1) + ng*nx*ny*(k - 1) end if - ! record value to function - get_matrix_idx = nidx + end subroutine indices_to_matrix - end function get_matrix_idx +!=============================================================================== +! MATRIX_TO_INDICES +!=============================================================================== + + subroutine matrix_to_indices(irow,g,i,j,k) + + integer :: i ! iteration counter for x + integer :: j ! iteration counter for y + integer :: k ! iteration counter for z + integer :: g ! iteration counter for groups + integer :: irow ! iteration counter over row (0 reference) + + ! compute indices + g = irow/(nx*ny*nz) + 1 + i = mod(irow,nx*ny*nz)/(ny*nz) + 1 + j = mod(irow,ny*nz)/nz + 1 + k = mod(irow,nz) + 1 + + end subroutine matrix_to_indices !=============================================================================== ! PRINT_M_OPERATOR @@ -317,11 +309,11 @@ contains PetscViewer :: viewer ! write out matrix in binary file (debugging) - call PetscViewerBinaryOpen(PETSC_COMM_SELF,trim(path_input)//'lossmat.bin' & + call PetscViewerBinaryOpen(PETSC_COMM_WORLD,trim(path_input)//'lossmat.bin' & & ,FILE_MODE_WRITE,viewer,ierr) call MatView(this%M,viewer,ierr) call PetscViewerDestroy(viewer,ierr) - +stop end subroutine print_M_operator !============================================================================== diff --git a/src/cmfd_slepc_solver.F90 b/src/cmfd_slepc_solver.F90 index 51c14ed6d7..dd0262caa8 100644 --- a/src/cmfd_slepc_solver.F90 +++ b/src/cmfd_slepc_solver.F90 @@ -30,10 +30,7 @@ contains !=============================================================================== subroutine cmfd_slepc_execute() - use global, only:time_cmfd - use timing -call timer_start(time_cmfd) ! initialize data call init_data()