Matrixfree modification

[Johann Schlamp]

Barry Smith wrote:
>    If you have the "matrix-based" version that you can run on the same
> problem then
> look at the residuals computed in the 0th and 1st step of the Krylov
> method
> to see how they are different in the two runs.
The residuals in the 0th and 1st step of the linear solver are
0.00363413
0.00189276
for the "matrix-based" version, and
0.00363413
1.27858e-17
for the matrix-free version. That's definitely smaller than epsilon, so
it converges. By the way, the "matrix-based" version doesn't converge
either, as I was not using a preconditioner for getting comparable results.

Simply thought: the residual is in the magnitude of machine accuracy, so
I would have concluded that the calculation of the residual (y-A*x)
results in zero with respect to some rounding errors. Unfortunately, I
don't completely understand the PETSc code for calculating the residual
and therfore cannot verify it for my new matrix structure.

>    Perhaps your matrix-free is corrupting memory? Run with -malloc_debug
> and put a CHKMEMQ; at the end of your matrix free multiply. Or better
> run through
> valgrind,. www.valgrind.org
Interesting thought! I will check this tomorrow.


Johann


> On Oct 28, 2008, at 3:39 PM, Johann Schlamp wrote:
>
>> 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
>>>>
>>>>
>>>
>>
>