Convergence (Laplace, 1st kind) needs stricter params

[labarba]

The new convergence tests with Laplace equation produced a line plot that curves upward at high N (i.e., losing convergence order) with the original parameters used: p = 10, tol = 1e-6.

After much head-scratching, we looked at the number of Gauss points. In the method section of Simon's thesis, it says the number of Gauss points is set to k = 3 or 4 for far integrals, and ~20 for near-singular integrals. The code has a default setting of k = 3, but you can modify it in args. We don't know the value of k used for the convergence plot in the thesis.

[labarba]

With p = 10, k = 4, tol = 1e-6, we get the following convergence plot:

[labarba]

With p = 12, k = 4, tol = 1e-6

[labarba]

Conclusion

The 1st-kind formulation shows the expected convergence with stricter parameters, as shown above: p = 12, k = 4, told = 1e-6.

we decided to include this plot in the paper

The 2nd-kind formulation appears to not be so sensitive, as the convergence line plot looked fine with the lower p and lower k. This formulation also requires very few iterations to converge.

[ncclementi]

Maybe this helps in clarifying more the "The 2nd-kind formulation appears to not be so sensitive".

• In the 1st-kind formulation the condition number DEPENDS ON the size of the matrix (number of elements) then the number of iterations in the linear solver increase with the number of elements.
• In the 2nd-kind formulation the condition number is INDEPENDENT of the size of the matrix (number of elements) then the linear solver will require less iterations.

I've read this in Chris thesis, and in Lu. et al 2006 - (4th page second column) though I couldn't find the 1st source that explain this. @cdcooper84 do you have any good source to recommend where they explain better de dependence and independence on the number of elements in the condition number for 1st and 2nd kind formulation respectively?

[cdcooper84]

Hey, this is a interesting discussion!

I think there are two different things here: one is the strange convergence behavior of the first kind equation, and the other is the condition number with respect to number of elements.

• Condition number: the eigenvalues of the matrix resulting from the double layer potential are bound between 0 and 2\pi (or 0 and 1/2 if you divide by 4pi). I can't find a reference for a derivation that shows this, but I did try it in pygbe. The matrix resulting from the single layer potential doesn't have bounded eigenvalues, and they just grow as N increases. This affects the iteration count, but we should still see convergence.
• Convergence: I would be more worried about the first kind equation converging as $1/\sqrt{N}$ than anything else. I know that Lu et al. 2006 report $1/\sqrt{N}$ convergence for their second kind formulation, but I don't think it should converge like that. For instance, check the convergence test for a sphere with Dirichlet conditions in PyGBe, this is a first kind equation and is converging as 1/N.

I remember seeing that it converged as 1/\sqrt{N} when I used the potential on one panel to test convergence, but this was not a smart move as the location of the panel wasn't exactly the same for every mesh. When I used an integral quantity (the energy) I did see 1/N convergence. @tingyu66, What are you guys using to test convergence? Are using the value of dphi/dn at a single point?

[tingyu66]

we used the l2 norm of the relative error of dphi / dn, so it is
considering every panel.

On Friday, March 11, 2016, cdcooper84 [email protected] wrote:

Hey, this is a interesting discussion!

I think there are two different things here: one is the strange
convergence behavior of the first kind equation, and the other is the
condition number with respect to number of elements.

Condition number: the eigenvalues of the matrix resulting from the
double layer potential are bound between 0 and 2\pi (or 0 and 1/2 if you
divide by 4pi). I can't find a reference for a derivation that shows this,
but I did try it in pygbe. The matrix resulting from the single layer
potential doesn't have bounded eigenvalues, and they just grow as N
increases. This affects the iteration count, but we should still see
convergence.

Convergence: I would be more worried about the first kind equation
converging as $1/\sqrt{N}$ than anything else. I know that Lu et al. 2006
report $1/\sqrt{N}$ convergence for their second kind formulation, but I
don't think it should converge like that. For instance, check the convergence
test for a sphere with Dirichlet conditions in PyGBe
https://github.com/barbagroup/pygbe/blob/master/bem_pycuda/regression_tests/figs/error_energy_dirichlet_surface.pdf,
this is a first kind equation and is converging as 1/N.

I remember seeing that it converged as 1/\sqrt{N} when I used the
potential on one panel to test convergence, but this was not a smart move
as the location of the panel wasn't exactly the same for every mesh. When I
used an integral quantity (the energy) I did see 1/N convergence.
@tingyu66 https://github.com/tingyu66, What are you guys using to test
convergence? Are using the value of dphi/dn at a single point?

—
Reply to this email directly or view it on GitHub
#15 (comment)
.

---

Tingyu Wang (Vince)

[labarba]

@cdcooper84 -- look what it says in Simon's thesis ... I'm sure he got this from you!

So, which is it?

[ncclementi]

@cdcooper84 Adding to Lorena's question another issue. You reported on your thesis (i.e First kind formulation) O(1/N) but I read in Geng & Krasny-2013 (i.e Second kind) that they got O(1/sqrt(N)). What is right the opposite of Simon's thesis says.

[cdcooper84]

Oh, I hadn't realized that they were the opposite!

I think that both should converge as first order. The problem is that, as far as I'm aware, nobody has been able to analytically show the convergence of collocation in BEM, even though experience shows that it does converge. This is the main reason why BEM people from the math side avoid collocation and use Galerkin. Hence, there is no way of saying "this should converge in such way", and people just trust equivalent codes.

This is something that I remember having discussed with Simon, and I did see the same behavior when I was not using an integral quantity to compare, but once I compared the energy (which is the quantity I was worried about), I did see 1st order. At that time, I blamed this behavior to the fact that the center of each panel does not exactly lie on the surface of the sphere, so when I looked at the norm of the error of dphi/dn, actually, I was not comparing at the same point. Perhaps you could try calculating phi in a point away from the surface, and compare that. In fact, this is what Liu does to look at errors in his Fast BEM book.

About the 1/sqrt(N) convergence in Geng & Krasny, I think this might be due to low accuracy in singular/near singular integrals, but I would need to test it to be sure. In fact, in the second kind formulation, they need to integrate the adjoint of the double layer, which is hard to do (check Tausch, Wang, White 2001). Again, this is just my impression.

[tingyu66]

Hey, Chris, thanks for the reply.
I searched the BEM books and numerical papers for the order of convergence for the collocation method, and indeed could not find the certain answer. Still lots of them claim that the solution converges well instead of reporting the order of convergence.

I went through Liu's FMM-BEM book, where he gave two examples of FMM-BEM for potentials problems: in the annual heat conducting plate example, phi & phi/dn on the surface are used as the quantity to compare with analytical solution; in the later example of electrostatic field around conducting spheres, charge densities on the spheres are used. So I didn't found the case where he picks an exterior point to calculate the phi.

I tried to use a high k for singular & near singular integrals, but still it yields to the same convergence curve (with a slope -0.75).

[cdcooper84]

Just re-checked Liu's book, and I see I was confused. The point is at the surface, however, they test the same point every time, which is on the surface for the mesh and the actual sphere.

[labarba]

An update: @tingyu66 did a new test with larger ncrit=700 (previous was 400) and was able to get almost linear convergence, with slope ~0.96 ... "close enough" ... My question to him was how many of those runs do any FMM at all, because it could be mostly P2P. He will check.

But the interesting thing is that mat-vec with FMM in this case is affecting the convergence (no relaxation here).

[cdcooper84]

Excellent! That is great news!
The sphere test is hard. The geometry is simple so the discretization error is low, and you need very fine parameters of the FMM and quadrature rules to get the convergence. In pygbe, we used the treecode with P=15 and GMRES tolerance 10^-9!

[tingyu66]

- red line: relative error of phi on an **external point** (p=12, k=4, ncrit=400, tol=1e-6) - black line: L2 norm of the relative error of phi on **surface panels** (p=12, k=7, ncrit=700, tol=1e-5) - solid line: 1st-kind - dotted line: 2nd-kind

result update, again.

• good news comes first. When I used the phi of an external point to compare with the analytical solution (red lines), I got the linear convergence rate for both 1st and 2nd kind problems (same as Chris's result).
• bad news is that, when I generated the black solid line using ncrit=700 yesterday, I mistakenly set solver's tolerance to 1e-5 in my auto scripts (I just ran too many of them recently), which is not enough. However, when I used a smaller tolerance, the slope of the L2norm relative error was 0.75 again.