mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
added gmres solver for non-petsc compilation
This commit is contained in:
Bryan Herman
parent
a726858186
commit
68a73a4a83
8 changed files with 1536 additions and 117 deletions
|
|
@ -6,13 +6,30 @@ energy_grid.o: global.o
|
|||
energy_grid.o: list_header.o
|
||||
energy_grid.o: output.o
|
||||
|
||||
ace_header.o: constants.o
|
||||
ace_header.o: endf_header.o
|
||||
|
||||
list_header.o: constants.o
|
||||
|
||||
petsc_solver.o: matrix_header.o
|
||||
petsc_solver.o: vector_header.o
|
||||
particle_header.o: constants.o
|
||||
|
||||
cmfd_slepc_solver.o: cmfd_loss_operator.o
|
||||
cmfd_slepc_solver.o: cmfd_prod_operator.o
|
||||
cmfd_slepc_solver.o: constants.o
|
||||
cmfd_slepc_solver.o: global.o
|
||||
|
||||
global.o: ace_header.o
|
||||
global.o: bank_header.o
|
||||
global.o: cmfd_header.o
|
||||
global.o: constants.o
|
||||
global.o: dict_header.o
|
||||
global.o: geometry_header.o
|
||||
global.o: hdf5_interface.o
|
||||
global.o: material_header.o
|
||||
global.o: mesh_header.o
|
||||
global.o: particle_header.o
|
||||
global.o: plot_header.o
|
||||
global.o: set_header.o
|
||||
global.o: source_header.o
|
||||
global.o: tally_header.o
|
||||
global.o: timer_header.o
|
||||
|
||||
cmfd_loss_operator.o: constants.o
|
||||
cmfd_loss_operator.o: global.o
|
||||
|
|
@ -29,8 +46,6 @@ particle_restart.o: physics.o
|
|||
particle_restart.o: random_lcg.o
|
||||
particle_restart.o: source.o
|
||||
|
||||
doppler.o: constants.o
|
||||
|
||||
tally_header.o: constants.o
|
||||
|
||||
cmfd_power_solver.o: cmfd_loss_operator.o
|
||||
|
|
@ -38,7 +53,7 @@ cmfd_power_solver.o: cmfd_prod_operator.o
|
|||
cmfd_power_solver.o: constants.o
|
||||
cmfd_power_solver.o: global.o
|
||||
cmfd_power_solver.o: matrix_header.o
|
||||
cmfd_power_solver.o: petsc_solver.o
|
||||
cmfd_power_solver.o: solver_interface.o
|
||||
cmfd_power_solver.o: vector_header.o
|
||||
|
||||
cmfd_data.o: cmfd_header.o
|
||||
|
|
@ -87,10 +102,8 @@ cmfd_prod_operator.o: constants.o
|
|||
cmfd_prod_operator.o: global.o
|
||||
cmfd_prod_operator.o: matrix_header.o
|
||||
|
||||
cmfd_slepc_solver.o: cmfd_loss_operator.o
|
||||
cmfd_slepc_solver.o: cmfd_prod_operator.o
|
||||
cmfd_slepc_solver.o: constants.o
|
||||
cmfd_slepc_solver.o: global.o
|
||||
ace_header.o: constants.o
|
||||
ace_header.o: endf_header.o
|
||||
|
||||
fission.o: ace_header.o
|
||||
fission.o: constants.o
|
||||
|
|
@ -140,7 +153,7 @@ cmfd_jfnk_solver.o: cmfd_prod_operator.o
|
|||
cmfd_jfnk_solver.o: constants.o
|
||||
cmfd_jfnk_solver.o: global.o
|
||||
cmfd_jfnk_solver.o: matrix_header.o
|
||||
cmfd_jfnk_solver.o: petsc_solver.o
|
||||
cmfd_jfnk_solver.o: solver_interface.o
|
||||
cmfd_jfnk_solver.o: vector_header.o
|
||||
|
||||
main.o: constants.o
|
||||
|
|
@ -170,21 +183,7 @@ interpolation.o: string.o
|
|||
|
||||
random_lcg.o: global.o
|
||||
|
||||
global.o: ace_header.o
|
||||
global.o: bank_header.o
|
||||
global.o: cmfd_header.o
|
||||
global.o: constants.o
|
||||
global.o: dict_header.o
|
||||
global.o: geometry_header.o
|
||||
global.o: hdf5_interface.o
|
||||
global.o: material_header.o
|
||||
global.o: mesh_header.o
|
||||
global.o: particle_header.o
|
||||
global.o: plot_header.o
|
||||
global.o: set_header.o
|
||||
global.o: source_header.o
|
||||
global.o: tally_header.o
|
||||
global.o: timer_header.o
|
||||
doppler.o: constants.o
|
||||
|
||||
string.o: constants.o
|
||||
string.o: error.o
|
||||
|
|
@ -212,7 +211,11 @@ tally.o: search.o
|
|||
tally.o: string.o
|
||||
tally.o: tally_header.o
|
||||
|
||||
particle_header.o: constants.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
|
||||
|
||||
output_interface.o: constants.o
|
||||
output_interface.o: error.o
|
||||
|
|
|
|||
|
|
@ -25,6 +25,7 @@ fixed_source.o \
|
|||
geometry.o \
|
||||
geometry_header.o \
|
||||
global.o \
|
||||
gmres_solver.o \
|
||||
hdf5_interface.o \
|
||||
hdf5_summary.o \
|
||||
initialize.o \
|
||||
|
|
@ -43,7 +44,6 @@ output_interface.o \
|
|||
particle_header.o \
|
||||
particle_restart.o \
|
||||
particle_restart_write.o \
|
||||
petsc_solver.o \
|
||||
physics.o \
|
||||
plot.o \
|
||||
plot_header.o \
|
||||
|
|
@ -51,6 +51,7 @@ ppmlib.o \
|
|||
random_lcg.o \
|
||||
search.o \
|
||||
set_header.o \
|
||||
solver_interface.o \
|
||||
source.o \
|
||||
source_header.o \
|
||||
state_point.o \
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@ module cmfd_jfnk_solver
|
|||
use cmfd_power_solver, only: cmfd_power_execute
|
||||
use cmfd_prod_operator, only: init_prod_matrix, build_prod_matrix
|
||||
use matrix_header, only: Matrix
|
||||
use petsc_solver, only: Petsc_jfnk, Jfnk_ctx
|
||||
use solver_interface, only: JFNKSolver, Jfnk_ctx
|
||||
use vector_header, only: Vector
|
||||
|
||||
implicit none
|
||||
|
|
@ -16,7 +16,7 @@ module cmfd_jfnk_solver
|
|||
type(Matrix) :: jac_prec ! Jacobian preconditioner object
|
||||
type(Matrix) :: loss ! CMFD loss matrix
|
||||
type(Matrix) :: prod ! CMFD production matrix
|
||||
type(Petsc_jfnk) :: jfnk ! JFNK solver object
|
||||
type(JFNKSolver) :: jfnk ! JFNK solver object
|
||||
type(Vector) :: resvec ! JFNK residual vector
|
||||
type(Vector) :: xvec ! JFNK solution vector
|
||||
|
||||
|
|
|
|||
|
|
@ -5,7 +5,7 @@ module cmfd_power_solver
|
|||
use cmfd_loss_operator, only: init_loss_matrix, build_loss_matrix
|
||||
use cmfd_prod_operator, only: init_prod_matrix, build_prod_matrix
|
||||
use matrix_header, only: Matrix
|
||||
use petsc_solver, only: Petsc_gmres
|
||||
use solver_interface, only: GMRESSolver
|
||||
use vector_header, only: Vector
|
||||
|
||||
implicit none
|
||||
|
|
@ -29,7 +29,7 @@ module cmfd_power_solver
|
|||
type(Vector) :: s_n ! new source vector
|
||||
type(Vector) :: s_o ! old flux vector
|
||||
type(Vector) :: serr_v ! error in source
|
||||
type(Petsc_gmres) :: gmres ! gmres solver
|
||||
type(GMRESSolver) :: gmres ! gmres solver
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -98,6 +98,7 @@ contains
|
|||
subroutine init_data(adjoint)
|
||||
|
||||
use constants, only: ONE, ZERO
|
||||
use global, only: cmfd_write_matrices
|
||||
|
||||
logical :: adjoint
|
||||
|
||||
|
|
@ -136,12 +137,16 @@ contains
|
|||
! Setup petsc for everything
|
||||
call loss % assemble()
|
||||
call prod % assemble()
|
||||
#ifdef PETSC
|
||||
call loss % setup_petsc()
|
||||
call prod % setup_petsc()
|
||||
call phi_n % setup_petsc()
|
||||
call phi_o % setup_petsc()
|
||||
call s_o % setup_petsc()
|
||||
call s_n % setup_petsc()
|
||||
if (cmfd_write_matrices) call loss % write_petsc_binary('loss.bin')
|
||||
if (cmfd_write_matrices) call prod % write_petsc_binary('prod.bin')
|
||||
#endif
|
||||
|
||||
! Set norms to 0
|
||||
norm_n = ZERO
|
||||
|
|
|
|||
1231
src/gmres_solver.F90
Normal file
1231
src/gmres_solver.F90
Normal file
File diff suppressed because it is too large
Load diff
|
|
@ -557,16 +557,10 @@ contains
|
|||
end do
|
||||
end if
|
||||
|
||||
! check for cmfd run
|
||||
! Check for cmfd run
|
||||
call lower_case(run_cmfd_)
|
||||
if (run_cmfd_ == 'true' .or. run_cmfd_ == '1') then
|
||||
cmfd_run = .true.
|
||||
#ifndef PETSC
|
||||
if (master) then
|
||||
message = 'CMFD is not available, compile OpenMC with PETSc'
|
||||
call fatal_error()
|
||||
end if
|
||||
#endif
|
||||
end if
|
||||
|
||||
end subroutine read_settings_xml
|
||||
|
|
|
|||
|
|
@ -12,8 +12,8 @@ module matrix_header
|
|||
integer :: nnz ! number of nonzeros in matrix
|
||||
integer :: n_count ! counter for length of matrix
|
||||
integer :: nz_count ! counter for number of non zeros
|
||||
integer, private, allocatable :: row(:) ! csr row vector
|
||||
integer, private, allocatable :: col(:) ! column vector
|
||||
integer, allocatable :: row(:) ! csr row vector
|
||||
integer, allocatable :: col(:) ! column vector
|
||||
real(8), allocatable :: val(:) ! matrix value vector
|
||||
# ifdef PETSC
|
||||
Mat :: petsc_mat
|
||||
|
|
|
|||
|
|
@ -1,8 +1,14 @@
|
|||
module petsc_solver
|
||||
module solver_interface
|
||||
|
||||
use error, only: fatal_error
|
||||
use global, only: message
|
||||
use matrix_header, only: Matrix
|
||||
use vector_header, only: Vector
|
||||
|
||||
#ifndef PETSC
|
||||
use gmres_solver, only: ilut, pgmres
|
||||
#endif
|
||||
|
||||
implicit none
|
||||
private
|
||||
|
||||
|
|
@ -10,18 +16,44 @@ module petsc_solver
|
|||
# include <finclude/petsc.h90>
|
||||
#endif
|
||||
|
||||
! PETSc GMRES solver type
|
||||
type, public :: Petsc_gmres
|
||||
! GMRES solver type
|
||||
type, public :: GMRESSolver
|
||||
#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
|
||||
procedure :: create => petsc_gmres_create
|
||||
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
|
||||
end type Petsc_gmres
|
||||
#endif
|
||||
end type GMRESSolver
|
||||
|
||||
! Derived type to contain list of data needed to be passed to procedures
|
||||
type, public :: Jfnk_ctx
|
||||
|
|
@ -29,8 +61,8 @@ module petsc_solver
|
|||
procedure (jac_interface), pointer, nopass :: jac_proc_ptr
|
||||
end type Jfnk_ctx
|
||||
|
||||
! Petsc SNES JFNK solver type
|
||||
type, public :: Petsc_jfnk
|
||||
! JFNK solver type
|
||||
type, public :: JFNKSolver
|
||||
#ifdef PETSC
|
||||
KSP :: ksp
|
||||
PC :: pc
|
||||
|
|
@ -39,13 +71,18 @@ module petsc_solver
|
|||
SNESLineSearch :: ls
|
||||
#endif
|
||||
contains
|
||||
#ifdef PETSC
|
||||
procedure :: create => petsc_jfnk_create
|
||||
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
|
||||
end type Petsc_jfnk
|
||||
|
||||
integer :: petsc_err
|
||||
#endif
|
||||
end type JFNKSolver
|
||||
|
||||
! Abstract interface stating how jacobian and residual routines look
|
||||
abstract interface
|
||||
|
|
@ -61,30 +98,28 @@ module petsc_solver
|
|||
end subroutine jac_interface
|
||||
end interface
|
||||
|
||||
#ifdef PETSC
|
||||
integer :: petsc_err ! petsc error code
|
||||
#endif
|
||||
|
||||
contains
|
||||
|
||||
#ifndef PETSC
|
||||
!===============================================================================
|
||||
! GMRES_CREATE sets up a PETSc GMRES solver
|
||||
!===============================================================================
|
||||
|
||||
subroutine gmres_create(self)
|
||||
|
||||
class(Petsc_gmres) :: self
|
||||
class(GMRESSolver) :: self
|
||||
|
||||
integer :: ilu_levels = 5
|
||||
real(8) :: rtol = 1.0e-10_8
|
||||
real(8) :: atol = 1.0e-10_8
|
||||
|
||||
#ifdef PETSC
|
||||
call KSPCreate(PETSC_COMM_WORLD, self % ksp, petsc_err)
|
||||
call KSPSetTolerances(self % ksp, rtol, atol, &
|
||||
PETSC_DEFAULT_DOUBLE_PRECISION, PETSC_DEFAULT_INTEGER, petsc_err)
|
||||
call KSPSetType(self % ksp, KSPGMRES, petsc_err)
|
||||
call KSPSetInitialGuessNonzero(self % ksp, PETSC_TRUE, petsc_err)
|
||||
call KSPGetPC(self % ksp, self % pc, petsc_err)
|
||||
call PCFactorSetLevels(self % pc, ilu_levels, petsc_err)
|
||||
call KSPSetFromOptions(self % ksp, petsc_err)
|
||||
#endif
|
||||
! 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
|
||||
|
||||
|
|
@ -94,15 +129,35 @@ contains
|
|||
|
||||
subroutine gmres_set_oper(self, prec_mat, mat_in)
|
||||
|
||||
class(Petsc_gmres) :: self
|
||||
class(GMRESSolver) :: self
|
||||
type(Matrix) :: prec_mat
|
||||
type(Matrix) :: mat_in
|
||||
type(Matrix), target :: mat_in
|
||||
|
||||
#ifdef PETSC
|
||||
call KSPSetOperators(self % ksp, mat_in % petsc_mat, prec_mat % petsc_mat, &
|
||||
SAME_NONZERO_PATTERN, petsc_err)
|
||||
call KSPSetUp(self % ksp, petsc_err)
|
||||
#endif
|
||||
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
|
||||
|
||||
|
|
@ -112,11 +167,19 @@ contains
|
|||
|
||||
subroutine gmres_destroy(self)
|
||||
|
||||
class(Petsc_gmres) :: self
|
||||
class(GMRESSolver) :: self
|
||||
|
||||
#ifdef PETSC
|
||||
call KSPDestroy(self % ksp, petsc_err)
|
||||
#endif
|
||||
! 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
|
||||
|
||||
|
|
@ -126,13 +189,20 @@ contains
|
|||
|
||||
subroutine gmres_solve(self, b, x)
|
||||
|
||||
class(Petsc_gmres) :: self
|
||||
class(GMRESSolver) :: self
|
||||
type(Vector) :: b
|
||||
type(Vector) :: x
|
||||
|
||||
#ifdef PETSC
|
||||
call KSPSolve(self % ksp, b % petsc_vec, x % petsc_vec, petsc_err)
|
||||
#endif
|
||||
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
|
||||
|
||||
|
|
@ -142,9 +212,130 @@ contains
|
|||
|
||||
subroutine jfnk_create(self)
|
||||
|
||||
class(Petsc_jfnk) :: 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
|
||||
!===============================================================================
|
||||
! PETSC_GMRES_CREATE sets up a PETSc GMRES solver
|
||||
!===============================================================================
|
||||
|
||||
subroutine petsc_gmres_create(self)
|
||||
|
||||
class(GMRESSolver) :: self
|
||||
|
||||
integer :: ilu_levels = 5
|
||||
real(8) :: rtol = 1.0e-10_8
|
||||
real(8) :: atol = 1.0e-10_8
|
||||
|
||||
call KSPCreate(PETSC_COMM_WORLD, self % ksp, petsc_err)
|
||||
call KSPSetTolerances(self % ksp, rtol, atol, &
|
||||
PETSC_DEFAULT_DOUBLE_PRECISION, PETSC_DEFAULT_INTEGER, petsc_err)
|
||||
call KSPSetType(self % ksp, KSPGMRES, petsc_err)
|
||||
call KSPSetInitialGuessNonzero(self % ksp, PETSC_TRUE, petsc_err)
|
||||
call KSPGetPC(self % ksp, self % pc, petsc_err)
|
||||
call PCFactorSetLevels(self % pc, ilu_levels, petsc_err)
|
||||
call KSPSetFromOptions(self % ksp, petsc_err)
|
||||
|
||||
end subroutine petsc_gmres_create
|
||||
|
||||
!===============================================================================
|
||||
! PETSC_GMRES_SET_OPER sets the matrix opetors for the GMRES solver
|
||||
!===============================================================================
|
||||
|
||||
subroutine petsc_gmres_set_oper(self, prec_mat, mat_in)
|
||||
|
||||
class(GMRESSolver) :: self
|
||||
type(Matrix) :: prec_mat
|
||||
type(Matrix) :: mat_in
|
||||
|
||||
call KSPSetOperators(self % ksp, mat_in % petsc_mat, prec_mat % petsc_mat, &
|
||||
SAME_NONZERO_PATTERN, petsc_err)
|
||||
call KSPSetUp(self % ksp, petsc_err)
|
||||
|
||||
end subroutine petsc_gmres_set_oper
|
||||
|
||||
!===============================================================================
|
||||
! PETSC_GMRES_DESTROY frees all memory associated with the GMRES solver
|
||||
!===============================================================================
|
||||
|
||||
subroutine petsc_gmres_destroy(self)
|
||||
|
||||
class(GMRESSolver) :: self
|
||||
|
||||
call KSPDestroy(self % ksp, petsc_err)
|
||||
|
||||
end subroutine petsc_gmres_destroy
|
||||
|
||||
!===============================================================================
|
||||
! PETSC_GMRES_SOLVE solves the linear system
|
||||
!===============================================================================
|
||||
|
||||
subroutine petsc_gmres_solve(self, b, x)
|
||||
|
||||
class(GMRESSolver) :: self
|
||||
type(Vector) :: b
|
||||
type(Vector) :: x
|
||||
|
||||
call KSPSolve(self % ksp, b % petsc_vec, x % petsc_vec, petsc_err)
|
||||
|
||||
end subroutine petsc_gmres_solve
|
||||
|
||||
!===============================================================================
|
||||
! PETSC_JFNK_CREATE sets up a JFNK solver using PETSc SNES
|
||||
!===============================================================================
|
||||
|
||||
subroutine petsc_jfnk_create(self)
|
||||
|
||||
class(JFNKSolver) :: self
|
||||
|
||||
#ifdef PETSC
|
||||
! Turn on mf_operator option for matrix free jacobian
|
||||
call PetscOptionsSetValue("-snes_mf_operator", "TRUE", petsc_err)
|
||||
|
||||
|
|
@ -166,75 +357,68 @@ contains
|
|||
call PCFactorSetLevels(self % pc, 5, petsc_err)
|
||||
call PCSetFromOptions(self % pc, petsc_err)
|
||||
call KSPSetFromOptions(self % ksp, petsc_err)
|
||||
#endif
|
||||
|
||||
end subroutine jfnk_create
|
||||
end subroutine petsc_jfnk_create
|
||||
|
||||
!===============================================================================
|
||||
! JFNK_DESTROY frees all memory associated with JFNK solver
|
||||
! PETSC_JFNK_DESTROY frees all memory associated with JFNK solver
|
||||
!===============================================================================
|
||||
|
||||
subroutine jfnk_destroy(self)
|
||||
subroutine petsc_jfnk_destroy(self)
|
||||
|
||||
class(Petsc_jfnk) :: self
|
||||
class(JFNKSolver) :: self
|
||||
|
||||
#ifdef PETSC
|
||||
call MatDestroy(self % jac_mf, petsc_err)
|
||||
call SNESDestroy(self % snes, petsc_err)
|
||||
#endif
|
||||
|
||||
end subroutine jfnk_destroy
|
||||
end subroutine petsc_jfnk_destroy
|
||||
|
||||
!===============================================================================
|
||||
! JFNK_SET_FUNCTIONS sets user functions and matrix free objects
|
||||
! PETSC_JFNK_SET_FUNCTIONS sets user functions and matrix free objects
|
||||
!===============================================================================
|
||||
|
||||
subroutine jfnk_set_functions(self, ctx, res, jac_prec)
|
||||
subroutine petsc_jfnk_set_functions(self, ctx, res, jac_prec)
|
||||
|
||||
class(Petsc_jfnk) :: self
|
||||
class(JFNKSolver) :: self
|
||||
type(Jfnk_ctx) :: ctx
|
||||
type(Vector) :: res
|
||||
type(Matrix) :: jac_prec
|
||||
|
||||
# ifdef PETSC
|
||||
! Set residual procedure
|
||||
call SNESSetFunction(self % snes, res % petsc_vec, jfnk_compute_residual, &
|
||||
ctx, petsc_err)
|
||||
call SNESSetFunction(self % snes, res % petsc_vec, &
|
||||
petsc_jfnk_compute_residual, ctx, petsc_err)
|
||||
|
||||
! Create the matrix free jacobian
|
||||
call MatCreateSNESMF(self % snes, self % jac_mf, petsc_err)
|
||||
|
||||
! Set Jacobian procedure
|
||||
call SNESSetJacobian(self % snes, self % jac_mf, jac_prec % petsc_mat, &
|
||||
jfnk_compute_jacobian, ctx, petsc_err)
|
||||
petsc_jfnk_compute_jacobian, ctx, petsc_err)
|
||||
|
||||
! Set up Jacobian Lags
|
||||
call SNESSetLagJacobian(self % snes, -2, petsc_err)
|
||||
call SNESSetLagPreconditioner(self % snes, -1, petsc_err)
|
||||
#endif
|
||||
|
||||
end subroutine jfnk_set_functions
|
||||
end subroutine petsc_jfnk_set_functions
|
||||
|
||||
!===============================================================================
|
||||
! JFNK_SOLVE solves the nonlinear system
|
||||
! PETSC_JFNK_SOLVE solves the nonlinear system
|
||||
!===============================================================================
|
||||
|
||||
subroutine jfnk_solve(self, xvec)
|
||||
subroutine petsc_jfnk_solve(self, xvec)
|
||||
|
||||
class(Petsc_jfnk) :: self
|
||||
class(JFNKSolver) :: self
|
||||
type(Vector) :: xvec
|
||||
|
||||
#ifdef PETSC
|
||||
call SNESSolve(self % snes, PETSC_NULL_DOUBLE, xvec % petsc_vec, petsc_err)
|
||||
#endif
|
||||
|
||||
end subroutine jfnk_solve
|
||||
end subroutine petsc_jfnk_solve
|
||||
|
||||
!===============================================================================
|
||||
! JFNK_COMPUTE_RESIDUAL buffer routine to user specifed residual routine
|
||||
! PETSC_JFNK_COMPUTE_RESIDUAL buffer routine to user specifed residual routine
|
||||
!===============================================================================
|
||||
# ifdef PETSC
|
||||
subroutine jfnk_compute_residual(snes, x, res, ctx, ierr)
|
||||
|
||||
subroutine petsc_jfnk_compute_residual(snes, x, res, ctx, ierr)
|
||||
|
||||
SNES :: snes
|
||||
Vec :: x
|
||||
|
|
@ -260,13 +444,13 @@ contains
|
|||
call VecRestoreArrayF90(x, xvec % val, ierr)
|
||||
call VecRestoreArrayF90(res, resvec % val, ierr)
|
||||
|
||||
end subroutine jfnk_compute_residual
|
||||
#endif
|
||||
end subroutine petsc_jfnk_compute_residual
|
||||
|
||||
!===============================================================================
|
||||
! JFNK_COMPUTE_JACOBIAN
|
||||
! PETSC_JFNK_COMPUTE_JACOBIAN buffer routine to user specified jacobian routine
|
||||
!===============================================================================
|
||||
#ifdef PETSC
|
||||
subroutine jfnk_compute_jacobian(snes, x, jac_mf, jac_prec, flag, ctx, &
|
||||
|
||||
subroutine petsc_jfnk_compute_jacobian(snes, x, jac_mf, jac_prec, flag, ctx, &
|
||||
ierr)
|
||||
|
||||
SNES :: snes
|
||||
|
|
@ -294,6 +478,7 @@ contains
|
|||
call MatAssemblyBegin(jac_mf, MAT_FINAL_ASSEMBLY, petsc_err)
|
||||
call MatAssemblyEnd(jac_mf, MAT_FINAL_ASSEMBLY, petsc_err)
|
||||
|
||||
end subroutine jfnk_compute_jacobian
|
||||
end subroutine petsc_jfnk_compute_jacobian
|
||||
#endif
|
||||
end module petsc_solver
|
||||
|
||||
end module solver_interface
|
||||
Loading…
Add table
Add a link
Reference in a new issue