Hi FiPy,
I've been solving transient solutions on a uniform 2D grid for a
convection-diffusion equation that initially centers around some arbitrary
position, e.g., (1,1), and moves towards the origin. However, I just
discovered that if I inverse the process and simulates a movement starting
from the origin towards (1,1), with everything else kept the same, the
solution blows up. Sample scripts are attached here. Note that only lines
31 & 62 are different between the two.
How could this be?
I also tried reducing the time step increment in the unstable example but
without any success..
Any suggestion is greatly appreciated.
Cheers,
Yun
--
Graduate Group of Ecology Doctoral Candidate
Department of Environmental Science and Policy
Center for Population Biology
University of California, Davis
#!/usr/bin/env python
# encoding: utf-8
"""
bivariate diffusion and advection
Created by Yun Tao on 2012-05-14.
Copyright (c) 2012 __MyCompanyName__. All rights reserved.
"""
import matplotlib.pyplot as plt
from matplotlib import colors
from numpy.random import normal
from matplotlib.mlab import bivariate_normal
from fipy import *
from fipy import numerix
fig = plt.figure
ax1 = plt.subplot((111))
# mesh config
l = 8.
nx = 201.
ny = nx
dx = l/nx
dy = dx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2]]
convection = VanLeerConvectionTerm
# initial condition
x, y = mesh.getCellCenters()
z = bivariate_normal(x, y, .5, .5, 1., 1.)
# variable config
phi = CellVariable(mesh=mesh, value=z)
# boundary conditions
facesTopRight = ((mesh.getFacesRight())
| (mesh.getFacesTop()))
facesBottomLeft = ((mesh.getFacesLeft())
| (mesh.getFacesBottom()))
BCs = (FixedFlux(faces=facesTopRight, value=0.),
FixedFlux(faces=facesBottomLeft, value=0.))
# viewer config
viewer1 = MatplotlibViewer(vars=phi,
datamin = 0.0,
limits={'xmin': -l/2, 'xmax': l/2, 'ymin': -l/2, 'ymax': l/2},
cmap = plt.cm.Spectral,
axes=ax1)
# parameters
D = 1.0
c = 1.0
b = ((c,), (c,))
alpha = 2./(dx * 0.5)
pos = FaceVariable(mesh=mesh, value=mesh.getFaceCenters(), rank=1)
xv, yv = mesh.getFaceCenters()
r = numerix.sqrt((xv)**2+(yv)**2)
# define equation
eq3 = TransientTerm() == DiffusionTerm(coeff=D) + convection(coeff=(b*numerix.tanh(alpha*r)*pos/r))
# solver config
timeStepDuration = 1. * 1./2 * dx / c
steps = 10#
for step in range(steps):
eq3.solve(var=phi,
boundaryConditions = BCs,
dt=timeStepDuration)
viewer1.plot()
#!/usr/bin/env python
# encoding: utf-8
"""
bivariate diffusion and advection
Created by Yun Tao on 2012-05-14.
Copyright (c) 2012 __MyCompanyName__. All rights reserved.
"""
import matplotlib.pyplot as plt
from matplotlib import colors
from numpy.random import normal
from matplotlib.mlab import bivariate_normal
from fipy import *
from fipy import numerix
fig = plt.figure
ax1 = plt.subplot((111))
# mesh config
l = 8.
nx = 201.
ny = nx
dx = l/nx
dy = dx
mesh = Grid2D(nx=nx, ny=ny, dx=dx, dy=dy) + [[-l/2]]
convection = VanLeerConvectionTerm
# initial condition
x, y = mesh.getCellCenters()
z = bivariate_normal(x, y, .5, .5, 0., 0.)
# variable config
phi = CellVariable(mesh=mesh, value=z)
# boundary conditions
facesTopRight = ((mesh.getFacesRight())
| (mesh.getFacesTop()))
facesBottomLeft = ((mesh.getFacesLeft())
| (mesh.getFacesBottom()))
BCs = (FixedFlux(faces=facesTopRight, value=0.),
FixedFlux(faces=facesBottomLeft, value=0.))
# viewer config
viewer1 = MatplotlibViewer(vars=phi,
datamin = 0.0,
limits={'xmin': -l/2, 'xmax': l/2, 'ymin': -l/2, 'ymax': l/2},
cmap = plt.cm.Spectral,
axes=ax1)
# parameters
D = 1.0
c = 1.0
b = ((c,), (c,))
alpha = 2./(dx * 0.5)
pos = FaceVariable(mesh=mesh, value=mesh.getFaceCenters(), rank=1)
xv, yv = mesh.getFaceCenters()
r = numerix.sqrt((xv-1.)**2+(yv-1.)**2)
# define equation
eq3 = TransientTerm() == DiffusionTerm(coeff=D) + convection(coeff=(b*numerix.tanh(alpha*r)*pos/r))
# solver config
timeStepDuration = 1. * 1./2 * dx / c
steps = 10
for step in range(steps):
eq3.solve(var=phi,
boundaryConditions = BCs,
dt=timeStepDuration)
viewer1.plot()
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
