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 <bsm...@petsc.dev> wrote:
> Saransh,
>
> I've add some code for SNESSetPicard() in the PETSc branch
> barry/2021-05-06/add-snes-picard-mf see also http
> 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 <knep...@gmail.com> wrote:
>
> On Thu, May 6, 2021 at 2:09 PM Saransh Saxena <
> saransh.saxena5...@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/>
>
>
>
