[petsc-users] Matrix free GMRES seems to ignore my initial guess

[Smith, Barry F.]

> On May 28, 2019, at 12:41 AM, Jan Izak Cornelius Vermaak <janicvermaak at gmail.com> 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. <bsmith at mcs.anl.gov> 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 <petsc-users at mcs.anl.gov> 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: janicvermaak at gmail.com
> > Cell:    +1-979-739-0789
> 
> 
> 
> -- 
> Jan Izak Cornelius Vermaak 
> (M.Eng Nuclear)
> Email: janicvermaak at gmail.com
> Cell:    +1-979-739-0789