Automated Efficient Solution of Nonlinear Partial Differential Equations

Solving nonlinear partial differential equations (PDEs) is hard. Solving nonlinear PDEs fast and accurately is even harder. Doing it all in an automated method from just a symbolic description is just plain fun. That's what we'd demonstrate here: how to solve a nonlinear PDE from a purely symbolic definition using the combination of ModelingToolkit, MethodOfLines, and DifferentialEquations.jl.

Note

This example is a combination of the Brusselator tutorial from MethodOfLines.jl and the Solving Large Stiff Equations tutorial from DifferentialEquations.jl.

Required Dependencies

The following parts of the SciML Ecosystem will be used in this tutorial:

Module	Description
ModelingToolkit.jl	The symbolic modeling environment
MethodOfLines.jl	The symbolic PDE discretization tooling
DifferentialEquations.jl	The numerical differential equation solvers
LinearSolve.jl	The numerical linear solvers

Problem Setup

The Brusselator PDE is defined as follows:

\[\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^2v + \alpha(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}) \end{align}\]

where

\[f(x, y, t) = \begin{cases} 5 & \quad \text{if } (x-0.3)^2+(y-0.6)^2 ≤ 0.1^2 \text{ and } t ≥ 1.1 \\ 0 & \quad \text{else} \end{cases}\]

and the initial conditions are

\[\begin{align} u(x, y, 0) &= 22\cdot (y(1-y))^{3/2} \\ v(x, y, 0) &= 27\cdot (x(1-x))^{3/2} \end{align}\]

with the periodic boundary condition

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

We wish to obtain the solution to this PDE on a timespan of $t \in [0,11.5]$.

Defining the symbolic PDEsystem with ModelingToolkit.jl

With ModelingToolkit.jl, we first symbolically define the system, see also the docs for PDESystem:

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, LinearSolve, DomainSets

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

brusselator_f(x, y, t) = (((x - 0.3)^2 + (y - 0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.0

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)

eq = [
    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 BCs
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(eq, bcs, domains, [x, y, t], [u(x, y, t), v(x, y, t)])

\[ \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} \]

Looks just like the LaTeX description, right? Now let's solve it.

Automated symbolic discretization with MethodOfLines.jl

Next we create the discretization. Here we will use the finite difference method via method of lines. Method of lines is a method of recognizing that a discretization of a partial differential equation transforms it into a new numerical problem. For example:

Discretization Form	Numerical Problem Type
Finite Difference, Finite Volume, Finite Element, discretizing all variables	`NonlinearProblem`
Finite Difference, Finite Volume, Finite Element, discretizing all variables except time	`ODEProblem`/`DAEProblem`
Physics-Informed Neural Network	`OptimizationProblem`
Feynman-Kac Formula	`SDEProblem`
Universal Stochastic Differential Equation (High dimensional PDEs)	`OptimizationProblem` inverse problem over `SDEProblem`

Thus the process of solving a PDE is fundamentally about transforming its symbolic form to a standard numerical problem and solving the standard numerical problem using one of the solvers in the SciML ecosystem! Here we will demonstrate one of the most classic methods: the finite difference method. Since the Brusselator is a time-dependent PDE with heavy stiffness in the time-domain, we will leave time undiscretized, which means that we will use the finite difference method in the x and y domains to obtain a representation of the equation at `u_i = u(x_i,y_i)grid point values, obtaining an ODEu_i' = \ldots that defines how the values at the grid points evolve over time.

To do this, we use the MOLFiniteDifference construct of MethodOfLines.jl as follows:

N = 32

dx = (x_max - x_min) / N
dy = (y_max - y_min) / N

order = 2

discretization = MOLFiniteDifference([x => dx, y => dy], t, approx_order = order,
    grid_align = center_align)

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

Next, we discretize the system, converting the PDESystem in to an ODEProblem:

prob = discretize(pdesys, discretization);

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 11.5)
u0: 2048-element Vector{Float64}:
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 ⋮
 2.1921268033293604
 1.9075279151581102
 1.605464573318729
 1.2924521756218903
 0.9766542925609525
 0.668640545523243
 0.3829487927768564
 0.1422185904446023
 0.0

Solving the PDE

Now your problem can be solved with an appropriate ODE solver. This is just your standard DifferentialEquations.jl usage, though we'll return to this point in a bit to talk about efficiency:

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

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.03125:1.0, 0.0:0.03125:1.0)
u: Dict{Symbolics.Num, Array{Float64, 3}} with 2 entries:
  u(x, y, t) => [0.0 0.115882 … 0.115882 0.0; 0.0 0.115882 … 0.115882 0.0; … ; …
  v(x, y, t) => [0.0 0.0 … 0.0 0.0; 0.142219 0.142219 … 0.142219 0.142219; … ; …

Examining Results via the Symbolic Solution Interface

Now that we have solved the ODE representation of the PDE, we have an PDETimeSeriesSolution that wraps an ODESolution, which we can get with sol.original_sol. If we look at the original sol, it represents $u_i' = \ldots$ at each of the grid points. If you check sol.original_sol.u inside the solution, that's those values... but that's not very helpful. How do you interpret original_sol[1]? How do you interpret original_sol[1,:]?

To make the handling of such cases a lot simpler, MethodOfLines.jl implements a symbolic interface for the solution object that allows for interpreting the computation through its original representation. For example, if we want to know how to interpret the values of the grid corresponding to the independent variables, we can just index using symbolic variables:

discrete_x = sol[x];
discrete_y = sol[y];
discrete_t = sol[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.5

What this tells us is that, for a solution at a given time point, say original_sol[1] for the solution at the initial time (the initial condition), the value original_sol[1][1] is the solution at the grid point (discrete_x[1], discrete_y[1]). For values that are not the initial time point, original_sol[i] corresponds to the solution at discrete_t[i].

But we also have two dependent variables, u and v. How do we interpret which of the results correspond to the different dependent variables? This is done by indexing the solution by the dependent variables! For example:

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

33×33×116 Array{Float64, 3}:
[:, :, 1] =
 0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
 0.142219  0.142219  0.142219  0.142219     0.142219  0.142219  0.142219
 0.382949  0.382949  0.382949  0.382949     0.382949  0.382949  0.382949
 0.668641  0.668641  0.668641  0.668641     0.668641  0.668641  0.668641
 0.976654  0.976654  0.976654  0.976654     0.976654  0.976654  0.976654
 1.29245   1.29245   1.29245   1.29245   …  1.29245   1.29245   1.29245
 1.60546   1.60546   1.60546   1.60546      1.60546   1.60546   1.60546
 1.90753   1.90753   1.90753   1.90753      1.90753   1.90753   1.90753
 2.19213   2.19213   2.19213   2.19213      2.19213   2.19213   2.19213
 2.45397   2.45397   2.45397   2.45397      2.45397   2.45397   2.45397
 ⋮                                       ⋱  ⋮                   
 2.19213   2.19213   2.19213   2.19213      2.19213   2.19213   2.19213
 1.90753   1.90753   1.90753   1.90753   …  1.90753   1.90753   1.90753
 1.60546   1.60546   1.60546   1.60546      1.60546   1.60546   1.60546
 1.29245   1.29245   1.29245   1.29245      1.29245   1.29245   1.29245
 0.976654  0.976654  0.976654  0.976654     0.976654  0.976654  0.976654
 0.668641  0.668641  0.668641  0.668641     0.668641  0.668641  0.668641
 0.382949  0.382949  0.382949  0.382949  …  0.382949  0.382949  0.382949
 0.142219  0.142219  0.142219  0.142219     0.142219  0.142219  0.142219
 0.0       0.0       0.0       0.0          0.0       0.0       0.0

[:, :, 2] =
 0.0      2.02429  2.02429  2.02429  …  2.02429  2.02429  2.02429  2.02429
 2.02429  2.02429  2.02429  2.02429     2.02429  2.02429  2.02429  2.02429
 2.02428  2.02428  2.02428  2.02428     2.02428  2.02428  2.02428  2.02428
 2.02427  2.02427  2.02427  2.02427     2.02427  2.02427  2.02427  2.02427
 2.02426  2.02426  2.02426  2.02426     2.02426  2.02426  2.02426  2.02426
 2.02424  2.02424  2.02424  2.02424  …  2.02424  2.02424  2.02424  2.02424
 2.02423  2.02423  2.02423  2.02423     2.02423  2.02423  2.02423  2.02423
 2.0242   2.02421  2.02421  2.02421     2.02421  2.02421  2.02421  2.0242
 2.02418  2.02418  2.02418  2.02418     2.02418  2.02418  2.02418  2.02418
 2.02416  2.02416  2.02416  2.02416     2.02416  2.02416  2.02416  2.02416
 ⋮                                   ⋱           ⋮                 
 2.02418  2.02418  2.02418  2.02418     2.02418  2.02418  2.02418  2.02418
 2.0242   2.02421  2.02421  2.02421  …  2.02421  2.02421  2.02421  2.0242
 2.02423  2.02423  2.02423  2.02423     2.02423  2.02423  2.02423  2.02423
 2.02424  2.02424  2.02424  2.02424     2.02424  2.02424  2.02424  2.02424
 2.02426  2.02426  2.02426  2.02426     2.02426  2.02426  2.02426  2.02426
 2.02427  2.02427  2.02427  2.02427     2.02427  2.02427  2.02427  2.02427
 2.02428  2.02428  2.02428  2.02428  …  2.02428  2.02428  2.02428  2.02428
 2.02429  2.02429  2.02429  2.02429     2.02429  2.02429  2.02429  2.02429
 2.02429  2.02429  2.02429  2.02429     2.02429  2.02429  2.02429  2.02429

[:, :, 3] =
 0.0      2.0791   2.0791   2.0791   …  2.0791   2.0791   2.0791   2.0791
 2.0791   2.0791   2.0791   2.0791      2.0791   2.0791   2.0791   2.0791
 2.07911  2.07911  2.07911  2.07911     2.07911  2.07911  2.07911  2.07911
 2.07912  2.07912  2.07912  2.07912     2.07912  2.07912  2.07912  2.07912
 2.07914  2.07914  2.07914  2.07914     2.07914  2.07914  2.07914  2.07914
 2.07916  2.07916  2.07916  2.07916  …  2.07916  2.07916  2.07916  2.07916
 2.07919  2.07919  2.07919  2.07919     2.07919  2.07919  2.07919  2.07919
 2.07922  2.07922  2.07922  2.07921     2.07921  2.07922  2.07922  2.07922
 2.07924  2.07924  2.07924  2.07924     2.07924  2.07924  2.07924  2.07924
 2.07927  2.07927  2.07927  2.07927     2.07927  2.07927  2.07927  2.07927
 ⋮                                   ⋱           ⋮                 
 2.07924  2.07924  2.07924  2.07924     2.07924  2.07924  2.07924  2.07924
 2.07922  2.07922  2.07922  2.07921  …  2.07921  2.07922  2.07922  2.07922
 2.07919  2.07919  2.07919  2.07919     2.07919  2.07919  2.07919  2.07919
 2.07916  2.07916  2.07916  2.07916     2.07916  2.07916  2.07916  2.07916
 2.07914  2.07914  2.07914  2.07914     2.07914  2.07914  2.07914  2.07914
 2.07912  2.07912  2.07912  2.07912     2.07912  2.07912  2.07912  2.07912
 2.07911  2.07911  2.07911  2.07911  …  2.07911  2.07911  2.07911  2.07911
 2.0791   2.0791   2.0791   2.0791      2.0791   2.0791   2.0791   2.0791
 2.0791   2.0791   2.0791   2.0791      2.0791   2.0791   2.0791   2.0791

;;; … 

[:, :, 114] =
 0.0      3.68145  3.68145  3.68145  …  3.68144  3.68145  3.68145  3.68145
 3.68145  3.68145  3.68145  3.68145     3.68144  3.68144  3.68144  3.68145
 3.68144  3.68145  3.68145  3.68145     3.68144  3.68144  3.68144  3.68144
 3.68144  3.68144  3.68145  3.68145     3.68143  3.68144  3.68144  3.68144
 3.68144  3.68144  3.68144  3.68144     3.68143  3.68143  3.68144  3.68144
 3.68144  3.68144  3.68144  3.68144  …  3.68143  3.68143  3.68143  3.68144
 3.68143  3.68144  3.68144  3.68144     3.68142  3.68143  3.68143  3.68143
 3.68143  3.68144  3.68144  3.68144     3.68142  3.68143  3.68143  3.68143
 3.68143  3.68143  3.68144  3.68144     3.68142  3.68142  3.68143  3.68143
 3.68143  3.68143  3.68144  3.68144     3.68142  3.68142  3.68143  3.68143
 ⋮                                   ⋱           ⋮                 
 3.68146  3.68146  3.68146  3.68146     3.68145  3.68146  3.68146  3.68146
 3.68146  3.68146  3.68146  3.68146  …  3.68145  3.68146  3.68146  3.68146
 3.68146  3.68146  3.68146  3.68146     3.68145  3.68146  3.68146  3.68146
 3.68146  3.68146  3.68146  3.68146     3.68145  3.68146  3.68146  3.68146
 3.68146  3.68146  3.68146  3.68146     3.68145  3.68145  3.68146  3.68146
 3.68146  3.68146  3.68146  3.68146     3.68145  3.68145  3.68145  3.68146
 3.68145  3.68146  3.68146  3.68146  …  3.68145  3.68145  3.68145  3.68145
 3.68145  3.68145  3.68146  3.68146     3.68145  3.68145  3.68145  3.68145
 3.68145  3.68145  3.68145  3.68145     3.68144  3.68145  3.68145  3.68145

[:, :, 115] =
 0.0      2.61478  2.61478  2.61478  …  2.61477  2.61478  2.61478  2.61478
 2.61478  2.61478  2.61478  2.61478     2.61477  2.61477  2.61478  2.61478
 2.61478  2.61478  2.61478  2.61478     2.61477  2.61477  2.61477  2.61478
 2.61477  2.61477  2.61478  2.61478     2.61476  2.61477  2.61477  2.61477
 2.61477  2.61477  2.61477  2.61477     2.61476  2.61476  2.61477  2.61477
 2.61477  2.61477  2.61477  2.61477  …  2.61476  2.61476  2.61476  2.61477
 2.61477  2.61477  2.61477  2.61477     2.61475  2.61476  2.61476  2.61477
 2.61476  2.61477  2.61477  2.61477     2.61475  2.61476  2.61476  2.61476
 2.61476  2.61476  2.61477  2.61477     2.61475  2.61475  2.61476  2.61476
 2.61476  2.61476  2.61477  2.61477     2.61475  2.61475  2.61476  2.61476
 ⋮                                   ⋱           ⋮                 
 2.61479  2.61479  2.61479  2.61479     2.61478  2.61479  2.61479  2.61479
 2.61479  2.61479  2.61479  2.61479  …  2.61478  2.61479  2.61479  2.61479
 2.61479  2.61479  2.61479  2.61479     2.61478  2.61479  2.61479  2.61479
 2.61479  2.61479  2.61479  2.61479     2.61478  2.61479  2.61479  2.61479
 2.61479  2.61479  2.61479  2.61479     2.61478  2.61479  2.61479  2.61479
 2.61479  2.61479  2.61479  2.61479     2.61478  2.61478  2.61479  2.61479
 2.61479  2.61479  2.61479  2.61479  …  2.61478  2.61478  2.61478  2.61479
 2.61478  2.61479  2.61479  2.61479     2.61478  2.61478  2.61478  2.61478
 2.61478  2.61478  2.61478  2.61478     2.61477  2.61478  2.61478  2.61478

[:, :, 116] =
 0.0      1.4548   1.4548   1.4548   …  1.45479  1.4548   1.4548   1.4548
 1.4548   1.4548   1.4548   1.4548      1.45479  1.45479  1.4548   1.4548
 1.4548   1.4548   1.4548   1.4548      1.45479  1.45479  1.45479  1.4548
 1.45479  1.4548   1.4548   1.4548      1.45479  1.45479  1.45479  1.45479
 1.45479  1.45479  1.45479  1.45479     1.45478  1.45479  1.45479  1.45479
 1.45479  1.45479  1.45479  1.45479  …  1.45478  1.45479  1.45479  1.45479
 1.45479  1.45479  1.45479  1.45479     1.45478  1.45478  1.45479  1.45479
 1.45479  1.45479  1.45479  1.45479     1.45478  1.45478  1.45479  1.45479
 1.45479  1.45479  1.45479  1.45479     1.45478  1.45478  1.45478  1.45479
 1.45479  1.45479  1.45479  1.45479     1.45478  1.45478  1.45478  1.45479
 ⋮                                   ⋱           ⋮                 
 1.45481  1.45481  1.45481  1.45481     1.4548   1.4548   1.4548   1.45481
 1.45481  1.45481  1.45481  1.45481  …  1.4548   1.4548   1.4548   1.45481
 1.45481  1.45481  1.45481  1.45481     1.4548   1.4548   1.4548   1.45481
 1.45481  1.45481  1.45481  1.45481     1.4548   1.4548   1.4548   1.45481
 1.4548   1.45481  1.45481  1.45481     1.4548   1.4548   1.4548   1.4548
 1.4548   1.4548   1.45481  1.45481     1.4548   1.4548   1.4548   1.4548
 1.4548   1.4548   1.4548   1.4548   …  1.4548   1.4548   1.4548   1.4548
 1.4548   1.4548   1.4548   1.4548      1.4548   1.4548   1.4548   1.4548
 1.4548   1.4548   1.4548   1.4548      1.45479  1.4548   1.4548   1.4548

This then gives an array of results for the u and v separately, each dimension corresponding to the discrete form of the independent variables.

Using this high-level indexing, we can create an animation of the solution of the Brusselator as follows. For u we receive:

using Plots
anim = @animate for k in 1:length(discrete_t)
    heatmap(solu[2:end, 2:end, k], title = "$(discrete_t[k])") # 2:end since end = 1, periodic condition
end
gif(anim, "plots/Brusselator2Dsol_u.gif", fps = 8)

Brusselator2Dsol_u

and for v:

anim = @animate for k in 1:length(discrete_t)
    heatmap(solv[2:end, 2:end, k], title = "$(discrete_t[k])")
end
gif(anim, "plots/Brusselator2Dsol_v.gif", fps = 8)

Brusselator2Dsol_v

Improving the Solution Process

Now, if all we needed was a single solution, then we're done. Budda bing budda boom, we got a solution, we're outta here. But if for example we're solving an inverse problem on a PDE, or we need to bump it up to higher accuracy, then we will need to make sure we solve this puppy more efficiently. So let's dive into how this can be done.

First of all, large PDEs generally are stiff and thus require an implicit solver. However, their stiffness is generally governed by a nonlinear system which as a sparse Jacobian. Handling that implicit system with sparsity is key to solving the system efficiently, so let's do that!

In order to enable such options, we simply need to pass the ModelingToolkit.jl problem construction options to the discretize call. This looks like:

# Analytical Jacobian expression and sparse Jacobian
prob_sparse = discretize(pdesys, discretization; jac = true, sparse = true)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 11.5)
u0: 2048-element Vector{Float64}:
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 0.11588181443634263
 ⋮
 2.1921268033293604
 1.9075279151581102
 1.605464573318729
 1.2924521756218903
 0.9766542925609525
 0.668640545523243
 0.3829487927768564
 0.1422185904446023
 0.0

Now when we solve the problem it will be a lot faster. We can use BenchmarkTools.jl to assess this performance difference:

using BenchmarkTools
@btime sol = solve(prob, TRBDF2(), saveat = 0.1);
@btime sol = solve(prob_sparse, TRBDF2(), saveat = 0.1);

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.03125:1.0, 0.0:0.03125:1.0)
u: Dict{Symbolics.Num, Array{Float64, 3}} with 2 entries:
  u(x, y, t) => [0.0 0.115882 … 0.115882 0.0; 0.0 0.115882 … 0.115882 0.0; … ; …
  v(x, y, t) => [0.0 0.0 … 0.0 0.0; 0.142219 0.142219 … 0.142219 0.142219; … ; …

But we can further improve this as well. Instead of just using the default linear solver, we can change this to a Newton-Krylov method by passing in the GMRES method:

@btime sol = solve(prob_sparse, TRBDF2(linsolve = KrylovJL_GMRES()), saveat = 0.1);

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.03125:1.0, 0.0:0.03125:1.0)
u: Dict{Symbolics.Num, Array{Float64, 3}} with 2 entries:
  u(x, y, t) => [0.0 0.115882 … 0.115882 0.0; 0.0 0.115882 … 0.115882 0.0; … ; …
  v(x, y, t) => [0.0 0.0 … 0.0 0.0; 0.142219 0.142219 … 0.142219 0.142219; … ; …

But to further improve performance, we can use an iLU preconditioner. This looks like as follows:

using IncompleteLU
function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
    if newW === nothing || newW
        Pl = ilu(convert(AbstractMatrix, W), τ = 50.0)
    else
        Pl = Plprev
    end
    Pl, nothing
end

@btime solve(prob_sparse,
    TRBDF2(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true),
    save_everystep = false);

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: 2-element Vector{Float64}:
  0.0
 11.5ivs: 3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 t
 x
 ydomain:([0.0, 11.5], 0.0:0.03125:1.0, 0.0:0.03125:1.0)
u: Dict{Symbolics.Num, Array{Float64, 3}} with 2 entries:
  u(x, y, t) => [0.0 0.115882 … 0.115882 0.0; 0.0 0.115882 … 0.115882 0.0; … ; …
  v(x, y, t) => [0.0 0.0 … 0.0 0.0; 0.142219 0.142219 … 0.142219 0.142219; … ; …

And now we're zooming! For more information on these performance improvements, check out the deeper dive in the DifferentialEquations.jl tutorials.

If you're interested in figuring out what's the fastest current solver for this kind of PDE, check out the Brusselator benchmark in SciMLBenchmarks.jl