7
$\begingroup$

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.

Improve this question
$\endgroup$
4
  • $\begingroup$ Out of curiosity, what is your nullspace looking like? The curl-curl Maxwell formulation has a very large nullspace, and sometimes very good approximations will faithfully represent this nullspace - and thus lead to a linear system which is singular. $\endgroup$
    – Reid.Atcheson
    Aug 15, 2013 at 6:41
  • $\begingroup$ I'm not sure what the nullspace looks like because the system is way too prohibitively large to actually compute the nullspace. I am discretizing quite finely and I assume the condition number is really bad (also can't really estimate this well due to the system size). The effective system matrix is size 23,360,400x23,360,400 , although I never explicitly formulate it as I perform the differencing operations directly on vectors. Are you suggesting that a full E/H formulation (rather than eliminating the H components) might have better convergence behavior for a system of this size? $\endgroup$
    – Costis
    Aug 15, 2013 at 8:04
  • $\begingroup$ It might, but without knowing the nullspace it's too hard to say for sure. The mixed formulation (E/H formulation as you put it) effectively deals with the nullspace issue, but without knowing if that is the real cause of the problem I can't say with confidence that it will help your situation. Unfortuntately though solving time harmonic wave propagation is a notoriously difficult problem, and there is no easy solution to your problem in the O(wavenumber) limits. Maxwell's equations make this situation even more difficult by having a massive nullspace. $\endgroup$
    – Reid.Atcheson
    Aug 15, 2013 at 8:09
  • $\begingroup$ As a test I'd suggest looking at the nullspace of smaller approximations, and see how this nullspace behaves under finer approximation. I suspect that as you increase the approximation power, you will find eigenvalues closer and closer to zero. This could account for some of your stagnation issues, though ultimately any iterative method you throw at this problem will fail in the high frequency limits (see paper in response). $\endgroup$
    – Reid.Atcheson
    Aug 15, 2013 at 8:20

2 Answers 2

Reset to default
4
$\begingroup$

You might be familiar with the following paper already:

http://link.springer.com/chapter/10.1007%2F978-3-642-22061-6_10

Problems which are highly indefinite and oscillatory are very difficult to design robust iterative methods for. The paper gives some suggestions which might be helpful to you though, many of them have been extended to the time-harmonic Maxwell case as well.

Edit: For the benefit of others, though you already have encountered this frustration, I'll mention that the difficulties you describe in proper convergence only increase with increasing wavenumber, as this parameter increases the oscillatory nature of the solution. The paper linked I think describes the nature of this problem far better than I could in a few paragraphs, and I highly recommend reading through it.

I should mention however for respect's sake that there have been developed some robust iterative methods which have even seen successful parallel implementation. Here's one such paper

http://arxiv.org/abs/1007.4291

My (primitive) understanding is this method is kind of like a multiplicative Schwarz preconditioner using PML boundary conditions, with some brilliant insights that allow the preconditioner to be applied with lower computational complexity than normal (using H-matrix algebra).

Improve this answer
$\endgroup$
4
  • $\begingroup$ My understanding is also fairly primitive, but the way the author has described the algorithm is that under a specific type of domain decomposition, the Schur complement system solved over each subdomain yields an subproblem with the structure of a lower-dimensional Helmholtz equation that can be handled efficiently by specific direct methods. @Costis, in case you're curious, the idea has been extended to Maxwell's as well. $\endgroup$
    – Jesse Chan
    Aug 15, 2013 at 15:38
  • $\begingroup$ @JLC: Are you referring to the sweeping PML preconditioner? If so, it indeed seems to work rather well (especially for large problems) based on the papers I read although I think it requires much programming "effort" to implement the whole system... $\endgroup$
    – Costis
    Aug 21, 2013 at 8:28
  • $\begingroup$ @Costis. Seems like it. A PhD thesis came out of it, so I guess it's not trivial :). $\endgroup$
    – Jesse Chan
    Aug 21, 2013 at 14:22
  • $\begingroup$ Yes it is not an easy thing to implement, but to my knowledge it is the only existing robust iterative method for time harmonic wave equations in the high frequency limits (with the possible exception of wave-ray preconditioners, which I have not seen developed far beyond the early papers of Brandt). There are some who frequent this exchange who know far more about it than I, but its dominating flaw seems to be difficult parallel implementation, but otherwise seems to be similar in nature to Schwarz preconditioning for which there are packages that may make the task easier. $\endgroup$
    – Reid.Atcheson
    Aug 22, 2013 at 0:08
3
$\begingroup$

Convergence of iterative methods are affected by the condition number of the matrix, which tends to increase as the mesh is refined. Benzi has done work showing that improved convergence can be obtained using ILU or approximate inverse preconditioners if one first performs permutations to maximize the diagonal entries of the matrix. Search for "Preconditioning symmetric and highly indefinite problems" by Bollhoffer and "Preconditioning highly indefinite and nonsymmetric matrices" by Benzi et al. They used either the Harwell MC64 algorithm, which I believe free for academic use, or their own variant of the underlying algorithm. Search for "The design and use of algorithms for permuting large entries to the diagonal of sparse matrices" for a description of MC64. Generally, one first reorders to either minimize fill or bandwidth and then permutes the matrix to maximize the diagonal. One can also regularize the matrix using row and column scaling matrices.

There is also a variant of BiCGStab that begins with multiple starting vectors. I am not sure if it will help in this situation. Seach for ML(n)BiCGStab. The author also supplies matlab code at:

http://www.uwyo.edu/mathmyeung/r19/index.html

The best bet may just be a better more robust preconditioner (e.g., ILUt code from Saad, ILUPACK code from Bollhoffer, Trilinos preconditioners developed at Sandia National Labs, etc). Jack Dongarra of ORNL maintains a nice site with a list of freely available direct and iterative solvers for sparse matrices at:

http://www.netlib.org/utk/people/JackDongarra/la-sw.html

Site includes links to the software's official website.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.