[petsc-users] matrix free non linear solver Fortran

[Smith, Barry F.]

  Also from the users manual for MatCreateSNESMF()

   Notes:
     You can call SNESSetJacobian() with MatMFFDComputeJacobian() if you are using matrix and not a different
     preconditioner matrix



> On Jan 16, 2018, at 12:33 PM, Matthew Knepley <knepley at gmail.com> wrote:
> 
> On Tue, Jan 16, 2018 at 1:26 PM, El Haddad, Fadi <f.el-haddad at imperial.ac.uk> wrote:
> Hello Petsc Users,
> 
> I am trying to set up a snes matrix free solver into an existing code.
> I obtained the following log error by petsc :
> 
> Setting it up is a little tricky, Why not just omit SNESSetJacobian() and it will automatically set it up?
> 
>   Thanks,
> 
>     Matt
>  
> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> [0]PETSC ERROR: No support for this operation for this object type
> [0]PETSC ERROR: Mat type mffd
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
> [0]PETSC ERROR: Petsc Development GIT revision: unknown  GIT Date: unknown
> [0]PETSC ERROR: ../Forse/JM62-FORSE/petsc-infconv/forse on a arch-linux2-c-debug named me-felhadda by felhadda Tue Jan 16 12:49:53 2018
> [0]PETSC ERROR: Configure options --with-cc=mpicc --with-cxx=mpic++ --with-fc=mpif90 --with-mpi-dir=/usr/lib/mpich --download-fblaslapack
> [0]PETSC ERROR: #1 MatZeroEntries() line 5617 in /opt/petsc/petsc-infconv/src/mat/interface/matrix.c
> [0]PETSC ERROR: #2 SNESComputeJacobianDefault() line 64 in /opt/petsc/petsc-infconv/src/snes/interface/snesj.c
> [0]PETSC ERROR: #3 oursnesjacobian() line 97 in /opt/petsc/petsc-infconv/src/snes/interface/ftn-custom/zsnesf.c
> [0]PETSC ERROR: #4 SNESComputeJacobian() line 2374 in /opt/petsc/petsc-infconv/src/snes/interface/snes.c
> [0]PETSC ERROR: #5 SNESSolve_NEWTONLS() line 222 in /opt/petsc/petsc-infconv/src/snes/impls/ls/ls.c
> [0]PETSC ERROR: #6 SNESSolve() line 4179 in /opt/petsc/petsc-infconv/src/snes/interface/snes.c
> 
> 
> Have anyone dealt with a similar error? I copied my snes matrix free set up here below : 
> 
> 
>                 external :: FormFunction
>                 external :: MySNESConverged
> 
> 
>                 call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
>                 if (ierr .ne. 0) then
>                         print*,'Unable to initialize PETSc'
>                         stop
>                 endif
>                 call MPI_Comm_size(PETSC_COMM_SELF,size,ierr)
>                 call MPI_Comm_rank(PETSC_COMM_SELF,rank,ierr)
> 
>                 
>                 call VecCreateSeq(PETSC_COMM_SELF,n,x,ierr)
>                 call VecCreateSeq(PETSC_COMM_SELF,n,r,ierr)
>                 
>                 call SNESCreate(PETSC_COMM_SELF,snes,ierr)
>                 call SNESSetApplicationContext(snes,userCtxptrs,ierr);CHKERRA(ierr)
>                 
>                 call SNESSetFunction(snes,r,FormFunction,userCtxptrs,ierr);CHKERRA(ierr)
>                 call MatCreateSnesMF(snes,JMF,ierr) 
>                 call MatMFFDWPSetComputeNormU(JMF,PETSC_FALSE,ierr)
> 
>                 call SNESSetJacobian(snes,Jmf,Jmf,SNESComputeJacobianDefault,userCtxptrs,ierr)
>                 call SNESSetUseMatrixFree(self%snes,PETSC_TRUE,PETSC_TRUE,ierr)
>                 call SNESSetConvergenceTest(snes,MySNESConverged,userCtxptrs,PETSC_NULL_FUNCTION,ierr)
>                 
>                 call SNESSetConvergedReason(snes,2,ierr)
>                 call SNESSetNormSchedule(snes,SNES_NORM_ALWAYS,ierr)
> 
>                 call SnesSetTolerances(snes,1.d-8,1.d-7,1.d-7,10,25,ierr)
>                 call SNESSetFromOptions(snes,ierr);CHKERRA(ierr)
>                 call SNESView(snes,PETSC_VIEWER_STDOUT_SELF,ierr)
>               
>                 call getpointersUserCtx(ds,pt,TypeY,TypeF,userCtxptrs)
>                 call FormInitialGuess(x,pt)
>                
>                 call SNESsolve(snes,PETSC_NULL_VEC,x,ierr)
>                
>                 call SNESGetConvergedReason(snes,reason,ierr)
>                 IF (reason .eq. 2 ) then 
>                     call VecGetArrayReadF90(x,xx,ierr)
>                     pt%Y = xx
>                     call VecRestoreArrayReadF90(x,xx,ierr)
>                     ConvFlag = True
>                 ELSE
>                     ConvFlag = False
>                 ENDIF
> 
> Am I missing a crucial subroutine call?  
> Thank you for the help
> 
> Fadi
> 
> 
> 
> 
> 
> 
> 
> 
> 
> -- 
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> -- Norbert Wiener
> 
> https://www.cse.buffalo.edu/~knepley/