From dd209cd1c47112d6a6b977cd0b2bbd98d4afc65c Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Fri, 27 Jan 2012 20:28:28 -0500 Subject: [PATCH] added separate modules for solver routines --- src/cmfd_power_solver.F90 | 278 ++++++++++++++++++++++++++++++++++++++ src/cmfd_slepc_solver.F90 | 175 ++++++++++++++++++++++++ src/cmfd_snes_solver.F90 | 251 ++++++++++++++++++++++++++++++++++ 3 files changed, 704 insertions(+) create mode 100644 src/cmfd_power_solver.F90 create mode 100644 src/cmfd_slepc_solver.F90 create mode 100644 src/cmfd_snes_solver.F90 diff --git a/src/cmfd_power_solver.F90 b/src/cmfd_power_solver.F90 new file mode 100644 index 0000000000..f50a6e9733 --- /dev/null +++ b/src/cmfd_power_solver.F90 @@ -0,0 +1,278 @@ +module cmfd_power_solver + + use cmfd_loss_operator, only: loss_operator,init_M_operator, & + & build_loss_matrix,print_M_operator,destroy_M_operator + use cmfd_prod_operator, only: prod_operator,init_F_operator, & + & build_prod_matrix,print_F_operator,destroy_F_operator + + + implicit none + private + public :: cmfd_power_execute + +#include + + type(loss_operator) :: loss + type(prod_operator) :: prod + + Vec :: phi_n ! new flux eigenvector + Vec :: phi_o ! old flux eigenvector + Vec :: S_n ! new source vector + Vec :: S_o ! old source vector + KSP :: krylov ! krylov solver + PC :: prec ! preconditioner for krylov + integer :: ierr ! error flag + real(8) :: k_n ! new k-eigenvalue + real(8) :: k_o ! old k-eigenvalue + logical :: iconv ! did the problem converged + +contains + +!=============================================================================== +! CMFD_POWER_EXECUTE +!=============================================================================== + + subroutine cmfd_power_execute() + + ! initialize solver + call init_solver() + + ! initialize matrices and vectors + call init_data() + + ! set up M loss matrix + call build_loss_matrix(loss) + + ! set up F production matrix + call build_prod_matrix(prod) + + ! set up krylov info + call KSPSetOperators(krylov, loss%M, loss%M, SAME_NONZERO_PATTERN, ierr) + call KSPSetUp(krylov,ierr) + + ! calculate preconditioner (ILU) + call PCFactorGetMatrix(prec,loss%M,ierr) + + ! begin power iteration + call execute_power_iter() + + ! extract results + call extract_results() + + ! deallocate petsc objects + call finalize() + + end subroutine cmfd_power_execute + +!=============================================================================== +! INIT_DATA allocates matrices vectors for CMFD solution +!=============================================================================== + + subroutine init_data() + + integer :: n ! problem size + real(8) :: guess ! initial guess + + ! set up matrices + call init_M_operator(loss) + call init_F_operator(prod) + + ! get problem size + n = loss%n + + ! 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 + guess = 1.0_8 + call VecSet(phi_n,guess,ierr) + call VecSet(phi_o,guess,ierr) + k_n = guess + k_o = guess + + end subroutine init_data + +!=============================================================================== +! INIT_SOLVER +!=============================================================================== + + subroutine init_solver() + + real(8) :: solvertol ! krylov tolerance + + ! set tolerance + solvertol = 1.0e-7_8 + + ! set up krylov solver + call KSPCreate(PETSC_COMM_WORLD,krylov,ierr) + call KSPSetTolerances(krylov,solvertol,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_solver + +!=============================================================================== +! EXECUTE_POWER_ITER in the main power iteration routine +! for the cmfd calculation +!=============================================================================== + + subroutine execute_power_iter() + + real(8) :: num ! numerator for eigenvalue update + real(8) :: den ! denominator for eigenvalue update + real(8) :: one=1.0_8 ! one + integer :: i ! iteration counter + + ! reset convergence flag + iconv = .FALSE. + + ! begin power iteration + do i = 1,10000 + + ! compute source vector + call MatMult(prod%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(prod%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() + + ! 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 + + end subroutine execute_power_iter + +!=============================================================================== +! CONVERGENCE checks the convergence of eigenvalue, eigenvector and source +!=============================================================================== + + subroutine convergence() + + real(8) :: ktol = 1.e-6_8 ! tolerance on keff + real(8) :: stol = 1.e-5_8 ! tolerance on source + real(8) :: kerr ! error in keff + real(8) :: serr ! error in source + real(8) :: one = -1.0_8 ! 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 out to user (TODO: make formatted) + print *,k_n,kerr,serr + + end subroutine convergence + +!============================================================================== +! EXTRACT_RESULTS +!============================================================================== + + subroutine extract_results() + + use global, only: cmfd + + PetscViewer :: viewer + PetscScalar, pointer :: phi_v(:) ! pointer to eigenvector info + integer :: n ! problem size + + ! get problem size + n = loss%n + + ! also allocate in cmfd object + if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n)) + + ! convert petsc phi_object to cmfd_obj + call VecGetArrayF90(phi_n,phi_v,ierr) + cmfd%phi = phi_v + + call VecRestoreArrayF90(phi_n,phi_v,ierr) + + ! save eigenvalue + cmfd%keff = k_n + + ! write out results + call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'fluxvec.bin',FILE_MODE_WRITE, & + & viewer,ierr) + call VecView(phi_n,viewer,ierr) + call PetscViewerDestroy(viewer,ierr) + + end subroutine extract_results + +!============================================================================== +! FINALIZE +!============================================================================== + + subroutine finalize() + + ! finalize data objects + call destroy_M_operator(loss) + call destroy_F_operator(prod) + + call VecDestroy(phi_n,ierr) + call VecDestroy(phi_o,ierr) + call VecDestroy(S_n,ierr) + call VecDestroy(S_o,ierr) + + ! finalize solver objects +! call KSPDestroy(krylov,ierr) + + end subroutine finalize + +end module cmfd_power_solver diff --git a/src/cmfd_slepc_solver.F90 b/src/cmfd_slepc_solver.F90 new file mode 100644 index 0000000000..2d4bf85c69 --- /dev/null +++ b/src/cmfd_slepc_solver.F90 @@ -0,0 +1,175 @@ +module cmfd_slepc_solver + + use cmfd_loss_operator, only: loss_operator,init_M_operator, & + & build_loss_matrix,print_M_operator,destroy_M_operator + use cmfd_prod_operator, only: prod_operator,init_F_operator, & + & build_prod_matrix,print_F_operator,destroy_F_operator + + implicit none + +#include +#include +#include + + type(loss_operator) :: loss + type(prod_operator) :: prod + + Vec :: phi ! eigenvector + EPS :: eps ! slepc eigenvalue object + ST :: st ! slepc spectral trans object + KSP :: ksp ! linear solver object + PC :: pc ! preconditioner object + integer :: ierr ! error flag + real(8) :: keff ! the converged eigenvalue + +contains + +!=============================================================================== +! CMFD_SLEPC_EXECUTE +!=============================================================================== + + subroutine cmfd_slepc_execute() + + ! initialize data + call init_data() + + ! initialize solver + call init_solver() + + ! build operators + call build_loss_matrix(loss) + call build_prod_matrix(prod) + + ! set operators to EPS object + call EPSSetOperators(eps,prod%F,loss%M,ierr) + + ! solve the system + call EPSSolve(eps,ierr) + + ! extracts results to cmfd object + call extract_results() + + ! deallocate all slepc data + call finalize() + + end subroutine cmfd_slepc_execute + +!=============================================================================== +! INIT_DATA allocates matrices vectors for CMFD solution +!=============================================================================== + + subroutine init_data() + + integer :: n ! problem size + + ! set up matrices + call init_M_operator(loss) + call init_F_operator(prod) + + ! get problem size + n = loss%n + + ! set up eigenvector + call VecCreate(PETSC_COMM_SELF,phi,ierr) + call VecSetSizes(phi,PETSC_DECIDE,n,ierr) + call VecSetFromOptions(phi,ierr) + + end subroutine init_data + +!=============================================================================== +! INIT_SOLVER +!=============================================================================== + + subroutine init_solver() + + 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) + + ! display information to user + write(*,*) 'EPS TYPE IS: ',epstype + write(*,*) 'ST TYPE IS: ',sttype + write(*,*) 'KSP TYPE IS: ',ksptype + write(*,*) 'PC TYPE IS: ',pctype + + end subroutine init_solver + +!============================================================================== +! EXTRACT_RESULTS +!============================================================================== + + subroutine extract_results() + + use global, only: cmfd + + integer :: n ! problem size + integer :: i_eig = 0 ! eigenvalue to extract + PetscViewer :: viewer ! petsc output object + PetscScalar, pointer :: phi_v(:) ! pointer to eigenvector info + + ! get problem size + n = loss%n + + ! also allocate in cmfd object + if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n)) + + ! extract run information + call EPSGetEigenpair(eps,i_eig,keff,PETSC_NULL,phi,PETSC_NULL_OBJECT,ierr) + + ! convert petsc phi_object to cmfd_obj + call VecGetArrayF90(phi,phi_v,ierr) + cmfd%phi = phi_v + call VecRestoreArrayF90(phi,phi_v,ierr) + + ! save eigenvalue + cmfd%keff = keff + + ! write out results + call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'fluxvec.bin',FILE_MODE_WRITE, & + & viewer,ierr) + call VecView(phi,viewer,ierr) + call PetscViewerDestroy(viewer,ierr) + + end subroutine extract_results + +!============================================================================== +! FINALIZE +!============================================================================== + + subroutine finalize() + + ! finalize data objects + call destroy_M_operator(loss) + call destroy_F_operator(prod) + call VecDestroy(phi,ierr) + + ! finalize solver objects + call EPSDestroy(eps,ierr) + + end subroutine finalize + +end module cmfd_slepc_solver diff --git a/src/cmfd_snes_solver.F90 b/src/cmfd_snes_solver.F90 new file mode 100644 index 0000000000..4a8d359044 --- /dev/null +++ b/src/cmfd_snes_solver.F90 @@ -0,0 +1,251 @@ +module cmfd_snes_solver + + use cmfd_loss_operator, only: loss_operator,init_M_operator, & + & build_loss_matrix,print_M_operator,destroy_M_operator + use cmfd_prod_operator, only: prod_operator,init_F_operator, & + & build_prod_matrix,print_F_operator,destroy_F_operator + use cmfd_slepc_solver, only: cmfd_slepc_execute + +implicit none + +#include + + + type(loss_operator) :: loss + type(prod_operator) :: prod + + Mat :: jac ! jacobian matrix + Vec :: resvec ! residual vector + Vec :: xvec ! results + KSP :: ksp ! linear solver context + PC :: pc ! preconditioner + SNES :: snes ! nonlinear solver context + integer :: ierr ! error flag + +contains + +!=============================================================================== +! CMFD_SNES_EXECUTE +!=============================================================================== + + subroutine cmfd_snes_execute() + + ! call slepc solver + call cmfd_slepc_execute() + + ! initialize data + call init_data() + + ! initialize solver + call init_solver() + + ! solve the system + call SNESSolve(snes,PETSC_NULL,xvec,ierr) + + ! extracts results to cmfd object + call extract_results() + + ! deallocate all slepc data + call finalize() + + end subroutine cmfd_snes_execute + +!=============================================================================== +! INIT_DATA allocates matrices vectors for CMFD solution +!=============================================================================== + + subroutine init_data() + + use global, only: cmfd + + integer :: k ! implied do counter + integer :: n ! problem size + real(8), pointer :: xptr(:) ! solution pointer + + ! set up matrices + call init_M_operator(loss) + call init_F_operator(prod) + + ! get problem size + n = loss%n + + ! create PETSc vectors + call VecCreate(PETSC_COMM_SELF,resvec,ierr) + call VecSetSizes(resvec,PETSC_DECIDE,n+1,ierr) + call VecSetFromOptions(resvec,ierr) + call VecCreate(PETSC_COMM_SELF,xvec,ierr) + call VecSetSizes(xvec,PETSC_DECIDE,n+1,ierr) + call VecSetFromOptions(xvec,ierr) + + ! set flux in guess + call VecSetValues(xvec,n,(/(k,k=0,n-1)/),cmfd%phi,INSERT_VALUES,ierr) + call VecAssemblyBegin(xvec,ierr) + call VecAssemblyEnd(xvec,ierr) + + ! set keff in guess + call VecGetArrayF90(xvec,xptr,ierr) + xptr(n+1) = 1.0_8/cmfd%keff + call VecRestoreArrayF90(xvec,xptr,ierr) + + end subroutine init_data + +!=============================================================================== +! INIT_SOLVER +!=============================================================================== + + subroutine init_solver() + + ! create SNES context + call SNESCreate(PETSC_COMM_SELF,snes,ierr) + + ! set the residual function + call SNESSetFunction(snes,resvec,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) + + end subroutine init_solver + +!=============================================================================== +! COMPUTE_NONLINEAR_RESIDUAL +!=============================================================================== + + subroutine compute_nonlinear_residual(snes,x,res,ierr) + + ! arguments + SNES :: snes ! nonlinear solver context + Vec :: x ! independent vector + Vec :: res ! residual vector + integer :: ierr ! error flag + + Vec :: phi ! flux vector + Vec :: rphi ! flux part of residual + Vec :: phiM ! M part of residual flux calc + integer :: n ! problem size + 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 problem size + n = loss%n + + ! 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) + + ! create operators + call build_loss_matrix(loss) + call build_prod_matrix(prod) + + ! 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(loss%M,phi,phiM,ierr) + call MatMult(prod%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 + +!============================================================================== +! EXTRACT_RESULTS +!============================================================================== + + subroutine extract_results() + + use global, only: cmfd + + integer :: n ! problem size + PetscViewer :: viewer + PetscScalar, pointer :: xptr(:) ! pointer to eigenvector info + + ! get problem size + n = loss%n + + ! also allocate in cmfd object + if (.not. allocated(cmfd%phi)) allocate(cmfd%phi(n)) + + ! convert petsc phi_object to cmfd_obj + call VecGetArrayF90(xvec,xptr,ierr) + cmfd%phi = xptr(1:n) + + ! save eigenvalue + cmfd%keff = 1.0_8 / xptr(n+1) + call VecRestoreArrayF90(xvec,xptr,ierr) + + ! write out results + call PetscViewerBinaryOpen(PETSC_COMM_WORLD,'residual.bin',FILE_MODE_WRITE, & + viewer,ierr) + call VecView(resvec,viewer,ierr) + call PetscViewerDestroy(viewer,ierr) + + end subroutine extract_results + +!============================================================================== +! FINALIZE +!============================================================================== + + subroutine finalize() + + ! finalize data objects + call destroy_M_operator(loss) + call destroy_F_operator(prod) + call VecDestroy(xvec,ierr) + call VecDestroy(resvec,ierr) + + ! finalize solver objects + call SNESDestroy(snes,ierr) + + end subroutine finalize + +end module cmfd_snes_solver