[petsc-users] Back tracking/preconditioning with matrix free solvers
Kenneth Hall, SC.D.
kenneth.c.hall at duke.edu
Tue Sep 4 09:07:53 CDT 2018
Hi All,
I am solving a nonlinear problem using a matrix free approach with Petsc. My code is written in Fortran90.
I have read the documentation, and have had some success, but I am having problems with two features.
1. The code runs fine with no preconditioning. But when I using the -user_precond option, my preconditioner never gets called by Petsc, even though I make the appropriate calls to PCSetType and PCShellSetApply.
2. I am using the GMRES solver, which defaults to the L2-norm minimization line search (SNESLINESEARCHL2). I would like to use backtracking (SNESLINESEARCHBT). I have specified backtracking using SNESLineSearchSetType, but I then get the following error:
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: MatMFFDSetBase() has not been called, this is often caused by forgetting to call
MatAssemblyBegin/End on the first Mat in the SNES compute function
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.9.3, Jul, 02, 2018
[0]PETSC ERROR: mustang on a arch-darwin-c-debug-double-real named Macintosh-11.local by hall Tue Sep 4 07:29:25 2018
[0]PETSC ERROR: Configure options --with-scalar-type=real
[0]PETSC ERROR: #1 MatMult_MFFD() line 314 in /Users/hall/Documents/Fortran_Codes/Packages/petsc-3.9.3/src/mat/impls/mffd/mffd.c
[0]PETSC ERROR: #2 MatMult() line 2311 in /Users/hall/Documents/Fortran_Codes/Packages/petsc-3.9.3/src/mat/interface/matrix.c
[0]PETSC ERROR: #3 SNESLineSearchApply_BT() line 130 in /Users/hall/Documents/Fortran_Codes/Packages/petsc-3.9.3/src/snes/linesearch/impls/bt/linesearchbt.c
[0]PETSC ERROR: #4 SNESLineSearchApply() line 645 in /Users/hall/Documents/Fortran_Codes/Packages/petsc-3.9.3/src/snes/linesearch/interface/linesearch.c
[0]PETSC ERROR: #5 SNESSolve_NGMRES() line 262 in /Users/hall/Documents/Fortran_Codes/Packages/petsc-3.9.3/src/snes/impls/ngmres/snesngmres.c
[0]PETSC ERROR: #6 SNESSolve() line 4314 in /Users/hall/Documents/Fortran_Codes/Packages/petsc-3.9.3/src/snes/interface/snes.c
So it appears that Petsc is trying to do a matrix multiply. But this is a matrix free method. How do I fix this?
The relevant code is attached below. Any help with these two issues is greatly appreciated.
Code follows signature. I have omitted non-relavent portions.
Kenneth Hall
=============================================================
module petsc_krylov
use petsc
use data_module
use solution_module
implicit none
#include <petsc/finclude/petscsys.h>
#include <petsc/finclude/petsc.h>
private
public :: gmres2
SNES :: snes ! nonlinear solver
PC :: pc ! preconditioner context
KSP :: ksp ! Krylov subspace method context
SNESLineSearch :: linesearch
Vec :: x,r ! solution, residual vectors
PetscErrorCode :: ierr
PetscInt :: n
integer :: neq
integer :: icount,level=0
contains
.
.
.
subroutine gmres2
implicit none
!
!.... declare local variables
PetscInt :: its ! iterations for convergence
PetscMPIInt :: rank
PetscBool :: snes_mf,snes_mf_set,user_precond
!
!.... set problem size
call get_equation_count
n = neq
!
!.... initialliize petsc
call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ; if (ierr .ne. 0) then ; print*,'Unable to initialize PETSc' ; endif
!
!.... Create nonlinear solver context
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRA(ierr)
call SNESSetType(snes,SNESNGMRES,ierr); CHKERRA(ierr)
call SNESSetTolerances(snes, 0.0, 0.0, 0.0, 20, 10000,ierr); CHKERRQ(ierr)
!
!.... Create vectors for solution and nonlinear function
call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr); CHKERRA(ierr)
call VecDuplicate(x,r,ierr); CHKERRA(ierr)
!
!.... Set function evaluation routine and vector
call SNESSetFunction(snes,r,FormFunction,0,ierr); CHKERRA(ierr)
!
!.... Customize nonlinear solver; set runtime options
snes_mf = PETSC_FALSE;
call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-snes_mf",snes_mf,snes_mf_set,ierr); CHKERRA(ierr)
if (snes_mf .eqv. PETSC_TRUE) then
call SNESGetKSP(snes,ksp,ierr); CHKERRA(ierr)
call KSPGetPC(ksp,pc,ierr); CHKERRA(ierr)
call SNESGetLineSearch(snes,linesearch,ierr); CHKERRA(ierr)
call SNESLineSearchSetType(linesearch,SNESLINESEARCHBT, ierr); CHKERRA(ierr)
user_precond = PETSC_FALSE;
call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-user_precond",user_precond,ierr); CHKERRA(ierr)
if (user_precond .eqv. PETSC_TRUE) then
call PCSetType(pc,PCSHELL,ierr); CHKERRA(ierr)
call PCShellSetApply(pc,MatrixFreePreconditioner,ierr); CHKERRA(ierr)
else
call PCSetType(pc,PCNONE,ierr); CHKERRA(ierr)
end if
end if
!
!.... set snes options
call SNESSetFromOptions(snes,ierr); CHKERRA(ierr)
!
!.... initialize solution
call load_x_from_u(x)
!
!.... solve
call SNESSolve(snes,PETSC_NULL_VEC,x,ierr)
!
!.... Free work space. All PETSc objects should be destroyed when they are no longer needed.
call VecDestroy(x,ierr)
call VecDestroy(r,ierr)
call SNESDestroy(snes,ierr)
call PetscFinalize(ierr)
end
!
!====================================================================================================
!====================================================================================================
!
subroutine MatrixFreePreconditioner(x1,y1)
implicit none
!
!.... declare passed variables
Vec :: x1,y1 .
.
.
.
end subroutine
!
!====================================================================================================
!====================================================================================================
!
subroutine FormFunction(snes,x1,f1,dummy,ierr)
implicit none
!
!.... declare passed variables
SNES :: snes
Vec :: x1,f1,x2
PetscErrorCode :: ierr
integer :: dummy(*)
.
.
.
end subroutine FormFunction
end module
More information about the petsc-users
mailing list