# [petsc-users] Integrating SNES in FEM code

Matthew Knepley knepley at gmail.com
Tue May 25 10:50:19 CDT 2021

```On Tue, May 25, 2021 at 11:29 AM Saransh Saxena <
saransh.saxena5571 at gmail.com> wrote:

> Hi Barry,
>
>  Mat Apetsc = A.getpetsc();
>  Vec bpetsc = b.getpetsc();
>
> 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.
>
> Vec output;
> VecDuplicate(x, &output);
> VecCopy(x, output);
>
> setdata(vec(b.getpointer()->getdofmanager(), output));
>
> 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.
>
> I was also browsing through the definition of SNESSetFunction and realized
> that it solves f'(x) x = -f(x), 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?
>

No. SNESSetFunction() provides the residual F, so that

F(x) = 0

You can use _many_ different algorithms to solve this system. One is Newton.

> Also in SNESSetPicard(), I need to pass a function to compute b. However,
> in my case b is constant. How do I use that? Also does Vec r in the
> definition refer to solution vector or residual vector?
>

Just keep returning that vector.

Thanks,

Matt

> Best regards,
> Saransh
>
>
>
>
>
> On Tue, May 25, 2021 at 10:15 AM Barry Smith <bsmith at petsc.dev> wrote:
>
>>
>>     VecNorm(F, NORM_2, &normvalres);
>>
>>
>>    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
>>
>> VecAXPY(F,-1.0, bpetsc);
>>
>>
>>
>>
>>     // Read pointers to A and b:
>>     Mat Apetsc = A.getpetsc();
>>     Vec bpetsc = b.getpetsc();
>>
>>
>>    Where are these things computed and are they both functions of output?
>>  or is b merely x (the current solution snes is working with)
>>
>> Vec output;
>>     VecDuplicate(x, &output);
>>     VecCopy(x, output);
>>
>>     setdata(vec(b.getpointer()->getdofmanager(), output));
>>
>>
>>    What is the line above doing?
>>
>>    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().
>>
>>    If you run with SNESSetPicard()  with no additional options it will
>> run the defect correction version of Picard
>>
>>    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
>>
>>    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.
>>
>>    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.
>>
>>   Barry
>>
>>
>> On May 25, 2021, at 2:48 AM, Saransh Saxena <saransh.saxena5571 at gmail.com>
>> wrote:
>>
>> Hi guys,
>>
>> 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.
>>
>> Residual code :-
>>
>> PetscErrorCode sl::myresidual(SNES snes, Vec x, Vec F, void *ctx)
>> {
>>     // Cast the application context:
>>     sl::formulCtx *user = (sl::formulCtx*)ctx;
>>
>>     formulation *thisformul = (*user).formul;
>>     thisformul->generate();
>>
>>     //vec *b = user->b;
>>     //mat *A = user->A;
>>
>>     vec b = thisformul->b();
>>     mat A = thisformul->A();
>>
>>     // Read pointers to A and b:
>>     Mat Apetsc = A.getpetsc();
>>     Vec bpetsc = b.getpetsc();
>>
>>     double normvalres, normvalsol;
>>     VecNorm(F, NORM_2, &normvalres);
>>     VecNorm(x, NORM_2, &normvalsol);
>>     std::cout <<
>> "----------------------------------------------------------------------------"
>> << std::endl;
>>     std::cout << "Entered residual function, norm of residual vector is :
>> " << normvalres << std::endl;
>>     std::cout << "Entered residual function, norm of solution vector is :
>> " << normvalsol << std::endl;
>>
>>     // Compute the residual as F = A*x - b
>>     MatMult(Apetsc, x, F);
>>     VecAXPY(F,-1.0, bpetsc);
>>
>>     Vec output;
>>     VecDuplicate(x, &output);
>>     VecCopy(x, output);
>>
>>     setdata(vec(b.getpointer()->getdofmanager(), output));
>>
>>     std::cout << "Writing the sol to fields \n";
>>
>>     return 0;
>> }
>>
>> SNES implementation :-
>>
>> void sl::solvenonlinear(formulation thisformul, double restol, int
>> maxitnum)
>> {
>>     // Make sure the problem is of the form Ax = b:
>>     if (thisformul.isdampingmatrixdefined() ||
>> thisformul.ismassmatrixdefined())
>>     {
>>         std::cout << "Error in 'sl' namespace: formulation to solve
>> cannot have a damping/mass matrix (use a time resolution algorithm)" <<
>> std::endl;
>>         abort();
>>     }
>>
>>     // Remove leftovers (if any):
>>     mat Atemp = thisformul.A(); vec btemp = thisformul.b();
>>
>>     // Create Application Context for formulation
>>     sl::formulCtx user;
>>     user.formul = &thisformul;
>>
>>     // Generate formulation to set PETSc SNES requirements:
>>     thisformul.generate();
>>
>>     mat A = thisformul.A();
>>     vec b = thisformul.b();
>>
>>     // SNES requirements:
>>     Vec bpetsc = b.getpetsc();
>>     Mat Apetsc = A.getpetsc();
>>
>>     vec residual(std::shared_ptr<rawvec>(new
>> rawvec(b.getpointer()->getdofmanager())));
>>     Vec residualpetsc = residual.getpetsc();
>>
>>     vec sol(std::shared_ptr<rawvec>(new
>> rawvec(b.getpointer()->getdofmanager())));
>>     Vec solpetsc = sol.getpetsc();
>>
>>     //Retrieve the SNES and KSP Context from A matrix:
>>     SNES* snes = A.getpointer()->getsnes();
>>     KSP* ksp = A.getpointer()->getksp();
>>
>>     // Create placeholder for preconditioner:
>>     PC pc;
>>
>>     // Create snes context:
>>     SNESCreate(PETSC_COMM_SELF, snes);
>>     SNESSetFunction(*snes, residualpetsc, sl::myresidual, &user);
>>     SNESSetTolerances(*snes, PETSC_DEFAULT, restol, PETSC_DEFAULT,
>> maxitnum, 5);
>>
>>     // Retrieve the KSP context automatically created:
>>     SNESGetKSP(*snes, ksp);
>>
>>     //Set KSP specific parameters/options:
>>     KSPSetOperators(*ksp, Apetsc, Apetsc);
>>     KSPSetFromOptions(*ksp);
>>     KSPGetPC(*ksp,&pc);
>>     PCSetType(pc,PCLU);
>>     PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
>>
>>     //Call SNES options to invoke changes from console:
>>     SNESSetFromOptions(*snes);
>>
>>     // Set SNES Monitor to retrieve convergence information:
>>     SNESMonitorSet(*snes, sl::mysnesmonitor, PETSC_NULL, PETSC_NULL);
>>     //SNESMonitorLGResidualNorm();
>>
>>     SNESSolve(*snes, PETSC_NULL, solpetsc);
>>
>>     // Print the norm of residual:
>>     double normres;
>>     VecNorm(residualpetsc, NORM_2, &normres);
>>     std::cout << "L2 norm of the residual is : " << normres << std::endl;
>>
>>     //Set the solution to all the fields:
>>     setdata(sol);
>>
>>     // Get the number of required iterations and the residual norm:
>>     //SNESGetIterationNumber(*snes, &maxitnum);
>>     //SNESGetResidualNorm(*snes, &restol);
>>
>>     // Destroy SNES context once done with computation:
>>     SNESDestroy(snes);
>>
>> }
>>
>> Output :-
>> <image.png>
>>
>> Am I doing something incorrect wrt SNES? When I use the linear solver
>> (KSP) and manually coded fixed point nonlinear iteration, it works fine.
>>
>> Best regards,
>> Saransh
>>
>>
>>
>> On Sun, May 9, 2021 at 10:43 PM Barry Smith <bsmith at petsc.dev> wrote:
>>
>>>
>>>   Saransh,
>>>
>>>      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.
>>>
>>>   Barry
>>>
>>>
>>> On May 9, 2021, at 2:07 PM, Saransh Saxena <saransh.saxena5571 at gmail.com>
>>> wrote:
>>>
>>> Thanks Barry and Matt,
>>>
>>> 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.
>>>
>>> Cheers,
>>> Saransh
>>>
>>> On Sat, May 8, 2021 at 5:18 AM Barry Smith <bsmith at petsc.dev> wrote:
>>>
>>>>   Saransh,
>>>>
>>>>   I've add some code for SNESSetPicard() in the PETSc branch
>>>> s://gitlab.com/petsc/petsc/-/merge_requests/3962 that will make your
>>>> coding much easier.
>>>>
>>>>    With this branch you can provide code that computes A(x), using
>>>> SNESSetPicard().
>>>>
>>>> 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
>>>>
>>>> 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 :-().
>>>>
>>>> 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.
>>>>
>>>> 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.
>>>>
>>>> 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().
>>>>
>>>>   Hope this helps,
>>>>
>>>>   Barry
>>>>
>>>>
>>>> On May 6, 2021, at 1:21 PM, Matthew Knepley <knepley at gmail.com> wrote:
>>>>
>>>> On Thu, May 6, 2021 at 2:09 PM Saransh Saxena <
>>>> saransh.saxena5571 at gmail.com> wrote:
>>>>
>>>>> Hi,
>>>>>
>>>>> 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.
>>>>>
>>>>> 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 :-
>>>>>
>>>>> 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
>>>>>
>>>>
>>>> You would give PETSc an outer function MyResidual() that looked like
>>>> this:
>>>>
>>>> PetscErrorCode MyResidual(SNES snes, Vec X, Vec F, void *ctx)
>>>> {
>>>>   <call your code to compute b, or pass it in using ctx>
>>>>   <call your code to compute A(X)>
>>>>   MatMult(A, X, F);
>>>>   VecAXPY(F, -1.0, b);
>>>> }
>>>>
>>>>
>>>>> 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?
>>>>>
>>>>
>>>> If you do nothing, we will compute it by default.
>>>>
>>>>   Thanks,
>>>>
>>>>     MAtt
>>>>
>>>>
>>>>> Best regards,
>>>>>
>>>>> Saransh
>>>>>
>>>>
>>>>
>>>> --
>>>> What most experimenters take for granted before they begin their
>>>> experiments is infinitely more interesting than any results to which their
>>>> -- Norbert Wiener
>>>>
>>>> https://www.cse.buffalo.edu/~knepley/
>>>> <http://www.cse.buffalo.edu/~knepley/>
>>>>
>>>>
>>>>
>>>
>>

--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their