Runtime error: Factor is exactly singular

Question

I am trying to implement 2 temperature models, the following equations:

1. C_e(∂T_e)/∂t=∇[k_e∇T_e ]-G(T_e-T_ph )+ A(r,t)

2. C_ph(∂T_ph)/∂t=∇[k_ph∇T_ph] + G(T_e-T_ph)

Code

from fipy.tools import numerix
import  scipy
import fipy
import numpy as np
from fipy import CylindricalGrid1D
from fipy import Variable, CellVariable, TransientTerm, DiffusionTerm, Viewer, LinearLUSolver, LinearPCGSolver, \
    LinearGMRESSolver, ImplicitDiffusionTerm, Grid1D

FIPY_SOLVERS = scipy

## Mesh


nr = 50
dr = 1e-7
# r = nr * dr
mesh = CylindricalGrid1D(nr=nr, dr=dr, origin=0)
x = mesh.cellCenters[0]

# Variables
T_e = CellVariable(name="electronTemp", mesh=mesh,hasOld=True)
T_e.setValue(300)
T_ph = CellVariable(name="phononTemp", mesh=mesh, hasOld=True)
T_ph.setValue(300)
G = CellVariable(name="EPC", mesh=mesh)
t = Variable()
# Material parameters
C_e = CellVariable(name="C_e", mesh=mesh)
k_e = CellVariable(name="k_e", mesh=mesh)

C_ph = CellVariable(name="C_ph", mesh=mesh)
k_ph = CellVariable(name="k_ph", mesh=mesh)

C_e = 4.15303 - (4.06897 * numerix.exp(T_e / -85120.8644))
C_ph = 4.10446 - 3.886 * numerix.exp(-T_ph / 373.8)
k_e = 0.1549 * T_e**-0.052
k_ph =1.24 + 16.29 * numerix.exp(-T_ph / 151.57)


G = numerix.exp(21.87 + 10.062 * numerix.log(numerix.log(T_e )- 5.4))

# Boundary conditions
T_e.constrain(300, where=x > 4.5e-6)
T_ph.constrain(300, where=x > 4.5e-6)

# Source  𝐴(𝑟,𝑡) = 𝑎𝐷(𝑟)𝜏−1 𝑒−𝑡/𝜏   , 𝐷(𝑟) = 𝑆𝑒 exp (−𝑟2/𝜎2)/√2𝜋𝜎2
sig = 1.0e-6
tau = 1e-15
S_e = 35

d_r = (S_e * 1.6e-9 * numerix.exp(-x**2 /sig**2)) / (numerix.sqrt(2. * 3.14 * sig**2))
A_t = numerix.exp(-t/tau)
a = (numerix.sqrt(2. * 3.14)) / (3.14 * sig)

A_r = a * d_r * tau**-1 * A_t 



eq0 = (TransientTerm(var=T_e, coeff=C_e) == DiffusionTerm(var=T_e, coeff=k_e) - G*(T_e - T_ph) + A_r

eq1 =(TransientTerm(var=T_ph, coeff=C_ph) == DiffusionTerm(var=T_ph, coeff=k_ph) + G*(T_e - T_ph)
eq = eq0 & eq1

dt = 1e-18
steps = 7000
elapsed = 0.
vi = Viewer((T_e, T_ph), datamin=0., datamax=2e4)


for step in range(steps):
    T_e.updateOld()
    T_ph.updateOld()
    vi.plot()
    res = 1e100
    dt *= 1.1
    while res > 1:
        res = eq.sweep(dt=dt)
        print(t, res)
    t.setValue(t + dt)

Problem

The code is working fine with very small dt = 1e-18, but I need to run it until e 1e-10.

With this time step is going to take very long time, and setting dt *= 1.1 the resduals at some point start to increase then

gives following runtime error:

factor is exactly singular

Even with very small increment dt*= 1.005 the same issue pop up. Using dt= 1.001 runs the code for quit long time then the residual get stuck at certain value.

Questions

1. I there any error in the fipy formalism of the equations?
2. What causes the error?
3. Is the error because of time step increase? If yes, how can I increase my time step?

Answer 1

In addition to @wd15's answer:

Your equations are extremely non-linear. You will likely benefit from Newton iterations to get decent convergence.

As @TimRoberts said, geometrically increasing the time step without bound is probably not a good idea.

I've recently posted a package called steppyngstounes that takes care of adapting timesteps. Although a standalone package, it's intended to work with FiPy. For example, you could change your solve loop to this:

from steppyngstounes import FixedStepper, PIDStepper

T_e.updateOld()
T_ph.updateOld()

for checkpoint in FixedStepper(start=0, stop=1e-10, size=1e-12):
    for step in PIDStepper(start=checkpoint.begin,
                           stop=checkpoint.end,
                           size=dt):
        res = 1e100
        for sweep in range(10):
            res = eq.sweep(dt=dt, underRelaxation=0.5)
            print(t, sweep, res)
        if step.succeeded(error=res / 1000):
            T_e.updateOld()
            T_ph.updateOld()
            t.value = step.end
        else:
            T_e.value = T_e.old
            T_ph.value = T_ph.old

        print('elapsed:', t.value)

    # the last step might have been smaller than possible,
    # if it was near the end of the checkpoint range
    dt = step.want
    _ = checkpoint.succeeded()
    vi.plot()

This code will update the viewer every 1e-12 time units, and adaptively make it's way between those checkpoints. There are other steppers in the package that would facilitate taking geometrically or exponentially increasing checkpoints, if that kept things more interesting.

You could probably get better overall performance by sweeping fewer times and letting the adapter take much smaller time steps in the beginning. I found that no time step was small enough to get the initial residual lower than 777.9. After the first couple of steps, the error metric could probably be much more aggressive, giving more accurate results.