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 >>>> >>>> >>> >> >
