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=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()
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=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=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=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)