Fipy gives incorrect results when applying Neumann’s conditions on one boundary

[z-soltani]

The results diverge after changing one boundary's condition from Dirichlet to Neumann.

We are solving a one-dimensional Coupled Continuity-Poisson equation, when applying Dirichlet's conditions on both boundaries, it gives the right result (which matches with the exact result), but when changing one of the boundary's conditions to Neumann, the result diverges after spending 0.2 of its given time. The code is as below (0<x<2.5). Any help is highly appreciated..

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) #Neumann valueright_phi=(1.53)e(-t0) # Dirichlet valueleft_ne=-6*e*(-t0) #Neumann valueright_ne=-9e**(-t0) # Dirichlet phi.faceGrad.constrain([valueleft_phi], m.facesLeft) phi.constrain([valueright_phi], m.facesRight) ne.faceGrad.constrain([valueleft_ne], m.facesLeft) ne.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)