[petsc-users] Problem with SNESsetFunction()
Barry Smith
bsmith at petsc.dev
Tue Nov 3 09:32:01 CST 2020
I recommend you upgrade to the latest version of PETSc, it makes it easier for us to help you. Anyone can install PETSc you don't need to be an administrator on your machine.
> On Nov 3, 2020, at 3:41 AM, baikadi pranay <pranayreddy865 at gmail.com> wrote:
>
> Thank you for the response Barry. I compiled the example ex1f.F90, which uses Newton method to solve a two variable system and this works fine. However, I have a couple of questions which I was hoping you could answer. This would help me in fixing the nonlinear solver I am trying to build to solve the poisson equation.
> 1) In the ex1f.F90 file (attached), MatAssembly is done only in the FormJacobian subroutine. When I do the same thing for my code, I get the following runtime error.
MatAssembly is always needed at the end of the FormJacobian>
>
> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> [0]PETSC ERROR: Object is in wrong state
> [0]PETSC ERROR: Not for unassembled matrix
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html <http://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.11.1, Apr, 12, 2019
> [0]PETSC ERROR: ./a.out on a linux-gnu-c-debug named cg17-7.agave.rc.asu.edu <http://cg17-7.agave.rc.asu.edu/> by pbaikadi Tue Nov 3 01:52:11 2020
> [0]PETSC ERROR: Configure options
> [0]PETSC ERROR: #1 MatGetOrdering() line 180 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/mat/order/sorder.c
> [0]PETSC ERROR: #2 PCSetUp_ILU() line 134 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/pc/impls/factor/ilu/ilu.c
> [0]PETSC ERROR: #3 PCSetUp() line 932 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/pc/interface/precon.c
> [0]PETSC ERROR: #4 KSPSetUp() line 391 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/ksp/interface/itfunc.c
> [0]PETSC ERROR: #5 KSPSolve() line 725 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/ksp/interface/itfunc.c
> [0]PETSC ERROR: #6 SNESSolve_NEWTONLS() line 225 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/snes/impls/ls/ls.c
> [0]PETSC ERROR: #7 SNESSolve() line 4560 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/snes/interface/snes.c
> [0]PETSC ERROR: #8 User provided function() line 0 in User file
>
> This error goes away when I do MatAssembly after calling SNESSetJacobian. Furthermore, I don't get any such error if I don't do VecAssembly after calling SNESSetFunction.
The Assembly routines are only needed if you use SetValues() for the vectors you used vecgetarray so do not need to call VecAssembly.
> I am not sure what I might be doing wrong here. Could you give me a hint as to why this might be the case.
>
> 2) I went ahead and added the MatAssembly routines after calling the SNESSetJacobian. Now, I get the following runtime error.
> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> [0]PETSC ERROR: Object is in wrong state
> [0]PETSC ERROR: Matrix is missing diagonal entry 0
> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html <http://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting.
> [0]PETSC ERROR: Petsc Release Version 3.11.1, Apr, 12, 2019
> [0]PETSC ERROR: ./a.out on a linux-gnu-c-debug named cg17-7.agave.rc.asu.edu <http://cg17-7.agave.rc.asu.edu/> by pbaikadi Tue Nov 3 02:00:05 2020
> [0]PETSC ERROR: Configure options
> [0]PETSC ERROR: #1 MatILUFactorSymbolic_SeqAIJ() line 1687 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/mat/impls/aij/seq/aijfact.c
> [0]PETSC ERROR: #2 MatILUFactorSymbolic() line 6737 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/mat/interface/matrix.c
> [0]PETSC ERROR: #3 PCSetUp_ILU() line 144 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/pc/impls/factor/ilu/ilu.c
> [0]PETSC ERROR: #4 PCSetUp() line 932 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/pc/interface/precon.c
> [0]PETSC ERROR: #5 KSPSetUp() line 391 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/ksp/interface/itfunc.c
> [0]PETSC ERROR: #6 KSPSolve() line 725 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/ksp/interface/itfunc.c
> [0]PETSC ERROR: #7 SNESSolve_NEWTONLS() line 225 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/snes/impls/ls/ls.c
> [0]PETSC ERROR: #8 SNESSolve() line 4560 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/snes/interface/snes.c
> [0]PETSC ERROR: #9 User provided function() line 0 in User file
Presumably your matrix has zeros on the diagonal. The ILU preconditioner cannot handle zeros on the diagonal but other preconditioners can.
>
> I figured the problem might be with the preconditioner, so I switched off preconditioning. The solver now runs without any runtime error, but the converged is -3. I looked up the meaning of the number -3 and found out that it refers to "linear solve failed" and the recommendation given is to use -pc_type lu. This unfortunately brings me back to the same runtime error.
"The same" meaning zero on diagonal? Yes LU can also sometimes fail with zero on the diagonal.
> Could you give some insight on how to proceed from here?
>
> 3) In the FormFunction routine of ex1f.F90, I see that the values are set in the preconditioning matrix B. And MatAssembly is done only for B. Does this mean that the Jacobian automatically becomes equal to the preconditioning matrix B?
This depends on how you call SNESSetJacobian() if you set both matrices as the same matrix then you only need to fill up one matrix in your FormJacobian (I would always start this way);
What you should do next.
1) determine if you matrix should have zeros on the diagonal.
You can run with a very small problem with the option -mat_view and it will print the Jacobian to the screen, does it look right structurally?
Run you program with -snes_test_jacobian -snes_test_jacobian_view
This will find any errors in your Jacobian
2) If your Jacobian is suppose to have zeros on the diagonal then generally you cannot use LU or ILU as preconditioner. You can start with -pc_type fieldsplit -pc_fieldsplit_detect_saddle_point
PCFIELDSPLIT has many options.
Barry
>
> Apologies for the long email. Thank you for the help.
> Best,
> Pranay.
>
>
> On Sat, Oct 31, 2020 at 7:14 AM Barry Smith <bsmith at petsc.dev <mailto:bsmith at petsc.dev>> wrote:
>
>
>> On Oct 30, 2020, at 10:11 PM, baikadi pranay <pranayreddy865 at gmail.com <mailto:pranayreddy865 at gmail.com>> wrote:
>>
>> Hello,
>> I have a couple of questions regarding SNESSetFunction usage, when programming in Fortran90.
>>
>> 1) I have the following usage paradigm.
>> call SNESSetFunction(snes,f_non,FormFunction,0,ierr)
>> subroutine FormFunction(snes,x,r,dummy,ierr)
>> In the FormFunction subroutine, the function values are stored in the vector r. I see that these values are formed correctly. But when I use FormFunction in SNESSetFunction(), the values are not getting populated into f_non and all of the values in f_non are zero.
>
>> Should the name of the variable used to store the function value be same in SNESSetFunction and FormFunction?
>
> It does not need to be the same, they are just the variables in each function
>
>> And should I be calling the SNESComputeFunction() after calling SNESSetFunction()?
>
> No, that is a developer function called in PETSc, one would not need to call that.
>
>>
>> 2) In the subroutine FormFunction, should the vector objects created be destroyed before ending the subroutine?
>
> What vectors? If you are creating any work vectors you need within the the FormFunction, yes those should be destroyed. But not the input and output functions.
>
> Here is any example from src/snes/tutorials/ex1f.F90 Note you call VecGetArrayF90() to access the arrays for the vectors, put the values into the arrays
>
>
> subroutine FormFunction(snes,X,F,user,ierr)
> implicit none
>
> ! Input/output variables:
> SNES snes
> Vec X,F
> PetscErrorCode ierr
> type (userctx) user
> DM da
>
> ! Declarations for use with local arrays:
> PetscScalar,pointer :: lx_v(:),lf_v(:)
> Vec localX
>
> ! Scatter ghost points to local vector, using the 2-step process
> ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
> ! By placing code between these two statements, computations can
> ! be done while messages are in transition.
> call SNESGetDM(snes,da,ierr);CHKERRQ(ierr)
> call DMGetLocalVector(da,localX,ierr);CHKERRQ(ierr)
> call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr)
> call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr)
>
> ! Get a pointer to vector data.
> ! - For default PETSc vectors, VecGetArray90() returns a pointer to
> ! the data array. Otherwise, the routine is implementation dependent.
> ! - You MUST call VecRestoreArrayF90() when you no longer need access to
> ! the array.
> ! - Note that the interface to VecGetArrayF90() differs from VecGetArray(),
> ! and is useable from Fortran-90 Only.
>
> call VecGetArrayF90(localX,lx_v,ierr);CHKERRQ(ierr)
> call VecGetArrayF90(F,lf_v,ierr);CHKERRQ(ierr)
>
> ! Compute function over the locally owned part of the grid
> call FormFunctionLocal(lx_v,lf_v,user,ierr);CHKERRQ(ierr)
>
> ! Restore vectors
> call VecRestoreArrayF90(localX,lx_v,ierr);CHKERRQ(ierr)
> call VecRestoreArrayF90(F,lf_v,ierr);CHKERRQ(ierr)
>
> ! Insert values into global vector
>
> call DMRestoreLocalVector(da,localX,ierr);CHKERRQ(ierr)
> call PetscLogFlops(11.0d0*user%ym*user%xm,ierr)
>
> ! call VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr)
> ! call VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr)
> return
> end subroutine formfunction
> end module f90module
>
>
>
>>
>> Please let me know if you need any further information. Thank you in advance.
>> Best regards,
>> Pranay.
>>
>> ᐧ
>
> ᐧ
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20201103/ff224b77/attachment.html>
More information about the petsc-users
mailing list