Singular matrix for non-linear conductivity in heat equation

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)

Hi Sindre,
by setting only the Dirichlet data of the temperature on the boundaries the initial guess uh for the Newton is zero in the interior. Evaluating (the first variation of) e.g. ln(uh) leads to NaNs in the matrix such that the matrix is not regular (functions like uh**2 work also for uh=0).Setting a non-zero function as initial guess including the appropriate boundary conditions solves the problem, in your example
#uh.Set(dirichletBCs, ng.BND)
cf = 2-ng.x
uh.Set(cf)
Another possibility is to first extend the boundary data to the domain by solving e.g. an auxiliary Poisson problem. 

Best,
Michael

Hi, Michael

Thank you for the clear and concise explanation!

For beginners like myself, I would like to add that writing
uh.Set(1)
uh.Set(dirichletBCs, ng.BND)
does not solve the problem because the second call to .Set somehow overwrites the first call.

Best,
Sindre

Hi, you can find some typical pitfalls here:

docu.ngsolve.org/latest/i-tutorials/unit.../subdomains.html#id2