Working with block matrices

General
1

Hi Everyone,

I am trying to solve the Maxwell eigenvalue problem in a 2.5D domain following the example in the documentation. Below is a mwe:

import numpy as np
from ngsolve import *
from netgen.geom2d import unit_square

mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))

r = y
m = 0

def a1(arg1, arg2):
    return r * curl(arg1)*curl(arg2)*dx
def a3(arg1, arg2):
    return m*(arg1[0]*grad(arg2)[0] + 1/r*arg1[0]*arg2[0])*dx
def a4(arg1, arg2):
    return m**2/r*(r**2*grad(arg1)[0]*grad(arg2)[0] + r*grad(arg1)[0]*arg2 + r*arg1*grad(arg2)[0] + arg1*arg2)*dx
def b1(arg1, arg2):
    return r*arg1*arg2*dx
def b2(arg1, arg2):
    return r*arg1*arg2*dx

order = 2
VrVz = HCurl(mesh, order=order, dirichlet=".*")
um, vm = VrVz.TnT()
G, Vphi = VrVz.CreateGradient()
V = VrVz * Vphi

(u, u_phi) = V.TrialFunction()
(v, v_phi) = V.TestFunction()

a = BilinearForm(a1(u, v) + a4(u_phi, v_phi) + a3(u, v_phi) + a3(v, u_phi))
b = BilinearForm(b1(u, v) + b2(u_phi, v_phi))
mmat = BilinearForm(b1(um, vm))
apre = BilinearForm(a1(u, v) + a4(u_phi, v_phi) + a3(u, v_phi) + a3(v, u_phi) + b1(u, v) + b2(u_phi, v_phi))
pre = Preconditioner(apre, type="direct", inverse="sparsecholesky")

with TaskManager():
    a.Assemble()
    b.Assemble()
    mmat.Assemble()
    apre.Assemble()

    I_phi = IdentityMatrix(Vphi.ndof, complex=True)
    GT = G.CreateTranspose()
    math1 = GT @ mmat.mat @ G
    math1[0, 0] += 1
    invh1 = math1.Inverse(inverse="sparsecholesky", freedofs=Vphi.FreeDofs())
    proj = IdentityMatrix(V.ndof, complex=True) - BlockMatrix([[G @ invh1 @ GT @ mmat.mat, None],
                                                               [None, I_phi]])
    projpre = proj @ pre.mat

    evals, evecs = solvers.PINVIT(a.mat, b.mat, pre=projpre, num=4, maxit=20, printrates=True)

When I try to run this, I get the error

---> 51 evals, evecs = solvers.PINVIT(a.mat, b.mat, pre=projpre, num=4, maxit=30, printrates=True)

File ~\AppData\Roaming\Python\Python311\site-packages\ngsolve\eigenvalues.py:87, in PINVIT(mata, matm, pre, num, maxit, printrates, GramSchmidt)
     85 for v in vecs[0:num]:
     86     v.SetRandom()
---> 87 uvecs[:] = pre * vecs[0:num]
     88 lams = Vector(num * [1])
     90 for i in range(maxit):

RuntimeError: Bad dynamic_cast!

which I think has something to do with the way that I am building the preconditioner BlockMatrix. How can I fix it?

Best regards

2

You cannot mix block-matrices with non-block mats, you do it that way:

# proj = IdentityMatrix(V.ndof, complex=True) - BlockMatrix([[G @ invh1 @ GT @ mmat.mat, None],
#                                                           [None, I_phi]])

embHc, embH1 = V.embeddings
proj = IdentityMatrix(V.ndof) - embHc @ G@invh1@GT@mmat.mat @ embHc.T - embH1 @ I_phi @ embH1.T

embHc and embH1 are embedding matrices.
the first is of shape (V.ndof, VrVz), with an identity matrix in the upper rows, and padded with 0, the second for the second space.

with

print (proj.GetOperatorInfo())

you can see the linear operator you have built up.

Joachim

3

Thanks a lot. This fixes the problem.

Best regards
Sos