Solve BlockMatrix system without a preconditioner

General
1

I am having issues solving a system given by a BlockMatrix. It seems I need a preconditioner (which I don’t have) as direct solvers are unavailable.

##### Geometry and mesh
h_max = 0.1

geo = CSG2d()
geo.Add(Rectangle(pmin=(0.0,0.0), pmax = (1.0,1.0), bc = "bdry"))
m = geo.GenerateMesh(maxh=h_max,quad_dominated=False)
mesh = Mesh(m)
mesh.Curve(1)

n = specialcf.normal(mesh.dim)
h = specialcf.mesh_size

# Physical parameters

b = 15.0
c = 0.5
m1 = 0.5
m2 = 0.5
m3 = 1.5

kk = CoefficientFunction((0.0,-1.0))

def tr(zeta): #For TnT
	return zeta[0,0] + zeta[1,1]

def devi(zeta): #For TnT
	return zeta - 0.5*tr(zeta)*CoefficientFunction( (1,0,0,1), (2,2))

def mu(fi):
	return (1-c*fi)**-2

def gamma(fi):
	return c*fi*(1-c*fi)**2

def theta(nfi):
	return m1 + m2*(1+ InnerProduct(nfi,nfi) )**(m3/2.0 - 1.0)

## FE spaces
Hhsig = HDiv(mesh, order=0, RT=True)**2 # Space for σₕ
Hhu = L2(mesh, order=0, dim=2) # Space for uₕ
Hhphi = H1(mesh, order=1, dirichlet="bdry") # Space for φₕ
Hhlag = NumberSpace(mesh)

# Trial and Test Functions
sigh, tauh = Hhsig.TnT()
uh, vh = Hhu.TnT()
phih, psih = Hhphi.TnT()
lamh, muh = Hhlag.TnT()


gfu = GridFunction(Hhu) # Fixed point function
gfphi = GridFunction(Hhphi) # Newton iterate
gfphi.Interpolate(phi_exact, definedon=mesh.Boundaries("bdry")) # Impose boundary conditions

# Bilinear forms for the fluid
a1F = BilinearForm(Hhsig)
a1F += ( InnerProduct(devi(sigh),devi(tauh))/(mu(gfphi) ) )*dx

a2F = BilinearForm(trialspace=Hhu,testspace=Hhsig)
a2F += ( InnerProduct(uh,div(tauh)) )*dx

lmF = BilinearForm(trialspace=Hhsig, testspace=Hhlag)
lmF += (tr(sigh)*muh)*dx

rhsF1 = LinearForm(Hhsig)
rhsF1 += (InnerProduct(u_exact,tauh.Trace()*n))*ds(definedon=mesh.Boundaries("bdry"))

rhsF2 = LinearForm(Hhu)
rhsF2 += (-gfphi*InnerProduct(ff,vh))*dx

# forms for the transport #We impose boundary conditions later
a1T = BilinearForm(Hhphi)
a1T += ( theta(grad(phih))*InnerProduct(grad(phih),grad(psih)) - phih*InnerProduct(gfu,grad(psih)) )*dx

rhsT1 = LinearForm(Hhphi)
rhsT1 += ( gamma(gfphi)*InnerProduct(kk,grad(psih)) + g*psih  )*dx

a1F.Assemble()
a2F.Assemble()
lmF.Assemble()
rhsF1.Assemble()
rhsF2.Assemble()
a1T.Assemble()
rhsT1.Assemble()

#We also need the following
gfsig = GridFunction(Hhsig)
gflag = GridFunction(Hhlag)

# For the fluid
MF = BlockMatrix([[a1F.mat, a2F.mat, lmF.mat], [a2F.mat.T, None, None], [lmF.mat.T, None, None]]) #System matrix
rhsF = BlockVector([rhsF1.vec, rhsF2.vec])  #RHS
solF = BlockVector([gfsig.vec, gfu.vec, gflag.vec]) #Solution vector for the problem

Then, the idea is to iterate between solving the “fluid” problem and the “temperature problem”

itt = 0;
tol = 10.0;

# Global solution
solTot = BlockVector([gfsig.vec,gfu.vec, gfphi.vec]) #We don't care about the Lagrange multiplier. We just need it to be 0
solTot_old = solTot.CreateVector() #To store previous iteration

res = rhsT1.vec.CreateVector() #Use to homogenize the temperature equations
print("Start of fixed point loop")

while ((tol > 1e-6) & (itt < 30)):
	itt += 1

	# calculating the solution for the fluid
	solF.data = krylovspace.GMResSolver(MF, rhsF, maxsteps=10000, printrates=True)

	# calculating the solution for the temperature. We homogenize to account for boundary conditions on φₕ
	res.data = rhsT1.vec - a1T.mat * gfphi.vec
	gfphi.vec.data = a1T.mat.Inverse(Hhphi.FreeDofs())*res

	# computing tol
	tol = sqrt( InnerProduct(solTot-solTot_old,solTot-solTot_old) / InnerProduct(solTot,solTot) )
	print(f"tol = {tol}")
	 
	# updating data for the next step
	solTot_old.data = solTot

Overall my code runs fine until I invoke the GMRes solver.

assert (freedofs is None) != (pre is None) # either pre or freedofs must be given

Any idea how I could fix this? Should I just proceed and assemble the system in just one big matrix instead of building a block system?

2

Based on BlockMatrix Direct Solver, it seems that I am doomed for the start and will have to deal with a huge sparse matrix. However, I am having issues with this approach too.

## FE spaces
Hhsig = VectorValued(HDiv(mesh, order=0, RT=True), dim=2) # Space for σₕ
Hhu = L2(mesh, order=0, dim=2) # Space for uₕ
Hhphi = H1(mesh, order=1, dirichlet="bdry") # Space for φₕ
Hhlag = NumberSpace(mesh)

fes = Hhsig * Hhu * Hhphi * Hhlag # Product FE space

This is, however, prompting the following error

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
netgen.libngpy._meshing.NgException: Product space needs same dimension for all components

Any ideas for how to solve it?

3

Hi, the dim=2 argument changes the block structure of the matrix which doesn’t work here
What you want is either

VectorValued(...)**2

or

MatrixValued(...)