mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
added separate modules for solver routines
This commit is contained in:
Bryan Herman
parent
5dde575415
commit
dd209cd1c4
3 changed files with 704 additions and 0 deletions
278
src/cmfd_power_solver.F90
Normal file
278
src/cmfd_power_solver.F90
Normal file
|
|
@ -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 <finclude/petsc.h90>
|
||||
|
||||
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
|
||||
175
src/cmfd_slepc_solver.F90
Normal file
175
src/cmfd_slepc_solver.F90
Normal file
|
|
@ -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 <finclude/petsc.h90>
|
||||
#include <finclude/slepcsys.h>
|
||||
#include <finclude/slepceps.h>
|
||||
|
||||
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
|
||||
251
src/cmfd_snes_solver.F90
Normal file
251
src/cmfd_snes_solver.F90
Normal file
|
|
@ -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 <finclude/petsc.h90>
|
||||
|
||||
|
||||
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
|
||||
Loading…
Add table
Add a link
Reference in a new issue