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.
