wiki:CookBook/BurgersEquation

The following examples were originally provided to us by fred2 <at> qnet.com in ​http://thread.gmane.org/gmane.comp.python.fipy/911

The linearized Burgers equation,

$$\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} - \mu \frac{\partial^2 u}{\partial x^2} = 0,$$
is a simple model equation for gas dynamics (Anderson, Tannehill, Pletcher [Anderson:1984]_, chapter 4). If $\mu = 0$, the wave equation is recovered. For the boundary conditions
$$u(0,t) = u_0$$
$$u(L,t) = 0,$$
the exact steady-state solution is
$$u = u_0 \left ( \frac{1-\exp{[R_L(x/L - 1)]}}{1-\exp{(-R_L)}} \right ),$$
where $R_L = cL/\mu$.

The following script integrates the equation (with $c=1/2$, $\mu=1/4$, and $L=1$) to large time. 

from fipy import *

mu = 0.25; c = 0.5; reyn = c/mu; u0 = 1.
nx = 50; dx = 1./nx

mesh = Grid1D(nx = nx, dx = dx)
x, = mesh.getCellCenters()

ux = CellVariable(name="ux", mesh=mesh, value=0.5*u0, hasOld=1, rank=0)
uxexact = CellVariable(name="uexact", mesh=mesh, value=0., rank=0)

expterm = numerix.exp(reyn*(x - 1.))
uxexact.setValue(u0*(1. - expterm) /(1. - numerix.exp(-reyn)))
ux.setValue(u0*(1. - x))
BC=(FixedValue(faces=mesh.getFacesLeft(), value=u0),FixedValue(faces=mesh.getFacesRight(), value=0))
cfl = 0.5; dt = cfl*dx/u0
nstep = 200
eqn =  TransientTerm() + VanLeerConvectionTerm(coeff=c) - ImplicitDiffusionTerm(coeff=mu)
rtol = 1.e-7

ttl = '''linearized unsteady burgers eqn
c=%g, mu=%3.2f, nx=%g''' % (c, mu, nx)

viewer = Viewer(vars=(ux,uxexact),title=ttl)
time = 0.
for step in range(nstep):
    time += dt
    ux.updateOld()
    resid = 1.e10
    while resid > rtol:
       resid = eqn.sweep(var = ux, dt = dt, 
                         solver = LinearLUSolver(),boundaryConditions=BC)
    viewer.plot()

For periodic boundary conditions, and the initial condition

$$u(x,0) = \sin{kx}$$
the exact unsteady solution is
$$u(x,t) = \exp{(-k^2\mu t)} \sin{[k(x-ct)]}.$$
The following script integrates the equation (with $c=1/2$, $\mu=1/4$, and $k=2\pi$) to moderate time on the interval $x\in [0,1]$. 

mu = 0.25; c = 0.5; reyn = c/mu
k = 2.*numerix.pi
nx = 200; dx = 1./nx
from fipy import *

mesh = PeriodicGrid1D(nx = nx, dx = dx)
x, = mesh.getCellCenters()

ux = CellVariable(name="ux", mesh=mesh, value=0., hasOld=1, rank=0)
uxexact = CellVariable(name="uexact", mesh=mesh, value=0., rank=0)

time = Variable()

ux.setValue(numerix.sin(k*x))
uxexact.setValue(numerix.exp(-k*k*mu*time) * numerix.sin(k*(x-c*time)))
cfl = 0.1; dt = cfl*dx/ux.value.max()
eqn =  TransientTerm() + VanLeerConvectionTerm(coeff=c) - ImplicitDiffusionTerm(coeff=mu)
rtol = 1.e-7

ttl = '''linearized unsteady burgers eqn
c=%g, mu=%3.2f, k/pi=%g, nx=%''' % (c, mu, k/numerix.pi, nx)

viewer = Viewer(vars=(ux,uxexact),
                title=ttl,
                limits={'datamin': -1.0, 'datamax': 1.0})

time.setValue(0.)
print 'time/timefinal, ux.allclose(uxexact, atol=1.e-3), max(|ux-uxexact|), max(|ux-uex|)/max(|ux|) '
timef = 0.2
while time() < timef:
    time.setValue(time() + dt)
    uxexact.setValue(numerix.exp(-k*k*mu*time) * numerix.sin(k*(x-c*time)))
    ux.updateOld()
    resid = 1.e10
    while resid > rtol:
       resid = eqn.sweep(var = ux, dt = dt, 
                   solver = LinearLUSolver())
    if (time() % (timef/20.)) < dt:
        #note, atol isn't too restrictive when ux decays to < 0.1!
        udif = abs(ux - uxexact).max()
        umax = abs(ux).max()
        print time.value, ux.allclose(uxexact, atol=1.e-3), udif, udif/umax
    viewer.plot()

The following script gives a solution to the steady-state, non-linear Burgers equation. 

from fipy import *

mu = 0.25; b = -1.; c = 0.5; u0 = 1.
#the exact stationary soln (4-169) requires these b,c ...
assert(b == -1.)
assert(c == 0.5)
#(& fig 4-37 requires this mu)
assert(mu == 0.25)
#fig 4-37 abscissa uses indep var x-x0
#sorry, cannot figure out origin in fipy numMesh/?
#make sure & use x-x0 instead of x as approp 
xrange = 7.; x0 = xrange/2
nx = 350; dx = xrange/nx
mesh = Grid1D(nx = nx, dx = dx)
x = mesh.getCellCenters()[0]
#
cfl = 0.5; dt = cfl*dx/u0
#arbitrary (large) steady state time
timef = 40.
rtol = 1.e-7
atol = 1.e-3

#fipy convterm var u
u = CellVariable(name="u", mesh=mesh, value=0.5*u0, hasOld=1, rank=1)
#depvar ux
ux = CellVariable(name="ux", mesh=mesh, value=0.5*u0, hasOld=1, rank=0)
uxexact = CellVariable(name="uexact", mesh=mesh, value=0., rank=0)
#
#eq (4-169) note x-x0
uxexact.setValue(-c / b * (1 + numerix.tanh(c*(x-x0)/(2.*mu))))
#
#conservation checks
cellvol = mesh.getCellVolumes()
utmp = uxexact.value
usumex = sum(map(lambda x,y: x*y, cellvol, utmp))
usqsumex = sum(map(lambda x,y: x*y*y, cellvol, utmp))
#
#simple linear gradient ic
ux.setValue(u0*(x/xrange))
#... or the exact soln 
#ux.setValue(uxexact.getValue())
#
#fixed bc at left & right ends, based on exact soln (4-169), applied at
#face
xl  = mesh.getFaceCenters()[0][0]
uleft = -c / b * (1 + numerix.tanh((c*(xl-x0))/(2.*mu)))
xr  = mesh.getFaceCenters()[0][-1]
uright = -c / b * (1 + numerix.tanh((c*(xr-x0))/(2.*mu)))

BC=(FixedValue(faces=mesh.getFacesLeft(), value=uleft), 
    FixedValue(faces=mesh.getFacesRight(), value=uright))

#generalized burgers eqn (4-168)
eqn =  TransientTerm() + VanLeerConvectionTerm(coeff=c+0.5*b*u) - \
       ImplicitDiffusionTerm(coeff=mu)

ttl = '''generalized burgers eqn
b=%g, c=%g, mu=%3.2f, nx=%g''' % (b, c, mu, nx)
viewer = viewers.make(vars=(ux,uxexact),title=ttl)
print '\nmomentum & energy conservation: \n' \
      'time\t sum(udV)/sumex\t sum(u^2dV)/sumex'
time = 0.
while time <= timef+dt:
    if (time % (timef/20.)) < dt: 
        utmp = ux.value
        usum = sum(map(lambda x,y: x*y, cellvol, utmp))
        usqsum = sum(map(lambda x,y: x*y*y, cellvol, utmp))
        print '%5.2f\t %8.5f\t %8.5f' % \
              (time/timef, usum/usumex, usqsum/usqsumex) 
        viewer.plot()
    time += dt
    ux.updateOld()
    resid = 1.e10
    #for sweep in range(sweeps):
    while resid > rtol:
       u[0] = ux
       resid = eqn.sweep(var = ux, dt = dt, 
                   solver = LinearLUSolver(),boundaryConditions=BC)

print 'step, ux.allclose(uxexact, atol=%g), ' \
      'max(|ux-uxexact|), max(|ux-uex|)/max(|ux|) ' %(atol)
udif = abs(ux - uxexact).max()
umax = abs(ux).max()
print int(time/dt), ux.allclose(uxexact, atol=atol), udif, udif/umax

Benton and Platzman [Benton:1972]_ describe 35 different exact solutions to the Burgers equation. 

@book{Anderson:1984,
   author =    {D. A. Anderson, J. C. Tannehill, and R. H. Pletcher}
   title =     {Computational Fluid Mechanics and Heat Transfer},
   year =      1984,
   publisher = {McGraw Hill},
}
@article{Benton:1972,
   author =  {E. R. Benton and G. W. Platzman},
   title =   {A table of solutions of the one-dimensional {B}urgers equation},
   journal = {Q. Appl. Math.},
   volume =  30,
   year =    1972,
   pages =   {195-212}
}
Last modified 5 years ago Last modified on 03/30/10 10:23:10

Attachments (3)

Download all attachments as: .zip