Dirichlet Laplace Indirect Method

Dirichlet Laplace Indirect Method#

keys: homogeneous Dirichlet bvp, single layer potential, ACA

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngbem 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\limits_\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 potential operator. \(\mathrm V\) is regular and symmetric.

Demo 1: Assemble the single layer operator \(V\) as dense matrix and solve for unknwon density \(j\):

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.185628444636173     
CG iteration 2, residual = 0.5511126197504223     
CG iteration 3, residual = 0.1311312543903476     
CG iteration 4, residual = 0.031075382341886653     
CG iteration 5, residual = 0.007368321718552715     
CG iteration 6, residual = 0.0017873102488004224     
CG iteration 7, residual = 0.0005185169856633695     
CG iteration 8, residual = 0.00025605265661500233     
CG iteration 9, residual = 0.0003206576934490297     
CG iteration 10, residual = 0.00011035885851815262     
CG iteration 11, residual = 0.00016752392864684382     
CG iteration 12, residual = 4.08958140444724e-05     
CG iteration 13, residual = 3.274472821469365e-05     
CG iteration 14, residual = 1.632572462927271e-05     
CG iteration 15, residual = 1.0240009168263005e-05     
CG iteration 16, residual = 9.235240212462255e-06     
CG iteration 17, residual = 7.881910845196794e-06     
CG iteration 18, residual = 4.84151330511191e-06     
CG iteration 19, residual = 3.082553945086529e-06     
CG iteration 20, residual = 3.862950278127602e-06     
CG iteration 21, residual = 1.8320936585767728e-06     
CG iteration 22, residual = 1.8105355171682958e-06     
CG iteration 23, residual = 1.4120630469467857e-06     
CG iteration 24, residual = 5.770021991178657e-07     
CG iteration 25, residual = 2.9626965296568e-07     
CG iteration 26, residual = 2.5728166005001245e-07     
CG iteration 27, residual = 4.359698034344633e-07     
CG iteration 28, residual = 1.7797189551278336e-07     
CG iteration 29, residual = 1.7317514759169652e-07     
CG iteration 30, residual = 1.7259325500745118e-07     
CG iteration 31, residual = 7.969862678472237e-08     
CG iteration 32, residual = 7.935418621481961e-08     
CG iteration 33, residual = 3.85136669005619e-08     
CG iteration 34, residual = 2.8360883014874435e-08     
CG iteration 35, residual = 1.6221553170952995e-08

Demo 2: Approximate the single layer operator \(V\) by ACA. There are three parameters that define the approximation:

  • leafsize is the minimal size of a cluster in the cluster tree

  • eta relates to the so-called admissibility condition and determines the type of block: either near or far-field

  • eps is the accuracy of low rank approximation of a far field block

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.00010372613791385414     
CG iteration 2, residual = 0.00011332303929202542     
CG iteration 3, residual = 4.8475130043928946e-05     
CG iteration 4, residual = 3.752744414069605e-05     
CG iteration 5, residual = 2.5063257847685102e-05     
CG iteration 6, residual = 2.201161658598163e-05     
CG iteration 7, residual = 4.0910741896836506e-05     
CG iteration 8, residual = 1.5722470629475495e-05     
CG iteration 9, residual = 1.1398173534547079e-05     
CG iteration 10, residual = 8.601114677452872e-06     
CG iteration 11, residual = 1.7985157762401756e-05     
CG iteration 12, residual = 6.24956732982615e-06     
CG iteration 13, residual = 6.840010024435814e-06     
CG iteration 14, residual = 3.1039996622566697e-06     
CG iteration 15, residual = 2.601687250773817e-06     
CG iteration 16, residual = 3.3905259452258448e-06     
CG iteration 17, residual = 2.3772167442172267e-06     
CG iteration 18, residual = 1.2484134584125856e-06     
CG iteration 19, residual = 1.0020809034945668e-06     
CG iteration 20, residual = 1.4679545198086889e-06     
CG iteration 21, residual = 1.345029797669753e-06     
CG iteration 22, residual = 6.233413646170433e-07     
CG iteration 23, residual = 5.316529772220773e-07     
CG iteration 24, residual = 3.8221074414510546e-07     
CG iteration 25, residual = 4.2254716868647137e-07     
CG iteration 26, residual = 1.9054567735665295e-07     
CG iteration 27, residual = 2.456064134128932e-07     
CG iteration 28, residual = 1.8998524316205275e-07     
CG iteration 29, residual = 2.173992286349492e-07     
CG iteration 30, residual = 1.2534619876156766e-07     
CG iteration 31, residual = 8.403493254757788e-08     
CG iteration 32, residual = 4.3332209270053146e-08     
CG iteration 33, residual = 1.0157994209403936e-07     
CG iteration 34, residual = 2.9761572643430306e-08     
CG iteration 35, residual = 4.0323202672285155e-08     
CG iteration 36, residual = 2.737427871216972e-08     
CG iteration 37, residual = 1.7371828631137093e-08     
CG iteration 38, residual = 3.452055776446457e-08     
CG iteration 39, residual = 1.3309568079797635e-08     
CG iteration 40, residual = 9.486019205138243e-09     
CG iteration 41, residual = 6.427829807220539e-09     
CG iteration 42, residual = 1.2848413223273628e-08     
CG iteration 43, residual = 3.761031417102926e-09     
CG iteration 44, residual = 4.679452921229483e-09     
CG iteration 45, residual = 2.190122037970437e-09     
CG iteration 46, residual = 1.996239646971565e-09     
CG iteration 47, residual = 2.3159899902351333e-09     
CG iteration 48, residual = 1.2723913819350062e-09     
CG iteration 49, residual = 9.31259051689337e-10     
CG iteration 50, residual = 1.1204428304854894e-09     
CG iteration 51, residual = 1.413388424930364e-09     
CG iteration 52, residual = 8.997634695803242e-10     
CG iteration 53, residual = 4.4531603486308826e-10     
CG iteration 54, residual = 4.5270100016281045e-10     
CG iteration 55, residual = 2.6152775675607896e-10     
CG iteration 56, residual = 4.617915988555627e-10     
CG iteration 57, residual = 1.5886672961052302e-10     
CG iteration 58, residual = 2.7846309529962185e-10     
CG iteration 59, residual = 8.723015476489554e-11     
CG iteration 60, residual = 2.5907217972154056e-10     
CG iteration 61, residual = 7.938117887022281e-11     
CG iteration 62, residual = 6.271285664807571e-11     
CG iteration 63, residual = 4.6670339271953476e-11     
CG iteration 64, residual = 5.1924137459351305e-11     
CG iteration 65, residual = 4.366007100324626e-11     
CG iteration 66, residual = 4.353862381900156e-11     
CG iteration 67, residual = 2.195362901022156e-11     
CG iteration 68, residual = 1.762485504806998e-11     
CG iteration 69, residual = 2.4145021673366737e-11     
CG iteration 70, residual = 9.719501974076612e-12     
CG iteration 71, residual = 9.661870415501173e-12     
CG iteration 72, residual = 1.0814296376953874e-11     
CG iteration 73, residual = 4.518647797308727e-12     
CG iteration 74, residual = 1.2093135661526256e-11     
CG iteration 75, residual = 5.420444183339184e-12     
CG iteration 76, residual = 2.952789092939208e-12     
CG iteration 77, residual = 1.7325571207483807e-12     
CG iteration 78, residual = 3.850113041607153e-12     
CG iteration 79, residual = 2.0103340936785782e-12     
CG iteration 80, residual = 1.8670382142067295e-12     
CG iteration 81, residual = 1.0746221727153705e-12     
CG iteration 82, residual = 7.385642929343038e-13     

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

{'name': 'ngbem FarFieldBlock - matvec2 LaplaceSL', 'time': 0.0076680862051428525, 'counts': 1831, 'flops': 12485592.0, 'Gflop/s': 1.6282539953223438}
{'name': 'ngbem FarFieldBlock - matvec1 LaplaceSL', 'time': 0.0016000993157137982, 'counts': 1831, 'flops': 12462408.0, 'Gflop/s': 7.788521548389369}
{'name': 'ngbem FarFieldBlock - norm LaplaceSL', 'time': 0.02413285989228843, 'counts': 1831, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock - kernel LaplaceSL', 'time': 0.005464391943065623, 'counts': 3662, 'flops': 5036040.0, 'Gflop/s': 0.9216103186724726}
{'name': 'ngbem FarFieldBlock - ACA LaplaceSL', 'time': 0.06518542551288707, 'counts': 174, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock LaplaceSL', 'time': 0.12686599146346794, 'counts': 189, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - HMatrix::MultAdd', 'time': 0.1696692086510814, 'counts': 147, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem disjoint LaplaceSL', 'time': 1.6634663368779543, 'counts': 94631, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem common edge LaplaceSL', 'time': 1.69745426576888, 'counts': 1146, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem common vertex LaplaceSL', 'time': 2.00345731113448, 'counts': 3371, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem identic panel LaplaceSL', 'time': 0.6705158937567653, 'counts': 382, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - loops LaplaceSL', 'time': 6.10330901470501, 'counts': 299, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - all LaplaceSL', 'time': 6.274055100880135, 'counts': 299, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - BaseClass::CalcHMatrix', 'time': 2.826205817393841, 'counts': 2, 'flops': 0.0, 'Gflop/s': 0.0}

Notes: