[petsc-users] Integrating SNES in FEM code

[Barry Smith]

  Saransh,
 
  I've add some code for SNESSetPicard() in the PETSc branch barry/2021-05-06/add-snes-picard-mf see also https://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 <mailto: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 experiments lead.
> -- Norbert Wiener
> 
> https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210507/361dbe17/attachment.html>