Matthew Knepley wrote: > On Tue, Oct 28, 2008 at 4:35 PM, Johann Schlamp <schlamp at in.tum.de> wrote: > >> 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 >> > > No this looks wrong. Shouldn't these be identical? It looks like you > are wiping out the input vector instead. > > Matt > Yes, they should be identical. That's excactly the point. My naive interpretation was that maybe only the residual's calculation is wrong.
After thinking again, I believe I have misunderstood Barry's hint on copying the 'u'. Apparently, the user-provided Jacobian calculation method gets a different input than my customized MatMult method called within KSP on my matrixfree context. But they are both needed for my approach, right? I haven't thought of that in advance, so it will take me some time to rewrite the code. I will report again tomorrow (it's 11 o'clock pm at my site). Thanks for your help! Johann >> 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 >>>>>> >>>>>> >>>>>> >> > > > >
