[petsc-users] matrix-free SNES with preconditioner

[Oleksandr Koshkarov]

Hello all,

I got lost a little with petsc and want to clarify something.

I want to do my nonlinear solves with matrix free method, but I want a 
real matrix preconditioner and invert it with LU.

without preconditioner it is simple:

Mat Jmf;

MatCreateSNESMF(snes,&Jmf);

SNESSetJacobian(snes,Jmf,Jmf,MatMFFDComputeJacobian,0);


But with preconditioner, I am lost:

Mat Jmf,J;

MatCreateSNESMF(snes,&Jmf);

DMCreateMatrix(da,&J);

SNESSetJacobian(snes,Jmf,J,FromJacobian,0);

//And here I have problems:

PetscErrorCode FormJacobian(SNES snes,Vec X,Mat Jmf,Mat J,void *ptr){

// I am constructing it only once in a while, not every nonlinear solve 
(this snes is inside time stepper).

if (sometimes) {  MatSetValuesStencil(J,...); ...    }

   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

// assembling Jmf
   MatAssemblyBegin(Jmf,MAT_FINAL_ASSEMBLY);
   MatAssemblyEnd(Jmf,MAT_FINAL_ASSEMBLY);

   return(0);
}

Here I have two concerns:

1) I am changing J only once in a while, because it is expansive to do 
otherwise. Am I doing it correctly? During iteration when I am not 
changing J, will petsc know this? I do not want petsc to perform LU 
every time when I am not changing J. And when I will change J, will it 
work? (My snes is inside time stepper, and I want to construct 
preconditioner only say every 100 time steps)

2) What about Jmf? how can I let petsc know that it should use matrix 
free finite difference? Will it work if I let it as it is? (I think so, 
but I am not sure - my information source is 
ts/examples/tutorials/ex15.c, but I do not understand it fully)

Thank you and best regards,

Alex.