[petsc-users] Integrating SNES in FEM code

Tue May 25 02:48:02 CDT 2021

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;

    // Read the formulation:
    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: 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
>> 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 <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
>> 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/20210525/446d17d0/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image.png
Type: image/png
Size: 160332 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210525/446d17d0/attachment-0001.png>