diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 47c2b50cd6..2c0a7381df 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -1,1643 +1,55 @@ module cmfd_execute - use cmfd_utils, only: get_matrix_idx,neutron_balance,set_coremap, & - & get_reflector_albedo,write_hdf5,read_hdf5, & - & allocate_cmfd,write_vtk + use cmfd_data + use cmfd_power_solver, only: cmfd_power_execute + use cmfd_slepc_solver, only: cmfd_slepc_execute + use cmfd_snes_solver, only: cmfd_snes_execute use global - use mesh, only: mesh_indices_to_bin - use mesh_header, only: StructuredMesh - use tally_header, only: TallyObject, TallyScore implicit none + +#include +#include +#include + contains -!=============================================================================== -! EXECUTE_CMFD is the highest level routine that controls CMFD calculation -!=============================================================================== +!============================================================================== +! EXECUTE_CMFD +!============================================================================== subroutine execute_cmfd() - ! allocate cmfd object - call allocate_cmfd() + integer :: ierr ! petsc error code - ! calculate all cross sections based on reaction rates from last batch - call compute_xs() + ! set up cmfd + call set_up_cmfd() - ! print out flux for debugginb -! write(100,*) cmfd % flux -! write(101,*) cmfd % totalxs -! write(102,*) cmfd % p1scattxs -! write(103,*) cmfd % scattxs -! write(104,*) cmfd % nfissxs -! write(105,*) cmfd % hxyz -! write(106,*) cmfd % current - - ! write out neutron balance - call neutron_balance() - - ! check for core map - if (allocated(cmfd % coremap)) then - call set_coremap() - end if - - ! compute dtilde terms - call compute_diffcoef() - - ! set dhats to zero - call compute_dhat() - - ! print dtilde and dhat -! write(107,*) cmfd % dtilde -! write(108,*) cmfd % dhat - - ! write cmfd object to hdf5 file - call write_hdf5() - - ! solve diffusion equation - call cmfd_solver() - - ! call vtk output - call write_vtk() - - end subroutine - -!=============================================================================== -! COMPUTE_XS takes tallies and computes macroscopic cross sections -!=============================================================================== - - subroutine compute_xs() - - 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 :: i ! iteration counter for x - integer :: j ! iteration counter for y - integer :: k ! iteration counter for z - integer :: g ! iteration counter for g - integer :: h ! iteration counter for outgoing groups - integer :: ijk(3) ! indices for mesh cell - integer :: score_index ! index to pull from tally object - integer :: bins(TALLY_TYPES) ! bins for filters - - real(8) :: flux ! temp variable for flux - - type(TallyObject), pointer :: t ! pointer for tally object - type(StructuredMesh), pointer :: m ! pointer for mesh object - - ! extract spatial and energy indices from object - nx = cmfd % indices(1) - ny = cmfd % indices(2) - nz = cmfd % indices(3) - ng = cmfd % indices(4) - - ! begin loop around space and energy groups - ZLOOP: do k = 1,nz - - YLOOP: do j = 1,ny - - XLOOP: do i = 1,nx - - OUTGROUP: do h = 1,ng - - ! begin with first tally - t => tallies(1) - m => meshes(t % mesh) - - ! set mesh widths - cmfd % hxyz(1,:,:,:) = m % width(1) ! set x width - cmfd % hxyz(2,:,:,:) = m % width(2) ! set y width - cmfd % hxyz(3,:,:,:) = m % width(3) ! set z width - - ! reset all bins to 1 - bins = 1 - - ! set ijk as mesh indices - ijk = (/ i, j, k /) - - ! get bin number for mesh indices - bins(T_MESH) = mesh_indices_to_bin(m,ijk) - - ! apply energy in filter - bins(T_ENERGYIN) = ng - h + 1 - - ! calculate score index from bins - score_index = sum((bins - 1) * t%stride) + 1 - - ! get flux - flux = t % scores(score_index,1) % val - cmfd % flux(h,i,j,k) = flux - - ! detect zero flux - if ((flux - 0.0D0) < 1.0D-10) then - if (.not. allocated(cmfd%coremap)) then - write(*,*) 'Fatal: detected zero flux without coremap' - stop - else - write(*,*) 'Warning: detected zero flux at:',i,j,k - flux = 99999.0D0 - if (.not. cmfd%coremap(i,j,k) == 99999) then - write(*,*) 'Fatal: need to check core map with zero flux' - stop - end if - end if - end if - - ! get total rr and convert to total xs - cmfd % totalxs(h,i,j,k) = t % scores(score_index,2) % val / flux - - ! get p1 scatter rr and convert to p1 scatter xs - cmfd % p1scattxs(h,i,j,k) = t % scores(score_index,3) % val / flux - - ! calculate diffusion coefficient - cmfd % diffcof(h,i,j,k) = 1/(3*(cmfd % totalxs(h,i,j,k) - & - & cmfd % p1scattxs(h,i,j,k))) - cmfd % diffcof(h,i,j,k) = 1/(3*(cmfd % totalxs(h,i,j,k))) - - ! begin loop to get energy out tallies - INGROUP: do g = 1,ng - - ! associate tally pointer to energy out tally object - t => tallies(2) - - ! set energy out bin - bins(T_ENERGYOUT) = ng - g + 1 - - ! calculate score index from bins - score_index = sum((bins - 1) * t%stride) + 1 - - ! get scattering - cmfd % scattxs(h,g,i,j,k) = t % scores(score_index,1) % val / flux - - ! get nu-fission - cmfd % nfissxs(h,g,i,j,k) = t % scores(score_index,2) % val / flux - - end do INGROUP - - ! extract surface currents - t => tallies(3) - - ! initialize and filter for energy - bins = 1 - bins(TS_ENERGYIN) = ng - h + 1 - - ! left surface - bins(1:3) = (/ i-1, j, k /) + 1 - bins(TS_SURFACE) = IN_RIGHT - score_index = sum((bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(1,h,i,j,k) = t % scores(score_index,1) % val - bins(TS_SURFACE) = OUT_RIGHT - score_index = sum((bins - 1) * t % stride) + 1 ! incoming - cmfd % current(2,h,i,j,k) = t % scores(score_index,1) % val - - ! right surface - bins(1:3) = (/ i, j, k /) + 1 - bins(TS_SURFACE) = IN_RIGHT - score_index = sum((bins - 1) * t % stride) + 1 ! incoming - cmfd % current(3,h,i,j,k) = t % scores(score_index,1) % val - bins(TS_SURFACE) = OUT_RIGHT - score_index = sum((bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(4,h,i,j,k) = t % scores(score_index,1) % val - - ! back surface - bins(1:3) = (/ i, j-1, k /) + 1 - bins(TS_SURFACE) = IN_FRONT - score_index = sum((bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(5,h,i,j,k) = t % scores(score_index,1) % val - bins(TS_SURFACE) = OUT_FRONT - score_index = sum((bins - 1) * t % stride) + 1 ! incoming - cmfd % current(6,h,i,j,k) = t % scores(score_index,1) % val - - ! front surface - bins(1:3) = (/ i, j, k /) + 1 - bins(TS_SURFACE) = IN_FRONT - score_index = sum((bins - 1) * t % stride) + 1 ! incoming - cmfd % current(7,h,i,j,k) = t % scores(score_index,1) % val - bins(TS_SURFACE) = OUT_FRONT - score_index = sum((bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(8,h,i,j,k) = t % scores(score_index,1) % val - - ! bottom surface - bins(1:3) = (/ i, j, k-1 /) + 1 - bins(TS_SURFACE) = IN_TOP - score_index = sum((bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(9,h,i,j,k) = t % scores(score_index,1) % val - bins(TS_SURFACE) = OUT_TOP - score_index = sum((bins - 1) * t % stride) + 1 ! incoming - cmfd % current(10,h,i,j,k) = t % scores(score_index,1) % val - - ! top surface - bins(1:3) = (/ i, j, k /) + 1 - bins(TS_SURFACE) = IN_TOP - score_index = sum((bins - 1) * t % stride) + 1 ! incoming - cmfd % current(11,h,i,j,k) = t % scores(score_index,1) % val - bins(TS_SURFACE) = OUT_TOP - score_index = sum((bins - 1) * t % stride) + 1 ! outgoing - cmfd % current(12,h,i,j,k) = t % scores(score_index,1) % val - - end do OUTGROUP - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - end subroutine compute_xs - -!=============================================================================== -! COMPUTE_DIFFCOEF computes the diffusion coupling coefficient -!=============================================================================== - - subroutine compute_diffcoef() - - ! local variables - integer :: nx ! maximum number of cells in x direction - integer :: ny ! maximum number of cells in y direction - integer :: nz ! maximum number of cells in z direction - integer :: ng ! maximum number of energy groups - integer :: nxyz(3,2) ! single vector containing boundary locations - 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 :: l ! iteration counter for leakages - integer :: xyz_idx ! index for determining if x,y or z leakage - integer :: dir_idx ! index for determining - or + face of cell - integer :: shift_idx ! parameter to shift index by +1 or -1 - integer :: neig_idx(3) ! spatial indices of neighbour - integer :: bound(6) ! vector containing indices for boudary check - real(8) :: albedo(6) ! albedo vector with global boundaries - real(8) :: cell_totxs ! total cross section of current ijk cell - real(8) :: cell_dc ! diffusion coef of current cell - real(8) :: cell_hxyz(3) ! cell dimensions of current ijk cell - real(8) :: neig_totxs ! total xs of neighbor cell - real(8) :: neig_dc ! diffusion coefficient of neighbor cell - real(8) :: neig_hxyz(3) ! cell dimensions of neighbor cell - real(8) :: dtilde ! finite difference coupling parameter - real(8) :: ref_albedo ! albedo to reflector - - ! get maximum of spatial and group indices - nx = cmfd%indices(1) - ny = cmfd%indices(2) - nz = cmfd%indices(3) - ng = cmfd%indices(4) - - ! create single vector of these indices for boundary calculation - nxyz(1,:) = (/1,nx/) - nxyz(2,:) = (/1,ny/) - nxyz(3,:) = (/1,nz/) - - ! get boundary condition information - albedo = cmfd%albedo - - ! geting loop over group and spatial indices - ZLOOP: do k = 1,nz - - YLOOP: do j = 1,ny - - XLOOP: do i = 1,nx - - GROUP: do g = 1,ng - - ! get cell data - cell_dc = cmfd%diffcof(g,i,j,k) - cell_hxyz = cmfd%hxyz(:,i,j,k) - - - ! setup of vector to identify boundary conditions - bound = (/i,i,j,j,k,k/) - - ! begin loop around sides of cell for leakage - LEAK: do l = 1,6 - - ! define xyz and +/- indices - xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3 - dir_idx = 2 - mod(l,2) ! -=1, +=2 - shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1 - - ! check if at a boundary - if (bound(l) == nxyz(xyz_idx,dir_idx)) then - - ! compute dtilde - dtilde = (2*cell_dc*(1-albedo(l)))/(4*cell_dc*(1+albedo(l)) + & - & (1-albedo(l))*cell_hxyz(xyz_idx)) - - else ! not a boundary - - ! compute neighboring cell indices - neig_idx = (/i,j,k/) ! begin with i,j,k - neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx) - - ! get neigbor cell data - neig_dc = cmfd%diffcof(g,neig_idx(1),neig_idx(2),neig_idx(3)) - neig_hxyz = cmfd%hxyz(:,neig_idx(1),neig_idx(2),neig_idx(3)) - - ! check for fuel-reflector interface - if (allocated(cmfd % coremap)) then - - if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) == & - & 99999 .and. cmfd % coremap(i,j,k) /= 99999) then - - ! get albedo - ref_albedo = get_reflector_albedo(l,g,i,j,k) - - ! compute dtilde - dtilde = (2*cell_dc*(1-ref_albedo))/(4*cell_dc*(1+ & - & ref_albedo)+(1-ref_albedo)*cell_hxyz(xyz_idx)) - - else ! not next to a reflector or no core map - - ! compute dtilde - dtilde = (2*cell_dc*neig_dc)/(neig_hxyz(xyz_idx)*cell_dc + & - & cell_hxyz(xyz_idx)*neig_dc) - - end if - - else ! no core map - - ! compute dtilde - dtilde = (2*cell_dc*neig_dc)/(neig_hxyz(xyz_idx)*cell_dc + & - & cell_hxyz(xyz_idx)*neig_dc) - - end if - - end if - - ! record dtilde in cmfd object - cmfd%dtilde(l,g,i,j,k) = dtilde - - end do LEAK - - end do GROUP - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - end subroutine compute_diffcoef - -!=============================================================================== -! COMPUTE_DHAT computes the nonlinear coupling coefficient -!=============================================================================== - - subroutine compute_dhat() - - ! local variables - integer :: nx ! maximum number of cells in x direction - integer :: ny ! maximum number of cells in y direction - integer :: nz ! maximum number of cells in z direction - integer :: ng ! maximum number of energy groups - integer :: nxyz(3,2) ! single vector containing boundary locations - 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 :: l ! iteration counter for leakages - integer :: xyz_idx ! index for determining if x,y or z leakage - integer :: dir_idx ! index for determining - or + face of cell - integer :: shift_idx ! parameter to shift index by +1 or -1 - integer :: neig_idx(3) ! spatial indices of neighbour - integer :: bound(6) ! vector containing indices for boudary check - real(8) :: cell_dtilde(6) ! cell dtilde for each face - real(8) :: cell_flux ! flux in current cell - real(8) :: current(12) ! area integrated cell current at each face - real(8) :: net_current ! net current on a face - real(8) :: neig_flux ! flux in neighbor cell - real(8) :: dhat ! dhat equivalence parameter - - ! get maximum of spatial and group indices - nx = cmfd%indices(1) - ny = cmfd%indices(2) - nz = cmfd%indices(3) - ng = cmfd%indices(4) - - ! create single vector of these indices for boundary calculation - nxyz(1,:) = (/1,nx/) - nxyz(2,:) = (/1,ny/) - nxyz(3,:) = (/1,nz/) - - ! geting loop over group and spatial indices - ZLOOP: do k = 1,nz - - YLOOP: do j = 1,ny - - XLOOP: do i = 1,nx - - GROUP: do g = 1,ng - - ! get cell data - cell_dtilde = cmfd%dtilde(:,g,i,j,k) - cell_flux = cmfd%flux(g,i,j,k)/product(cmfd%hxyz(:,i,j,k)) - current = cmfd%current(:,g,i,j,k) - - ! setup of vector to identify boundary conditions - bound = (/i,i,j,j,k,k/) - - ! begin loop around sides of cell for leakage - LEAK: do l = 1,6 - - ! define xyz and +/- indices - xyz_idx = int(ceiling(real(l)/real(2))) ! x=1, y=2, z=3 - dir_idx = 2 - mod(l,2) ! -=1, +=2 - shift_idx = -2*mod(l,2) +1 ! shift neig by -1 or +1 - - ! calculate net current on l face (divided by surf area) - net_current = (current(2*l) - current(2*l-1)) / & - & product(cmfd%hxyz(:,i,j,k)) * & - & cmfd%hxyz(xyz_idx,i,j,k) - - ! check if at a boundary - if (bound(l) == nxyz(xyz_idx,dir_idx)) then - - ! compute dhat - dhat = (net_current - shift_idx*cell_dtilde(l)*cell_flux) / & - & cell_flux - - else ! not a boundary - - ! compute neighboring cell indices - neig_idx = (/i,j,k/) ! begin with i,j,k - neig_idx(xyz_idx) = shift_idx + neig_idx(xyz_idx) - - ! get neigbor flux - neig_flux = cmfd%flux(g,neig_idx(1),neig_idx(2),neig_idx(3)) / & - product(cmfd%hxyz(:,neig_idx(1),neig_idx(2),neig_idx(3))) - - ! check for fuel-reflector interface - if (allocated(cmfd % coremap)) then - - if (cmfd % coremap(neig_idx(1),neig_idx(2),neig_idx(3)) == & - & 99999 .and. cmfd % coremap(i,j,k) /= 99999) then - - ! compute dhat - dhat = (net_current - shift_idx*cell_dtilde(l)*cell_flux) /& - & cell_flux - - else ! not a fuel-reflector interface - - ! compute dhat - dhat = (net_current + shift_idx*cell_dtilde(l)* & - & (neig_flux - cell_flux))/(neig_flux + cell_flux) - - end if - - else ! not for fuel-reflector case - - ! compute dhat - dhat = (net_current + shift_idx*cell_dtilde(l)* & - & (neig_flux - cell_flux))/(neig_flux + cell_flux) - - end if - - end if - - ! record dtilde in cmfd object - cmfd%dhat(l,g,i,j,k) = dhat - - end do LEAK - - end do GROUP - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - end subroutine compute_dhat - -!=============================================================================== -! ACCUMULATE accumulates cycle to cycle diffusion parameters -!=============================================================================== - - subroutine accumulate() - - end subroutine accumulate - -!=============================================================================== -! CMFD_SOLVER in the main power iteration routine for the cmfd calculation -!=============================================================================== - - subroutine cmfd_solver() - -use output, only: write_message -use timing, only: timer_start, timer_stop - -#include -#include -#include - - Mat :: M ! loss matrix - Mat :: F ! production matrix - Vec :: phi ! eigenvector - EPS :: eps ! slepc eigenvalue object - ST :: st ! slepc spectral trans object - KSP :: ksp ! linear solver object - PC :: pc ! preconditioner object - PetscViewer :: viewer ! viewer for eigenvector - PetscScalar, pointer :: phi_v(:) ! pointer to eigenvector info - - integer :: ierr ! error flag - integer :: i_eig=0 ! eigenvalue to extract - integer :: its ! number of iterations to eigenvalue solve - real(8) :: keff ! eigenvalue - - ! initialize PETSc + ! initialize slepc/petsc call SlepcInitialize(PETSC_NULL_CHARACTER,ierr) - ! initialize matrices and vectors - message = "Initializing and building matrices..." - call write_message(1) - call timer_start(time_mat) - call init_data(M,F,phi) - call init_solver(eps,st,ksp,pc) + ! execute diffusion solver + call cmfd_power_execute() - ! set up M loss matrix - call loss_matrix(M) + ! print results + print *,'Power Eigenvalue is:',cmfd%keff - ! set up F production matrix - call prod_matrix(F) - call timer_stop(time_mat) + ! execute slepc solver + call cmfd_slepc_execute() - ! set eigenvalue operators - call EPSSetOperators(eps,F,M,ierr) - - ! begin timer for eigenvalue solution - message = "Beginning Eigenvalue Calculation..." - call write_message(1) - call timer_start(time_eigen) + ! print results + print *,'SLEPC Eigenvalue is:',cmfd%keff - ! solve system - call EPSSolve(eps,ierr) + ! execute snes solver + call cmfd_snes_execute() - ! end eigenvalue timer - call timer_stop(time_eigen) + ! print results + print *,'SNES Eigenvalue is:',cmfd%keff - ! extract run information - call EPSGetIterationNumber(eps,its,ierr) - call EPSGetEigenpair(eps,i_eig,keff,PETSC_NULL,phi,PETSC_NULL_OBJECT,ierr) - - ! print output to user - print *, "Eigenvalue is: ",keff - print *, "Number of iterations: ",its - print *, "Matrix building time (s): ",time_mat%elapsed - print *, "Eigenvalue solution time (s): ",time_eigen%elapsed - - ! convert petsc phi object to cmfd obj - call VecGetArrayF90(phi,phi_v,ierr) - cmfd%phi = phi_v - call VecRestoreArrayF90(phi,phi_v,ierr) - - ! compute source pdf and record in cmfd object - ! call source_pdf(S_n) - - ! output answers - call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'fluxvec.bin',FILE_MODE_WRITE, & - viewer,ierr) - call VecView(phi,viewer,ierr) - - ! solve nonlinear system - call solve_cmfd_nonlinear(phi,keff) - - ! destroy all PETSc objects - call destroy_data(M,F,phi,eps,viewer) - - ! finalize SLEPc + ! finalize slepc call SlepcFinalize(ierr) - end subroutine cmfd_solver - -!=============================================================================== -! INIT_DATA allocates matrices vectors for CMFD solution -!=============================================================================== - - subroutine init_data(M,F,phi) - -#include - - ! arguments - Mat :: M ! loss matrix - Mat :: F ! production matrix - Vec :: phi ! eigenvector - - ! local variables - integer :: n ! dimensions of matrix - integer :: nz_M ! number of nonzeros in loss matrix - integer :: nz_F ! number of nonzeros in prod matrix - integer :: ierr ! error flag - integer :: nx ! maximum number of x cells - integer :: ny ! maximum number of y cells - integer :: nz ! maximum number of z cells - integer :: ng ! maximum number of groups - integer :: nzM ! max number of nonzeros in a row for M - integer :: nzF ! max number of nonzeros in a row for F - - ! get maximum number of cells in each direction - nx = cmfd%indices(1) - ny = cmfd%indices(2) - nz = cmfd%indices(3) - ng = cmfd%indices(4) - - ! calculate dimensions of matrix - if (allocated(cmfd % coremap)) then - n = cmfd % mat_dim * ng - else - n = nx*ny*nz*ng - end if - - ! maximum number of nonzeros in each matrix - nzM = 7+ng-1 - nzF = ng - - ! set up loss matrix - call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,nzM,PETSC_NULL_INTEGER,M,ierr) - call MatSetOption(M,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE,ierr) - call MatSetOption(M,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr) - call MatSetOption(M,MAT_USE_HASH_TABLE,PETSC_TRUE,ierr) - - ! set up production matrix - call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,nzF,PETSC_NULL_INTEGER,F,ierr) - call MatSetOption(F,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE,ierr) - call MatSetOption(F,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr) - call MatSetOption(F,MAT_USE_HASH_TABLE,PETSC_TRUE,ierr) - - ! set up eigenvector - call VecCreate(PETSC_COMM_SELF,phi,ierr) - call VecSetSizes(phi,PETSC_DECIDE,n,ierr) - call VecSetFromOptions(phi,ierr) - - ! also allocate in cmfd object - if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n)) - - end subroutine init_data - -!=============================================================================== -! INIT_SOLVER setups the eigenvalue problem solvers -!=============================================================================== - - subroutine init_solver(eps,st,ksp,pc) - -#include -#include -#include - - integer :: ierr ! error flag - EPS :: eps ! slepc eigenvalue object - ST :: st ! slepc spetral trans object - KSP :: ksp ! linear solver object - PC :: pc ! preconditioner object - PetscViewer :: viewer ! viewer for answer - character(LEN=20) :: epstype,sttype,ksptype,pctype - - ! create EPS Object - call EPSCreate(PETSC_COMM_WORLD,eps,ierr) - call EPSSetProblemType(eps,EPS_GNHEP,ierr) - call EPSSetType(eps,EPSPOWER,ierr) - call EPSSetFromOptions(eps,ierr) - call EPSSetWhichEigenpairs(eps,EPS_LARGEST_MAGNITUDE,ierr) - - ! get ST, KSP and PC objects - call EPSGetST(eps,st,ierr) - call STGetKSP(st,ksp,ierr) - call KSPGetPC(ksp,pc,ierr) - - ! set GMRES default - call KSPSetType(ksp,KSPGMRES,ierr) - - ! set precursor type - call PCSetType(pc,PCILU,ierr) - call PCFactorSetLevels(pc,4,ierr) - call PCSetFromOptions(pc,ierr) - - ! get all types and print - call EPSGetType(eps,epstype,ierr) - call STGetType(st,sttype,ierr) - call KSPGetType(ksp,ksptype,ierr) - call PCGetType(pc,pctype,ierr) - print *,'EPS TYPE IS: ',epstype - print *,'ST TYPE IS: ',sttype - print *,'KSP TYPE IS: ',ksptype - print *,'PC TYPE IS: ',pctype - - end subroutine init_solver - -!=============================================================================== -! SOLVE_CMFD_NONLINEAR solves nonlinear system with JFNK -!=============================================================================== - - subroutine solve_cmfd_nonlinear(phi,keff) - -#include - - Vec :: phi ! flux vector from slepc - real(8) :: keff ! eigenvalue - - Mat :: jac ! jacobian matrix - Vec :: res ! residual vector - Vec :: x ! results - KSP :: ksp ! linear solver context - PC :: pc ! preconditioner - SNES :: snes ! nonlinear solver context - integer :: ierr ! error flag - integer :: n ! dimensions of matrix - integer :: nx ! maximum number of x cells - integer :: ny ! maximum number of y cells - integer :: nz ! maximum number of z cells - integer :: ng ! maximum number of groups - PetscViewer :: viewer ! viewer for output - - real(8), pointer :: fptr(:) ! flux pointer - real(8), pointer :: xptr(:) ! solution pointer - - ! get maximum number of cells in each direction - nx = cmfd%indices(1) - ny = cmfd%indices(2) - nz = cmfd%indices(3) - ng = cmfd%indices(4) - - ! calculate dimensions of matrix - if (allocated(cmfd % coremap)) then - n = cmfd % mat_dim * ng - else - n = nx*ny*nz*ng - end if - - ! create PETSc vectors - call VecCreate(PETSC_COMM_SELF,res,ierr) - call VecSetSizes(res,PETSC_DECIDE,n+1,ierr) - call VecSetFromOptions(res,ierr) - call VecCreate(PETSC_COMM_SELF,x,ierr) - call VecSetSizes(x,PETSC_DECIDE,n+1,ierr) - call VecSetFromOptions(x,ierr) - - ! get pointers - call VecGetArrayF90(x,xptr,ierr) - call VecGetArrayF90(phi,fptr,ierr) - - ! set new xptr - xptr(1:n) = fptr(1:n) - xptr(n+1) = 1.0_8/keff - - ! restore array - call VecRestoreArrayF90(x,xptr,ierr) - call VecRestoreArrayF90(phi,fptr,ierr) - - ! create SNES context - call SNESCreate(PETSC_COMM_SELF,snes,ierr) - - ! set the residual function - call SNESSetFunction(snes,res,compute_nonlinear_residual,PETSC_NULL,ierr) - - ! set GMRES solver - call SNESGetKSP(snes,ksp,ierr) - call KSPSetType(ksp,KSPGMRES,ierr) - - ! set preconditioner - call KSPGetPC(ksp,pc,ierr) - call PCSetType(pc,PCNONE,ierr) - - ! set matrix free finite difference - call MatCreateSNESMF(snes,jac,ierr) - call SNESSetJacobian(snes,jac,jac,MatMFFDComputeJacobian,PETSC_NULL,ierr) - - ! set SNES options - call SNESSetFromOptions(snes,ierr) - - ! solve nonlinear system - call SNESSolve(snes,PETSC_NULL,x,ierr) - - ! output answers - call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'nonlinear.bin',FILE_MODE_WRITE, & - viewer,ierr) - call VecView(x,viewer,ierr) - call PetscViewerDestroy(viewer,ierr) - - end subroutine solve_cmfd_nonlinear - -!=============================================================================== -! COMPUTE_NONLINEAR_RESIDUAL -!=============================================================================== - - subroutine compute_nonlinear_residual(snes,x,res,ierr) - -#include - - ! arguments - SNES :: snes ! nonlinear solver context - Vec :: x ! independent vector - Vec :: res ! residual vector - integer :: ierr ! error flag - - Mat :: M ! loss matrix - Mat :: F ! production matrix - Vec :: phi ! flux vector - Vec :: rphi ! flux part of residual - Vec :: phiM ! M part of residual flux calc - integer :: n ! dimensions of matrix - integer :: nx ! maximum number of x cells - integer :: ny ! maximum number of y cells - integer :: nz ! maximum number of z cells - integer :: ng ! maximum number of groups - integer :: nzM ! max number of nonzeros in a row for M - integer :: nzF ! max number of nonzeros in a row for F - 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 - - ! get maximum number of cells in each direction - nx = cmfd%indices(1) - ny = cmfd%indices(2) - nz = cmfd%indices(3) - ng = cmfd%indices(4) - - ! calculate dimensions of matrix - if (allocated(cmfd % coremap)) then - n = cmfd % mat_dim * ng - else - n = nx*ny*nz*ng - end if - - ! get pointers to vectors - call VecGetArrayF90(x,xptr,ierr) - call VecGetArrayF90(res,rptr,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 VecCreate(PETSC_COMM_SELF,rphi,ierr) - call VecSetSizes(rphi,PETSC_DECIDE,n,ierr) - call VecSetFromOptions(rphi,ierr) - - ! extract flux and place in petsc vector - call VecPlaceArray(phi,xptr,ierr) - call VecPlaceArray(rphi,rptr,ierr) - - ! extract eigenvalue - lambda = xptr(n+1) - - ! maximum number of nonzeros in each matrix - nzM = 7+ng-1 - nzF = ng - - ! set up loss matrix - call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,nzM,PETSC_NULL_INTEGER,M,ierr) - call MatSetOption(M,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE,ierr) - call MatSetOption(M,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr) - call MatSetOption(M,MAT_USE_HASH_TABLE,PETSC_TRUE,ierr) - - ! set up production matrix - call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,nzF,PETSC_NULL_INTEGER,F,ierr) - call MatSetOption(F,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE,ierr) - call MatSetOption(F,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr) - call MatSetOption(F,MAT_USE_HASH_TABLE,PETSC_TRUE,ierr) - - ! create operators - call loss_matrix(M) - call prod_matrix(F) - - ! create new petsc vectors to perform math - call VecCreate(PETSC_COMM_SELF,phiM,ierr) - call VecSetSizes(phiM,PETSC_DECIDE,n,ierr) - call VecSetFromOptions(phiM,ierr) - - ! calculate flux part of residual vector - call MatMult(M,phi,phiM,ierr) - call MatMult(F,phi,rphi,ierr) - call VecAYPX(rphi,-1.0_8*lambda,phiM,ierr) - - ! set eigenvalue part of residual vector - call VecDot(phi,phi,reslamb,ierr) - - ! map to ptr - rptr(n+1) = 0.5_8 - 0.5_8*reslamb - - ! reset arrays that are not used - call VecResetArray(phi,ierr) - call VecResetArray(rphi,ierr) - - ! restore arrays for residual and solution - call VecRestoreArrayF90(x,xptr,ierr) - call VecRestoreArrayF90(res,rptr,ierr) - - ! destroy all temp vectors - call VecDestroy(phi,ierr) - call VecDestroy(phiM,ierr) - call VecDestroy(rphi,ierr) - - end subroutine compute_nonlinear_residual - -!=============================================================================== -! DESTROY_DATA deletes all PETSc and SLEPc objects -!=============================================================================== - - subroutine destroy_data(M,F,phi,eps,viewer) - -#include -#include -#include - - Mat :: M ! loss matrix - Mat :: F ! production matrix - Vec :: phi ! eigenvector - EPS :: eps ! slepc eigenvalue object - PetscViewer :: viewer ! viewer for answer - - integer :: ierr ! petsc error code - - ! destroy all - call PetscViewerDestroy(viewer,ierr) - call EPSDestroy(eps,ierr) - call VecDestroy(phi,ierr) - call MatDestroy(M,ierr) - call MatDestroy(F,ierr) - - end subroutine destroy_data - -!=============================================================================== -! LOSS_MATRIX creates the matrix representing loss of neutrons -!=============================================================================== - - subroutine loss_matrix(M) - -#include - - ! arguments - Mat :: M ! loss matrix - - ! local variables - integer :: nx ! maximum number of cells in x direction - integer :: ny ! maximum number of cells in y direction - integer :: nz ! maximum number of cells in z direction - integer :: ng ! maximum number of energy groups - integer :: nxyz(3,2) ! single vector containing bound. locations - 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 :: 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 - integer :: xyz_idx ! index for determining if x,y or z leakage - integer :: dir_idx ! index for determining - or + face of cell - integer :: neig_idx(3) ! spatial indices of neighbour - integer :: shift_idx ! parameter to shift index by +1 or -1 - integer :: ierr ! Persc error code - integer :: kount ! integer for counting values in vector - 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 - real(8) :: dtilde(6) ! finite difference coupling parameter - real(8) :: dhat(6) ! nonlinear coupling parameter - real(8) :: hxyz(3) ! cell lengths in each direction - real(8) :: jn ! direction dependent leakage coeff to neig - real(8) :: jo(6) ! leakage coeff in front of cell flux - real(8) :: jnet ! net leakage from jo - real(8) :: val ! temporary variable before saving to matrix - PetscViewer :: viewer ! viewer to write out matrix to binary file - - ! initialize matrix for building - call MatAssemblyBegin(M,MAT_FLUSH_ASSEMBLY,ierr) - - ! get maximum indices - nx = cmfd%indices(1) - ny = cmfd%indices(2) - nz = cmfd%indices(3) - ng = cmfd%indices(4) - - ! create single vector of these indices for boundary calculation - nxyz(1,:) = (/1,nx/) - nxyz(2,:) = (/1,ny/) - nxyz(3,:) = (/1,nz/) - - ! begin iteration loops - ZLOOP: do k = 1,nz - - YLOOP: do j = 1,ny - - XLOOP: do i = 1,nx - - ! check if not including reflector - if (allocated(cmfd % coremap)) then - - ! check if at a reflector - if (cmfd % coremap(i,j,k) == 99999) then - cycle - end if - - end if - - GROUP: do g = 1,ng - - ! get matrix index of cell - cell_mat_idx = get_matrix_idx(g,i,j,k,ng,nx,ny) - - ! 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 and get dhat - if (allocated(cmfd%dhat)) then - dhat = cmfd%dhat(:,g,i,j,k) - else - dhat = 0.0 - 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 (allocated(cmfd % coremap)) then - - ! 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 - neig_mat_idx = get_matrix_idx(g,neig_idx(1),neig_idx(2), & - & neig_idx(3),ng,nx,ny) - - ! compute value and record to bank - val = jn/hxyz(xyz_idx) - - ! record value in matrix - call MatSetValue(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),ng,nx,ny) - - ! compute value and record to bank - val = jn/hxyz(xyz_idx) - - ! record value in matrix - call MatSetValue(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(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,ng,nx,ny) - - ! get scattering macro xs - scattxshg = cmfd%scattxs(h,g,i,j,k) - - ! record value in matrix (negate it) - val = -scattxshg - - call MatSetValue(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 - call MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY,ierr) - - ! write out matrix in binary file (debugging) - call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'lossmat.bin',FILE_MODE_WRITE, & - viewer,ierr) - call MatView(M,viewer,ierr) - call PetscViewerDestroy(viewer,ierr) - - end subroutine loss_matrix - -!=============================================================================== -! PROD_MATRIX creates the matrix representing production of neutrons -!=============================================================================== - - subroutine prod_matrix(F) - - use cmfd_utils, only: get_matrix_idx - -#include - - ! arguments - Mat :: F ! production matrix - - ! local variables - integer :: nx ! maximum number of cells in x direction - integer :: ny ! maximum number of cells in y direction - integer :: nz ! maximum number of cells in z direction - integer :: ng ! maximum number of energy groups - 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 :: l ! iteration counter for leakages - integer :: h ! energy group when doing scattering - integer :: gmat_idx ! index in matrix for energy group g - integer :: hmat_idx ! index in matrix for energy group h - integer :: ierr ! Petsc error code - real(8) :: nfissxs ! nufission cross section h-->g - real(8) :: val ! temporary variable for nfissxs - PetscViewer :: viewer ! viewer to print out matrix - - ! initialize matrix for building - call MatAssemblyBegin(F,MAT_FLUSH_ASSEMBLY,ierr) - - ! get maximum indices - nx = cmfd%indices(1) - ny = cmfd%indices(2) - nz = cmfd%indices(3) - ng = cmfd%indices(4) - - ! begin loop around energy groups and spatial indices - ZLOOP: do k = 1,nz - - YLOOP: do j = 1,ny - - XLOOP: do i = 1,nx - - ! check if not including reflector - if (allocated(cmfd % coremap)) then - - ! check if at a reflector - if (cmfd % coremap(i,j,k) == 99999) then - cycle - end if - - end if - - GROUP: do g = 1,ng - - NFISS: do h = 1,ng - - ! get cell data - nfissxs = cmfd%nfissxs(h,g,i,j,k) - - ! get matrix location - gmat_idx = get_matrix_idx(g,i,j,k,ng,nx,ny) - hmat_idx = get_matrix_idx(h,i,j,k,ng,nx,ny) - - ! reocrd value in matrix - val = nfissxs - call MatSetValue(F,gmat_idx-1,hmat_idx-1,val,INSERT_VALUES,ierr) - - end do NFISS - - end do GROUP - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - ! finalize matrix assembly - call MatAssemblyEnd(F,MAT_FINAL_ASSEMBLY,ierr) - - ! write out matrix in binary file (debugging) - call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'prodmat.bin',FILE_MODE_WRITE, & - viewer,ierr) - call MatView(F,viewer,ierr) - call PetscViewerDestroy(viewer,ierr) - - end subroutine prod_matrix - -!=============================================================================== -! SOURCE_PDF calculates the probability distribution of the cmfd fission source -!=============================================================================== - - subroutine source_pdf(source) - - use cmfd_utils, only: get_matrix_idx - -#include - - ! arguments - Vec :: source ! new source vector - - ! local variables - integer :: nx ! maximum number of cells in x direction - integer :: ny ! maximum number of cells in y direction - integer :: nz ! maximum number of cells in z direction - integer :: ng ! maximum number of energy groups - 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 :: ierr ! PETSC error code - integer :: idx ! index in vector - real(8) :: total ! sum of source vector - real(8) :: hxyz(3) ! cell dimensions of current ijk cell - PetscScalar, pointer :: source_ptr(:) ! pointer to petsc vector for fortran - - ! get maximum of spatial and group indices - nx = cmfd%indices(1) - ny = cmfd%indices(2) - nz = cmfd%indices(3) - ng = cmfd%indices(4) - - ! get source vector from petsc - call VecGetArrayF90(source,source_ptr,ierr) - - ! loop around indices to map to cmfd object - ZLOOP: do k = 1,nz - - YLOOP: do j = 1,ny - - XLOOP: do i = 1,nx - - GROUP: do g = 1,ng - - ! get dimensions of cell - hxyz = cmfd%hxyz(:,i,j,k) - - ! get index - idx = get_matrix_idx(g,i,j,k,ng,nx,ny) - - ! multiply source density by volume and record in object - cmfd%sourcepdf(g,i,j,k) = source_ptr(idx)*hxyz(1)*hxyz(2)*hxyz(3) - - end do GROUP - - end do XLOOP - - end do YLOOP - - end do ZLOOP - - ! normalize source such that it sums to 1.0 - cmfd%sourcepdf = cmfd%sourcepdf/sum(cmfd%sourcepdf) - - ! restore petsc vector - call VecRestoreArrayF90(source,source_ptr,ierr) - - end subroutine source_pdf - -!=============================================================================== -! COUNT_SOURCE determines the number of source sites in each mesh box -!=============================================================================== - - subroutine count_source() - - end subroutine count_source - -!=============================================================================== -! WEIGHT_FACTORS calculates the weight adjustment factors for next MC cycle -!=============================================================================== - - subroutine weight_factors() - - end subroutine weight_factors - -!=============================================================================== -! ADJUST_WEIGHT adjusts the initial weight of the particle -!=============================================================================== - - subroutine adjust_weight() - -! should this do all source particles at once or called each time before a -! neutron is born in the MC cycle? MC21 is the latter - - end subroutine adjust_weight - -!=============================================================================== -! CMFD_SOLVER in the main power iteration routine for the cmfd calculation -!=============================================================================== - - subroutine cmfd_solver_power() - -use timing, only: timer_start, timer_stop - -#include - - Mat :: M ! loss matrix - Mat :: F ! production matrix - Vec :: phi_n ! new flux eigenvector - Vec :: phi_o ! old flux eigenvector - Vec :: S_n ! new source vector - Vec :: S_o ! old source vector - real(8) :: k_n ! new k-eigenvalue - real(8) :: k_o ! old k-eigenvlaue - real(8) :: num ! numerator for eigenvalue update - real(8) :: den ! denominator for eigenvalue update - real(8) :: one=1.0 ! one - integer :: ierr ! error flag - KSP :: krylov ! krylov solver - PC :: prec ! preconditioner for krylov - PetscViewer :: viewer ! viewer for answer - real(8) :: info(MAT_INFO_SIZE) - real(8) :: mall - real(8) :: nza,nzu,nzun - integer :: i ! iteration counter - logical :: iconv ! is problem converged - - ! reset convergence flag - iconv = .FALSE. - - ! initialize PETSc - call PetscInitialize(PETSC_NULL_CHARACTER,ierr) - - ! initialize matrices and vectors - print *,"Initializing and building matrices" - call timer_start(time_mat) - call init_data_power(M,F,phi_n,phi_o,S_n,S_o,k_n,k_o,krylov,prec) - - ! set up M loss matrix - call loss_matrix(M) -! call MatGetInfo(M,MAT_LOCAL,info,ierr) -! mall = info(MAT_INFO_MEMORY) -! nza = info(MAT_INFO_NZ_ALLOCATED) -! nzu = info(MAT_INFO_NZ_USED) -! nzun = info(MAT_INFO_NZ_UNNEEDED) - - ! set up F production matrix - call prod_matrix(F) - call timer_stop(time_mat) - - ! set up krylov info - call KSPSetOperators(krylov, M, M, SAME_NONZERO_PATTERN, ierr) - call KSPSetUp(krylov,ierr) - - ! calculate preconditioner (ILU) - call PCFactorGetMatrix(prec,M,ierr) - - ! begin timer for power iteration - print *,"Beginning power iteration" - call timer_start(time_eigen) - - ! begin power iteration - do i = 1,10000 - - ! compute source vector - call MatMult(F,phi_o,S_o,ierr) - - ! normalize source vector - call VecScale(S_o,one/k_o,ierr) - - ! compute new flux vector - call KSPSolve(krylov,S_o,phi_n,ierr) - - ! compute new source vector - call MatMult(F,phi_n,S_n,ierr) - - ! compute new k-eigenvalue - call VecSum(S_n,num,ierr) - call VecSum(S_o,den,ierr) - k_n = num/den - - ! renormalize the old source - call VecScale(S_o,k_o,ierr) - - ! check convergence - call convergence(S_n,S_o,k_o,k_n,iconv) - - ! to break or not to break - if (iconv) exit - - ! record old values - call VecCopy(phi_n,phi_o,ierr) - k_o = k_n - - end do - - ! print out keff - print *,'k-effective:',k_n - - ! end power iteration timer - call timer_stop(time_eigen) - print *,"Matrix building time (s):",time_mat%elapsed - print *,"Power iteration time (s):",time_eigen%elapsed - print *,"Power iteration time per iteration (s):",time_eigen%elapsed/i - print *,"Number of Power iterations:",i - - ! compute source pdf and record in cmfd object -! call source_pdf(S_n) - - ! output answers - call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'fluxvec.bin',FILE_MODE_WRITE, & - viewer,ierr) - call VecView(phi_n,viewer,ierr) - call PetscViewerDestroy(viewer,ierr) - - ! finalize PETSc - call PetscFinalize(ierr) - - end subroutine cmfd_solver_power - -!=============================================================================== -! INIT_DATA allocates matrices vectors for CMFD solution -!=============================================================================== - - subroutine init_data_power(M,F,phi_n,phi_o,S_n,S_o,k_n,k_o,krylov,prec) - -#include - - ! arguments - Mat :: M ! loss matrix - Mat :: F ! production matrix - Vec :: phi_n ! new flux eigenvector - Vec :: phi_o ! old flux eigenvector - Vec :: S_n ! new source vector - Vec :: S_o ! old source vector - real(8) :: k_n ! new k-eigenvalue - real(8) :: k_o ! old k-eigenvalue - KSP :: krylov ! krylov solver - PC :: prec ! preconditioner for krylov - - ! local variables - integer :: n ! dimensions of matrix - integer :: nz_M ! number of nonzeros in loss matrix - integer :: nz_F ! number of nonzeros in prod matrix - integer :: ierr ! error flag - integer :: nx ! maximum number of x cells - integer :: ny ! maximum number of y cells - integer :: nz ! maximum number of z cells - integer :: ng ! maximum number of groups - integer :: nzM ! max number of nonzeros in a row for M - integer :: nzF ! max number of nonzeros in a row for F - real(8) :: guess=1.0 ! initial guess - real(8) :: ktol=1.e-7 ! krylov tolerance - real(8) :: mem - - - ! get maximum number of cells in each direction - nx = cmfd%indices(1) - ny = cmfd%indices(2) - nz = cmfd%indices(3) - ng = cmfd%indices(4) - - ! calculate dimensions of matrix - if (allocated(cmfd % coremap)) then - n = cmfd % mat_dim * ng - else - n = nx*ny*nz*ng - end if - - ! maximum number of nonzeros in each matrix - nzM = 7+ng-1 - nzF = ng - - ! set up loss matrix - call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,nzM,PETSC_NULL_INTEGER,M,ierr) - call MatSetOption(M,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE,ierr) - call MatSetOption(M,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr) - call MatSetOption(M,MAT_USE_HASH_TABLE,PETSC_TRUE,ierr) - - ! set up production matrix - call MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,nzF,PETSC_NULL_INTEGER,F,ierr) - call MatSetOption(F,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE,ierr) - call MatSetOption(F,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE,ierr) - call MatSetOption(F,MAT_USE_HASH_TABLE,PETSC_TRUE,ierr) - - - ! set up flux vectors - call VecCreate(PETSC_COMM_WORLD,phi_n,ierr) - call VecSetSizes(phi_n,PETSC_DECIDE,n,ierr) - call VecSetFromOptions(phi_n,ierr) - call VecCreate(PETSC_COMM_WORLD,phi_o,ierr) - call VecSetSizes(phi_o,PETSC_DECIDE,n,ierr) - call VecSetFromOptions(phi_o,ierr) - - - ! set up source vectors - call VecCreate(PETSC_COMM_WORLD,S_n,ierr) - call VecSetSizes(S_n,PETSC_DECIDE,n,ierr) - call VecSetFromOptions(S_n,ierr) - call VecCreate(PETSC_COMM_WORLD,S_o,ierr) - call VecSetSizes(S_o,PETSC_DECIDE,n,ierr) - call VecSetFromOptions(S_o,ierr) - - ! set initial guess - call VecSet(phi_n,guess,ierr) - call VecSet(phi_o,guess,ierr) - k_n = guess - k_o = guess - - ! set up krylov solver - call KSPCreate(PETSC_COMM_WORLD,krylov,ierr) - call KSPSetTolerances(krylov,ktol,PETSC_DEFAULT_DOUBLE_PRECISION, & - & PETSC_DEFAULT_DOUBLE_PRECISION, & - & PETSC_DEFAULT_INTEGER,ierr) - call KSPSetType(krylov,KSPGMRES,ierr) - call KSPSetInitialGuessNonzero(krylov,PETSC_TRUE,ierr) - call KSPSetInitialGuessNonzero(krylov,PETSC_TRUE,ierr) - call KSPGetPC(krylov,prec,ierr) - call PCSetType(prec,PCILU,ierr) - call PCFactorSetLevels(prec,5,ierr) - call KSPSetFromOptions(krylov,ierr) - - end subroutine init_data_power - -!=============================================================================== -! CONVERGENCE checks the convergence of eigenvalue, eigenvector and source -!=============================================================================== - - subroutine convergence(S_n,S_o,k_o,k_n,iconv) - -#include - - ! arguments - Vec :: phi_n ! new flux eigenvector - Vec :: phi_o ! old flux eigenvector - Vec :: S_n ! new source vector - Vec :: S_o ! old source vector - real(8) :: k_n ! new k-eigenvalue - real(8) :: k_o ! old k-eigenvalue - logical :: iconv ! is the problem converged - - ! local variables - real(8) :: ktol = 1.e-6 ! tolerance on keff - real(8) :: stol = 1.e-5 ! tolerance on source - real(8) :: kerr ! error in keff - real(8) :: serr ! error in source - real(8) :: one = -1.0 ! one - real(8) :: norm_n ! L2 norm of new source - real(8) :: norm_o ! L2 norm of old source - integer :: floc ! location of max error in flux - integer :: sloc ! location of max error in source - integer :: ierr ! petsc error code - integer :: n ! vector size - - ! reset convergence flag - iconv = .FALSE. - - ! calculate error in keff - kerr = abs(k_o - k_n)/k_n - - ! calculate max error in source - call VecNorm(S_n,NORM_2,norm_n,ierr) - call VecNorm(S_o,NORM_2,norm_o,ierr) - serr = abs(norm_n-norm_o)/norm_n - - ! check for convergence - if(kerr < ktol .and. serr < stol) iconv = .TRUE. - - print *,k_n,kerr,serr - - end subroutine convergence + end subroutine execute_cmfd end module cmfd_execute