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=nx*dx
e=math.e
m = Grid1D(nx=nx, dx=dx)
print(np.log(e))
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()
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=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=-6e(-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)
+ImplicitSourceTerm(ne,var=ne) , DiffusionTerm(coeff=D, var=ne)
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)