Laplace solver using indirect formulation (Dirichlet)#

keys: Laplace single layer potential, electrostatic potential, homogeneous Dirichlet bvp

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from libbem import *
from ngsolve import Projector, Preconditioner
from ngsolve.krylovspace import CG

Loading ngbem library

We consider the Dirichlet boundary value problem

\[\begin{split} \left\{ \begin{array}{rcl l} \Delta u &=& 0, \quad &\Omega \subset \mathbb R^3\,,\\ \gamma_0 u&=& u_0, \quad &\Gamma = \partial \Omega\,.\end{array} \right. \end{split}\]

Let us choose the following ansatz for the solution \(u\in H^1(\Omega)\) (indirect ansatz)

\[ u(x) = \underbrace{ \int_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{1}{\| x-y\|} } \, j(y)\, \mathrm{d}\sigma_y }_{\displaystyle{ \mathrm{SL}(j) } }\]

and solve for the density \(j\in H^{-\frac12}(\Gamma)\) by the boundary element method, i.e. the numerical solution of the variational formulation

\[ \forall \, v\in H^{-\frac12}(\Gamma): \quad \left\langle \gamma_0 \left(\mathrm{SL}(j)\right), v \right\rangle_{-\frac12} = \left\langle u_0, v\right\rangle_{-\frac12} \,. \]

Define the geometry \(\Omega \subset \mathbb R^3\) and create a mesh:

sp = Sphere( (0,0,0), 1)
mesh = Mesh( OCCGeometry(sp).GenerateMesh(maxh=0.3)).Curve(3)

Create test and trial function finite element spaces for \(H^{-\frac12}(\Gamma)\) according to the given mesh:

fesL2 = SurfaceL2(mesh, order=3, dual_mapping=True)
u,v = fesL2.TnT()

Define Dirichlet data \(u_0\) and compute the right hand side vector:

u0 = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
rhs = LinearForm (u0*v.Trace()*ds(bonus_intorder=3)).Assemble()

The discretisation of the above variational formulation leads to a system of linear equations, ie

\[ \mathrm{V} \, \mathrm{j} = \mathrm{rhs} \,, \]

where \(\mathrm{V}\) is the single layer operator. \(\mathrm V\) is regular and symmetric.

Demo 1: Assemble the single layer operator \(V\) as dense matrix and solve for \(j\) with an iterative solver:

j = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager(pajetrace=1000*1000*1000):
    V = SingleLayerPotentialOperator(fesL2, intorder=10, method="dense")
    CG(mat = V.mat, pre=pre, rhs = rhs.vec, sol=j.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=True)

Draw (j);

CG iteration 1, residual = 2.1856284446361722     
CG iteration 2, residual = 0.5511126197504219     
CG iteration 3, residual = 0.13113125439034773     
CG iteration 4, residual = 0.031075382341886598     
CG iteration 5, residual = 0.007368321718552593     
CG iteration 6, residual = 0.0017873102488002246     
CG iteration 7, residual = 0.0005185169856628255     
CG iteration 8, residual = 0.0002560526565371585     
CG iteration 9, residual = 0.0003206576807033507     
CG iteration 10, residual = 0.00011035833531823466     
CG iteration 11, residual = 0.0001674770829883903     
CG iteration 12, residual = 4.0896513121476776e-05     
CG iteration 13, residual = 3.274473300635716e-05     
CG iteration 14, residual = 1.632572463283343e-05     
CG iteration 15, residual = 1.0240009168202185e-05     
CG iteration 16, residual = 9.235240211917132e-06     
CG iteration 17, residual = 7.881910820038137e-06     
CG iteration 18, residual = 4.841512597141705e-06     
CG iteration 19, residual = 3.0817964457634693e-06     

CG iteration 20, residual = 3.864433067165785e-06     
CG iteration 21, residual = 1.8320946929674603e-06     
CG iteration 22, residual = 1.810535521215225e-06     
CG iteration 23, residual = 1.4120630469529938e-06     
CG iteration 24, residual = 5.77002199113252e-07     
CG iteration 25, residual = 2.962696505682998e-07     
CG iteration 26, residual = 2.572788414147209e-07     
CG iteration 27, residual = 4.357077550694016e-07     
CG iteration 28, residual = 1.779905987387955e-07     
CG iteration 29, residual = 1.7317505071728763e-07     
CG iteration 30, residual = 1.725920301877434e-07     
CG iteration 31, residual = 7.969876086897504e-08     
CG iteration 32, residual = 7.935418732541848e-08     
CG iteration 33, residual = 3.851366708217703e-08     
CG iteration 34, residual = 2.8360881195691114e-08     
CG iteration 35, residual = 1.621976036918233e-08

Demo 2: Compte an approximation of the single layer operator \(V\) and solve for \(j\) with iterative solver:

pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager(pajetrace=1000*1000*1000):
    V = SingleLayerPotentialOperator(fesL2, intorder=10, leafsize=40, eta=3., eps=1e-4, method="aca", testhmatrix=True)    
    CG(mat = V.mat, pre=pre, rhs = rhs.vec, sol=j.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=True)

Draw (j);

dense: 3080 x 3080
compression rate 0.349931  2-norm error 0.000589323
CG iteration 1, residual = 0.00010372613361416254     
CG iteration 2, residual = 0.00011332304421634448     
CG iteration 3, residual = 4.8475130100534105e-05     
CG iteration 4, residual = 3.752744413900531e-05     
CG iteration 5, residual = 2.5063257844046756e-05     
CG iteration 6, residual = 2.201161649282029e-05     
CG iteration 7, residual = 4.091074210729724e-05     
CG iteration 8, residual = 1.5722470650406946e-05     
CG iteration 9, residual = 1.1398173533714403e-05     
CG iteration 10, residual = 8.601114682987693e-06     
CG iteration 11, residual = 1.7985157935834744e-05     
CG iteration 12, residual = 6.2495673291056765e-06     
CG iteration 13, residual = 6.840010014190462e-06     
CG iteration 14, residual = 3.1039996619119485e-06     
CG iteration 15, residual = 2.6016873278899223e-06     
CG iteration 16, residual = 3.390525782088221e-06     
CG iteration 17, residual = 2.377216742960555e-06     
CG iteration 18, residual = 1.2484134581618446e-06     
CG iteration 19, residual = 1.0020809167544435e-06     
CG iteration 20, residual = 1.467954530079095e-06     
CG iteration 21, residual = 1.3450297598977127e-06     
CG iteration 22, residual = 6.233413652711596e-07     
CG iteration 23, residual = 5.316529768772103e-07     
CG iteration 24, residual = 3.822107798886752e-07     
CG iteration 25, residual = 4.225471204391166e-07     
CG iteration 26, residual = 1.9054567751438192e-07     
CG iteration 27, residual = 2.4560641325559065e-07     
CG iteration 28, residual = 1.899852454813338e-07     
CG iteration 29, residual = 2.1739925181351957e-07     
CG iteration 30, residual = 1.2534619390055752e-07     
CG iteration 31, residual = 8.403493245788256e-08     
CG iteration 32, residual = 4.333220932454454e-08     
CG iteration 33, residual = 1.015799643685005e-07     
CG iteration 34, residual = 2.9761572078100446e-08     
CG iteration 35, residual = 4.032320271101176e-08     
CG iteration 36, residual = 2.7374278722324998e-08     
CG iteration 37, residual = 1.737182905755099e-08     
CG iteration 38, residual = 3.4520557190232194e-08     
CG iteration 39, residual = 1.3309567923488516e-08     
CG iteration 40, residual = 9.486019266672398e-09     
CG iteration 41, residual = 6.427829801110218e-09     
CG iteration 42, residual = 1.2848416741338577e-08     
CG iteration 43, residual = 3.761031330423225e-09     
CG iteration 44, residual = 4.679452910898918e-09     
CG iteration 45, residual = 2.190122036559631e-09     
CG iteration 46, residual = 1.9962400043626533e-09     
CG iteration 47, residual = 2.3159894368478814e-09     
CG iteration 48, residual = 1.2723913823070666e-09     
CG iteration 49, residual = 9.312590538032909e-10     
CG iteration 50, residual = 1.1204428392127928e-09     
CG iteration 51, residual = 1.4133887218374485e-09     
CG iteration 52, residual = 8.997633944422024e-10     
CG iteration 53, residual = 4.4531603485629705e-10     
CG iteration 54, residual = 4.527009997484055e-10     
CG iteration 55, residual = 2.615277828888713e-10     
CG iteration 56, residual = 4.6179145984743047e-10     
CG iteration 57, residual = 1.588667294015795e-10     
CG iteration 58, residual = 2.784630954652936e-10     
CG iteration 59, residual = 8.723015565024802e-11     
CG iteration 60, residual = 2.590722046322888e-10     
CG iteration 61, residual = 7.938117752372998e-11     
CG iteration 62, residual = 6.271285706220372e-11     
CG iteration 63, residual = 4.6670339142744413e-11     
CG iteration 64, residual = 5.192415044188918e-11     
CG iteration 65, residual = 4.36600639174352e-11     
CG iteration 66, residual = 4.353862325618007e-11     
CG iteration 67, residual = 2.1953629100643763e-11     
CG iteration 68, residual = 1.7624855025373323e-11     
CG iteration 69, residual = 2.414503450699543e-11     
CG iteration 70, residual = 9.719501134026967e-12     
CG iteration 71, residual = 9.661870417526893e-12     
CG iteration 72, residual = 1.08142964186749e-11     

CG iteration 73, residual = 4.518647929051921e-12     
CG iteration 74, residual = 1.209313546344474e-11     
CG iteration 75, residual = 5.420444020705122e-12     
CG iteration 76, residual = 2.9527890842994155e-12     
CG iteration 77, residual = 1.7325571305513662e-12     
CG iteration 78, residual = 3.850114886249477e-12     
CG iteration 79, residual = 2.0103338433322948e-12     
CG iteration 80, residual = 1.8670381963343142e-12     
CG iteration 81, residual = 1.0746221761820988e-12     
CG iteration 82, residual = 7.385642983609332e-13     

for t in Timers():
    if "ngbem" in t["name"]:
        print (t)

{'name': 'ngbem FarFieldBlock - matvec2 LaplaceSL', 'time': 0.02510996370511808, 'counts': 2123, 'flops': 12485592.0, 'Gflop/s': 0.497236561017215}
{'name': 'ngbem FarFieldBlock - matvec1 LaplaceSL', 'time': 0.00559765081292412, 'counts': 2123, 'flops': 12473136.0, 'Gflop/s': 2.228280472801454}
{'name': 'ngbem FarFieldBlock - norm LaplaceSL', 'time': 0.04782915262406372, 'counts': 2123, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock - kernel LaplaceSL', 'time': 0.07387245444163758, 'counts': 4246, 'flops': 5039028.0, 'Gflop/s': 0.06821254333685432}
{'name': 'ngbem FarFieldBlock - ACA LaplaceSL', 'time': 0.18769478386642635, 'counts': 204, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock LaplaceSL', 'time': 0.4301981212551047, 'counts': 223, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - HMatrix::MultAdd', 'time': 0.605981293988976, 'counts': 147, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem disjoint LaplaceSL', 'time': 8.139544149046872, 'counts': 94750, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem common edge LaplaceSL', 'time': 4.3827823884261825, 'counts': 1124, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem common vertex LaplaceSL', 'time': 5.183806652877994, 'counts': 3304, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem identic panel LaplaceSL', 'time': 1.675121835609008, 'counts': 371, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - loops LaplaceSL', 'time': 19.622057012616107, 'counts': 269, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - all LaplaceSL', 'time': 20.462699723234323, 'counts': 269, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - BaseClass::CalcHMatrix', 'time': 9.06872795152745, 'counts': 2, 'flops': 0.0, 'Gflop/s': 0.0}

Notes: