2D Brusselator PDE

2D Brusselator PDE#

The Brusselator PDE:

\[\begin{split} \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} \end{split}\]

where

\[\begin{split} 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} \end{split}\]

and the initial conditions are

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

with the periodic boundary condition

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

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

using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets
using Plots

Setup parameters, variables, and differential operators

@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)

∇² (generic function with 1 method)

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

x_min = y_min = t_min = 0.0
x_max = y_max = 1.0
t_max = 11.5
α = 10.0

u0(x, y, t) = 22 * (y * (1 - y))^(3 / 2)
v0(x, y, t) = 27 * (x * (1 - x))^(3 / 2)

v0 (generic function with 1 method)

PDEs

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))
]

\[\begin{split} \begin{align} \frac{\mathrm{d}}{\mathrm{d}t} u\left( x, y, t \right) &= 1 + 10 \left( \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( x, y, t \right) + \frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} u\left( x, y, t \right) \right) - 4.4 u\left( x, y, t \right) + 5 \left( \left( -0.3 + x \right)^{2} + \left( -0.6 + y \right)^{2} \leq 0.01 \right) \left( t \geq 1.1 \right) + \left( u\left( x, y, t \right) \right)^{2} v\left( x, y, t \right) \\ \frac{\mathrm{d}}{\mathrm{d}t} v\left( x, y, t \right) &= 3.4 u\left( x, y, t \right) + 10 \left( \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} v\left( x, y, t \right) + \frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} v\left( x, y, t \right) \right) - \left( u\left( x, y, t \right) \right)^{2} v\left( x, y, t \right) \end{align} \end{split}\]

Space and time domains

domains = [
    x ∈ Interval(x_min, x_max),
    y ∈ Interval(y_min, y_max),
    t ∈ Interval(t_min, t_max)
]

3-element Vector{Symbolics.VarDomainPairing}:
 Symbolics.VarDomainPairing(x, 0.0 .. 1.0)
 Symbolics.VarDomainPairing(y, 0.0 .. 1.0)
 Symbolics.VarDomainPairing(t, 0.0 .. 11.5)

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)
]

\[\begin{split} \begin{align} u\left( x, y, 0 \right) &= 22 \left( 1 - y \right)^{1.5} y^{1.5} \\ u\left( 0, y, t \right) &= u\left( 1, y, t \right) \\ u\left( x, 0, t \right) &= u\left( x, 1, t \right) \\ v\left( x, y, 0 \right) &= 27 x^{1.5} \left( 1 - x \right)^{1.5} \\ v\left( 0, y, t \right) &= v\left( 1, y, t \right) \\ v\left( x, 0, t \right) &= v\left( x, 1, t \right) \end{align} \end{split}\]

PDE system

@named pdesys = PDESystem(eqs, bcs, domains, [x, y, t], [u(x, y, t), v(x, y, t)])

\[\begin{split} \begin{align} \frac{\mathrm{d}}{\mathrm{d}t} u\left( x, y, t \right) &= 1 + 10 \left( \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} u\left( x, y, t \right) + \frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} u\left( x, y, t \right) \right) - 4.4 u\left( x, y, t \right) + 5 \left( \left( -0.3 + x \right)^{2} + \left( -0.6 + y \right)^{2} \leq 0.01 \right) \left( t \geq 1.1 \right) + \left( u\left( x, y, t \right) \right)^{2} v\left( x, y, t \right) \\ \frac{\mathrm{d}}{\mathrm{d}t} v\left( x, y, t \right) &= 3.4 u\left( x, y, t \right) + 10 \left( \frac{\mathrm{d}}{\mathrm{d}x} \frac{\mathrm{d}}{\mathrm{d}x} v\left( x, y, t \right) + \frac{\mathrm{d}}{\mathrm{d}y} \frac{\mathrm{d}}{\mathrm{d}y} v\left( x, y, t \right) \right) - \left( u\left( x, y, t \right) \right)^{2} v\left( x, y, t \right) \end{align} \end{split}\]

Discretization to an ODE system

disc = let N = 32, order = 2
    MOLFiniteDifference([x=>N, y=>N], t, approx_order=order)
end

prob = discretize(pdesys, disc)

┌ 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/EiZyI/src/system_parsing/pde_system_transformation.jl:43

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 11.5)
u0: 1922-element Vector{Float64}:
12134432813715874
12134432813715874
12134432813715874
12134432813715874
12134432813715874
12134432813715874
12134432813715874
12134432813715874
12134432813715874
12134432813715874
 ⋮
262066755474226
9735248771761977
6650051119921914
3433640166822105
0172186542655526
6977464117458192
4003323380813969
14892258453196736
0

Solvers: https://diffeq.sciml.ai/stable/solvers/ode_solve/

alg = TRBDF2()
@time sol = solve(prob, alg, saveat=0.1)

262.535893 seconds (508.55 M allocations: 25.780 GiB, 4.33% gc time, 92.56% compilation time)

retcode: Success
Interpolation: Dict{Symbolics.Num, Interpolations.GriddedInterpolation{Float64, 3, Array{Float64, 3}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}}}}
t: 116-element Vector{Float64}:
  0.0
  0.1
  0.2
  0.3
  0.4
  0.5
  0.6
  0.7
  0.8
  0.9
  ⋮
 10.7
 10.8
 10.9
 11.0
 11.1
 11.2
 11.3
 11.4
 11.5ivs: 3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 t
 x
 ydomain:([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9  …  10.6, 10.7, 10.8, 10.9, 11.0, 11.1, 11.2, 11.3, 11.4, 11.5], 0.0:0.03225806451612903:1.0, 0.0:0.03225806451612903:1.0)
u: Dict{Symbolics.Num, Array{Float64, 3}} with 2 entries:
  u(x, y, t) => [0.0 0.121344 … 0.121344 0.0; 0.0 0.121344 … 0.121344 0.0; … ; …
  v(x, y, t) => [0.0 0.0 … 0.0 0.0; 0.148923 0.148923 … 0.148923 0.148923; … ; …

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.966509417157691

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 /home/runner/work/jl-pde/jl-pde/.cache/docs/tmp.mp4

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 /home/runner/work/jl-pde/jl-pde/.cache/docs/tmp.mp4

This notebook was generated using Literate.jl.

2D Brusselator PDE

Contents

2D Brusselator PDE#

Visualization#