Diffusion in 2D

$$\frac{\partial u}{\partial t} = \nabla\cdot D\nabla x$$
In [1]:
%matplotlib inline
from pylab import *
from numpy import *
In [2]:
#dimensions of the square computational domain:
maxx = 10.
maxt = 1.

#difffusivity
D = 0.5
In [3]:
#discretization parameters
nx = 100   # number of unknown grid points in spatial direction
SC = 1.0  # stability criterion SC = 2*D*dt/dx^2, for FTCS should be SC <= 1

def diffusion_init(maxx, maxt, D, nx, SC, stride=1):
    #choose time step according to CFL condition
    dx = maxx/(nx+1)
    dt = SC*dx**2/(4*D)
    nt = int(maxt/dt)+1
    
    #define array for storing the solution
    U = zeros((int(nt/stride), nx+2, nx+2))
    
    x = arange(nx+2)*dx
    t = arange(nt)*dt
    return U, dx, dt, x, t
In [4]:
def diffusion_solve(U, dx, dt, nx, nt, D, stride=1):
    Ui = U[0,:,:]
    for it in range(0,nt-1):
        Ui[1:-1, 1:-1] = Ui[1:-1, 1:-1] + D*dt/dx**2*(Ui[2:, 1:-1] - 2*Ui[1:-1, 1:-1] + Ui[0:-2, 1:-1])\
            + D*dt/dx**2*(Ui[1:-1,2:] - 2*Ui[1:-1, 1:-1] + Ui[1:-1, 0:-2])
        if it%stride == 0:
            U[(it+1)/stride, :, :] = Ui

Let's try to propagate a square initial condition

In [5]:
U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=0.9)
X, Y = meshgrid(x, x)
U[0, :, :] = 0.
U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1.
pcolormesh(U[0,:,:], rasterized=True, vmin=-1, vmax=1)
Out[5]:
<matplotlib.collections.QuadMesh at 0x7fd398ed5128>
In [7]:
diffusion_solve(U, dx, dt, nx, len(t), D)
pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)
Out[7]:
<matplotlib.collections.QuadMesh at 0x7fd398dd44e0>

The time step limit seems to be to strict, let's try increasing the time step twice

In [8]:
U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=1.1)
X, Y = meshgrid(x, x)
U[0, :, :] = 0.
U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1.
In [9]:
diffusion_solve(U, dx, dt, nx, len(t), D)
pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)
Out[9]:
<matplotlib.collections.QuadMesh at 0x7fd398db1b38>
In [ ]:
from scipy.sparse import dia_matrix, eye, kron
from scipy.sparse.linalg import factorized
def d2matrix(nelem):
    elements = ones((3,nelem))
    elements[1,:] *= -2
    return dia_matrix((elements, [-1,0,1]), shape=(nelem,nelem)).tocsc()

def diffusion_solve_implicit(U, dx, dt, nx, nt, D):
    alpha = D*dt/dx**2
    mat2D = kron(d2matrix(nx), eye(nx)) + kron(eye(nx), d2matrix(nx))
    M = eye(nx*nx)-mat2D*alpha
    LU = factorized(M.tocsc())
    for it in range(0,nt-1):
        U[it+1,1:-1, 1:-1] = LU(U[it,1:-1, 1:-1].ravel()).reshape(nx, nx)
In [ ]:
U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=10)
X, Y = meshgrid(x, x)
U[0, :, :] = 0.
U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1.
In [ ]:
diffusion_solve_implicit(U, dx, dt, nx, len(t), D)
pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)

Again, Crank-Nicolson follows:

In [ ]:
def diffusion_solve_CN(U, dx, dt, nx, nt, D):
    alpha2 = D*dt/dx**2*0.5
    mat2D = kron(d2matrix(nx), eye(nx)) + kron(eye(nx), d2matrix(nx))
    M1 = eye(nx*nx)-mat2D*alpha2
    M2 = eye(nx*nx)+mat2D*alpha2
    LU = factorized(M1.tocsc())
    for it in range(0,nt-1):
        U[it+1,1:-1, 1:-1] = LU(M2.dot(U[it,1:-1, 1:-1].ravel())).reshape(nx, nx)
In [ ]:
diffusion_solve_CN(U, dx, dt, nx, len(t), D)
pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)

Neumann BC

In [ ]:
def diffusion_solve_CN_neumann_zero(U, dx, dt, nx, nt, D):
    alpha = D*dt/dx**2
    d2x = d2matrix(nx+2)
    d2x[0, 1] += 1
    d2x[-1, -2] += 1
    mat2D = kron(d2x, eye(nx+2)) + kron(eye(nx+2), d2x)
    M1 = eye((nx+2)**2)-mat2D*alpha/2
    M2 = eye((nx+2)**2)+mat2D*alpha/2
    LU = factorized(M1.tocsc())
    for it in range(0,nt-1):
        # U[it, boundary] -= alpha*dx*(NBC[it]+NBC[it+1])
        U[it+1,:, :] = LU(M2.dot(U[it,:, :].ravel())).reshape(nx+2, nx+2)
In [ ]:
U, dx, dt, x, t = diffusion_init(maxx, maxt*50, D, nx, SC=50)
X, Y = meshgrid(x, x)
U[0, :, :] = 0.
U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1.
In [ ]:
diffusion_solve_CN_neumann_zero(U, dx, dt, nx, len(t), D)
In [ ]:
pcolormesh(U[:, 30, :], rasterized=True, vmin=0, vmax=.1)
In [ ]:
# constants taken from http://en.wikipedia.org/wiki/File:Brusselator_space.gif
#initial field values
U0 = 1.
V0 = 1.
noise = 2.

#diffusion coefficients
DU = 0.2
DV = 0.02

#constant concentrations
A = 1
B = 3

def RHS(U, V):
    fV = B*U - U**2*V
    fU = A - fV - U
    return fU, fV

maxt = 100.
maxx = 100.

U, dx, dt, x, t = diffusion_init(maxx, maxt, max((DU, DV)), nx, SC=.05)
V, dx, dt, x, t = diffusion_init(maxx, maxt, max((DU, DV)), nx, SC=.05)
X, Y = meshgrid(x, x)
U[0, :, :] = U0 + noise*randn(nx+2, nx+2)
V[0, :, :] = V0 + noise*randn(nx+2, nx+2)
#for i in range(10):
#    U[0, randint(0, nx), randint(0, nx)] = 10
#    V[0, randint(0, nx), randint(0, nx)] = 10
In [ ]:
def brusselator(U, V, dx, dt, nx, nt, DU, DV):
    alphaU = DU*dt/dx**2
    alphaV = DV*dt/dx**2
    d2x = d2matrix(nx+2)
    d2x[0, 1] += 1
    d2x[-1, -2] += 1
    mat2D = kron(d2x, eye(nx+2)) + kron(eye(nx+2), d2x)
    M1U = eye((nx+2)**2)-mat2D*alphaU/2
    M2U = eye((nx+2)**2)+mat2D*alphaU/2
    M1V = eye((nx+2)**2)-mat2D*alphaV/2
    M2V = eye((nx+2)**2)+mat2D*alphaV/2
    LUU = factorized(M1U.tocsc())
    LUV = factorized(M1V.tocsc())
    for it in range(10):
        U[0] = M2V.dot(U[it,:, :].ravel()).reshape(nx+2, nx+2)
    for it in range(0,nt-1):
        # U[it, boundary] -= alpha*dx*(NBC[it]+NBC[it+1])
        # diffusion operator
        U[it+1,:, :] = LUU(M2U.dot(U[it,:, :].ravel())).reshape(nx+2, nx+2)
        V[it+1,:, :] = LUV(M2V.dot(V[it,:, :].ravel())).reshape(nx+2, nx+2)
        
        fU, fV = RHS(U[it+1], V[it+1])
        U[it+1] += fU*dt
        V[it+1] += fV*dt
In [ ]:
brusselator(U, V, dx, dt, nx, len(t), DU, DV)
In [ ]:
figure(figsize=(12,12))
subplot(221)
pcolormesh(U[:, 30, :], rasterized=True, vmin=0, vmax=5)
subplot(223)
pcolormesh(V[:, 30, :], rasterized=True, vmin=0, vmax=5)
subplot(222)
pcolormesh(U[1400, :-2, :-2], rasterized=True, vmin=0, vmax=5)
subplot(224)
pcolormesh(V[1400, :-2, :-2], rasterized=True, vmin=0, vmax=5)
In [ ]: