From e180d09d23fecdc12e85aa4cba6801ad4f7bfbc1 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 10:06:29 +0200 Subject: [PATCH 01/15] add heat equation tutorial --- docs/make.jl | 1 + docs/src/heat_equation.md | 217 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 218 insertions(+) create mode 100644 docs/src/heat_equation.md diff --git a/docs/make.jl b/docs/make.jl index e4307ce1..9ae54b06 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -75,6 +75,7 @@ makedocs(modules = [PositiveIntegrators], "Home" => "index.md", "Tutorials" => [ "Linear Advection" => "linear_advection.md", + "Heat Equation" => "heat_equation.md", ], "API reference" => "api_reference.md", "Contributing" => "contributing.md", diff --git a/docs/src/heat_equation.md b/docs/src/heat_equation.md new file mode 100644 index 00000000..7fb4b652 --- /dev/null +++ b/docs/src/heat_equation.md @@ -0,0 +1,217 @@ +# [Tutorial: Solution of the heat equation](@id tutorial-heat-equation) + +Similar to the +[tutorial on the linear advection equation](@ref tutorial-linear-advection), +we will demonstrate how to solve a conservative production-destruction +system (PDS) resulting from a PDE discretization and means to improve +the performance. + + +## Definition of the production-destruction system + +Consider the heat equation + +```math +\partial_t u(t,x) = \mu \partial_x^2 u(t,x),\quad u(0,x)=u_0(x), +``` + +with ``μ ≥ 0``, ``t≥ 0``, ``x\in[0,1]``, and homogeneous Neumann boundary conditions. +We choose the classical central finite difference discretization of the +Laplacian, resulting in the ODE + +```math +\partial_t u(t) = L u(t), +\quad +L = \frac{\mu}{\Delta x^2} \begin{pmatrix} + -1 & 1 \\ + 1 & -2 & 1 \\ + & \ddots & \ddots & \ddots \\ + && 1 & -2 & 1 \\ + &&& 1 & -1 +\end{pmatrix}. +``` + +The system can be written as a conservative PDS with production terms + +```math +\begin{aligned} +&p_{i,i-1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=2,\dots,N, +&p_{i,i+1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=1,\dots,N-1, +\end{aligned} +``` + +and destruction terms ``d_{i,j} = p_{j,i}``. + + +## Solution of the production-destruction system + +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 ``N=100`` nodes and the time domain ``t\in[0,1]``. +Moreover, we choose the initial condition + +```math +u_0(x) = \cos(\pi x). +``` + +```@example HeatEquation +x_boundaries = range(0, 1, length = 101) +x = x_boundaries[1:end-1] .+ step(x_boundaries) / 2 +u0 = @. cospi(x) # initial solution +tspan = (0.0, 1.0) # time domain + +nothing #hide +``` + +We will choose three different matrix types for the production terms and +the resulting linear systems: + +1. standard dense matrices (default) +2. sparse matrices (from SparseArrays.jl) +3. tridiagonal matrices (from LinearAlgebra.jl) + + +### Standard dense matrices + +```@example HeatEquation +using PositiveIntegrators # load ConservativePDSProblem + +function heat_eq_P!(P, u, μ, t) + fill!(P, 0) + N = length(u) + Δx = 1 / N + μ_Δx2 = μ / Δx^2 + + let i = 1 + # Neumann boundary condition + P[i, i + 1] = u[i + 1] * μ_Δx2 + end + + for i in 2:(length(u) - 1) + # interior stencil + P[i, i - 1] = u[i - 1] * μ_Δx2 + P[i, i + 1] = u[i + 1] * μ_Δx2 + end + + let i = length(u) + # Neumann boundary condition + P[i, i - 1] = u[i - 1] * μ_Δx2 + end + + return nothing +end + +μ = 1.0e-2 +prob = ConservativePDSProblem(heat_eq_P!, u0, tspan, μ) # create the PDS + +sol = solve(prob, MPRK22(1.0); save_everystep = false) + +nothing #hide +``` + +```@example HeatEquation +using Plots + +plot(x, u0; label = "u0", xguide = "x", yguide = "u") +plot!(x, last(sol.u); label = "u") +``` + + +### Sparse matrices + +To use different matrix types for the production terms and linear systems, +you can use the keyword argument `p_prototype` of +[`ConservativePDSProblem`](@ref) and [`PDSProblem`](@ref). + +```@example HeatEquation +using SparseArrays +p_prototype = spdiagm(-1 => ones(eltype(u0), length(u0) - 1), + +1 => ones(eltype(u0), length(u0) - 1)) +prob_sparse = ConservativePDSProblem(heat_eq_P!, u0, tspan, μ; + p_prototype = p_prototype) + +sol_sparse = solve(prob_sparse, MPRK22(1.0); save_everystep = false) + +nothing #hide +``` + +```@example HeatEquation +plot(x,u0; label = "u0", xguide = "x", yguide = "u") +plot!(x, last(sol_sparse.u); label = "u") +``` + + +### Tridiagonal matrices + +The sparse matrices used in this case have a very special structure +since they are in fact tridiagonal matrices. Thus, we can also use +the special matrix type `Tridiagonal` from the standard library +`LinearAlgebra`. + +```@example HeatEquation +using LinearAlgebra +p_prototype = Tridiagonal(ones(eltype(u0), length(u0) - 1), + ones(eltype(u0), length(u0)), + ones(eltype(u0), length(u0) - 1)) +prob_tridiagonal = ConservativePDSProblem(heat_eq_P!, u0, tspan, μ; + p_prototype = p_prototype) + +sol_tridiagonal = solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false) + +nothing #hide +``` + +```@example HeatEquation +plot(x,u0; label = "u0", xguide = "x", yguide = "u") +plot!(x, last(sol_tridiagonal.u); label = "u") +``` + + + +### Performance comparison + +Finally, we use [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) +to compare the performance of the different implementations. + +```@example HeatEquation +using BenchmarkTools +@benchmark solve(prob, MPRK22(1.0); save_everystep = false) +``` + +```@example HeatEquation +@benchmark solve(prob_sparse, MPRK22(1.0); save_everystep = false) +``` + +By default, we use an LU factorization for the linear systems. At the time of +writing, Julia uses +[SparseArrays.jl](https://github.com/JuliaSparse/SparseArrays.jl) +defaulting to UMFPACK from SuiteSparse in this case. However, the linear +systems do not necessarily have the structure for which UMFPACK is optimized + for. Thus, it is often possible to gain performance by switching to KLU + instead. + +```@example LinearAdvection +using LinearSolve +@benchmark solve(prob_sparse, MPRK22(1.0; linsolve = KLUFactorization()); save_everystep = false) +``` + +```@example HeatEquation +@benchmark solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false) +``` + + +## Package versions + +These results were obtained using the following versions. +```@example LinearAdvection +using InteractiveUtils +versioninfo() +println() + +using Pkg +Pkg.status(["PositiveIntegrators", "SparseArrays", "KLU", "LinearSolve", "OrdinaryDiffEq"], + mode=PKGMODE_MANIFEST) +nothing # hide +``` From 37cdb0d472cbd100796e7f3a920e956def75c3f0 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 10:54:25 +0200 Subject: [PATCH 02/15] add LinearAlgebra to Project.toml --- docs/Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index 286d97ce..5ad7099b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,7 @@ [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -8,6 +9,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] BenchmarkTools = "1" Documenter = "1" +LinearAlgebra = "1.7" OrdinaryDiffEq = "6.59" Plots = "1" SparseArrays = "1.7" From 12ee81b6b99c5f99bd57b98a3219955d664b6018 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 10:55:51 +0200 Subject: [PATCH 03/15] update ref --- docs/src/heat_equation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/heat_equation.md b/docs/src/heat_equation.md index 7fb4b652..398793cc 100644 --- a/docs/src/heat_equation.md +++ b/docs/src/heat_equation.md @@ -1,7 +1,7 @@ # [Tutorial: Solution of the heat equation](@id tutorial-heat-equation) Similar to the -[tutorial on the linear advection equation](@ref tutorial-linear-advection), +[tutorial on linear advection](@ref tutorial-linear-advection), we will demonstrate how to solve a conservative production-destruction system (PDS) resulting from a PDE discretization and means to improve the performance. From ce42f3386fdd5cf4f96c7b58062b27caa6fa2c5e Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 12:09:27 +0200 Subject: [PATCH 04/15] add missing ID --- docs/src/linear_advection.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md index 6694d441..95c1d5a7 100644 --- a/docs/src/linear_advection.md +++ b/docs/src/linear_advection.md @@ -1,4 +1,4 @@ -# Tutorial: Solution of the linear advection equation +# [Tutorial: Solution of the linear advection equation](@id tutorial-linear-advection) This tutorial is about the efficient solution of production-destruction systems (PDS) with a large number of differential equations. We will explore several ways to represent such large systems and assess their efficiency. From c101f0da3d554307b5eb6f68760b7881cc1ae014 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 12:22:23 +0200 Subject: [PATCH 05/15] fix typo --- docs/src/heat_equation.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/heat_equation.md b/docs/src/heat_equation.md index 398793cc..a7ec4596 100644 --- a/docs/src/heat_equation.md +++ b/docs/src/heat_equation.md @@ -192,7 +192,7 @@ systems do not necessarily have the structure for which UMFPACK is optimized for. Thus, it is often possible to gain performance by switching to KLU instead. -```@example LinearAdvection +```@example HeatEquation using LinearSolve @benchmark solve(prob_sparse, MPRK22(1.0; linsolve = KLUFactorization()); save_everystep = false) ``` @@ -205,7 +205,7 @@ using LinearSolve ## Package versions These results were obtained using the following versions. -```@example LinearAdvection +```@example HeatEquation using InteractiveUtils versioninfo() println() From ec83dc23385992aeb445dfed8d6e8fda8d528a0e Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 12:36:16 +0200 Subject: [PATCH 06/15] add missing packages to docs/Project.toml --- docs/Project.toml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index 5ad7099b..cc100db0 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,15 +1,21 @@ [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] BenchmarkTools = "1" Documenter = "1" +InteractiveUtils = "1" LinearAlgebra = "1.7" +LinearSolve = "2.21" OrdinaryDiffEq = "6.59" +Pkg = "1" Plots = "1" SparseArrays = "1.7" From 446074d92ef243c344393dce876ddb0b62f10f57 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 13:28:12 +0200 Subject: [PATCH 07/15] Apply suggestions from code review --- docs/src/heat_equation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/heat_equation.md b/docs/src/heat_equation.md index a7ec4596..7283fa9b 100644 --- a/docs/src/heat_equation.md +++ b/docs/src/heat_equation.md @@ -35,7 +35,7 @@ The system can be written as a conservative PDS with production terms ```math \begin{aligned} -&p_{i,i-1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=2,\dots,N, +&p_{i,i-1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=2,\dots,N, \\ &p_{i,i+1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=1,\dots,N-1, \end{aligned} ``` From c01afd8712e4fffb322f590665de60680e675f88 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 15:43:48 +0200 Subject: [PATCH 08/15] improve tutorial --- docs/src/heat_equation.md | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/docs/src/heat_equation.md b/docs/src/heat_equation.md index 7283fa9b..b1585e09 100644 --- a/docs/src/heat_equation.md +++ b/docs/src/heat_equation.md @@ -7,7 +7,7 @@ system (PDS) resulting from a PDE discretization and means to improve the performance. -## Definition of the production-destruction system +## Definition of the conservative production-destruction system Consider the heat equation @@ -16,8 +16,12 @@ Consider the heat equation ``` with ``μ ≥ 0``, ``t≥ 0``, ``x\in[0,1]``, and homogeneous Neumann boundary conditions. -We choose the classical central finite difference discretization of the -Laplacian, resulting in the ODE +We use a finite volume discretization, i.e., we split the domain ``[0, 1]`` into +``N`` uniform cells of width ``\Delta x = 1 / N``. As degrees of freedom, we use +the mean values of ``u(t)`` in each cell approximated by the point value ``u_i(t)`` +in the center of cell ``i``. Finally, we use the classical central finite difference +discretization of the Laplacian with homogeneous Neumann boundary conditions, +resulting in the ODE ```math \partial_t u(t) = L u(t), @@ -36,30 +40,30 @@ The system can be written as a conservative PDS with production terms ```math \begin{aligned} &p_{i,i-1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=2,\dots,N, \\ -&p_{i,i+1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=1,\dots,N-1, +&p_{i,i+1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i+1}(t),\quad i=1,\dots,N-1, \end{aligned} ``` and destruction terms ``d_{i,j} = p_{j,i}``. -## Solution of the production-destruction system +## Solution of the conservative production-destruction system 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 ``N=100`` nodes and the time domain ``t\in[0,1]``. +In the following we use ``N = 100`` nodes and the time domain ``t \in [0,1]``. Moreover, we choose the initial condition ```math -u_0(x) = \cos(\pi x). +u_0(x) = \cos(\pi x)^2. ``` ```@example HeatEquation x_boundaries = range(0, 1, length = 101) x = x_boundaries[1:end-1] .+ step(x_boundaries) / 2 -u0 = @. cospi(x) # initial solution +u0 = @. cospi(x)^2 # initial solution tspan = (0.0, 1.0) # time domain nothing #hide From dc4d6add9702dffa47e3a9ff7251cebbaf511c86 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 16:01:28 +0200 Subject: [PATCH 09/15] another tutorial with Dirichlet BCs --- docs/make.jl | 3 +- docs/src/heat_equation_dirichlet.md | 239 ++++++++++++++++++ ...t_equation.md => heat_equation_neumann.md} | 28 +- 3 files changed, 255 insertions(+), 15 deletions(-) create mode 100644 docs/src/heat_equation_dirichlet.md rename docs/src/{heat_equation.md => heat_equation_neumann.md} (90%) diff --git a/docs/make.jl b/docs/make.jl index 94659170..bc400550 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -75,7 +75,8 @@ makedocs(modules = [PositiveIntegrators], "Home" => "index.md", "Tutorials" => [ "Linear Advection" => "linear_advection.md", - "Heat Equation" => "heat_equation.md", + "Heat Equation, Neumann BCs" => "heat_equation_neumann.md", + "Heat Equation, Dirichlet BCs" => "heat_equation_dirichlet.md", ], "Troubleshooting, FAQ" => "faq.md", "API reference" => "api_reference.md", diff --git a/docs/src/heat_equation_dirichlet.md b/docs/src/heat_equation_dirichlet.md new file mode 100644 index 00000000..d7bbd393 --- /dev/null +++ b/docs/src/heat_equation_dirichlet.md @@ -0,0 +1,239 @@ +# [Tutorial: Solution of the heat equation with Dirichlet boundary conditions](@id tutorial-heat-equation-dirichlet) + +We continue the previous tutorial on +[solving the heat equation with Neumann boundary conditions](@ref tutorial-heat-equation-neumann) +by looking at Dirichlet boundary conditions instead, resulting in a non-conservative +production-destruction system. + + +## Definition of the (non-conservative) production-destruction system + +Consider the heat equation + +```math +\partial_t u(t,x) = \mu \partial_x^2 u(t,x),\quad u(0,x)=u_0(x), +``` + +with ``μ ≥ 0``, ``t≥ 0``, ``x\in[0,1]``, and homogeneous Dirichlet boundary conditions. +We use again a finite volume discretization, i.e., we split the domain ``[0, 1]`` into +``N`` uniform cells of width ``\Delta x = 1 / N``. As degrees of freedom, we use +the mean values of ``u(t)`` in each cell approximated by the point value ``u_i(t)`` +in the center of cell ``i``. Finally, we use the classical central finite difference +discretization of the Laplacian with homogeneous Dirichlet boundary conditions, +resulting in the ODE + +```math +\partial_t u(t) = L u(t), +\quad +L = \frac{\mu}{\Delta x^2} \begin{pmatrix} + -2 & 1 \\ + 1 & -2 & 1 \\ + & \ddots & \ddots & \ddots \\ + && 1 & -2 & 1 \\ + &&& 1 & -2 +\end{pmatrix}. +``` + +The system can be written as a non-conservative PDS with production terms + +```math +\begin{aligned} +&p_{i,i-1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i-1}(t),\quad i=2,\dots,N, \\ +&p_{i,i+1}(t,\mathbf u(t)) = \frac{\mu}{\Delta x^2} u_{i+1}(t),\quad i=1,\dots,N-1, +\end{aligned} +``` + +and destruction terms ``d_{i,j} = p_{j,i}`` for ``i \ne j`` as well as the +non-conservative destruction terms + +```math +d_{1,1}(t,\mathbf u(t)) &= \frac{\mu}{\Delta x^2} u_{1}(t), \\ +d_{N,N}(t,\mathbf u(t)) &= \frac{\mu}{\Delta x^2} u_{N}(t). +``` + + +## Solution of the non-conservative production-destruction system + +Now we are ready to define a [`PDSProblem`](@ref) 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 ``N = 100`` nodes and the time domain ``t \in [0,1]``. +Moreover, we choose the initial condition + +```math +u_0(x) = \sin(\pi x)^2. +``` + +```@example HeatEquationDirichlet +x_boundaries = range(0, 1, length = 101) +x = x_boundaries[1:end-1] .+ step(x_boundaries) / 2 +u0 = @. sinpi(x)^2 # initial solution +tspan = (0.0, 1.0) # time domain + +nothing #hide +``` + +We will choose three different matrix types for the production terms and +the resulting linear systems: + +1. standard dense matrices (default) +2. sparse matrices (from SparseArrays.jl) +3. tridiagonal matrices (from LinearAlgebra.jl) + + +### Standard dense matrices + +```@example HeatEquationDirichlet +using PositiveIntegrators # load ConservativePDSProblem + +function heat_eq_P!(P, u, μ, t) + fill!(P, 0) + N = length(u) + Δx = 1 / N + μ_Δx2 = μ / Δx^2 + + let i = 1 + # Dirichlet boundary condition + P[i, i + 1] = u[i + 1] * μ_Δx2 + end + + for i in 2:(length(u) - 1) + # interior stencil + P[i, i - 1] = u[i - 1] * μ_Δx2 + P[i, i + 1] = u[i + 1] * μ_Δx2 + end + + let i = length(u) + # Dirichlet boundary condition + P[i, i - 1] = u[i - 1] * μ_Δx2 + end + + return nothing +end + +function heat_eq_D!(D, u, μ, t) + fill!(D, 0) + N = length(u) + Δx = 1 / N + μ_Δx2 = μ / Δx^2 + + # Dirichlet boundary condition + D[begin] = u[begin] * μ_Δx2 + D[end] = u[end] * μ_Δx2 + + return nothing +end + +μ = 1.0e-2 +prob = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ) # create the PDS + +sol = solve(prob, MPRK22(1.0); save_everystep = false) + +nothing #hide +``` + +```@example HeatEquationDirichlet +using Plots + +plot(x, u0; label = "u0", xguide = "x", yguide = "u") +plot!(x, last(sol.u); label = "u") +``` + + +### Sparse matrices + +To use different matrix types for the production terms and linear systems, +you can use the keyword argument `p_prototype` of +[`ConservativePDSProblem`](@ref) and [`PDSProblem`](@ref). + +```@example HeatEquationDirichlet +using SparseArrays +p_prototype = spdiagm(-1 => ones(eltype(u0), length(u0) - 1), + +1 => ones(eltype(u0), length(u0) - 1)) +prob_sparse = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ; + p_prototype = p_prototype) + +sol_sparse = solve(prob_sparse, MPRK22(1.0); save_everystep = false) + +nothing #hide +``` + +```@example HeatEquationDirichlet +plot(x,u0; label = "u0", xguide = "x", yguide = "u") +plot!(x, last(sol_sparse.u); label = "u") +``` + + +### Tridiagonal matrices + +The sparse matrices used in this case have a very special structure +since they are in fact tridiagonal matrices. Thus, we can also use +the special matrix type `Tridiagonal` from the standard library +`LinearAlgebra`. + +```@example HeatEquationDirichlet +using LinearAlgebra +p_prototype = Tridiagonal(ones(eltype(u0), length(u0) - 1), + ones(eltype(u0), length(u0)), + ones(eltype(u0), length(u0) - 1)) +prob_tridiagonal = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ; + p_prototype = p_prototype) + +sol_tridiagonal = solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false) + +nothing #hide +``` + +```@example HeatEquationDirichlet +plot(x,u0; label = "u0", xguide = "x", yguide = "u") +plot!(x, last(sol_tridiagonal.u); label = "u") +``` + + + +### Performance comparison + +Finally, we use [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) +to compare the performance of the different implementations. + +```@example HeatEquationDirichlet +using BenchmarkTools +@benchmark solve(prob, MPRK22(1.0); save_everystep = false) +``` + +```@example HeatEquationDirichlet +@benchmark solve(prob_sparse, MPRK22(1.0); save_everystep = false) +``` + +By default, we use an LU factorization for the linear systems. At the time of +writing, Julia uses +[SparseArrays.jl](https://github.com/JuliaSparse/SparseArrays.jl) +defaulting to UMFPACK from SuiteSparse in this case. However, the linear +systems do not necessarily have the structure for which UMFPACK is optimized + for. Thus, it is often possible to gain performance by switching to KLU + instead. + +```@example HeatEquationDirichlet +using LinearSolve +@benchmark solve(prob_sparse, MPRK22(1.0; linsolve = KLUFactorization()); save_everystep = false) +``` + +```@example HeatEquationDirichlet +@benchmark solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false) +``` + + +## Package versions + +These results were obtained using the following versions. +```@example HeatEquationDirichlet +using InteractiveUtils +versioninfo() +println() + +using Pkg +Pkg.status(["PositiveIntegrators", "SparseArrays", "KLU", "LinearSolve", "OrdinaryDiffEq"], + mode=PKGMODE_MANIFEST) +nothing # hide +``` diff --git a/docs/src/heat_equation.md b/docs/src/heat_equation_neumann.md similarity index 90% rename from docs/src/heat_equation.md rename to docs/src/heat_equation_neumann.md index b1585e09..0497a90a 100644 --- a/docs/src/heat_equation.md +++ b/docs/src/heat_equation_neumann.md @@ -1,4 +1,4 @@ -# [Tutorial: Solution of the heat equation](@id tutorial-heat-equation) +# [Tutorial: Solution of the heat equation with Neumann boundary conditions](@id tutorial-heat-equation-neumann) Similar to the [tutorial on linear advection](@ref tutorial-linear-advection), @@ -49,7 +49,7 @@ and destruction terms ``d_{i,j} = p_{j,i}``. ## Solution of the conservative production-destruction system -Now we are ready to define a `ConservativePDSProblem` and to solve this +Now we are ready to define a [`ConservativePDSProblem`](@ref) 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/). @@ -60,7 +60,7 @@ Moreover, we choose the initial condition u_0(x) = \cos(\pi x)^2. ``` -```@example HeatEquation +```@example HeatEquationNeumann x_boundaries = range(0, 1, length = 101) x = x_boundaries[1:end-1] .+ step(x_boundaries) / 2 u0 = @. cospi(x)^2 # initial solution @@ -79,7 +79,7 @@ the resulting linear systems: ### Standard dense matrices -```@example HeatEquation +```@example HeatEquationNeumann using PositiveIntegrators # load ConservativePDSProblem function heat_eq_P!(P, u, μ, t) @@ -115,7 +115,7 @@ sol = solve(prob, MPRK22(1.0); save_everystep = false) nothing #hide ``` -```@example HeatEquation +```@example HeatEquationNeumann using Plots plot(x, u0; label = "u0", xguide = "x", yguide = "u") @@ -129,7 +129,7 @@ To use different matrix types for the production terms and linear systems, you can use the keyword argument `p_prototype` of [`ConservativePDSProblem`](@ref) and [`PDSProblem`](@ref). -```@example HeatEquation +```@example HeatEquationNeumann using SparseArrays p_prototype = spdiagm(-1 => ones(eltype(u0), length(u0) - 1), +1 => ones(eltype(u0), length(u0) - 1)) @@ -141,7 +141,7 @@ sol_sparse = solve(prob_sparse, MPRK22(1.0); save_everystep = false) nothing #hide ``` -```@example HeatEquation +```@example HeatEquationNeumann plot(x,u0; label = "u0", xguide = "x", yguide = "u") plot!(x, last(sol_sparse.u); label = "u") ``` @@ -154,7 +154,7 @@ since they are in fact tridiagonal matrices. Thus, we can also use the special matrix type `Tridiagonal` from the standard library `LinearAlgebra`. -```@example HeatEquation +```@example HeatEquationNeumann using LinearAlgebra p_prototype = Tridiagonal(ones(eltype(u0), length(u0) - 1), ones(eltype(u0), length(u0)), @@ -167,7 +167,7 @@ sol_tridiagonal = solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false) nothing #hide ``` -```@example HeatEquation +```@example HeatEquationNeumann plot(x,u0; label = "u0", xguide = "x", yguide = "u") plot!(x, last(sol_tridiagonal.u); label = "u") ``` @@ -179,12 +179,12 @@ plot!(x, last(sol_tridiagonal.u); label = "u") Finally, we use [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) to compare the performance of the different implementations. -```@example HeatEquation +```@example HeatEquationNeumann using BenchmarkTools @benchmark solve(prob, MPRK22(1.0); save_everystep = false) ``` -```@example HeatEquation +```@example HeatEquationNeumann @benchmark solve(prob_sparse, MPRK22(1.0); save_everystep = false) ``` @@ -196,12 +196,12 @@ systems do not necessarily have the structure for which UMFPACK is optimized for. Thus, it is often possible to gain performance by switching to KLU instead. -```@example HeatEquation +```@example HeatEquationNeumann using LinearSolve @benchmark solve(prob_sparse, MPRK22(1.0; linsolve = KLUFactorization()); save_everystep = false) ``` -```@example HeatEquation +```@example HeatEquationNeumann @benchmark solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false) ``` @@ -209,7 +209,7 @@ using LinearSolve ## Package versions These results were obtained using the following versions. -```@example HeatEquation +```@example HeatEquationNeumann using InteractiveUtils versioninfo() println() From a5c22b9ec970b78ccfcc6587c731199474ad68da Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 16:01:52 +0200 Subject: [PATCH 10/15] bump version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5bb53446..fb03dd67 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PositiveIntegrators" uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5" authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"] -version = "0.1.15" +version = "0.1.16" [deps] FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" From de14473255d26bf48c42f9c5465b0564649248e0 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 19:27:01 +0200 Subject: [PATCH 11/15] Update docs/src/heat_equation_neumann.md Co-authored-by: Stefan Kopecz --- docs/src/heat_equation_neumann.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/heat_equation_neumann.md b/docs/src/heat_equation_neumann.md index 0497a90a..76b3f71e 100644 --- a/docs/src/heat_equation_neumann.md +++ b/docs/src/heat_equation_neumann.md @@ -142,7 +142,7 @@ nothing #hide ``` ```@example HeatEquationNeumann -plot(x,u0; label = "u0", xguide = "x", yguide = "u") +plot(x, u0; label = "u0", xguide = "x", yguide = "u") plot!(x, last(sol_sparse.u); label = "u") ``` From c58923738678dc2ac2a233ff71bb0daf7260dd93 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 19:27:29 +0200 Subject: [PATCH 12/15] Update docs/src/heat_equation_neumann.md Co-authored-by: Stefan Kopecz --- docs/src/heat_equation_neumann.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/heat_equation_neumann.md b/docs/src/heat_equation_neumann.md index 76b3f71e..3ae00ce7 100644 --- a/docs/src/heat_equation_neumann.md +++ b/docs/src/heat_equation_neumann.md @@ -168,7 +168,7 @@ nothing #hide ``` ```@example HeatEquationNeumann -plot(x,u0; label = "u0", xguide = "x", yguide = "u") +plot(x, u0; label = "u0", xguide = "x", yguide = "u") plot!(x, last(sol_tridiagonal.u); label = "u") ``` From 0b08d550e915f5f53755d4ea6608d09b75f05ede Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 19:27:54 +0200 Subject: [PATCH 13/15] Update docs/src/heat_equation_dirichlet.md Co-authored-by: Stefan Kopecz --- docs/src/heat_equation_dirichlet.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/src/heat_equation_dirichlet.md b/docs/src/heat_equation_dirichlet.md index d7bbd393..6148af40 100644 --- a/docs/src/heat_equation_dirichlet.md +++ b/docs/src/heat_equation_dirichlet.md @@ -47,8 +47,10 @@ and destruction terms ``d_{i,j} = p_{j,i}`` for ``i \ne j`` as well as the non-conservative destruction terms ```math +\begin{aligned} d_{1,1}(t,\mathbf u(t)) &= \frac{\mu}{\Delta x^2} u_{1}(t), \\ d_{N,N}(t,\mathbf u(t)) &= \frac{\mu}{\Delta x^2} u_{N}(t). +\end{aligned} ``` From 86fe176d10a9270cfd2e121b20fc631f405f6a93 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 19:28:05 +0200 Subject: [PATCH 14/15] Update docs/src/heat_equation_dirichlet.md Co-authored-by: Stefan Kopecz --- docs/src/heat_equation_dirichlet.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/heat_equation_dirichlet.md b/docs/src/heat_equation_dirichlet.md index 6148af40..89c69d09 100644 --- a/docs/src/heat_equation_dirichlet.md +++ b/docs/src/heat_equation_dirichlet.md @@ -188,7 +188,7 @@ nothing #hide ``` ```@example HeatEquationDirichlet -plot(x,u0; label = "u0", xguide = "x", yguide = "u") +plot(x, u0; label = "u0", xguide = "x", yguide = "u") plot!(x, last(sol_tridiagonal.u); label = "u") ``` From 5117ee0e01e077d65ce52eb6638453fb6bf935fb Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 19:28:11 +0200 Subject: [PATCH 15/15] Update docs/src/heat_equation_dirichlet.md Co-authored-by: Stefan Kopecz --- docs/src/heat_equation_dirichlet.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/heat_equation_dirichlet.md b/docs/src/heat_equation_dirichlet.md index 89c69d09..b2848d49 100644 --- a/docs/src/heat_equation_dirichlet.md +++ b/docs/src/heat_equation_dirichlet.md @@ -162,7 +162,7 @@ nothing #hide ``` ```@example HeatEquationDirichlet -plot(x,u0; label = "u0", xguide = "x", yguide = "u") +plot(x, u0; label = "u0", xguide = "x", yguide = "u") plot!(x, last(sol_sparse.u); label = "u") ```