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

@time disc = let N = 20, order = 2
    dx = 1 / N
    MOLFiniteDifference([x=>dx, y=>dx], t; approx_order=order)
end

  0.031708 seconds (26.74 k allocations: 1.546 MiB, 99.09% compilation time)

MethodOfLines.MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}(Dict{Symbolics.Num, Float64}(y => 0.05, x => 0.05), t, 2, MethodOfLines.UpwindScheme(1), MethodOfLines.CenterAlignedGrid(), true, false, MethodOfLines.ScalarizedDiscretization(), true, Any[], Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())

convert the PDESystem in to an ODEProblem or a NonlinearProblem

@time 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/JDu9X/src/system_parsing/pde_system_transformation.jl:43
131.603393 seconds (814.85 M allocations: 25.907 GiB, 2.62% gc time, 3.66% compilation time: 19% of which was recompilation)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 11.5)
u0: 800-element Vector{Float64}:
0
27951439475454604
7289999999999999
229218368262938
7279999999999998
1921268033293604
598320419039961
9297857723517944
174538706646999
324501774232494
 ⋮
22775246980000022
22775246980000022
22775246980000022
22775246980000022
22775246980000022
22775246980000022
22775246980000022
22775246980000022
22775246980000022

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

@time sol = solve(prob, KenCarp47(), saveat=0.1)

167.519230 seconds (182.59 M allocations: 7.458 GiB, 2.23% gc time, 88.99% 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.05:1.0, 0.0:0.05:1.0)
u: Dict{Symbolics.Num, Array{Float64, 3}} with 2 entries:
  u(x, y, t) => [0.0 0.227752 … 0.227752 0.0; 0.0 0.227752 … 0.227752 0.0; … ; …
  v(x, y, t) => [0.0 0.0 … 0.0 0.0; 0.279514 0.279514 … 0.279514 0.279514; … ; …

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)

5.003356858277019

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#