<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class=""><div class=""><br class=""></div> 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.<div class=""><br class=""><div class=""><div><br class=""><blockquote type="cite" class=""><div class="">On Nov 3, 2020, at 3:41 AM, baikadi pranay <<a href="mailto:pranayreddy865@gmail.com" class="">pranayreddy865@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class=""><div dir="ltr" class="">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.<br class=""><b class="">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. </b><br class=""></div></div></div></blockquote><div><br class=""></div> MatAssembly is always needed at the end of the <b class="">FormJacobian></b></div><div><b class=""><br class=""></b><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><br class="">[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br class="">[0]PETSC ERROR: Object is in wrong state<br class="">[0]PETSC ERROR: Not for unassembled matrix<br class="">[0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank" class="">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br class="">[0]PETSC ERROR: Petsc Release Version 3.11.1, Apr, 12, 2019<br class="">[0]PETSC ERROR: ./a.out on a linux-gnu-c-debug named <a href="http://cg17-7.agave.rc.asu.edu/" target="_blank" class="">cg17-7.agave.rc.asu.edu</a> by pbaikadi Tue Nov 3 01:52:11 2020<br class="">[0]PETSC ERROR: Configure options<br class="">[0]PETSC ERROR: #1 MatGetOrdering() line 180 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/mat/order/sorder.c<br class="">[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<br class="">[0]PETSC ERROR: #3 PCSetUp() line 932 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/pc/interface/precon.c<br class="">[0]PETSC ERROR: #4 KSPSetUp() line 391 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/ksp/interface/itfunc.c<br class="">[0]PETSC ERROR: #5 KSPSolve() line 725 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/ksp/interface/itfunc.c<br class="">[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<br class="">[0]PETSC ERROR: #7 SNESSolve() line 4560 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/snes/interface/snes.c<br class="">[0]PETSC ERROR: #8 User provided function() line 0 in User file<br class=""></div></div></div></blockquote><div><br class=""></div> <br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><br class=""><b class="">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. </b></div></div></div></blockquote><div><br class=""></div> The Assembly routines are only needed if you use SetValues() for the vectors you used vecgetarray so do not need to call VecAssembly.</div><div><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><b class="">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.<br class=""></b><br class=""><b class="">2) I went ahead and added the MatAssembly routines after calling the SNESSetJacobian. Now, I get the following runtime error. <br class=""></b>[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br class="">[0]PETSC ERROR: Object is in wrong state<br class="">[0]PETSC ERROR: Matrix is missing diagonal entry 0<br class="">[0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank" class="">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br class="">[0]PETSC ERROR: Petsc Release Version 3.11.1, Apr, 12, 2019<br class="">[0]PETSC ERROR: ./a.out on a linux-gnu-c-debug named <a href="http://cg17-7.agave.rc.asu.edu/" target="_blank" class="">cg17-7.agave.rc.asu.edu</a> by pbaikadi Tue Nov 3 02:00:05 2020<br class="">[0]PETSC ERROR: Configure options<br class="">[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<br class="">[0]PETSC ERROR: #2 MatILUFactorSymbolic() line 6737 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/mat/interface/matrix.c<br class="">[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<br class="">[0]PETSC ERROR: #4 PCSetUp() line 932 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/pc/interface/precon.c<br class="">[0]PETSC ERROR: #5 KSPSetUp() line 391 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/ksp/interface/itfunc.c<br class="">[0]PETSC ERROR: #6 KSPSolve() line 725 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/ksp/ksp/interface/itfunc.c<br class="">[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<br class="">[0]PETSC ERROR: #8 SNESSolve() line 4560 in /packages/7x/petsc/3.11.1/petsc-3.11.1/src/snes/interface/snes.c<br class="">[0]PETSC ERROR: #9 User provided function() line 0 in User file<br class=""></div></div></div></blockquote><div><br class=""></div> Presumably your matrix has zeros on the diagonal. The ILU preconditioner cannot handle zeros on the diagonal but other preconditioners can.</div><div><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><br class=""><b class="">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.</b></div></div></div></blockquote><div><br class=""></div> "The same" meaning zero on diagonal? Yes LU can also sometimes fail with zero on the diagonal.</div><div><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><b class=""> Could you give some insight on how to proceed from here?<br class=""></b><br class=""><b class="">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? <br class=""></b></div></div></div></blockquote><div><br class=""></div> 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);</div><div><br class=""></div><div> What you should do next.</div><div><br class=""></div><div>1) determine if you matrix should have zeros on the diagonal. </div><div><br class=""></div><div> 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?</div><div><br class=""></div><div> Run you program with -snes_test_jacobian -snes_test_jacobian_view</div><div><br class=""></div><div> This will find any errors in your Jacobian</div><div><br class=""></div><div>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 </div><div><br class=""></div><div> PCFIELDSPLIT has many options.</div><div><br class=""></div><div> Barry</div><div><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class=""><br class="">Apologies for the long email. Thank you for the help.<br class="">Best,<br class="">Pranay.<br class=""><br class=""></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sat, Oct 31, 2020 at 7:14 AM Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank" class="">bsmith@petsc.dev</a>> wrote:<br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class=""><br class=""><div class=""><br class=""><blockquote type="cite" class=""><div class="">On Oct 30, 2020, at 10:11 PM, baikadi pranay <<a href="mailto:pranayreddy865@gmail.com" target="_blank" class="">pranayreddy865@gmail.com</a>> wrote:</div><br class=""><div class=""><div dir="ltr" class="">Hello,<br class="">I have a couple of questions regarding SNESSetFunction usage, when programming in Fortran90.<br class=""><br class="">1) I have the following usage paradigm.<br class=""> call SNESSetFunction(snes,f_non,FormFunction,0,ierr)<br class=""> subroutine FormFunction(snes,x,r,dummy,ierr)<br class="">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. <br class=""></div></div></blockquote><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class="">Should the name of the variable used to store the function value be same in
SNESSetFunction and FormFunction?</div></div></blockquote><div class=""><br class=""></div> It does not need to be the same, they are just the variables in each function</div><div class=""><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""> And should I be calling the SNESComputeFunction() after calling SNESSetFunction()?<br class=""></div></div></blockquote><div class=""><br class=""></div> No, that is a developer function called in PETSc, one would not need to call that.</div><div class=""><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><br class="">2) In the subroutine FormFunction, should the vector objects created be destroyed before ending the subroutine?<br class=""></div></div></blockquote><div class=""><br class=""></div> 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.</div><div class=""><br class=""></div><div class=""> 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</div><div class=""><br class=""></div><div class=""><br class=""></div><div class=""><div class="">subroutine FormFunction(snes,X,F,user,ierr)</div><div class=""> implicit none</div><div class=""><br class=""></div><div class="">! Input/output variables:</div><div class=""> SNES snes</div><div class=""> Vec X,F</div><div class=""> PetscErrorCode ierr</div><div class=""> type (userctx) user</div><div class=""> DM da</div><div class=""><br class=""></div><div class="">! Declarations for use with local arrays:</div><div class=""> PetscScalar,pointer :: lx_v(:),lf_v(:)</div><div class=""> Vec localX</div><div class=""><br class=""></div><div class="">! Scatter ghost points to local vector, using the 2-step process</div><div class="">! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().</div><div class="">! By placing code between these two statements, computations can</div><div class="">! be done while messages are in transition.</div><div class=""> call SNESGetDM(snes,da,ierr);CHKERRQ(ierr)</div><div class=""> call DMGetLocalVector(da,localX,ierr);CHKERRQ(ierr)</div><div class=""> call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr)</div><div class=""> call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr);CHKERRQ(ierr)</div><div class=""><br class=""></div><div class="">! Get a pointer to vector data.</div><div class="">! - For default PETSc vectors, VecGetArray90() returns a pointer to</div><div class="">! the data array. Otherwise, the routine is implementation dependent.</div><div class="">! - You MUST call VecRestoreArrayF90() when you no longer need access to</div><div class="">! the array.</div><div class="">! - Note that the interface to VecGetArrayF90() differs from VecGetArray(),</div><div class="">! and is useable from Fortran-90 Only.</div><div class=""><br class=""></div><div class=""> call VecGetArrayF90(localX,lx_v,ierr);CHKERRQ(ierr)</div><div class=""> call VecGetArrayF90(F,lf_v,ierr);CHKERRQ(ierr)</div><div class=""><br class=""></div><div class="">! Compute function over the locally owned part of the grid</div><div class=""> call FormFunctionLocal(lx_v,lf_v,user,ierr);CHKERRQ(ierr)</div><div class=""><br class=""></div><div class="">! Restore vectors</div><div class=""> call VecRestoreArrayF90(localX,lx_v,ierr);CHKERRQ(ierr)</div><div class=""> call VecRestoreArrayF90(F,lf_v,ierr);CHKERRQ(ierr)</div><div class=""><br class=""></div><div class="">! Insert values into global vector</div><div class=""><br class=""></div><div class=""> call DMRestoreLocalVector(da,localX,ierr);CHKERRQ(ierr)</div><div class=""> call PetscLogFlops(11.0d0*user%ym*user%xm,ierr)</div><div class=""><br class=""></div><div class="">! call VecView(X,PETSC_VIEWER_STDOUT_WORLD,ierr)</div><div class="">! call VecView(F,PETSC_VIEWER_STDOUT_WORLD,ierr)</div><div class=""> return</div><div class=""> end subroutine formfunction</div><div class=""> end module f90module</div><div class=""><br class=""></div><div class=""><br class=""></div><div class=""><br class=""></div><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><br class="">Please let me know if you need any further information. Thank you in advance.<br class="">Best regards,<br class="">Pranay.<br class=""><br class=""></div><div hspace="streak-pt-mark" style="max-height:1px" class=""><img alt="" style="width:0px;max-height:0px;overflow:hidden" src="https://mailfoogae.appspot.com/t?sender=acHJhbmF5cmVkZHk4NjVAZ21haWwuY29t&type=zerocontent&guid=72a73203-6731-4e13-aa3e-75b68ac30ecf" class=""><font color="#ffffff" size="1" class="">ᐧ</font></div>
</div></blockquote></div><br class=""></div></blockquote></div></div><div hspace="streak-pt-mark" style="max-height:1px" class=""><img alt="" style="width:0px;max-height:0px;overflow:hidden" src="https://mailfoogae.appspot.com/t?sender=acHJhbmF5cmVkZHk4NjVAZ21haWwuY29t&type=zerocontent&guid=f7e72ede-7007-4677-ba5f-3bf3f377773c" class=""><font color="#ffffff" size="1" class="">ᐧ</font></div>
</div></blockquote></div><br class=""></div></div></body></html>