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

[z-soltani]

Code not working when applying Neumann boundary conditions on both sides

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=nxdx 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=3e*(-t0) valueright_phi=6.75e**(-t0) valueleft_ne=-6e**(-t0) valueright_ne=-6e**(-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) #+ImplicitSourceTerm(ne,var=ne) , DiffusionTerm(coeff=D, var=ne) eqn = eqn0 & eqn1 #____________________________________________________ steps = 1.e4 dt=1.e-4 T=dtsteps 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)