On Tue, May 25, 2021 at 11:29 AM Saransh Saxena <
saransh.saxena5...@gmail.com> wrote:
> 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?
>
No. SNESSetFunction() provides the residual F, so that
F(x) = 0
You can use _many_ different algorithms to solve this system. One is Newton.
> 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?
>
Just keep returning that vector.
Thanks,
Matt
> Best regards,
> Saransh
>
>
>
>
>
> On Tue, May 25, 2021 at 10:15 AM Barry Smith <bsm...@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.saxena5...@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 <bsm...@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.saxena5...@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 <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/>
>>>>
>>>>
>>>>
>>>
>>
--
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/>
