[petsc-users] BiCGSTAB for general use

[Paul Anton Letnes]

On 14. aug. 2011, at 16.42, Jed Brown wrote:

> On Sun, Aug 14, 2011 at 06:18, Paul Anton Letnes <paul.anton.letnes at gmail.com> wrote:
> I believe of the first kind - there is. Our approach is to discretize the integral equation. The equations we are "really solving" are the Maxwell equations.
> 
> Is this the sort of system you're working with?
> 
> http://link.aps.org/doi/10.1103/PhysRevLett.104.223904

Ingve is my supervisor - so it's a good guess! But, no. BiCGSTAB works well for his code, by the way. He gets a massive speedup compared to direct solvers.

I am working on a different equation called the Reduced Rayleigh Equation (RRE) for 2 dimensions, which has not been solved directly before (AFAIK). Have a look here:
http://www.tandfonline.com/doi/abs/10.1088/0959-7174/6/3/006
The form of the solution is given in Eq. 5, and the equation itself is given in Eq. 8 with some definitions in Eq. 9-11. As far as I can tell there is only 2 terms, so it's a first kind equation. But by all means, feel free to correct me on that one. The cited paper uses a completely different approach (perturbation theory where 0. order is the flat surface) to solving the RRE, by the way. Our approach is more "brute force" - discretize the integral on a 2D grid in \vec{q}_\parallel space and solve the linear equation system.

The idea is that the RRE is less memory intensive than the rigorous approach in the paper you cited. However, if we can't avoid the LU factorization, the RRE approach will be much more CPU intensive, and of lesser interest than we would have hoped...

We have submitted a first paper to ArXiv; it will appear on monday, I think.

> Note that the system has the form
> 
> J_H(x_\parallel | \omega) = J_H(x_\parallel | \omega)_{inc} + \int (...) G(x | x') J_H(x_\parallel' | \omega) + \int (...) G(x | x') J_E(x_\parallel' | \omega)
> 
> which is the form of a second order integral equation. I assume the incident field J_H(...)_{inc} is known in this equation. If you dropped the term on the left hand side in this equation, you would have a Fredholm integral equation of the first kind to "solve", which is problematic at a mathematical level due to ill-posedness.

But since it's of the second kind, it's better posed, right?

> I have downloaded and attempted to use a different BiCGSTAB code. It converges, but only after several hundred (about 400) for a very small (not physically interesting) problem. It would appear that if we are to get good performance, some form of preconditioning is necessary.
> 
> Do the eigenvalues decay quickly? Can you plot some eigenvalues? They should decay rapidly to a positive value like (with appropriate scaling) 1.

I'll have a look at the eigenvalues. Is your suggestion to plot the magnitude of the complex eigenvectors as a function of their index, after sorting them? I guess that must be what you mean by "plotting".

cheers
Paul