2D Brusselator PDE

The Brusselator PDE:

\begin{align} \frac{\partial u}{\partial t} &= 1 + u^2v - 4.4u + \alpha (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) + f(x, y, t) \\ \frac{\partial v}{\partial t} &= 3.4u - u^2 v + \alpha (\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}) \end{align}

(1)

where

f(x, y, t) = \begin{cases} 5 \qquad \text{if} (x - 0.3)^2 + (y - 0.6)^2 \leq 0.1^2 \ \text{and} \ t \geq 1.1 \\ 0 \qquad \text{otherwise} \end{cases}

(2)

and the initial conditions are

\begin{align} u(x, y, 0) &= 22(y(1-y))^{1.5} \\ v(x, y, 0) &= 27(x(1-x))^{1.5} \end{align}

(3)

with the periodic boundary condition

\begin{align} u(x+1, y, 0) &= u(x, y, t) \\ u(x, y+1, 0) &= u(x, y, t) \end{align}

(4)

on a time span of $t \in [0, 11.5]$ .

using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets
using Plots
using OrdinaryDiffEqBDF

Setup parameters, variables, and differential operators

function model_bru(; x_min=0.0, x_max=1.0, y_min=0.0, y_max=1.0, t_min=0.0, t_max=11.5, α=10.0, N::Int = 20, order::Int = 2)
    @independent_variables x y t
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dy = Differential(y)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    ∇²(u) = Dxx(u) + Dyy(u)

    # Dynamics on each grid point
    brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5
    u0(x, y, t) = 22 * (y * (1 - y))^(3 / 2)
    v0(x, y, t) = 27 * (x * (1 - x))^(3 / 2)

    eqs = [
        Dt(u(x, y, t)) ~ 1.0 + v(x, y, t) * u(x, y, t)^2 - 4.4 * u(x, y, t) + α * ∇²(u(x, y, t)) + brusselator_f(x, y, t),
        Dt(v(x, y, t)) ~ 3.4 * u(x, y, t) - v(x, y, t) * u(x, y, t)^2 + α * ∇²(v(x, y, t))
    ]
    domains = [
        x ∈ Interval(x_min, x_max),
        y ∈ Interval(y_min, y_max),
        t ∈ Interval(t_min, t_max)
    ]
    # Periodic boundary conditions
    bcs = [
        u(x, y, 0) ~ u0(x, y, 0),
        u(0, y, t) ~ u(1, y, t),
        u(x, 0, t) ~ u(x, 1, t),
        v(x, y, 0) ~ v0(x, y, 0),
        v(0, y, t) ~ v(1, y, t),
        v(x, 0, t) ~ v(x, 1, t)
    ]
    @named pdesys = PDESystem(eqs, bcs, domains, [x, y, t], [u(x, y, t), v(x, y, t)])
    disc = MOLFiniteDifference([x=>N, y=>N], t; approx_order=order)
    prob = discretize(pdesys, disc)
    return (; prob, x, y, t, u, v)
end

@time "Build problem" prob, x, y, t, u, v = model_bru()
@time "Solve problem" sol = solve(prob, DFBDF(), saveat=0.1);

┌ Warning: The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature.
└ @ MethodOfLines ~/.julia/packages/MethodOfLines/0m4Gg/src/system_parsing/pde_system_transformation.jl:55
Build problem: 81.661626 seconds (147.41 M allocations: 7.873 GiB, 2.68% gc time, 99.35% compilation time: 7% of which was recompilation)
Solve problem: 66.777245 seconds (39.00 M allocations: 2.688 GiB, 1.84% gc time, 96.40% compilation time: <1% of which was recompilation)

Extract data

discrete_x = sol[x]
discrete_y = sol[y]
discrete_t = sol[t]

solu = sol[u(x, y, t)]
solv = sol[v(x, y, t)]

umax = maximum(maximum, solu)
vmax = maximum(maximum, solv)

4.930747443605161

Visualization¶

Interval == 2:end since in periodic condition, end == 1.

anim = @animate for k in eachindex(discrete_t)
    heatmap(solu[2:end, 2:end, k], title="u @ t=$(discrete_t[k])", clims = (0.0, umax))
end

mp4(anim, fps = 8)

[ Info: Saved animation to /tmp/jl_n0LEbAdIlS.mp4

Loading...

anim = @animate for k in eachindex(discrete_t)
    heatmap(solv[2:end, 2:end, k], title="v @ t=$(discrete_t[k])", clims = (0.0, vmax))
end

mp4(anim, fps = 8)

[ Info: Saved animation to /tmp/jl_twDOn8RCNy.mp4

Loading...

This notebook was generated using Literate.jl.