[petsc-users] Integrating SNES in FEM code
Saransh Saxena
saransh.saxena5571 at gmail.com
Tue May 25 10:28:47 CDT 2021
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?
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?
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;
>
> // 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.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/12d921f1/attachment-0001.html>
More information about the petsc-users
mailing list