Just some feedback. I found the problem. For reference my solve was called as follows
KSPSolve(ksp,b,phi_new) Inside my matrix operation (the "Matrix-Action" or MAT_OP_MULT) I was using phi_new for a computation and that overwrite my initial guess everytime. Looks like the solver still holds on to phi_new and uses it internally and therefore when I modify it it basically changes the entire behavior. Lesson learned: Internal to the custom MAT_OP_MULT, do not modify, b or phi_new. Thanks for the help. Regards, Jan Vermaak On Tue, May 28, 2019 at 1:16 AM Smith, Barry F. <[email protected]> wrote: > > > > On May 28, 2019, at 12:41 AM, Jan Izak Cornelius Vermaak < > [email protected]> wrote: > > > > Checking the matrix would be hard as I have a really big operator. Its > transport sweeps. > > > > When I increase the restart interval the solution converges to the right > one. > > Run with -ksp_monitor_true_residual what are the true residuals being > printed? > > The GMRES code has been in continuous use for 25 years, it would > stunning if you suddenly found a bug in it. > > How it works, within a restart, the GMRES algorithm uses a simple > recursive formula to compute an "estimate" for the residual norm. At > restart it actually computes the current solution and then uses that to > compute an accurate residual norm via the formula b - A x. When the > residual computed by the b - A x is different than that computed by the > recursive formula it means the recursive formula has run into some > difficulty (bad operator, bad preconditioner, null space in the operator) > and is not computing correct values. Now if you increase the restart to > past the point when it "converges" you are hiding the incorrectly computed > values computed via the recursive formula. > > I urge you to check the residual norm by b - A x at the end of the solve > and double check that it is small. It seems unlikely GMRES is providing the > correct answer for your problem. > > Barry > > > > Checked against a reference solution and Classic Richardson. It is > really as if the initial guess is completely ignored. > > > > [0] Computing b > > [0] Iteration 0 Residual 169.302 > > [0] Iteration 1 Residual 47.582 > > [0] Iteration 2 Residual 13.2614 > > [0] Iteration 3 Residual 4.46795 > > [0] Iteration 4 Residual 1.03038 > > [0] Iteration 5 Residual 0.246807 > > [0] Iteration 6 Residual 0.0828341 > > [0] Iteration 7 Residual 0.0410627 > > [0] Iteration 8 Residual 0.0243749 > > [0] Iteration 9 Residual 0.0136067 > > [0] Iteration 10 Residual 0.00769078 > > [0] Iteration 11 Residual 0.00441658 > > [0] Iteration 12 Residual 0.00240794 > > [0] Iteration 13 Residual 0.00132048 > > [0] Iteration 14 Residual 0.00073003 > > [0] Iteration 15 Residual 0.000399504 > > [0] Iteration 16 Residual 0.000217677 > > [0] Iteration 17 Residual 0.000120408 > > [0] Iteration 18 Residual 6.49719e-05 > > [0] Iteration 19 Residual 3.44523e-05 > > [0] Iteration 20 Residual 1.87909e-05 > > [0] Iteration 21 Residual 1.02385e-05 > > [0] Iteration 22 Residual 5.57859e-06 > > [0] Iteration 23 Residual 3.03431e-06 > > [0] Iteration 24 Residual 1.63696e-06 > > [0] Iteration 25 Residual 8.78202e-07 > > > > On Mon, May 27, 2019 at 11:55 PM Smith, Barry F. <[email protected]> > wrote: > > > > This behavior where the residual norm jumps at restart indicates > something is very very wrong. Run with the option > -ksp_monitor_true_residual and I think you'll see the true residual is not > decreasing as is the preconditioned residual. My guess is that your "action > of the matrix" is incorrect and not actually a linear operator. Try using > MatComputeExplicitOperator() and see what explicit matrix it produces, is > it what you expect? > > > > Barry > > > > > > > > > > > On May 27, 2019, at 11:33 PM, Jan Izak Cornelius Vermaak via > petsc-users <[email protected]> wrote: > > > > > > Hi all, > > > > > > So I am faced with this debacle. I have a neutron transport solver > with a sweep code that can compute the action of the matrix on a vector. > > > > > > I use a matrix shell to set up the action of the matrix. The method > works but only if I can get the solution converged before GMRES restarts. > It gets the right answer. Now my first problem is (and I only saw this when > I hit the first restart) is that it looks like the solver completely resets > after the GMRES-restart. Below is an iteration log with restart interval > set to 10. At first I thought it wasn't updating the initial guess but it > became clear that it initial guess always had no effect. I do set > KSPSetInitialGuessNonZero but it has no effect. > > > > > > Is the matrix-free business defaulting my initial guess to zero > everytime? What can I do to actually supply an initial guess? I've used > PETSc for diffusion many times and the initial guess always works, just not > now. > > > > > > [0] Computing b > > > [0] Iteration 0 Residual 169.302 > > > [0] Iteration 1 Residual 47.582 > > > [0] Iteration 2 Residual 13.2614 > > > [0] Iteration 3 Residual 4.46795 > > > [0] Iteration 4 Residual 1.03038 > > > [0] Iteration 5 Residual 0.246807 > > > [0] Iteration 6 Residual 0.0828341 > > > [0] Iteration 7 Residual 0.0410627 > > > [0] Iteration 8 Residual 0.0243749 > > > [0] Iteration 9 Residual 0.0136067 > > > [0] Iteration 10 Residual 169.302 > > > [0] Iteration 11 Residual 47.582 > > > [0] Iteration 12 Residual 13.2614 > > > [0] Iteration 13 Residual 4.46795 > > > [0] Iteration 14 Residual 1.03038 > > > [0] Iteration 15 Residual 0.246807 > > > [0] Iteration 16 Residual 0.0828341 > > > [0] Iteration 17 Residual 0.0410627 > > > [0] Iteration 18 Residual 0.0243749 > > > [0] Iteration 19 Residual 0.0136067 > > > [0] Iteration 20 Residual 169.302 > > > > > > -- > > > Jan Izak Cornelius Vermaak > > > (M.Eng Nuclear) > > > Email: [email protected] > > > Cell: +1-979-739-0789 > > > > > > > > -- > > Jan Izak Cornelius Vermaak > > (M.Eng Nuclear) > > Email: [email protected] > > Cell: +1-979-739-0789 > > -- Jan Izak Cornelius Vermaak (M.Eng Nuclear) Email: [email protected] Cell: +1-979-739-0789
