I am solving 3D time-harmonic Maxwell FDFD problems (which result in huge sparse linear systems) using BiCGStab(l). I have tried out a bunch of different methods and for my specific use case, it seems like this algorithm outperforms the other common ones such as BiCG and restarted GMRES in terms of total solution time. However, I have noticed that as I have increased my problem size as well as increased the meshing accuracy (I am using a non-uniform grid, but the smallest $dx$ seems to be important for the convergence behavior) that I am running into an issue where often after many iterations for a specific initial solution, BiCGStab(l) stagnates and cannot reduce the residual further. This seems to happen about 10% of the time and so far I have always been able to solve it by restarting from a different initial solution. Each problem takes a long time to solve, however, so it would be ideal if I could be guaranteed convergence on the first go. I have tried setting l=2,4,5, and even 8 and the same thing happens.
I am wondering if doing something simple like detecting stagnation, taking the current solution and randomly perturbing it a little bit, followed by restarting the solver with that as the initial guess would help. Sounds like a complete hack, but I'm not sure how else I can resolve this without having to completely restart.
In case it is important for anything, I am using a simple preconditioner which is pretty much a diagonal preconditioner. I found that for my specific problems it leads to faster convergence than an ILU(0) preconditioner. I have not tried any other preconditioners than these two.