Solution stability depends on convectional target?

[Yun Tao]

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
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
  [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]