Johann Schlamp wrote: > 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 >> >>
After some extended testing, I think I really messed up the input vector somehow. I dare say I can take it from here. If not - I know where to get competent support. :) Thanks for all help! Johann > 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 >>>>>>> >>>>>>> >>>>>>> >>>>>>> >>> >>> >> >> >> > >
