Fipy error:’’The Factor is exactly singular’’, when applying Neumann boundary conditions

Question

We’re trying to solve a one-dimensional Coupled Continuity-Poisson problem in Fipy. When applying Dirichlet’s conditions, it gives the correct results, but when we change the boundaries conditions to Neumann’s which is closer to our problem, it gives “The Factor is exactly singular’’ error. Any help is highly appreciated. The code is as follows (0<x<2.5):

from fipy import *
from fipy import Grid1D, CellVariable, TransientTerm, DiffusionTerm, Viewer
import numpy as np
import math
import matplotlib.pyplot as plt
from matplotlib import cm
from cachetools import cached, TTLCache #caching to increase the speed of python
cache = TTLCache(maxsize=100, ttl=86400) #creating the cache object: the 
#first argument= the number of objects we store in the cache.
#____________________________________________________
nx=50
dx=0.05
L=nx*dx
e=math.e
m = Grid1D(nx=nx, dx=dx)
print(np.log(e))
#____________________________________________________
phi = CellVariable(mesh=m, hasOld=True, value=0.)
ne = CellVariable(mesh=m, hasOld=True, value=0.)
phi_face = phi.faceValue
ne_face = ne.faceValue
x = m.cellCenters[0]
t0 = Variable()
phi.setValue((x-1)**3)
ne.setValue(-6*(x-1))
#____________________________________________________
@cached(cache)
def S(x,t):
    f=6*(x-1)*e**(-t)+54*((x-1)**2)*e**(-2.*t)
    return f
#____________________________________________________
#Boundary Condition:
valueleft_phi=3*e**(-t0)
valueright_phi=6.75*e**(-t0)
valueleft_ne=-6*e**(-t0)
valueright_ne=-6*e**(-t0)
phi.faceGrad.constrain([valueleft_phi], m.facesLeft)
phi.faceGrad.constrain([valueright_phi], m.facesRight)
ne.faceGrad.constrain([valueleft_ne], m.facesLeft)
ne.faceGrad.constrain([valueright_ne], m.facesRight)
#____________________________________________________
eqn0 = DiffusionTerm(1.,var=phi)==ImplicitSourceTerm(-1.,var=ne)
eqn1 = TransientTerm(1.,var=ne) == 
VanLeerConvectionTerm(phi.faceGrad,var=ne)+S(x,t0) 
eqn = eqn0 & eqn1
#____________________________________________________
steps = 1.e4
dt=1.e-4
T=dt*steps
F=dt/(dx**2)
print('F=',F)
#____________________________________________________
vi = Viewer(phi)
with open('out2.txt', 'w') as output:
    while t0()<T:
        print(t0)
        phi.updateOld()
        ne.updateOld()
        res=1.e30
        #for sweep in range(steps):
        while res > 1.e-4:
            res = eqn.sweep(dt=dt)
        t0.setValue(t0()+dt)
        for m in range(nx):
            output.write(str(phi[m])+' ') #+ os.linesep
        output.write('\n')
        if __name__ == '__main__':
            vi.plot()
#____________________________________________________
data = np.loadtxt('out2.txt')
X, T = np.meshgrid(np.linspace(0, L, len(data[0,:])), np.linspace(0, T, 
len(data[:,0])))
fig = plt.figure(3)
ax = fig.add_subplot(111,projection='3d')
ax.plot_surface(X, T, Z=data)
plt.show(block=True)

Answer 1

The issue with these equations, particularly eqn0, is that they admit an infinite number of solutions when Neumann boundary conditions are applied on both boundaries. You can fix this by pinning a value somewhere with an internal fixed value. E.g., based on the analytical solution given in the comments, phi = (x-1)**3 * exp(-t), we can pin phi = 0 at x = 1 with

mask = (m.x > 1-dx/2) & (m.x < 1+dx/2)
largeValue = 1e6
value = 0.
#____________________________________________________
eqn0 = (DiffusionTerm(1.,var=phi)==ImplicitSourceTerm(-1.,var=ne) 
        + ImplicitSourceTerm(mask * largeValue, var=phi) - mask * largeValue * value)

At this point, the solutions still do not agree with the expected solutions. This is because, while you have called ne.faceGrad.constrain() for the left and right boundaries, does not appear in the discretized equations. You can see this if you plot ne; the gradient is zero at both boundaries despite the constraint because FiPy never "sees" the constraint.

What does appear is the flux . By applying fixed flux boundary conditions, I obtain the expected solutions:

ne_left = 6 * numerix.exp(-t0)
ne_right = -9 * numerix.exp(-t0)

eqn1 = (TransientTerm(1.,var=ne) 
        == VanLeerConvectionTerm(phi.faceGrad * m.interiorFaces,var=ne)
        + S(x,t0) 
        + (m.facesLeft * ne_left * phi.faceGrad).divergence
        + (m.facesRight * ne_right * phi.faceGrad).divergence)

You can probably get better convergence properties with

eqn1 = (TransientTerm(1.,var=ne) 
        == DiffusionTerm(coeff=ne.faceValue * m.interiorFaces, var=phi)
        + S(x,t0) 
        + (m.facesLeft * ne_left * phi.faceGrad).divergence
        + (m.facesRight * ne_right * phi.faceGrad).divergence)

but either formulation seems to work.

Note: phi.faceGrad.constrain() is fine, because the flux does appear in DiffusionTerm(coeff=1., var=phi).

Separately, it appears (based on "The Factor is exactly singular") that you are solving with the SciPy LinearLUSolver. The PETSc LinearLUSolver does better, but the baseline value of the solution wanders all over the place. Calling

res = eqn.sweep(dt=dt, solver=LinearGMRESSolver())

also seems to produce stable results (without pinning an internal value). This behavior probably shouldn't be relied on; pinning a value is the right thing to do.