BlockMatrix Direct Solver

NGSolve
1

Dear NGSolve,

Is there a possibility to use a direct solver for a BlockMatrix?
I read through the Documentation and the NGS build-in solvers but I didn’t figure out how to use/adapt things.
Here is a small example

from netgen.geom2d import unit_square
from ngsolve import *


mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
V = H1(mesh, order= 1)
u,v=V.TnT()
a=BilinearForm(V,symmetric=True)
a+=grad(u)*grad(v)*dx
a.Assemble()
f=LinearForm(V)
f.Assemble()
f+=v*dx

N=NumberSpace(mesh)
p,q=N.TnT()
b=BilinearForm(trialspace=V,testspace=N)
b+=u*q*dx
b.Assemble()
g=LinearForm(N)
g.Assemble()

K=BlockMatrix( [[a.mat,b.mat.T],[b.mat,None]] )
rhs=BlockVector( [f.vec,g.vec] )
gfu=GridFunction(V)
gfN=GridFunction(N)
sol=BlockVector( [gfu.vec,gfN.vec] )

sol.data=K.Inverse(inverse="pardiso")*rhs #???[/code]

Best regards,
Nepomuk

2

Hi Nepomuk,

as you have found: you cannot directly call a direct solver from a BlockMatrix
BlockMatrix entries can be everything, not only sparse matrices.

Recommendation:
use a product-space:

fes = H1(...) * NumberSpace(...)
(u,p), (v,q) = fes.TnT()

then NGSolve assembles one big sparse matrix.

If you want to build a big sparse matrix from blocks, you have to go via numpy, e.g. using the coordinate format, see attached example.

Joachim

blocksparse.ipynb (3.2 KB)