|
| 1 | +# Tutorial: Solution of the linear advection equation |
| 2 | + |
| 3 | +This tutorial is about the efficient solution of production-destruction systems (PDS) with a large number of differential equations. |
| 4 | +We will explore several ways to represent such large systems and assess their efficiency. |
| 5 | + |
| 6 | +## Definition of the production-destruction system |
| 7 | + |
| 8 | +One example of the occurrence of a PDS with a large number of equations is the space discretization of a partial differential equation. In this tutorial we want to solve the linear advection equation |
| 9 | + |
| 10 | +```math |
| 11 | +\partial_t u(t,x)=-a\partial_x u(t,x),\quad u(0,x)=u_0(x) |
| 12 | +``` |
| 13 | + |
| 14 | +with ``a>0``, ``t≥ 0``, ``x\in[0,1]`` and periodic boundary conditions. To keep things as simple as possible, we |
| 15 | +discretize the space domain as ``0=x_0<x_1\dots <x_{N-1}<x_N=1`` with ``x_i = i Δ x`` for ``i=0,\dots,N`` and ``Δx=1/N``. An upwind discretization of the spatial derivative yields the ODE system |
| 16 | + |
| 17 | +```math |
| 18 | +\begin{aligned} |
| 19 | +&\partial_t u_1(t) =-\frac{a}{Δx}\bigl(u_1(t)-u_{N}(t)\bigr),\\ |
| 20 | +&\partial_t u_i(t) =-\frac{a}{Δx}\bigl(u_i(t)-u_{i-1}(t)\bigr),\quad i=2,\dots,N, |
| 21 | +\end{aligned} |
| 22 | +``` |
| 23 | + |
| 24 | +where ``u_i(t)`` is an approximation of ``u(t,x_i)`` for ``i=1,\dots, N``. |
| 25 | +This system can also be written as ``\partial_t \mathbf u(t)=\mathbf A\mathbf u(t)`` with ``\mathbf u(t)=(u_1(t),\dots,u_N(t))`` and |
| 26 | + |
| 27 | +```math |
| 28 | +\mathbf A= \frac{a}{Δ x}\begin{bmatrix}-1&0&\dots&0&1\\1&-1&\ddots&&0\\0&\ddots&\ddots&\ddots&\vdots\\ \vdots&\ddots&\ddots&\ddots&0\\0&\dots&0&1&-1\end{bmatrix}. |
| 29 | +``` |
| 30 | + |
| 31 | +In particular the matrix ``\mathbf A`` shows that there is a single production term and a single destruction term per equation. |
| 32 | +Furthermore, the system is conservative as ``\mathbf A`` has column sum zero. |
| 33 | +To be precise, the production matrix ``\mathbf P = (p_{i,j})`` of this conservative PDS is given by |
| 34 | + |
| 35 | +```math |
| 36 | +\begin{aligned} |
| 37 | +&p_{1,N}(t,\mathbf u(t)) = \frac{a}{Δ x}u_N(t),\\ |
| 38 | +&p_{i,i-1}(t,\mathbf u(t)) = \frac{a}{Δ x}u_{i-1}(t),\quad i=2,\dots,N. |
| 39 | +\end{aligned} |
| 40 | +``` |
| 41 | + |
| 42 | +Since the PDS is conservative, we have ``d_{i,j}=p_{j,i}`` and the system is fully determined by the production matrix ``\mathbf P``. |
| 43 | + |
| 44 | +## Solution of the production-destruction system |
| 45 | + |
| 46 | +Now we are ready to define a `ConservativePDSProblem` and to solve this problem with a method of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) or [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). In the following we use ``a=1``, ``N=1000`` and the time domain ``t\in[0,1]``. Moreover, we choose the step function |
| 47 | + |
| 48 | +```math |
| 49 | +u_0(x)=\begin{cases}1, & 0.4 ≤ x ≤ 0.6,\\ 0,& \text{elsewhere}\end{cases} |
| 50 | +``` |
| 51 | + |
| 52 | +as initial condition. Due to the periodic boundary conditions and the transport velocity ``a=1``, the solution at time ``t=1`` is identical to the initial distribution, i.e. ``u(1,x) = u_0(x)``. |
| 53 | + |
| 54 | +```@example LinearAdvection |
| 55 | +N = 1000 # number of subintervals |
| 56 | +dx = 1/N # mesh width |
| 57 | +x = LinRange(dx, 1.0, N) # discretization points x_1,...,x_N = x_0 |
| 58 | +u0 = @. 0.0 + (0.4 ≤ x ≤ 0.6) * 1.0 # initial solution |
| 59 | +tspan = (0.0, 1.0) # time domain |
| 60 | +
|
| 61 | +nothing #hide |
| 62 | +``` |
| 63 | + |
| 64 | +As mentioned above, we will try different approaches to solve this PDS and compare their efficiency. These are |
| 65 | +1. an in-place implementation with a dense matrix, |
| 66 | +2. an in-place implementation with a sparse matrix. |
| 67 | + |
| 68 | +### Standard in-place implementation |
| 69 | + |
| 70 | +```@example LinearAdvection |
| 71 | +using PositiveIntegrators # load ConservativePDSProblem |
| 72 | +
|
| 73 | +function lin_adv_P!(P, u, p, t) |
| 74 | + P .= 0.0 |
| 75 | + N = length(u) |
| 76 | + dx = 1 / N |
| 77 | + P[1, N] = u[N] / dx |
| 78 | + for i in 2:N |
| 79 | + P[i, i - 1] = u[i - 1] / dx |
| 80 | + end |
| 81 | + return nothing |
| 82 | +end |
| 83 | +
|
| 84 | +prob = ConservativePDSProblem(lin_adv_P!, u0, tspan) # create the PDS |
| 85 | +
|
| 86 | +sol = solve(prob, MPRK43I(1.0, 0.5); save_everystep = false) |
| 87 | +
|
| 88 | +nothing #hide |
| 89 | +``` |
| 90 | + |
| 91 | +```@example LinearAdvection |
| 92 | +using Plots |
| 93 | +
|
| 94 | +plot(x, u0; label = "u0", xguide = "x", yguide = "u") |
| 95 | +plot!(x, last(sol.u); label = "u") |
| 96 | +``` |
| 97 | + |
| 98 | +### Using sparse matrices |
| 99 | + |
| 100 | +TODO: Some text |
| 101 | + |
| 102 | +```@example LinearAdvection |
| 103 | +using SparseArrays |
| 104 | +p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1), |
| 105 | + N - 1 => ones(eltype(u0), 1)) |
| 106 | +prob_sparse = ConservativePDSProblem(lin_adv_P!, u0, tspan; p_prototype=p_prototype) |
| 107 | +
|
| 108 | +sol_sparse = solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false) |
| 109 | +
|
| 110 | +nothing #hide |
| 111 | +``` |
| 112 | + |
| 113 | +```@example LinearAdvection |
| 114 | +plot(x,u0; label = "u0", xguide = "x", yguide = "u") |
| 115 | +plot!(x, last(sol_sparse.u); label = "u") |
| 116 | +``` |
| 117 | + |
| 118 | +### Performance comparison |
| 119 | + |
| 120 | +Finally, we use [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) |
| 121 | +to compare the performance of the different implementations. |
| 122 | + |
| 123 | +```@example LinearAdvection |
| 124 | +using BenchmarkTools |
| 125 | +@benchmark solve(prob, MPRK43I(1.0, 0.5); save_everystep = false) |
| 126 | +``` |
| 127 | + |
| 128 | +```@example LinearAdvection |
| 129 | +@benchmark solve(prob_sparse, MPRK43I(1.0, 0.5); save_everystep = false) |
| 130 | +``` |
0 commit comments