removed built-in GMRES solver only PETSc available again

This commit is contained in:
Bryan Herman 2013-08-27 16:11:27 -04:00
parent 4b6b1a7b98
commit e2d82f35a2
5 changed files with 19 additions and 1426 deletions

View file

@ -327,7 +327,6 @@ set_header.o: list_header.o
solver_interface.o: error.o
solver_interface.o: global.o
solver_interface.o: gmres_solver.o
solver_interface.o: matrix_header.o
solver_interface.o: vector_header.o

View file

@ -46,19 +46,25 @@ contains
call init_data()
! Initialize solver
#ifdef PETSC
call jfnk % create()
#endif
! Set up residual and jacobian routines
jfnk_data % res_proc_ptr => compute_nonlinear_residual
jfnk_data % jac_proc_ptr => build_jacobian_matrix
#ifdef PETSC
call jfnk % set_functions(jfnk_data, resvec, jac_prec)
#endif
! Stop timer for build
call time_cmfdbuild % stop()
! Solve the system
call time_cmfdsolve % start()
#ifdef PETSC
call jfnk % solve(xvec)
#endif
call time_cmfdsolve % stop()
! Extracts results to cmfd object
@ -378,7 +384,9 @@ contains
call jac_prec % destroy()
call xvec % destroy()
call resvec % destroy()
call jfnk % destroy()
#ifdef PETSC
call jfnk % destroy()
#endif
nullify(jfnk_data % res_proc_ptr)
nullify(jfnk_data % jac_proc_ptr)

View file

@ -63,7 +63,9 @@ contains
call time_cmfdbuild % start()
! Initialize solver
#ifdef PETSC
call gmres % create()
#endif
! Initialize matrices and vectors
call init_data(physical_adjoint)
@ -73,7 +75,9 @@ contains
call compute_adjoint()
! Set up krylov info
#ifdef PETSC
call gmres % set_oper(loss, loss)
#endif
! Stop timer for build
call time_cmfdbuild % stop()
@ -196,7 +200,9 @@ contains
s_o % val = s_o % val / k_o
! Compute new flux vector
#ifdef PETSC
call gmres % solve(s_o, phi_n)
#endif
! Compute new source vector
call prod % vector_multiply(phi_n, s_n)
@ -324,7 +330,9 @@ contains
subroutine finalize()
! Destroy all objects
call gmres % destroy()
#ifdef PETSC
call gmres % destroy()
#endif
call loss % destroy()
call prod % destroy()
call phi_n % destroy()

File diff suppressed because it is too large Load diff

View file

@ -5,10 +5,6 @@ module solver_interface
use matrix_header, only: Matrix
use vector_header, only: Vector
#ifndef PETSC
use gmres_solver, only: ilut, pgmres
#endif
implicit none
private
@ -21,25 +17,6 @@ module solver_interface
#ifdef PETSC
KSP :: ksp
PC :: pc
#else
integer :: im ! iterations before restart
integer :: iout ! unit number for printing gmres iters
integer :: maxits ! maximum iterations
integer :: lfil ! fill in parameter
integer :: iwk ! size of MSR storage
integer, allocatable :: ju(:) ! row pointers for ILUT result
integer, allocatable :: jlu(:) ! MSR storage of j
integer, allocatable :: jr(:) ! temp work array
integer, allocatable :: jwu(:) ! temp work array
integer, allocatable :: jwl(:) ! temp work array
real(8) :: tol ! gmres tolerance
real(8) :: droptol ! ILU drop tolerance
real(8), allocatable :: alu(:) ! MSR storage of ILU A
real(8), allocatable :: vv(:) ! GMRES work array
real(8), allocatable :: wu(:) ! temp work array
real(8), allocatable :: wl(:) ! temp work array
real(8), allocatable :: rhs(:) ! copy of rhs array
type(Matrix), pointer :: A ! coefficient matrix
#endif
contains
#ifdef PETSC
@ -47,11 +24,6 @@ module solver_interface
procedure :: set_oper => petsc_gmres_set_oper
procedure :: destroy => petsc_gmres_destroy
procedure :: solve => petsc_gmres_solve
#else
procedure :: create => gmres_create
procedure :: set_oper => gmres_set_oper
procedure :: destroy => gmres_destroy
procedure :: solve => gmres_solve
#endif
end type GMRESSolver
@ -76,11 +48,6 @@ module solver_interface
procedure :: destroy => petsc_jfnk_destroy
procedure :: set_functions => petsc_jfnk_set_functions
procedure :: solve => petsc_jfnk_solve
#else
procedure :: create => jfnk_create
procedure :: destroy => jfnk_destroy
procedure :: set_functions => jfnk_set_functions
procedure :: solve => jfnk_solve
#endif
end type JFNKSolver
@ -104,165 +71,7 @@ module solver_interface
contains
#ifndef PETSC
!===============================================================================
! GMRES_CREATE sets up a PETSc GMRES solver
!===============================================================================
subroutine gmres_create(self)
class(GMRESSolver) :: self
! Set up default values
self % im = 40
self % iout = 0
self % maxits = 1000
self % lfil = 0
self % tol = 1.e-10_8
self % droptol = 0.0_8
end subroutine gmres_create
!===============================================================================
! GMRES_SET_OPER sets the matrix opetors for the GMRES solver
!===============================================================================
subroutine gmres_set_oper(self, prec_mat, mat_in)
class(GMRESSolver) :: self
type(Matrix) :: prec_mat
type(Matrix), target :: mat_in
integer :: ierr
! Set pointer to matrix
self % A => mat_in
! Set up data
self % iwk = self % A % nnz + 1
! Set up arrays
allocate(self % ju(self % A % n))
allocate(self % jlu(self % iwk))
allocate(self % alu(self % iwk))
allocate(self % jr(self % A % n))
allocate(self % jwu(self % A % n))
allocate(self % jwl(self % A % n))
allocate(self % vv(self % A % n * (self % im + 1)))
allocate(self % wu(self % A % n + 1))
allocate(self % wl(self % A % n))
allocate(self % rhs(self % A % n))
! Perform ILU Preconditioning
call ilut(self % A % n, self % A % val, self % A % col, self % A % row, &
self % lfil, self % droptol, self % alu, self % jlu, self % ju, &
self % iwk, self % wu, self % wl, self % jr, self % jwl, &
self % jwu, ierr)
end subroutine gmres_set_oper
!===============================================================================
! GMRES_DESTROY frees all memory associated with the GMRES solver
!===============================================================================
subroutine gmres_destroy(self)
class(GMRESSolver) :: self
! Destroy all arrays associated with GMRES solver
deallocate(self % ju)
deallocate(self % jlu)
deallocate(self % alu)
deallocate(self % jr)
deallocate(self % jwu)
deallocate(self % jwl)
deallocate(self % vv)
deallocate(self % wu)
deallocate(self % wl)
deallocate(self % rhs)
end subroutine gmres_destroy
!===============================================================================
! GMRES_SOLVE solves the linear system
!===============================================================================
subroutine gmres_solve(self, b, x)
class(GMRESSolver) :: self
type(Vector) :: b
type(Vector) :: x
integer :: ierr ! error code
! Copy rhs
self % rhs = b % val
! Perform GMRES calculation
call pgmres(self % A % n, self % im, self % rhs, x % val, self % vv, &
self % tol, self % maxits, self % iout, self % A % val, &
self % A % col, self % A % row, self % alu, self % jlu, &
self % ju, ierr)
end subroutine gmres_solve
!===============================================================================
! JFNK_CREATE sets up a JFNK solver using PETSc SNES
!===============================================================================
subroutine jfnk_create(self)
class(JFNKSolver) :: self
message = 'JFNK cannot be used without PETSc.'
call fatal_error()
end subroutine jfnk_create
!===============================================================================
! JFNK_DESTROY frees all memory associated with JFNK solver
!===============================================================================
subroutine jfnk_destroy(self)
class(JFNKSolver) :: self
message = 'JFNK cannot be used without PETSc.'
call fatal_error()
end subroutine jfnk_destroy
!===============================================================================
! JFNK_SET_FUNCTIONS sets user functions and matrix free objects
!===============================================================================
subroutine jfnk_set_functions(self, ctx, res, jac_prec)
class(JFNKSolver) :: self
type(Jfnk_ctx) :: ctx
type(Vector) :: res
type(Matrix) :: jac_prec
message = 'JFNK cannot be used without PETSc.'
call fatal_error()
end subroutine jfnk_set_functions
!===============================================================================
! JFNK_SOLVE solves the nonlinear system
!===============================================================================
subroutine jfnk_solve(self, xvec)
class(JFNKSolver) :: self
type(Vector) :: xvec
message = 'JFNK cannot be used without PETSc.'
call fatal_error()
end subroutine jfnk_solve
#else
#ifdef PETSC
!===============================================================================
! PETSC_GMRES_CREATE sets up a PETSc GMRES solver
!===============================================================================