Matrixfree modification

[Johann Schlamp]

Thanks for your reply, Barry!

Barry Smith wrote:
>    Your general approach seems fine. I would put a break point in
> your MatMult routine for a tiny problem and verify that the input
> vector for u and
> x are what you expect (and then at the end of the function make sure the
> output y is what you expect).
I have already done this, everything is as expected.

>    Here is my guess: The matrix vector product is J(u)*x; when your
> calculation
> is done you need to use the correct u value. This is the vector that is
> passed into your "empty method as user-provided Jacobian method".
> In your "empty method as user-provided Jacobian method" you should make
> a copy of this vector so that you have it for each of the matrix
> vector products.
> At each Newton step your "empty method as user-provided Jacobian method"
> will be called and you will copy the new value of u over.
It took me some time to understand what you mean. But after that I got
excited trying it out. :)
It definitely had some effect on the nonlinear solution, but the linear
solver still finishes after one iteration with the way too small
residual norm.

Anyway, thanks for the anticipated bugfix!

Do you have any further suggestions?


Best regards,
Johann


> On Oct 28, 2008, at 10:08 AM, Johann Schlamp <schlamp at in.tum.de> wrote:
>
>> 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
>>
>>
>