|
| 1 | +# [Tutorial: Solution of the heat equation with Dirichlet boundary conditions](@id tutorial-heat-equation-dirichlet) |
| 2 | + |
| 3 | +We continue the previous tutorial on |
| 4 | +[solving the heat equation with Neumann boundary conditions](@ref tutorial-heat-equation-neumann) |
| 5 | +by looking at Dirichlet boundary conditions instead, resulting in a non-conservative |
| 6 | +production-destruction system. |
| 7 | + |
| 8 | + |
| 9 | +## Definition of the (non-conservative) production-destruction system |
| 10 | + |
| 11 | +Consider the heat equation |
| 12 | + |
| 13 | +```math |
| 14 | +\partial_t u(t,x) = \mu \partial_x^2 u(t,x),\quad u(0,x)=u_0(x), |
| 15 | +``` |
| 16 | + |
| 17 | +with ``μ ≥ 0``, ``t≥ 0``, ``x\in[0,1]``, and homogeneous Dirichlet boundary conditions. |
| 18 | +We use again a finite volume discretization, i.e., we split the domain ``[0, 1]`` into |
| 19 | +``N`` uniform cells of width ``\Delta x = 1 / N``. As degrees of freedom, we use |
| 20 | +the mean values of ``u(t)`` in each cell approximated by the point value ``u_i(t)`` |
| 21 | +in the center of cell ``i``. Finally, we use the classical central finite difference |
| 22 | +discretization of the Laplacian with homogeneous Dirichlet boundary conditions, |
| 23 | +resulting in the ODE |
| 24 | + |
| 25 | +```math |
| 26 | +\partial_t u(t) = L u(t), |
| 27 | +\quad |
| 28 | +L = \frac{\mu}{\Delta x^2} \begin{pmatrix} |
| 29 | + -2 & 1 \\ |
| 30 | + 1 & -2 & 1 \\ |
| 31 | + & \ddots & \ddots & \ddots \\ |
| 32 | + && 1 & -2 & 1 \\ |
| 33 | + &&& 1 & -2 |
| 34 | +\end{pmatrix}. |
| 35 | +``` |
| 36 | + |
| 37 | +The system can be written as a non-conservative PDS with production terms |
| 38 | + |
| 39 | +```math |
| 40 | +\begin{aligned} |
| 41 | +&p_{i,i-1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=2,\dots,N, \\ |
| 42 | +&p_{i,i+1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i+1}(t),\quad i=1,\dots,N-1, |
| 43 | +\end{aligned} |
| 44 | +``` |
| 45 | + |
| 46 | +and destruction terms ``d_{i,j} = p_{j,i}`` for ``i \ne j`` as well as the |
| 47 | +non-conservative destruction terms |
| 48 | + |
| 49 | +```math |
| 50 | +\begin{aligned} |
| 51 | +d_{1,1}(t,\mathbf u(t)) &= \frac{\mu}{\Delta x^2} u_{1}(t), \\ |
| 52 | +d_{N,N}(t,\mathbf u(t)) &= \frac{\mu}{\Delta x^2} u_{N}(t). |
| 53 | +\end{aligned} |
| 54 | +``` |
| 55 | + |
| 56 | + |
| 57 | +## Solution of the non-conservative production-destruction system |
| 58 | + |
| 59 | +Now we are ready to define a [`PDSProblem`](@ref) and to solve this |
| 60 | +problem with a method of |
| 61 | +[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) or |
| 62 | +[OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). |
| 63 | +In the following we use ``N = 100`` nodes and the time domain ``t \in [0,1]``. |
| 64 | +Moreover, we choose the initial condition |
| 65 | + |
| 66 | +```math |
| 67 | +u_0(x) = \sin(\pi x)^2. |
| 68 | +``` |
| 69 | + |
| 70 | +```@example HeatEquationDirichlet |
| 71 | +x_boundaries = range(0, 1, length = 101) |
| 72 | +x = x_boundaries[1:end-1] .+ step(x_boundaries) / 2 |
| 73 | +u0 = @. sinpi(x)^2 # initial solution |
| 74 | +tspan = (0.0, 1.0) # time domain |
| 75 | +
|
| 76 | +nothing #hide |
| 77 | +``` |
| 78 | + |
| 79 | +We will choose three different matrix types for the production terms and |
| 80 | +the resulting linear systems: |
| 81 | + |
| 82 | +1. standard dense matrices (default) |
| 83 | +2. sparse matrices (from SparseArrays.jl) |
| 84 | +3. tridiagonal matrices (from LinearAlgebra.jl) |
| 85 | + |
| 86 | + |
| 87 | +### Standard dense matrices |
| 88 | + |
| 89 | +```@example HeatEquationDirichlet |
| 90 | +using PositiveIntegrators # load ConservativePDSProblem |
| 91 | +
|
| 92 | +function heat_eq_P!(P, u, μ, t) |
| 93 | + fill!(P, 0) |
| 94 | + N = length(u) |
| 95 | + Δx = 1 / N |
| 96 | + μ_Δx2 = μ / Δx^2 |
| 97 | +
|
| 98 | + let i = 1 |
| 99 | + # Dirichlet boundary condition |
| 100 | + P[i, i + 1] = u[i + 1] * μ_Δx2 |
| 101 | + end |
| 102 | +
|
| 103 | + for i in 2:(length(u) - 1) |
| 104 | + # interior stencil |
| 105 | + P[i, i - 1] = u[i - 1] * μ_Δx2 |
| 106 | + P[i, i + 1] = u[i + 1] * μ_Δx2 |
| 107 | + end |
| 108 | +
|
| 109 | + let i = length(u) |
| 110 | + # Dirichlet boundary condition |
| 111 | + P[i, i - 1] = u[i - 1] * μ_Δx2 |
| 112 | + end |
| 113 | +
|
| 114 | + return nothing |
| 115 | +end |
| 116 | +
|
| 117 | +function heat_eq_D!(D, u, μ, t) |
| 118 | + fill!(D, 0) |
| 119 | + N = length(u) |
| 120 | + Δx = 1 / N |
| 121 | + μ_Δx2 = μ / Δx^2 |
| 122 | +
|
| 123 | + # Dirichlet boundary condition |
| 124 | + D[begin] = u[begin] * μ_Δx2 |
| 125 | + D[end] = u[end] * μ_Δx2 |
| 126 | +
|
| 127 | + return nothing |
| 128 | +end |
| 129 | +
|
| 130 | +μ = 1.0e-2 |
| 131 | +prob = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ) # create the PDS |
| 132 | +
|
| 133 | +sol = solve(prob, MPRK22(1.0); save_everystep = false) |
| 134 | +
|
| 135 | +nothing #hide |
| 136 | +``` |
| 137 | + |
| 138 | +```@example HeatEquationDirichlet |
| 139 | +using Plots |
| 140 | +
|
| 141 | +plot(x, u0; label = "u0", xguide = "x", yguide = "u") |
| 142 | +plot!(x, last(sol.u); label = "u") |
| 143 | +``` |
| 144 | + |
| 145 | + |
| 146 | +### Sparse matrices |
| 147 | + |
| 148 | +To use different matrix types for the production terms and linear systems, |
| 149 | +you can use the keyword argument `p_prototype` of |
| 150 | +[`ConservativePDSProblem`](@ref) and [`PDSProblem`](@ref). |
| 151 | + |
| 152 | +```@example HeatEquationDirichlet |
| 153 | +using SparseArrays |
| 154 | +p_prototype = spdiagm(-1 => ones(eltype(u0), length(u0) - 1), |
| 155 | + +1 => ones(eltype(u0), length(u0) - 1)) |
| 156 | +prob_sparse = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ; |
| 157 | + p_prototype = p_prototype) |
| 158 | +
|
| 159 | +sol_sparse = solve(prob_sparse, MPRK22(1.0); save_everystep = false) |
| 160 | +
|
| 161 | +nothing #hide |
| 162 | +``` |
| 163 | + |
| 164 | +```@example HeatEquationDirichlet |
| 165 | +plot(x, u0; label = "u0", xguide = "x", yguide = "u") |
| 166 | +plot!(x, last(sol_sparse.u); label = "u") |
| 167 | +``` |
| 168 | + |
| 169 | + |
| 170 | +### Tridiagonal matrices |
| 171 | + |
| 172 | +The sparse matrices used in this case have a very special structure |
| 173 | +since they are in fact tridiagonal matrices. Thus, we can also use |
| 174 | +the special matrix type `Tridiagonal` from the standard library |
| 175 | +`LinearAlgebra`. |
| 176 | + |
| 177 | +```@example HeatEquationDirichlet |
| 178 | +using LinearAlgebra |
| 179 | +p_prototype = Tridiagonal(ones(eltype(u0), length(u0) - 1), |
| 180 | + ones(eltype(u0), length(u0)), |
| 181 | + ones(eltype(u0), length(u0) - 1)) |
| 182 | +prob_tridiagonal = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ; |
| 183 | + p_prototype = p_prototype) |
| 184 | +
|
| 185 | +sol_tridiagonal = solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false) |
| 186 | +
|
| 187 | +nothing #hide |
| 188 | +``` |
| 189 | + |
| 190 | +```@example HeatEquationDirichlet |
| 191 | +plot(x, u0; label = "u0", xguide = "x", yguide = "u") |
| 192 | +plot!(x, last(sol_tridiagonal.u); label = "u") |
| 193 | +``` |
| 194 | + |
| 195 | + |
| 196 | + |
| 197 | +### Performance comparison |
| 198 | + |
| 199 | +Finally, we use [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) |
| 200 | +to compare the performance of the different implementations. |
| 201 | + |
| 202 | +```@example HeatEquationDirichlet |
| 203 | +using BenchmarkTools |
| 204 | +@benchmark solve(prob, MPRK22(1.0); save_everystep = false) |
| 205 | +``` |
| 206 | + |
| 207 | +```@example HeatEquationDirichlet |
| 208 | +@benchmark solve(prob_sparse, MPRK22(1.0); save_everystep = false) |
| 209 | +``` |
| 210 | + |
| 211 | +By default, we use an LU factorization for the linear systems. At the time of |
| 212 | +writing, Julia uses |
| 213 | +[SparseArrays.jl](https://github.com/JuliaSparse/SparseArrays.jl) |
| 214 | +defaulting to UMFPACK from SuiteSparse in this case. However, the linear |
| 215 | +systems do not necessarily have the structure for which UMFPACK is optimized |
| 216 | + for. Thus, it is often possible to gain performance by switching to KLU |
| 217 | + instead. |
| 218 | + |
| 219 | +```@example HeatEquationDirichlet |
| 220 | +using LinearSolve |
| 221 | +@benchmark solve(prob_sparse, MPRK22(1.0; linsolve = KLUFactorization()); save_everystep = false) |
| 222 | +``` |
| 223 | + |
| 224 | +```@example HeatEquationDirichlet |
| 225 | +@benchmark solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false) |
| 226 | +``` |
| 227 | + |
| 228 | + |
| 229 | +## Package versions |
| 230 | + |
| 231 | +These results were obtained using the following versions. |
| 232 | +```@example HeatEquationDirichlet |
| 233 | +using InteractiveUtils |
| 234 | +versioninfo() |
| 235 | +println() |
| 236 | +
|
| 237 | +using Pkg |
| 238 | +Pkg.status(["PositiveIntegrators", "SparseArrays", "KLU", "LinearSolve", "OrdinaryDiffEq"], |
| 239 | + mode=PKGMODE_MANIFEST) |
| 240 | +nothing # hide |
| 241 | +``` |
0 commit comments