Laplace solver using direct formulation (Dirichlet)

Laplace solver using direct formulation (Dirichlet)#

keys: Laplace double layer potential operator, electrostatics, homogeneous Dirichlet bvp, Neumann data

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

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)\) (direct ansatz)

\[ u(x) = \underbrace{ \int_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{1}{\| x-y\|} } \, u_1(y)\, \mathrm{d}\sigma_y}_{\displaystyle{ \mathrm{ SL}(u_1) }} - \underbrace{ \int_\Gamma \displaystyle{\frac{1}{4\,\pi}\, \frac{\langle n(y) , x-y\rangle }{\| x-y\|^3} } \, u_0(y)\, \mathrm{d}\sigma_y}_{\displaystyle{ \mathrm{DL}(u_0) }}\]

and solve for the Neumann data \(u_1 \in H^{-\frac12}(\Gamma)\) by the boundary element method, i.e.,

\[ \forall \, v\in H^{-\frac12}(\Gamma): \quad \left\langle \gamma_0 \left(\mathrm{SL}(u_1)\right), v \right\rangle_{-\frac12}= \left\langle u_0, v\right\rangle_{-\frac12} + \left\langle \gamma_0 \left(\mathrm{DL}(u_0)\right), 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.2)).Curve(4)

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

fesL2 = SurfaceL2(mesh, order=3, dual_mapping=True)
u,v = fesL2.TnT()
fesH1 = H1(mesh, order=4)
uH1,vH1 = fesH1.TnT()
print ("ndofL2 = ", fesL2.ndof, "ndof H1 = ", fesH1.ndof)

ndofL2 =  7220 ndof H1 =  22641

Compute the interpolation of exact Dirichlet data \(u_0\) in finite element space:

uexa = 1/ sqrt( (x-1)**2 + (y-1)**2 + (z-1)**2 )
u0 = GridFunction(fesH1)
u0.Interpolate (uexa)
Draw (u0, mesh, draw_vol=False, order=3)

BaseWebGuiScene

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

\[ \mathrm{V} \, \mathrm{u}_1 = \left( \frac12 \,\mathrm{M} + \mathrm{K} \right) \, \mathrm{u}_0\,, \]

where the linear operators are as follows

\(\mathrm{V}\) is the single layer operator. \(\mathrm V\) is regular and symmetric.
\(\mathrm{M}\) is a mass matrix.
\(\mathrm{K}\) is the double layer operator.

We approximate the linear operators as hmatrices and solve for the Neumann data \(u_1\) with an iterative solver:

u1 = GridFunction(fesL2)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager(): 
    V = SingleLayerPotentialOperator(fesL2, intorder=12, leafsize=40, eta=3., eps=1e-4, method="aca", testhmatrix=False)

    M = BilinearForm( uH1 * v.Trace() * ds(bonus_intorder=3)).Assemble()
    K = DoubleLayerPotentialOperator(fesH1, fesL2, intorder=12, leafsize=40, eta=3., eps=1e-4, method="aca")
    
    rhs = ( (0.5 * M.mat + K.mat)*u0.vec).Evaluate()
    CG(mat = V.mat, pre=pre, rhs = rhs, sol=u1.vec, tol=1e-8, maxsteps=200, initialize=False, printrates=True)

Draw (u1, mesh, draw_vol=False, order=3);

CG iteration 1, residual = 0.2688253801037296     
CG iteration 2, residual = 0.06553558458073674     
CG iteration 3, residual = 0.015544599077948142     
CG iteration 4, residual = 0.004632058606436214     
CG iteration 5, residual = 0.004593036294366361     
CG iteration 6, residual = 0.0008022189950998422     
CG iteration 7, residual = 0.00019766046284216818     
CG iteration 8, residual = 7.15892463816215e-05     
CG iteration 9, residual = 0.0001030320559235765     
CG iteration 10, residual = 7.699909626524041e-05     
CG iteration 11, residual = 3.1502554443343895e-05     
CG iteration 12, residual = 1.8228831474569672e-05     
CG iteration 13, residual = 1.5343515082076535e-05     
CG iteration 14, residual = 1.7402177351037414e-05     
CG iteration 15, residual = 6.589075458534609e-06     
CG iteration 16, residual = 5.802384353983555e-06     
CG iteration 17, residual = 1.1118866472261058e-05     
CG iteration 18, residual = 4.963671995941951e-06     
CG iteration 19, residual = 6.6925470414869674e-06     
CG iteration 20, residual = 4.66702470235276e-06     
CG iteration 21, residual = 2.1838341886573426e-06     
CG iteration 22, residual = 1.4625200656887405e-06     
CG iteration 23, residual = 1.7967828880434654e-06     
CG iteration 24, residual = 1.7350325927730044e-06     
CG iteration 25, residual = 1.589123284127319e-06     
CG iteration 26, residual = 1.0526815231904489e-06     
CG iteration 27, residual = 8.878557830071291e-07     
CG iteration 28, residual = 4.2594350522941694e-07     
CG iteration 29, residual = 3.804753677237142e-07     
CG iteration 30, residual = 3.831943457883365e-07     
CG iteration 31, residual = 3.537955230879237e-07     

CG iteration 32, residual = 3.497449268176541e-07     
CG iteration 33, residual = 1.95000222228805e-07     
CG iteration 34, residual = 2.095693642993343e-07     
CG iteration 35, residual = 3.858478903961728e-07     
CG iteration 36, residual = 1.4045168299810223e-07     
CG iteration 37, residual = 1.1182182403058054e-07     
CG iteration 38, residual = 5.517831213831011e-08     
CG iteration 39, residual = 9.929492891920444e-08     
CG iteration 40, residual = 6.762790390359661e-08     
CG iteration 41, residual = 5.3668453774058395e-08     
CG iteration 42, residual = 3.1798784297157776e-08     
CG iteration 43, residual = 2.574529250963146e-08     
CG iteration 44, residual = 2.2837643129845963e-08     
CG iteration 45, residual = 3.621838882517188e-08     
CG iteration 46, residual = 2.893687540235177e-08     
CG iteration 47, residual = 2.1352335237316694e-08     
CG iteration 48, residual = 1.040593074767152e-08     
CG iteration 49, residual = 9.194248692424211e-09     
CG iteration 50, residual = 8.885071352902476e-09     
CG iteration 51, residual = 1.2971281649348865e-08     
CG iteration 52, residual = 5.514226106326673e-09     
CG iteration 53, residual = 6.345431815835099e-09     
CG iteration 54, residual = 3.1571079296941915e-09     
CG iteration 55, residual = 3.3915831149934337e-09     
CG iteration 56, residual = 4.1965363423052695e-09     
CG iteration 57, residual = 2.364657159649791e-09     

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

{'name': 'ngbem - HMatrix::MultAdd', 'time': 0.3538058838800736, 'counts': 58, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem disjoint LaplaceDL', 'time': 22.12617002633463, 'counts': 38897, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem common edge LaplaceDL', 'time': 34.428677623375286, 'counts': 1592, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem common vertex LaplaceDL', 'time': 41.603784699149855, 'counts': 4794, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem identic panel LaplaceDL', 'time': 13.598112534351337, 'counts': 533, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock - matvec2 LaplaceDL', 'time': 0.17181667322479216, 'counts': 4615, 'flops': 73607065.0, 'Gflop/s': 0.42840466887458584}
{'name': 'ngbem FarFieldBlock - matvec1 LaplaceDL', 'time': 0.063394616502485, 'counts': 4615, 'flops': 143971849.0, 'Gflop/s': 2.271042194795144}
{'name': 'ngbem FarFieldBlock - norm LaplaceDL', 'time': 0.39298400624405583, 'counts': 4615, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock - kernel LaplaceDL', 'time': 0.7901241040355118, 'counts': 9230, 'flops': 38862145.0, 'Gflop/s': 0.049184861974864345}
{'name': 'ngbem FarFieldBlock - ACA LaplaceDL', 'time': 1.750626997102634, 'counts': 382, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - loops LaplaceDL', 'time': 112.2799173875817, 'counts': 804, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock LaplaceDL', 'time': 4.243423060776662, 'counts': 494, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - all LaplaceDL', 'time': 119.2479335451965, 'counts': 804, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock - matvec2 LaplaceSL', 'time': 0.1320554461924254, 'counts': 7960, 'flops': 60723201.0, 'Gflop/s': 0.45983109936652533}
{'name': 'ngbem FarFieldBlock - matvec1 LaplaceSL', 'time': 0.027091304741423405, 'counts': 7960, 'flops': 60632551.0, 'Gflop/s': 2.2380816124847263}
{'name': 'ngbem FarFieldBlock - norm LaplaceSL', 'time': 0.2473252484690905, 'counts': 7960, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock - kernel LaplaceSL', 'time': 0.3803503632251338, 'counts': 15920, 'flops': 25839954.0, 'Gflop/s': 0.06793724023527495}
{'name': 'ngbem FarFieldBlock - ACA LaplaceSL', 'time': 0.9927976793333007, 'counts': 793, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem FarFieldBlock LaplaceSL', 'time': 2.64535421986481, 'counts': 877, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem disjoint LaplaceSL', 'time': 4.062457333895154, 'counts': 9396, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem common edge LaplaceSL', 'time': 9.924644570033442, 'counts': 467, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem common vertex LaplaceSL', 'time': 11.96600206818119, 'counts': 1469, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem identic panel LaplaceSL', 'time': 3.8488542258607916, 'counts': 135, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - loops LaplaceSL', 'time': 29.889819128635487, 'counts': 728, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - all LaplaceSL', 'time': 36.35680725366288, 'counts': 728, 'flops': 0.0, 'Gflop/s': 0.0}
{'name': 'ngbem - BaseClass::CalcHMatrix', 'time': 45.988357938299124, 'counts': 2, 'flops': 0.0, 'Gflop/s': 0.0}

Let’s have a look at the exact Neumann data and compute the error of the numerical solution:

graduexa = CF( (uexa.Diff(x), uexa.Diff(y), uexa.Diff(z)) )
n = specialcf.normal(3)
u1exa = graduexa*n
Draw (u1exa, mesh, draw_vol=False, order=3);
print ("L2-error =", sqrt (Integrate ( (u1exa-u1)**2, mesh, BND)))

L2-error = 0.0021893572396393434