I’m attempting to implement a stray field solver within ngsolve, using ngbem to compute the double layer potential operator. The algorithm is from Redirecting, specifically Algorithm 12 =
I have some code that works for step (i) and (iii) using the standard methods of NGSolve. I am not convinced that step (ii) is functioning properly, but I am unsure of what the error is, or if my discrete formulation is correct.
Attached are four pieces of code:
- The
strayfield.pyfile contains all of the actual code for the computations. - A test on a unit sphere with m=(0,0,1), for which the
stray_fieldfunction should return (0,0,-1/3), and an energy of 2pi/9. - A test on a unit cube with m=(0,0,1). The stray field has no analytical form, but the energy should be exactly 1/6.
- An auxiliary file.
General_Functions.py (2.2 KB)
strayfield.py (8.4 KB)
strayfield_cube_test.ipynb (2.1 KB)
strayfield_sphere_test.ipynb (2.3 KB)
When I run these tests, I don’t see the convergence I would expect. It is “close”, so it cannot be too wrong. In the sphere test, we generally see convergence as h->0 (generating meshes with maxh=1/n, for n=1,2,4,...,32)
but it doesn’t seem to be fast enough.
|hmax | energy | L^2 error|
|1.986897527 | 0.6244154164 | 0.0006382110957|
|1.993270202 | 0.6343580189 | 0.0002835307397|
|0.6891751754 | 0.6780127155 | 2.15E-05|
|0.3123171932 | 0.6934280217 | 1.67E-06|
|0.1694168532 | 0.6969098422 | 7.87E-07|
|0.09040233825 | 0.6978044914 | 4.29E-07|
As for the cube test, it is even worse. Again, it gets closer to the expected value of 1/6, but eventually misses and goes past the target.
|hmax | energy|
|1.732050808 | 0.1350023107|
|1 | 0.1465266853|
|0.5097136556 | 0.1642313063|
|0.2726903446 | 0.1664041535|
|0.0688505752 | 0.1669602132|
|0.0372839163 | 0.1675898603|
Within the strayfield.py file, step (i) is dealt with in compute_internal. We feed the output of this into DL_boundary_projection to achieve step (ii). Then for step (iii), we use the gfu.Set method to put this onto the boundary of an H1 function and then apply the BVP function. This is all combined in stray_field.