[petsc-users] advice about preconditioning

[Smith, Barry F.]

> On Aug 31, 2018, at 8:31 PM, Oleksandr Koshkarov <olk548 at mail.usask.ca> wrote:
> 
> Dear All,
> 
> I hope you can help me with some advice on how to proceed with constructing a preconditioner for my petsc based simulation.
> 
> I have a large nonlinear system dU/dt = L(U) + N(U,U), with U being a state vector, L linear operator and N nonlinear one. It comes from discretizing hyperbolic PDE with finite volume. L is advection and N is some nonlinear source which does not have spatial derivatives. My U is PETSc 3D DMDA with dof = from 1e2 up to 2e3. I solve it with petsc TS object and with Crank Nicolson method. In order to save space, I do use finite difference Jacobian, which petsc generates itself. Something like this
> 
>     MatCreateSNESMF(snes,&J);
>     SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,0);
> 
> Now I want to precondition my linear solves inside Newton iterations with linearized Jacobian, i.e, with something like this
> 
>  P_n = 1 - dt/2 * ( L + N(U^n,*) + N(*,U^n) )        // where U^n is state from previous time step
> 
> I first checked, that this works on small test problem. I used matrix free preconditioner (PCSHELL) inside which I inverted P_n with KSP.
> 
> Where  P_n is matrix free matrix (shell matrix). It reduces the number of KSP dramatically, which is not a surprise - it was just a sanity check.
> 
> Now I am trying to understand how to proceed... and I am not sure what to do. In my head I have two options:
> 
> 1) Try to construct matrix-free preconditioner. For example, I can keep inverting P_n every iteration inside the preconditioner with KSP (or some other method) but now precondition it with something else, to make it cheap. So I would have double preconditioner...  I am not sure if this is a right direction, and it is even possible to make effective matrix-free preconditioner in my case
> 
> 2) Construct P_n as real matrix, and feed it to petsc (SNESSetJacobian(snes,J,P,MatMFFDComputeJacobian,0);) and ask to perform something like LU for preconditioning. It is very expensive (as I have huge system) and I would need to probably contract P_n only every 100 steps or something like this... It is probably the right direction I have to go, but it seems it will take too much memory.
> 
> Is there other better way? :)
> 
> Small question. If I will need to construct P_n element by element, does petsc provide some help with DMDA indexing. I usually access state U with x-index (xi), y-index (yi), z-index(zi), and dof like u[zi][yi][xi][dof]. Can I use similar indexing to access matrix elements?

   MatSetValuesStencil()


   I would construct the P matrix and use LU to start, make sure you have good Newton convergence then after everything is working well is time to investigate preconditioners for P.

     Barry

> 
> 
> Thank you and best regards,
> 
> Oleksandr.
>