# [petsc-users] Integrating SNES in FEM code

Matthew Knepley knepley at gmail.com
Thu May 6 13:21:57 CDT 2021

```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