FETI methods

80. FETI methods#

Finite Element Tearing and Interconnection.

As the name says, we break the global (finite element) system apart, and then enforce continuity.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import CSG2d, Rectangle

geo = CSG2d()

mx, my = 3,3
for i in range(mx): 
    for j in range(my):
        rect = Rectangle(pmin=(i/mx,j/my), \
                         pmax=((i+1)/mx,(j+1)/mx), \
                         mat='mat'+str(i)+str(j), \
                         bc = 'default', \
                         bottom = 'bot' if j == 0 else None)
                  
        geo.Add(rect)
        
mesh = Mesh(geo.GenerateMesh(maxh=0.04))
print (mesh.GetMaterials())
print (mesh.GetBoundaries())
Draw (mesh);

('mat00', 'mat01', 'mat02', 'mat10', 'mat11', 'mat12', 'mat20', 'mat21', 'mat22')
('bot', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'bot', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'bot', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default')

We define \(H^1\)-spaces on sub-domains \(\Omega_i\). The product spaces is

\[ V = \Pi_i H^1(\Omega_i) \]

fes = None
for dom in mesh.Materials('.*').Split():
    fesi = Compress(H1(mesh, definedon=dom, dirichlet="bot"))
    fes = fes * fesi if fes else fesi

print ("ndof =", fes.ndof)

u, v = fes.TnT()

domtrial = {} 
domtest = {}
for nr,dom in enumerate (mesh.Materials('.*').Split()):
    domtrial[dom] = u[nr]
    domtest[dom] = v[nr]

ndof = 817

We identify the interfaces (i.e. internal boundaries) between two sub-domains. On these \(\gamma_{ij}\) we define spaces for the Lagrange parameter. Although the spaces are \(H^{-1/2}\), we use \(H^1(\gamma_{ij})\) to obtain the same number of constraints as we have basis functions on the interface.

feslam = None
for inter in mesh.Boundaries('.*').Split():
    doms = inter.Neighbours(VOL).Split()
    if len(doms) == 2:
        feslami = Compress(H1(mesh, definedon=inter))
        feslam = feslam * feslami if feslam else feslami 
    
print ("ndof-lam:", feslam.ndof)

lam, mu = feslam.TnT()

intertrial = {} 
intertest = {}
nr = 0
for inter in mesh.Boundaries('.*').Split():
    doms = inter.Neighbours(VOL).Split()
    if len(doms) == 2:
        intertrial[inter] = lam[nr]
        intertest[inter] = mu[nr]
        nr = nr+1

ndof-lam: 108

We define bilinear-forms on the sub-domains:

\[ a(u,v) = \sum_{\Omega_i} \int_{\Omega_i} \nabla u_i \nabla v_i + u_i v_i \]

and the constraint equations

\[ b(u,\mu) = \sum_{\gamma_{ij}} \int_{\gamma_{ij}} (u_i - u_j) \mu_{ij} \]

a = BilinearForm(fes)
f = LinearForm(fes)
b = BilinearForm(trialspace=fes, testspace=feslam)

for (ui,vi) in zip(u,v):
    a += grad(ui)*grad(vi)*dx + 1*ui*vi*dx
    f += y*x*vi*dx
    

for inter in mesh.Boundaries('.*').Split():
    doms = inter.Neighbours(VOL).Split()
    if len(doms) == 2:
        dom1,dom2 = doms
        # a += 1*(domtrial[dom1]-domtrial[dom2])*(domtest[dom1]-domtest[dom2])*ds(inter)
        b += (domtrial[dom1]-domtrial[dom2]) * intertest[inter] * ds(inter)
        
a.Assemble()
b.Assemble()
f.Assemble()

<ngsolve.comp.LinearForm at 0x10e1f1570>

Obviously, if we only solve the decomposed sub-problems, we don’t get the correct solution:

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(inverse="sparsecholesky", freedofs=fes.FreeDofs())*f.vec
gftot = CF ( list(gfu.components) )
Draw (gftot, mesh);

Next we solve the saddle-point problem

\[\begin{split} \left( \begin{array}{cc} A & B^T \\ B & 0 \end{array} \right) \left( \begin{array}{c} u \\ \lambda \end{array} \right) = \left( \begin{array}{c} f \\ 0 \end{array} \right) \end{split}\]

We explicitly build the Schur-complement \(S = B A^{-1} B^T\), and use conjugate gradients to solve for the Lagrange parameter \(\lambda\). Then we recover the primal variable \(u\) from the first equation:

ainv = a.mat.Inverse(inverse="sparsecholesky", freedofs=fes.FreeDofs())
S = b.mat @ ainv @ b.mat.T
g = (b.mat @ ainv * f.vec).Evaluate()

from ngsolve.krylovspace import CGSolver
invS = CGSolver(S, pre=IdentityMatrix(feslam.ndof), printrates="\r", maxiter=500)

lam = g.CreateVector()
lam.data = invS * g
gfu.vec.data = ainv * (f.vec - b.mat.T * lam)

CG iteration 1, residual = 0.08403599336703109     
CG iteration 2, residual = 0.03971919917200752     
CG iteration 3, residual = 0.01229555425486588     
CG iteration 4, residual = 0.00209015543474033     
CG iteration 5, residual = 0.0015921163492495247     
CG iteration 6, residual = 0.00032163144527917325     
CG iteration 7, residual = 0.001486650169981125     
CG iteration 8, residual = 0.00021211111708816004     
CG iteration 9, residual = 0.0008992713275613841     
CG iteration 10, residual = 7.97870424274299e-05     
CG iteration 11, residual = 0.00010989556910275802     
CG iteration 12, residual = 6.807312175357524e-05     
CG iteration 13, residual = 4.4893134619348585e-05     
CG iteration 14, residual = 4.221168193883369e-05     
CG iteration 15, residual = 4.772794989088839e-05     
CG iteration 16, residual = 0.0003902376619710067     
CG iteration 17, residual = 6.116080611146976e-05     
CG iteration 18, residual = 4.7907524159514695e-05     
CG iteration 19, residual = 0.00015104356593841114     
CG iteration 20, residual = 0.00011455858334166559     
CG iteration 21, residual = 2.8842911389966815e-05     
CG iteration 22, residual = 1.1308682699054299e-05     
CG iteration 23, residual = 6.720277923249615e-05     
CG iteration 24, residual = 1.0728344173676783e-05     
CG iteration 25, residual = 1.2562185071306156e-05     
CG iteration 26, residual = 2.843722051876385e-05     
CG iteration 27, residual = 1.575570123894986e-05     
CG iteration 28, residual = 2.9851038052631544e-05     
CG iteration 29, residual = 9.262409501838687e-06     
CG iteration 30, residual = 5.278453604751753e-06     
CG iteration 31, residual = 2.48165354867617e-05     
CG iteration 32, residual = 4.286505224127795e-06     
CG iteration 33, residual = 1.8670169881917116e-05     
CG iteration 34, residual = 2.897947407783534e-06     
CG iteration 35, residual = 1.0063391913992392e-05     
CG iteration 36, residual = 2.8810774163262462e-06     
CG iteration 37, residual = 9.514952075563536e-06     
CG iteration 38, residual = 1.3353005724913072e-06     
CG iteration 39, residual = 4.608585554903282e-06     
CG iteration 40, residual = 8.411348617119802e-07     
CG iteration 41, residual = 5.881460064130038e-07     
CG iteration 42, residual = 6.421373232085135e-07     
CG iteration 43, residual = 4.954697536053879e-07     
CG iteration 44, residual = 4.173145243523698e-07     
CG iteration 45, residual = 1.203579372682444e-06     
CG iteration 46, residual = 2.82748412619732e-07     
CG iteration 47, residual = 1.6831414237428172e-06     
CG iteration 48, residual = 4.407631613365212e-07     
CG iteration 49, residual = 8.235698848966272e-07     
CG iteration 50, residual = 2.582460201324559e-07     
CG iteration 51, residual = 1.4014811948944164e-06     
CG iteration 52, residual = 1.89537770278461e-07     
CG iteration 53, residual = 7.18239326239495e-07     
CG iteration 54, residual = 2.1354489295710067e-07     
CG iteration 55, residual = 1.2794268561723303e-07     
CG iteration 56, residual = 6.939175519254065e-08     
CG iteration 57, residual = 3.785161810542102e-08     
CG iteration 58, residual = 2.697838018857469e-07     
CG iteration 59, residual = 9.007767462483043e-08     
CG iteration 60, residual = 2.6741148658528377e-08     
CG iteration 61, residual = 6.171242169527921e-08     
CG iteration 62, residual = 2.252011089157917e-08     
CG iteration 63, residual = 8.2570709463492e-08     
CG iteration 64, residual = 1.845707984343893e-08     
CG iteration 65, residual = 1.6491627752406698e-08     
CG iteration 66, residual = 1.3715758714675145e-08     
CG iteration 67, residual = 1.5284325663857127e-08     
CG iteration 68, residual = 3.3509906761220354e-08     
CG iteration 69, residual = 3.500116704668778e-08     
CG iteration 70, residual = 6.269432811736806e-09     
CG iteration 71, residual = 7.499121459722843e-09     
CG iteration 72, residual = 5.842549500246563e-09     
CG iteration 73, residual = 1.9544491028129525e-08     
CG iteration 74, residual = 5.4846707043150605e-09     
CG iteration 75, residual = 1.0493228166340138e-08     
CG iteration 76, residual = 3.7213460179708344e-09     
CG iteration 77, residual = 1.7642322774630083e-08     
CG iteration 78, residual = 6.1457055730195295e-09     
CG iteration 79, residual = 1.6612253938470774e-09     
CG iteration 80, residual = 1.6581012527089735e-08     
CG iteration 81, residual = 1.3606471893290851e-09     
CG iteration 82, residual = 1.775540546546009e-09     
CG iteration 83, residual = 5.436024997236365e-09     
CG iteration 84, residual = 9.209841802955364e-10     
CG iteration 85, residual = 9.189004510910192e-10     
CG iteration 86, residual = 6.543414380099021e-09     
CG iteration 87, residual = 6.39376057770696e-10     
CG iteration 88, residual = 6.690793878626996e-10     
CG iteration 89, residual = 4.954524988040401e-10     
CG iteration 90, residual = 4.512365097809827e-10     
CG iteration 91, residual = 4.184518646616133e-10     
CG iteration 92, residual = 2.376294787821936e-09     
CG iteration 93, residual = 1.5802375801172456e-10     
CG iteration 94, residual = 1.1631690327272827e-10     
CG iteration 95, residual = 8.210516495096637e-11     
CG iteration 96, residual = 1.1166378431252478e-10     
CG iteration 97, residual = 2.2799027620592605e-10     
CG iteration 98, residual = 1.1160794010521349e-10     
CG iteration 99, residual = 7.164002978391902e-11     
CG iteration 100, residual = 8.642007949458749e-11     
CG iteration 101, residual = 7.228111988549506e-11     
CG iteration 102, residual = 2.1495557599322867e-11     
CG iteration 103, residual = 1.2482396528288556e-10     
CG iteration 104, residual = 6.354153128306859e-11     
CG iteration 105, residual = 2.3759257345796876e-11     
CG iteration 106, residual = 2.203340764549582e-11     
CG iteration 107, residual = 2.0502118520135045e-11     
CG iteration 108, residual = 1.0024127234859912e-11     
CG iteration 109, residual = 1.989051753270956e-11     
CG iteration 110, residual = 3.9583025195272006e-11     
CG iteration 111, residual = 1.2444781742724771e-11     
CG iteration 112, residual = 5.2926088093011324e-11     
CG iteration 113, residual = 6.3767551669114815e-12     
CG iteration 114, residual = 6.064378099172241e-12     
CG iteration 115, residual = 1.2346168202002402e-11     
CG iteration 116, residual = 3.1731567271102746e-12     
CG iteration 117, residual = 1.0034274924070966e-11     
CG iteration 118, residual = 3.868009654151403e-12     
CG iteration 119, residual = 3.4895926840424588e-12     
CG iteration 120, residual = 1.1085370483308783e-12     
CG iteration 121, residual = 1.1222150109879732e-11     
CG iteration 122, residual = 1.5419012858384775e-12     
CG iteration 123, residual = 7.298152195730669e-13     
CG iteration 124, residual = 1.4903482436642743e-12     
CG iteration 125, residual = 9.128385842495598e-13     
CG iteration 126, residual = 4.895117949331075e-13     
CG iteration 127, residual = 3.786976191748678e-13     
CG iteration 128, residual = 3.9236962476294173e-13     
CG iteration 129, residual = 4.0718173613822063e-13     
CG iteration 130, residual = 1.9800158109928926e-12     
CG iteration 131, residual = 7.472437150273745e-13     
CG iteration 132, residual = 2.27921192240865e-13     
CG iteration 133, residual = 2.8600961405495873e-13     
CG iteration 134, residual = 3.093847080520956e-13     
CG iteration 135, residual = 2.625844366509997e-13     
CG iteration 136, residual = 2.0708562851366104e-13     
CG iteration 137, residual = 1.0321664556812178e-13     
CG iteration 138, residual = 1.0855068742571685e-12     
CG iteration 139, residual = 1.4074313271603042e-13     
CG iteration 140, residual = 3.318629762870249e-14     
CG converged in 140 iterations to residual 3.318629762870249e-14

gftot = CF ( list(gfu.components) )
Draw(gftot, mesh);

80.1. Preconditioner for \(S\)#

The function space for the Lagrange parameter is \(\Pi_{ij} H^{-1/2} (\gamma_{ij})\). The Schur complement matrix scales like a bilinear-form in \(H^{-1/2}\). As a preconditioner, we need a matrix which scales like \(\Pi_{ij} H^{1/2}\). We cheat a bit with the non-additivity of the \(H^{1/2}\)-norm, and use \(\sum_{\Omega_i} \| \operatorname{tr} u \|_{H^{1/2}(\partial \Omega_i)}^2\)

We have to map functions from the skeleton onto the domain:

bnddofs = fes.GetDofs(mesh.Boundaries(".*"))
innerdofs = ~bnddofs

massbnd = BilinearForm(fes)
for (ui,vi) in zip(u,v):
    massbnd += ui*vi*ds
massbnd.Assemble()
invmassbnd = massbnd.mat.Inverse(inverse="sparsecholesky", freedofs=bnddofs)

massinter = BilinearForm(feslam)
for inter in mesh.Boundaries('.*').Split():
    doms = inter.Neighbours(VOL).Split()
    if len(doms) == 2:
        massinter += intertrial[inter]*intertest[inter]*ds(inter)
massinter.Assemble()

emb = invmassbnd@b.mat.T@massinter.mat.Inverse(inverse="sparsecholesky")

used dof inconsistency
(silence this warning by setting BilinearForm(...check_unused=False) )

The \(H^{1/2}(\partial \Omega_i)\)-norms are obtained by forming Schur-complements of the sub-domain matrices with respect to boundary dofs:

SchurDir = a.mat - a.mat@a.mat.Inverse(inverse="sparsecholesky", freedofs=innerdofs)@a.mat

pre = emb.T @ SchurDir @ emb

invS = CGSolver(S, pre=pre, printrates="\r", maxiter=500)

lam = g.CreateVector()
lam.data = invS * g
gfu.vec.data = ainv * (f.vec - b.mat.T * lam)

CG iteration 1, residual = 2.00544510021893     
CG iteration 2, residual = 0.8601107411343858     
CG iteration 3, residual = 0.2607019225332258     
CG iteration 4, residual = 0.07092390115693398     
CG iteration 5, residual = 0.18854148259682096     
CG iteration 6, residual = 0.04148923094504725     
CG iteration 7, residual = 0.02005492668175809     
CG iteration 8, residual = 0.0049385274075988955     
CG iteration 9, residual = 0.005065786308053268     
CG iteration 10, residual = 0.00274976526375405     
CG iteration 11, residual = 0.0007957605690849792     
CG iteration 12, residual = 0.00019849103261233912     
CG iteration 13, residual = 5.3657345829248144e-05     
CG iteration 14, residual = 3.470109807340424e-05     
CG iteration 15, residual = 2.039923161314892e-05     
CG iteration 16, residual = 4.923210728679797e-06     
CG iteration 17, residual = 2.537520206911699e-06     
CG iteration 18, residual = 3.625231949488012e-06     
CG iteration 19, residual = 9.88978359878323e-07     
CG iteration 20, residual = 6.517507837203933e-07     
CG iteration 21, residual = 1.472796970933967e-06     
CG iteration 22, residual = 3.138056451376175e-07     
CG iteration 23, residual = 2.565176484800696e-07     
CG iteration 24, residual = 1.7684222523134668e-07     
CG iteration 25, residual = 1.7948968585386554e-07     
CG iteration 26, residual = 1.5336307237815367e-08     
CG iteration 27, residual = 4.06672284005037e-09     
CG iteration 28, residual = 1.0658429693010248e-09     
CG iteration 29, residual = 5.794336920486616e-10     
CG iteration 30, residual = 1.7521849481247322e-10     
CG iteration 31, residual = 3.5087570130423114e-10     
CG iteration 32, residual = 1.450846077211376e-11     
CG iteration 33, residual = 3.2950048008498856e-11     
CG iteration 34, residual = 2.7805639122786753e-11     
CG iteration 35, residual = 1.2981497810886503e-11     
CG iteration 36, residual = 8.433399195958605e-12     
CG iteration 37, residual = 1.130024942605001e-11     
CG iteration 38, residual = 7.905832846232014e-13     
CG converged in 38 iterations to residual 7.905832846232014e-13

Draw(gftot, mesh);

FETI methods

Contents

80. FETI methods#

80.1. Preconditioner for \(S\)#