Hello,

I am looking to solve the steady heat equation for various systems with various temperature-dependent thermal conductivity functions. For certain conductivity functions, I encounter a warning from UMFPACK stating that my assembled matrix is singular. Subsequently, the iterative solver crashes and the program exits. Unfortunately, my theoretical understanding of FEM is rather limited, so I do not understand why some of the functions I have explored yield a singular matrix while others do not.

For example, it appears to be problematic to have float exponents. E.g. k(T) = 1 + T^1.1 (T=temperature, k=conductivity). k(T) = 1 + T^2 does not yield the same problem. k(T) = 1 + 0.9*sin(pi*T) also works. However, k(T) = 1 + ln(T) does not, even if the boundary conditions ensure k>0 throughout.

If anyone could shed some light on why I sometimes get singular matrices, that would be greatly apprectiated. I have attached an example code below. If it can be modified in a simple way to provide non-singular matrices for e.g. k(T) = 1 + T^1.1 or k(T) = 1 + ln(T), I would love to know how.

Best regards,

Sindre

import ngsolve as ng

from ngsolve.webgui import Draw

from netgen.geom2d import unit_square

from ngsolve.solvers import Newton

mesh = ng.Mesh(unit_square.GenerateMesh(maxh=0.1))

V = ng.H1(mesh, order=1, dirichlet=“left|right”)

u = V.TrialFunction()

v = V.TestFunction()

uh = ng.GridFunction(V)

dirichletBCs = mesh.BoundaryCF({‘left’:2, ‘right’:1})

uh.Set(dirichletBCs, ng.BND)

k = lambda T: 1 + T**1.1

a = ng.BilinearForm(V)

a += k(u) * ng.grad(u) * ng.grad(v) * ng.dx

Newton(a, uh)