Bug with initial guess for Krylov solvers

[prisae]

This is a bug which has always been an issue, in one or another way. The Krylov methods in its current implementation have a problem with initial guesses, either they do not converge, or the field is not updated in place. I have currently no idea whatsoever how to fix this.

No, the bug was fixed in #32 by simply resetting initial guesses to zero if an sslsolver is used.

However, I leave this issue open as an 'Enhancement'-request to re-enable that feature again.

Original initial message:

There is a bug in solver.krylov for very rare cases. I thought I fixed it with
d255e2d
but that was wrong (overwriting a potentially provided efield). So that commit was undone again in
f34dcb6

[prisae]

import emg3d
import discretize
import numpy as np
freqs = np.array([1, 2])
skind = 503./np.sqrt(3*freqs)
xx = emg3d.utils.get_stretched_h(h_min=160, domain=np.array([-5000, 5000]), nx=2**5)

for fi, frq in enumerate(freqs):

    xmin = -4000 - 5*skind[fi]

    grid = discretize.TensorMesh([xx, xx, xx], x0=np.array([xmin, xmin, xmin]))
    model = emg3d.utils.Model(grid, res_x=1., freq=frq)
    sfield = emg3d.utils.get_source_field(grid, [0, 0, 0, 0, 0], frq)

    if fi== 0:
        efield = emg3d.utils.Field(grid)
    else:
        efield = emg3d.utils.grid2grid(old_grid, efield, grid, method='cubic', extrapolate=False)
        
    emg3d.solver.solver(
        grid, model, sfield, efield=efield, sslsolver=True, semicoarsening=True,
        linerelaxation=True, verb=3, maxit=3, tol=1e-15)

    old_grid = grid

:: emg3d START :: 15:06:47 ::

   MG-cycle       : 'F'                 sslsolver : 'bicgstab'
   semicoarsening : True [1 2 3]        tol       : 1e-15
   linerelaxation : True [4 5 6]        maxit     : 3 (3)
   nu_{i,1,c,2}   : 0, 2, 1, 2          verb      : 3
   Original grid  :  32 x  32 x  32     => 32,768 cells
   Coarsest grid  :   2 x   2 x   2     => 8 cells
   Coarsest level :   4 ;   4 ;   4   

   [hh:mm:ss]  rel. error            solver              MG          l s

       h_
      2h_ \                  /
      4h_  \          /\    / 
      8h_   \    /\  /  \  /  
     16h_    \/\/  \/    \/   

   [15:06:48]   2.817e-03  after                       1 F-cycles    4 1
   [15:06:49]   5.250e-05  after                       2 F-cycles    5 2
   [15:06:49]   6.931e-07  after                       3 F-cycles    6 3
   [15:06:50]   1.158e-08  after                       4 F-cycles    4 1
   [15:06:51]   2.283e-10  after                       5 F-cycles    5 2
   [15:06:51]   6.743e-12  after                       6 F-cycles    6 3
   [15:06:51]   6.201e-12  after   1 bicgstab-cycles
   [15:06:52]   1.404e-13  after                       7 F-cycles    4 1
   [15:06:52]   7.189e-15  after                       8 F-cycles    5 2
   [15:06:53]   4.983e-16  after                       9 F-cycles    6 3
   [15:06:53]   1.306e-15  after   2 bicgstab-cycles

   > CONVERGED
   > Solver steps     : 2
   > MG prec. steps   : 9
   > Final rel. error : 1.306e-15

:: emg3d END   :: 15:06:53 :: runtime = 0:00:05


:: emg3d START :: 15:06:53 ::

   MG-cycle       : 'F'                 sslsolver : 'bicgstab'
   semicoarsening : True [1 2 3]        tol       : 1e-15
   linerelaxation : True [4 5 6]        maxit     : 3 (3)
   nu_{i,1,c,2}   : 0, 2, 1, 2          verb      : 3
   Original grid  :  32 x  32 x  32     => 32,768 cells
   Coarsest grid  :   2 x   2 x   2     => 8 cells
   Coarsest level :   4 ;   4 ;   4   

   [hh:mm:ss]  rel. error            solver              MG          l s

       h_
      2h_ \                  /
      4h_  \          /\    / 
      8h_   \    /\  /  \  /  
     16h_    \/\/  \/    \/   

   [15:06:54]   5.018e-04  after                       1 F-cycles    4 1
   [15:06:54]   6.633e-06  after                       2 F-cycles    5 2
   [15:06:55]   3.520e-06  after                       3 F-cycles    6 3
   [15:06:56]   3.519e-06  after                       4 F-cycles    4 1
   [15:06:56]   3.519e-06  after                       5 F-cycles    5 2
   [15:06:57]   3.519e-06  after                       6 F-cycles    6 3
   [15:06:57]   3.519e-06  after   1 bicgstab-cycles
   [15:06:57]   3.519e-06  after                       7 F-cycles    4 1
   [15:06:58]   3.519e-06  after                       8 F-cycles    5 2
   [15:06:58]   3.519e-06  after                       9 F-cycles    6 3
   [15:06:59]   3.519e-06  after                      10 F-cycles    4 1
   [15:06:59]   3.519e-06  after                      11 F-cycles    5 2
   [15:07:00]   3.519e-06  after                      12 F-cycles    6 3
   [15:07:00]   3.519e-06  after   2 bicgstab-cycles
   [15:07:00]   2.711e-02  after                      13 F-cycles    4 1
   [15:07:01]   2.711e-02  after                      14 F-cycles    5 2
   [15:07:01]   2.711e-02  after                      15 F-cycles    6 3
   [15:07:02]   3.519e-06  after                      16 F-cycles    4 1
   [15:07:02]   3.519e-06  after                      17 F-cycles    5 2
   [15:07:03]   3.519e-06  after                      18 F-cycles    6 3
   [15:07:03]   3.519e-06  after   3 bicgstab-cycles

   > MAX. ITERATION REACHED, NOT CONVERGED
   > Solver steps     : 3
   > MG prec. steps   : 18
   > Final rel. error : 3.519e-06

:: emg3d END   :: 15:07:03 :: runtime = 0:00:10

Solver gets stuck, sometimes the error after BICGSTAB becomes massive (1e+28 and more), is corrected afterwards again, just to jump up again. And, the l2-norm should never be EXACTLY the same...

Any ideas?

[prisae]

So far it is only confirmed to affect ssl.bicgstab with an initial field provided, in certain cases. No idea what is causing it.

[prisae]

it could be resolved by always returning the efield.

[prisae]

A few important points for the future:

• This has always been an issue, in one or another way.
• The Krylov methods in its current implementation have a problem with initial guesses.
• Either they do not converge, or the field is not updated in place.
• No idea whatsoever currently how to fix this.
• BUG: Easy-fix for initial guess with sslsolver (#31) #32 is an easy fix for the moment: an initial guess is reset to zero if sslsolver is not False. This does not resolve the problem, but it avoids stagnation/divergence issues. So no initial guesses possible for sslsolvers. It might still be worth providing an initial guess. If it is good enough it will return immediately, and the result is still put in place.