Symbolic Brusselator PDE

Symbolic Brusselator PDE#

Source

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

[ Info: Precompiling IJuliaExt [2f4121a4-3b3a-5ce6-9c5e-1f2673ce168a]

@parameters 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) =& 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) + 3.4 u\left( x, y, t \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) =& 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) + 3.4 u\left( x, y, t \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

discretization = let N = 16, order = 2
    MOLFiniteDifference([x=>N, y=>N], t, approx_order=order)
end
prob = discretize(pdesys, discretization)

┌ 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 /srv/juliapkg/packages/MethodOfLines/xyn4D/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: 450-element Vector{Float64}:
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 0.3414608815189256
 ⋮
 2.3349038524102017
 2.8284271247461903
 3.174538706646999
 3.352525018549451
 3.3525250185494513
 3.174538706646999
 2.8284271247461907
 2.3349038524102017
 1.7279999999999998
 1.060596058827299
 0.41906562731868124
 0.0

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

sol = solve(prob, TRBDF2(), saveat=0.1)

retcode: Success
Interpolation: Dict{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
  1.0
  1.1
  1.2
  ⋮
 10.4
 10.5
 10.6
 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.06666666666666667:1.0, 0.0:0.06666666666666667:1.0)
u: Dict{Num, Array{Float64, 3}} with 2 entries:
  u(x, y, t) => [0.0 0.341461 … 0.341461 0.0; 0.0 0.341461 … 0.341461 0.0; … ; …
  v(x, y, t) => [0.0 0.0 … 0.0 0.0; 0.419066 0.419066 … 0.419066 0.419066; … ; …

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, solu)

4.1049465699855725

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, 4.2))
end

mp4(anim, fps = 8)

[ Info: Saved animation to /tmp/docs/pde/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, 4.2))
end

mp4(anim, fps = 8)

[ Info: Saved animation to /tmp/docs/pde/tmp.mp4

Runtime environment#

using Pkg
Pkg.status()

Status `/tmp/docs/pde/Project.toml`
⌅ [5b8099bc] DomainSets v0.6.7
  [94925ecb] MethodOfLines v0.11.0
  [961ee093] ModelingToolkit v9.15.0
  [8913a72c] NonlinearSolve v3.12.4
  [1dea7af3] OrdinaryDiffEq v6.80.1
  [91a5bcdd] Plots v1.40.4
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

using InteractiveUtils
InteractiveUtils.versioninfo()

Julia Version 1.10.3
Commit 0b4590a5507 (2024-04-30 10:59 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 4 × AMD EPYC 7763 64-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 4 default, 0 interactive, 2 GC (on 4 virtual cores)
Environment:
  JULIA_CI = true
  JULIA_NUM_THREADS = auto
  JULIA_CONDAPKG_BACKEND = Null
  JULIA_PATH = /usr/local/julia/
  JULIA_DEPOT_PATH = /srv/juliapkg/
Contents