brusselator.jl

using OrdinaryDiffEq, Plots
plotly()
import GR: meshgrid

const N = 64
const xd, yd = linspace(0,1,N), linspace(0,1,N)
function brusselator_loop(du, u, p, t)
    @inbounds begin
        A, B, α = p
        # Interior
        for i in 2:N-1, j in 2:N-1
            du[i,j,1] = α*(u[i-1,j,1] + u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
            B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        end
        for i in 2:N-1, j in 2:N-1
            du[i,j,2] = α*(u[i-1,j,2] + u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
            A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
        end

        # Boundary @ edges
        for j in 2:N-1
            i = 1
            du[1,j,1] = α*(2u[i+1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
            B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        end
        for j in 2:N-1
            i = 1
            du[1,j,2] = α*(2u[i+1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
            A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
        end
        for j in 2:N-1
            i = N
            du[end,j,1] = α*(2u[i-1,j,1] + u[i,j+1,1] + u[i,j-1,1] - 4u[i,j,1]) +
            B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        end
        for j in 2:N-1
            i = N
            du[end,j,2] = α*(2u[i-1,j,2] + u[i,j+1,2] + u[i,j-1,2] - 4u[i,j,2]) +
            A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
        end
        for i in 2:N-1
            j = 1
            du[i,1,1] = α*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
            B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        end
        for i in 2:N-1
            j = 1
            du[i,1,2] = α*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
            A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
        end
        for i in 2:N-1
            j = N
            du[i,end,1] = α*(u[i-1,j,1] + u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
            B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        end
        for i in 2:N-1
            j = N
            du[i,end,2] = α*(u[i-1,j,2] + u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
            A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
        end

        # Boundary @ four vertexes
        i = 1; j = 1
        du[1,1,1] = α*(2u[i+1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
        B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        du[1,1,2] = α*(2u[i+1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
        A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]

        i = 1; j = N
        du[1,N,1] = α*(2u[i+1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
        B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        du[1,N,2] = α*(2u[i+1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
        A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]

        i = N; j = 1
        du[N,1,1] = α*(2u[i-1,j,1] + 2u[i,j+1,1] - 4u[i,j,1]) +
        B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        du[N,1,2] = α*(2u[i-1,j,2] + 2u[i,j+1,2] - 4u[i,j,2]) +
        A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]

        i = N; j = N
        du[end,end,1] = α*(2u[i-1,j,1] + 2u[i,j-1,1] - 4u[i,j,1]) +
        B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        du[end,end,2] = α*(2u[i-1,j,2] + 2u[i,j-1,2] - 4u[i,j,2]) +
        A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
    end
end

function uv0_fun()
    u = zeros(N, N, 2)
    x, y = meshgrid(xd, yd)
    for i in 1:N, j in 1:N
        u[i,j,1] = (2 + 0.25y[i,j])
        u[i,j,2] = (1 + 0.8x[i,j])
    end
    u
end

const u0 = uv0_fun()
const dx = step(xd)
function laplace_operator(N)
    e1 = ones(N-1)
    e1[end] = 2
    T = spdiagm((e1, -2ones(N), reverse(e1)), (-1,0,1)) ./ dx^2
    # A ⊗ B * x == B * x * A'
    T
end

const D = laplace_operator(N)
const Dt= D'
const buffer = similar(u0)
function brusselator_oper(du, u, p, t)
    @inbounds begin
        A, B, α = p
        A_mul_B!(@view(buffer[:,:,1]), D, @view(u[:,:,1]))
        A_mul_B!(@view(buffer[:,:,2]), D, @view(u[:,:,2]))
        A_mul_B!(@view(du[:,:,1]), @view(u[:,:,1]), Dt)
        A_mul_B!(@view(du[:,:,2]), @view(u[:,:,2]), Dt)
        for i in 1:N, j in 1:N
            du[i,j,1] = α*(du[i,j,1] + buffer[i,j,1]) +
            B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1]
        end
        for i in 1:N, j in 1:N
            du[i,j,2] = α*(du[i,j,2] + buffer[i,j,2]) +
            A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
        end
    end
end

function benchloop(t)
    prob_loop = ODEProblem(brusselator_loop, u0, (0.,t), (3.4, 1., 0.002/dx^2))
    sol_loop = solve(prob_loop, BS3(), save_at=5.)
end

function benchoper(t)
    prob_oper = ODEProblem(brusselator_oper, u0, (0.,t), (3.4, 1., 0.002))
    sol_oper = solve(prob_oper, BS3(), save_at=t)
end

sol_loop = benchloop(7)
sol_oper = benchoper(7)
surface(xd, yd, sol_loop(7)[:,:,1], color=:darktest)
surface(xd, yd, sol_oper(7)[:,:,1], color=:darktest)




#=

plots = [surface(xd, yd, sol(i)[:,:,1], color=:darktest, zlims = (0,6)) for i in 0:.05:15]
film = roll(plots, fps=40)
write("output.mp4", film)

=#

########################################################
# 1D
########################################################
using DiffEqOperators, OrdinaryDiffEq, Plots
N = 64
D = DerivativeOperator{Float64}(2,2,1/(N-1),N,:Neumann0,:Neumann0)
function brusselator_f2(du, u, p, t)
    A, B = p
    @. du[:, 1] = B + u[:,1]^2*u[:, 2] - (A+1)*u[:,1]
    @. du[:, 2] = A*u[:,1] - u[:,1]^2*u[:, 2]
end
x = linspace(0, 1, N)
u0 = @. sech(x) - sech(5)
v0 = @. sech(x) - sech(5)
init = hcat(u0, v0)
prob = SplitODEProblem(DiffEqArrayOperator((0.002full(D))), brusselator_f2, init, (0,10.), (3.4, 1.))
sol = solve(prob, LawsonEuler(), dt=0.01)
plot([sol(i)[:,1] for i in 0:1:10])

function brusselator(du, u, p, t)
    A, B, α = p
    Du = D*u[:,1]
    Dv = D*u[:,2]
    @. du[:, 1] = B + u[:,1]^2*u[:, 2] - (A+1)*u[:,1] + α*Du
    @. du[:, 2] = A*u[:,1] - u[:,1]^2*u[:, 2] + α*Dv
end
prob = ODEProblem(brusselator, init, (0,10.), (3.4, 1., 0.002))
sol = solve(prob, BS3())
plot!([sol(i)[:,1] for i in 0:1:10])