> On May 25, 2021, at 10:28 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?
>
> Also in SNESSetPicard(), I need to pass a function to compute b. However, in
> my case b is constant. How do I use that?
You have two choices.
* As Matt says you can just provide a function that copies over your
constant b vector each time.
* Or you can pass a NULL for the function and call SNESSolve(snes,b,x)
where b is your constant b vector (this will be slightly more efficient)
> Also does Vec r in the definition refer to solution vector or residual vector?
r is a vector that SNES will use to compute the residual in. It is just a
work vector you can pass in. You can pass in NULL and PETSc will create the
work vector it needs internally.
>
> Best regards,
> Saransh
>
>
>
>
>
> On Tue, May 25, 2021 at 10:15 AM Barry Smith <bsm...@petsc.dev
> <mailto: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
>> <mailto: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
>> <mailto: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
>>> <mailto: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
>>> <mailto: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
>>> 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 <knep...@gmail.com
>>>> <mailto:knep...@gmail.com>> wrote:
>>>>
>>>> On Thu, May 6, 2021 at 2:09 PM Saransh Saxena
>>>> <saransh.saxena5...@gmail.com <mailto: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/>
>>>
>>
>
