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.9sin(piT) 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)