Matrixfree modification

[Johann Schlamp <schlamp@in.tum.de>]

Hello folks,

I have to implement some advanced matrixfree calculation.

Usually, we use our own analytical method for calculating the Jacobian
matrix and provide it via SNESSetJacobian(). For complex problems, the
global matrix takes too much memory. So here's the idea:

First, I set some empty method as user-provided Jacobian method. After a
corresponding call, SNES hopefully thinks the Jacobian matrix got
calculated right. Then it will use it in some multiplication like y=A*x.
For that I created the matrix A with MatCreateShell() and set up my own
MatMult method. This new method iterates over my grid, calculates local
Jacobians, multiplies them with the corresponding part of the vector x
and writes them into y. After this full iteration, it should look like I
had the full global Jacobian and multiplied it with the full vector x.
In y, the result will be the same.
The matrix A and the empty Jacobian method are set through
SNESSetJacobian().

I have implemented some unittests for generally proofing the idea, seems
to work.

But if I run a complete simulation, the KSP converges after the first
iteration with a residual norm of about 1.0e-18, which is definitely not
right.


Now my question: does this procedure make sense at all, and if so - is
it possible that just the calculation of the residual norm goes wrong
due to my new matrix structure? I searched the PETSc code, but I wasn't
able to find a proof or solution for that.

Any help would be appreciated.


Best regards,
Johann Schlamp