Need help with a coding section - Convergence and UmfpackInverse Error

NGSolve
1

Hello,

I’m seeking assistance with a specific section of my code. Here’s the relevant snippet:

V = VectorL2(mesh, order=order)
Qi = H1(mesh, order=order, dirichlet=“gamma1|gamma2”)
QI = Compress(Qi, active_dofs=Qi.GetDofs(mesh.Boundaries(“gammaI”)))
X = FESpace([V, Q, Qbar, QI])
(u, p, pbar, pI), (v, q, qbar, qI) = X.TnT()

dx1 = dx(definedon=mesh.Materials(“Domain1”))
ds1 = dx(element_boundary=True, definedon=mesh.Materials(“Domain1”))
dx2 = dx(definedon=mesh.Materials(“Domain2”))
ds2 = dx(element_boundary=True, definedon=mesh.Materials(“Domain2”))
dsI = ds(definedon=mesh.Boundaries(“gammaI”))
dsI2 = ds(element_boundary=True, definedon=mesh.Boundaries(“gammaI”))
dsI3 = ds(skeleton=True, definedon=mesh.Boundaries(“gammaI”))

Bilinear form terms with pressure in the interface

…

a += pI.Deriv().Trace() * qI.Deriv().Trace() * dsI2
a += jump_u * v * n * dsI2
a += av_u * v * n * dsI2
a += jump_u * v.Other() * n * dsI2
a += -av_u * v.Other() * n * dsI3
a += pI * jump_v * dsI3
a += qI * jump_u * dsI3
…

Order=1 I am not getting convergence, and for order>1, I encounter the following error: “UmfpackInverse: Numeric factorization failed.”

I suspect the issue lies within the aforementioned terms of the bilinear form, specifically regarding the integration differences between dsI, dsI2, and dsI3. Could someone please shed some light on this problem?

I would greatly appreciate any assistance provided. Thank you!

2

It is always easier to help with a running minimal example.
Best
Christopher

4

Umfpack gives you the reason:

UMFPACK V6.1.0 (Jan 17, 2023): WARNING: matrix is singular

My first guess is, you defined oz1 and oz2, but only used oz1. Maybe a typo?

6

Hi Dori,

dsI = ds(definedon=mesh.Boundaries(“gammaI”)), or more compactly ds(“gammaI”) yields a boundary integrator over the boundary gammaI.
dsI2 = ds(element_boundary=True, definedon=mesh.Boundaries(“gammaI”)) will go over all boundary elements of gammaI and then integrate over the boundaries of these elements. I.e. for each segment, the corresponding endpoints get evaluated.

dsI3 = ds(skeleton=True, definedon=mesh.Boundaries(“gammaI”)) is similar to the first one (dsI). However, it includes the information of the attached volume element such that functions can be evaluated, which do not have a well-defined trace like L2 test/trial functions.

For example, the line
a+=IfPos(oz1, 1.0, 0.0)*pbar*v*n*ds
leads to a segmentation fault as you try to evaluate an L2 testfunction, v, at the boundary, which does not have a well-defined boundary trace.

Best,
Michael