[petsc-users] advice about preconditioning

[Oleksandr Koshkarov]

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?


Thank you and best regards,

Oleksandr.