<html><head><meta http-equiv="Content-Type" content="text/html; charset=us-ascii"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class=""><br class=""><div><br class=""><blockquote type="cite" class=""><div class="">On May 25, 2021, at 10:28 AM, Saransh Saxena <<a href="mailto:saransh.saxena5571@gmail.com" class="">saransh.saxena5571@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class="">Hi Barry,<br class=""><br class=""><blockquote type="cite" style="color:rgb(80,0,80)" class=""><div dir="ltr" class=""> Mat Apetsc = A.getpetsc();<br class=""> Vec bpetsc = b.getpetsc();</div></blockquote><font class="">Apetsc and bpetsc are Matrix A and vector b (in petsc format), A and b are using different class structure (as per the FEM code) in solving the nonlinear equation A(x).x = b. b is the RHS vector (applied forces in my case) and A is global stiffness matrix (K for static simulations in FEM terms). x is the solution vector (displacements in my case for FEM simulation). r is the residual vector of the form r = b - A(x).x. Only Matrix A is a function of the output in the current case, but I am implementing for a general case where b might also depend on the output.</font><div class=""><font class=""><br class=""></font></div><div class=""><blockquote type="cite" style="color:rgb(80,0,80)" class=""><div dir="ltr" class="">Vec output;<br class="">VecDuplicate(x, &output);<br class="">VecCopy(x, output);<br class=""><br class="">setdata(vec(b.getpointer()->getdofmanager(), output));</div></blockquote><font class="">The above lines store the solution for the current iteration so that when the .generate() function is called, updated A matrix is obtained (and updated b as well for a general case where both A and b vary with x, the output). I have to do it by copying the x vector to output because setdata() destroys the vector when called.</font></div><div class=""><font class=""><br class=""></font></div><div class=""><font class="">I was also browsing through the definition of SNESSetFunction and realized that it solves </font><span style="" class="">f'(x) x = -f(x), <font face="arial, sans-serif" class="">however, in newton raphson x_(n+1) = x_(n) - f(x_(n))/f'(x_(n)). So am I solving for delta_x here with SNESSetFunction?</font></span></div><div class=""><font face="arial, sans-serif" class=""><br class=""></font></div><div class=""><font face="arial, sans-serif" class="">Also in SNESSetPicard(), I need to pass a function to compute b. However, in my case b is constant. How do I use that?</font></div></div></div></blockquote><div><br class=""></div><div> You have two choices. </div><div> * As Matt says you can just provide a function that copies over your constant b vector each time. </div><div> * Or you can pass a NULL for the function and call SNESSolve(snes,b,x) where b is your constant b vector (this will be slightly more efficient)</div><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div class=""><font face="arial, sans-serif" class=""> Also does Vec r in the definition refer to solution vector or residual vector?<br class=""></font></div></div></div></blockquote><div><br class=""></div> r is a vector that SNES will use to compute the residual in. It is just a work vector you can pass in. You can pass in NULL and PETSc will create the work vector it needs internally.</div><div><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div class=""><font face="arial, sans-serif" class=""><br class="">Best regards,<br class="">Saransh<br class=""></font><div class=""><font class=""><br class=""> </font><div class=""><font class=""><br class=""></font></div><div class=""><font class=""><br class=""></font></div></div></div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, May 25, 2021 at 10:15 AM Barry Smith <<a href="mailto:bsmith@petsc.dev" 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 style="overflow-wrap: break-word;" class=""><div class=""><br class=""></div><blockquote type="cite" class=""><div dir="ltr" class=""> VecNorm(F, NORM_2, &normvalres);</div></blockquote><div class=""><br class=""></div> The F has not yet been computed by you so you shouldn't take the norm here. F could have anything in it. You should take the norm after the line<div class=""><br class=""></div><div class=""><blockquote type="cite" class=""><div dir="ltr" class="">VecAXPY(F,-1.0, bpetsc);</div></blockquote><div class=""><br class=""></div><br class=""><div class=""><br class=""></div><div class=""><blockquote type="cite" class=""><div dir="ltr" class=""> // Read pointers to A and b:<br class=""> Mat Apetsc = A.getpetsc();<br class=""> Vec bpetsc = b.getpetsc();</div></blockquote></div><div class=""><br class=""></div> Where are these things computed and are they both functions of output? or is b merely x (the current solution snes is working with)</div><div class=""><br class=""></div><div class=""><blockquote type="cite" class=""><div dir="ltr" class="">Vec output;<br class=""> VecDuplicate(x, &output);<br class=""> VecCopy(x, output);<br class=""><br class=""> setdata(vec(b.getpointer()->getdofmanager(), output));</div></blockquote><div class=""><br class=""></div><div class=""> What is the line above doing?</div><div class=""><br class=""></div> I think you using using Picard iteration A(x^n) x^{n+1} = b(x^n). (Sometimes people call this a fixed-point iteration) If so you should use SNESSetPicard() and not SNESSetFunction(). </div><div class=""><br class=""></div><div class=""> If you run with SNESSetPicard() with no additional options it will run the defect correction version of Picard</div><div class=""><br class=""></div><div class=""> If you run with SNESSetPicard() and use -snes_mf_operator then SNES will run matrix-free Newton's method using your A as the preconditioner for the Jacobian</div><div class=""><br class=""></div><div class=""> If you run with SNESSetPicard() and use -snes_fd then SNES will form explicitly the Jacobian and run Newton's method with it. This will be very slow but you gives you an idea of how Newton's method works on your problem. </div><div class=""><br class=""></div><div class=""> If you call SNESSetFromOptions() before SNESSolve() then you can use -snes_monitor -ksp_monitor -snes_converged_reason and many other options to monitor the convergence, then you will not have to compute the norms yourself and put print statements in your code for the norms.<br class=""><div class=""><br class=""></div><div class=""> Barry</div><div class=""><br class=""><div class=""><br class=""><blockquote type="cite" class=""><div class="">On May 25, 2021, at 2:48 AM, Saransh Saxena <<a href="mailto:saransh.saxena5571@gmail.com" target="_blank" class="">saransh.saxena5571@gmail.com</a>> wrote:</div><br class=""><div class=""><div dir="ltr" class="">Hi guys,<br class=""><br class="">I've written an implementation of SNES within my code to use the petsc nonlinear solvers but for some reason, I am getting results I can't make sense of. To summarize, I've written a function to calculate residual using Matthew's suggestion. However, when I run the code, the behaviour is odd, the solver seems to enter the myresidual function initially. However, after that it never updates the iteration counter and the solution vector remains unchanged (and a really small value) while the residual vector explodes in value. <br class=""><br class="">Residual code :-<br class=""><br class="">PetscErrorCode sl::myresidual(SNES snes, Vec x, Vec F, void *ctx)<br class="">{<br class=""> // Cast the application context:<br class=""> sl::formulCtx *user = (sl::formulCtx*)ctx;<br class=""><br class=""> // Read the formulation:<br class=""> formulation *thisformul = (*user).formul;<br class=""> thisformul->generate();<br class=""><br class=""> //vec *b = user->b;<br class=""> //mat *A = user->A;<br class=""><br class=""> vec b = thisformul->b();<br class=""> mat A = thisformul->A();<br class=""><br class=""> // Read pointers to A and b:<br class=""> Mat Apetsc = A.getpetsc();<br class=""> Vec bpetsc = b.getpetsc();<br class=""><br class=""> double normvalres, normvalsol;<br class=""> VecNorm(F, NORM_2, &normvalres);<br class=""> VecNorm(x, NORM_2, &normvalsol);<br class=""> std::cout << "----------------------------------------------------------------------------" << std::endl;<br class=""> std::cout << "Entered residual function, norm of residual vector is : " << normvalres << std::endl;<br class=""> std::cout << "Entered residual function, norm of solution vector is : " << normvalsol << std::endl;<br class=""><br class=""> // Compute the residual as F = A*x - b<br class=""> MatMult(Apetsc, x, F);<br class=""> VecAXPY(F,-1.0, bpetsc);<br class=""><br class=""> Vec output;<br class=""> VecDuplicate(x, &output);<br class=""> VecCopy(x, output);<br class=""><br class=""> setdata(vec(b.getpointer()->getdofmanager(), output));<br class=""><br class=""> std::cout << "Writing the sol to fields \n";<br class=""><br class=""> return 0; <br class="">}<br class=""><div class=""><br class=""></div><div class="">SNES implementation :-<br class=""><br class="">void sl::solvenonlinear(formulation thisformul, double restol, int maxitnum)<br class="">{ <br class=""> // Make sure the problem is of the form Ax = b:<br class=""> if (thisformul.isdampingmatrixdefined() || thisformul.ismassmatrixdefined())<br class=""> {<br class=""> std::cout << "Error in 'sl' namespace: formulation to solve cannot have a damping/mass matrix (use a time resolution algorithm)" << std::endl;<br class=""> abort(); <br class=""> }<br class=""><br class=""> // Remove leftovers (if any):<br class=""> mat Atemp = thisformul.A(); vec btemp = thisformul.b();<br class=""><br class=""> // Create Application Context for formulation<br class=""> sl::formulCtx user;<br class=""> user.formul = &thisformul;<br class=""><br class=""> // Generate formulation to set PETSc SNES requirements:<br class=""> thisformul.generate();<br class=""><br class=""> mat A = thisformul.A();<br class=""> vec b = thisformul.b();<br class=""><br class=""> // SNES requirements:<br class=""> Vec bpetsc = b.getpetsc();<br class=""> Mat Apetsc = A.getpetsc();<br class=""><br class=""> vec residual(std::shared_ptr<rawvec>(new rawvec(b.getpointer()->getdofmanager())));<br class=""> Vec residualpetsc = residual.getpetsc();<br class=""><br class=""> vec sol(std::shared_ptr<rawvec>(new rawvec(b.getpointer()->getdofmanager())));<br class=""> Vec solpetsc = sol.getpetsc();<br class=""><br class=""> //Retrieve the SNES and KSP Context from A matrix:<br class=""> SNES* snes = A.getpointer()->getsnes();<br class=""> KSP* ksp = A.getpointer()->getksp();<br class=""> <br class=""> // Create placeholder for preconditioner:<br class=""> PC pc;<br class=""><br class=""> // Create snes context:<br class=""> SNESCreate(PETSC_COMM_SELF, snes);<br class=""> SNESSetFunction(*snes, residualpetsc, sl::myresidual, &user);<br class=""> SNESSetTolerances(*snes, PETSC_DEFAULT, restol, PETSC_DEFAULT, maxitnum, 5);<br class=""><br class=""> // Retrieve the KSP context automatically created:<br class=""> SNESGetKSP(*snes, ksp);<br class=""><br class=""> //Set KSP specific parameters/options:<br class=""> KSPSetOperators(*ksp, Apetsc, Apetsc);<br class=""> KSPSetFromOptions(*ksp);<br class=""> KSPGetPC(*ksp,&pc);<br class=""> PCSetType(pc,PCLU);<br class=""> PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);<br class=""><br class=""> //Call SNES options to invoke changes from console:<br class=""> SNESSetFromOptions(*snes);<br class=""><br class=""> // Set SNES Monitor to retrieve convergence information:<br class=""> SNESMonitorSet(*snes, sl::mysnesmonitor, PETSC_NULL, PETSC_NULL);<br class=""> //SNESMonitorLGResidualNorm();<br class=""><br class=""> SNESSolve(*snes, PETSC_NULL, solpetsc);<br class=""><br class=""> // Print the norm of residual:<br class=""> double normres;<br class=""> VecNorm(residualpetsc, NORM_2, &normres);<br class=""> std::cout << "L2 norm of the residual is : " << normres << std::endl; <br class=""><br class=""> //Set the solution to all the fields:<br class=""> setdata(sol);<br class=""><br class=""> // Get the number of required iterations and the residual norm:<br class=""> //SNESGetIterationNumber(*snes, &maxitnum);<br class=""> //SNESGetResidualNorm(*snes, &restol);<br class=""><br class=""> // Destroy SNES context once done with computation:<br class=""> SNESDestroy(snes);<br class=""><br class="">}<br class=""></div><div class=""><br class=""></div><div class="">Output :-<br class=""><span id="gmail-m_6268095669457834418cid:ii_kp3qilht0" class=""><image.png></span><br class=""><br class=""></div><div class="">Am I doing something incorrect wrt SNES? When I use the linear solver (KSP) and manually coded fixed point nonlinear iteration, it works fine.<br class=""><br class="">Best regards,<br class="">Saransh<br class=""><br class=""></div><div class=""><br class=""></div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sun, May 9, 2021 at 10:43 PM 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=""><div class=""><br class=""></div> Saransh,<div class=""><br class=""></div><div class=""> If Picard or Newton's method does not converge, you can consider adding pseudo-transient and/or other continuation methods. For example, if the problem is made difficult by certain physical parameters you can start with "easier" values of the parameters, solve the nonlinear system, then use its solution as the initial guess for slightly more "difficult" parameters, etc. Or, depending on the problem grid sequencing may be appropriate. We have some tools to help with all these approaches.</div><div class=""><br class=""></div><div class=""> Barry</div><div class=""><br class=""><div class=""><br class=""><blockquote type="cite" class=""><div class="">On May 9, 2021, at 2:07 PM, Saransh Saxena <<a href="mailto:saransh.saxena5571@gmail.com" target="_blank" class="">saransh.saxena5571@gmail.com</a>> wrote:</div><br class=""><div class=""><div dir="ltr" class=""><div class=""><div class=""><div class="">Thanks Barry and Matt,<br class=""><br class=""></div>Till now I was only using a simple fixed point nonlinear solver manually coded instead of ones provided by PETSc. However, the problem I am trying to solve is highly nonlinear so I suppose I'll need at least a newton based solver to start with. I'll get back to you guys if I have any questions.<br class=""><br class=""></div>Cheers,<br class=""></div>Saransh<br class=""></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Sat, May 8, 2021 at 5:18 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=""><div class=""> Saransh,</div><div class=""> </div> I've add some code for SNESSetPicard() in the PETSc branch barry/2021-05-06/add-snes-picard-mf see also http<a class="">s://gitlab.com/petsc/petsc/-/merge_requests/3962</a> that will make your coding much easier.<div class=""><br class=""></div><div class=""> With this branch you can provide code that computes A(x), using SNESSetPicard(). </div><div class=""><br class=""></div><div class="">1) by default it will use the defection-correction form of Picard iteration A(x^n)(x^{n+1} - x^{n}) = b - A(x^n) to solve, which can be slower than Newton </div><div class=""><br class=""></div><div class="">2) with -snes_fd_color it will compute the Jacobian via coloring using SNESComputeJacobianDefaultColor() (assuming the true Jacobian has the same sparsity structure as A). The true Jacobian is J(x^n) = A'(x^n)[x^n] - A(x^n) where A'(x^n) is the third order tensor of the derivatives of A() and A'(x^n)[x^n] is a matrix, I do not know if, in general, it has the same nonzero structure as A. (I'm lost beyond matrices :-().</div><div class=""></div><div class=""><br class=""></div><div class="">3) with -snes_mf_operator it will apply the true Jacobian matrix-free and precondition it with a preconditioner built from A(x^n) matrix, for some problems this works well. </div><div class=""><br class=""></div><div class="">4) with -snes_fd it uses SNESComputeJacobianDefault() and computes the Jacobian by finite differencing one column at a time, thus it is very slow and not useful for large problems. But useful for testing with small problems.</div><div class=""><br class=""></div><div class="">So you can provide A() and need not worrying about providing the Jacobian or even the function evaluation code. It is all taken care of by SNESSetPicard().</div><div class=""><br class=""></div><div class=""> Hope this helps,</div><div class=""><br class=""></div><div class=""> Barry</div><div class=""><br class=""></div><div class=""><div class=""><div class=""><br class=""><blockquote type="cite" class=""><div class="">On May 6, 2021, at 1:21 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank" class="">knepley@gmail.com</a>> wrote:</div><br class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class="">On Thu, May 6, 2021 at 2:09 PM Saransh Saxena <<a href="mailto:saransh.saxena5571@gmail.com" target="_blank" class="">saransh.saxena5571@gmail.com</a>> wrote:<br class=""></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr" class="">Hi,<br class=""><br class="">I am trying to incorporate newton method in solving a nonlinear FEM equation using SNES from PETSc. The overall equation is of the type A(x).x = b, where b is a vector of external loads, x is the solution field (say displacements for e.g.) and A is the combined LHS matrix derived from the discretization of weak formulation of the governing finite element equation. <br class=""><br class="">While going through the manual and examples of snes, I found that I need to define the function of residual using SNESSetFunction() and jacobian using SNESSetJacobian(). In that context I had a couple of questions :-<div class=""><br class=""></div><div class="">1. In the snes tutorials I've browsed through, the functions for computing residual passed had arguments only for x, the solution vector and f, the residual vector. Is there a way a user can pass an additional vector (b) and matrix (A) for computing the residual as well? as in my case, f = b - A(x).x<br class=""></div></div></blockquote><div class=""><br class=""></div><div class="">You would give PETSc an outer function MyResidual() that looked like this:</div><div class=""><br class=""></div><div class="">PetscErrorCode MyResidual(SNES snes, Vec X, Vec F, void *ctx)</div><div class="">{</div><div class=""> <call your code to compute b, or pass it in using ctx></div><div class=""> <call your code to compute A(X)></div><div class=""> MatMult(A, X, F);</div><div class=""> VecAXPY(F, -1.0, b);</div><div class="">}</div><div 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 dir="ltr" class=""><div class="">2. Since computing jacobian is not that trivial, I would like to use one of the pre-built jacobian methods. Is there any other step other than setting the 3rd argument in SNESSetJacobian to SNESComputeJacobianDefault?<br class=""></div></div></blockquote><div class=""><br class=""></div><div class="">If you do nothing, we will compute it by default.</div><div class=""><br class=""></div><div class=""> Thanks,</div><div class=""><br class=""></div><div class=""> MAtt</div><div 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 dir="ltr" class=""><div class="">Best regards,<br class=""><br class="">Saransh</div></div>
</blockquote></div><br clear="all" class=""><div class=""><br class=""></div>-- <br class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class="">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">-- Norbert Wiener</div><div class=""><br class=""></div><div class=""><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank" class="">https://www.cse.buffalo.edu/~knepley/</a><br class=""></div></div></div></div></div></div></div></div>
</div></blockquote></div><br class=""></div></div></div></blockquote></div>
</div></blockquote></div><br class=""></div></div></blockquote></div>
</div></blockquote></div><br class=""></div></div></div></blockquote></div>
</div></blockquote></div><br class=""></body></html>