|
| 1 | +# [Tutorial: Solution of an NPZD model](@id tutorial-npzd) |
| 2 | + |
| 3 | +This tutorial is about the efficient solution of production-destruction systems (PDS) with a small number of differential equations. |
| 4 | +We will compare the use of standard arrays and static arrays from [StaticArrays.jl](https://juliaarrays.github.io/StaticArrays.jl/stable/) and assess their efficiency. |
| 5 | + |
| 6 | +## Definition of the production-destruction system |
| 7 | + |
| 8 | +The NPZD model we want to solve was described by Burchard, Deleersnijder and Meister in [Application of modified Patankar schemes to stiff biogeochemical models for the water column](https://doi.org/10.1007/s10236-005-0001-x). The model reads |
| 9 | +```math |
| 10 | +\begin{aligned} |
| 11 | +N' &= 0.01P + 0.01Z + 0.003D - \frac{NP}{0.01 + N},\\ |
| 12 | +P' &= \frac{NP}{0.01 + N}- 0.01P - 0.5( 1 - e^{-1.21P^2})Z - 0.05P,\\ |
| 13 | +Z' &= 0.5(1 - e^{-1.21P^2})Z - 0.01Z - 0.02Z,\\ |
| 14 | +D' &= 0.05P + 0.02Z - 0.003D, |
| 15 | +\end{aligned} |
| 16 | +``` |
| 17 | +and we consider the initial conditions ``N=8``, ``P=2``, ``Z=1`` and ``D=4``. The time domain of interest is ``t\in[0,10]``. |
| 18 | + |
| 19 | +The model can be represented as a conservative PDS with production terms |
| 20 | +```math |
| 21 | +\begin{aligned} |
| 22 | +p_{12} &= 0.01 P, & p_{13} &= 0.01 Z, & p_{14} &= 0.003 D,\\ |
| 23 | +p_{21} &= \frac{NP}{0.01 + N}, & p_{32} &= 0.5 (1.0 - e^{-1.21 P^2}) Z,& p_{42} &= 0.05 P,\\ |
| 24 | +p_{43} &= 0.02 Z, |
| 25 | +\end{aligned} |
| 26 | +``` |
| 27 | +whereby production terms not listed have the value zero. Since the PDS is conservative, we have ``d_{i,j}=p_{j,i}`` and the system is fully determined by the production matrix ``(p_{ij})_{i,j=1}^4``. |
| 28 | + |
| 29 | +## Solution of the production-destruction system |
| 30 | + |
| 31 | +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/). |
| 32 | + |
| 33 | +As mentioned above, we will try different approaches to solve this PDS and compare their efficiency. These are |
| 34 | +1. an out-of-place implementation with standard (dynamic) matrices and vectors, |
| 35 | +2. an in-place implementation with standard (dynamic) matrices and vectors, |
| 36 | +3. an out-of-place implementation with static matrices and vectors from [StaticArrays.jl](https://juliaarrays.github.io/StaticArrays.jl/stable/). |
| 37 | + |
| 38 | +### Standard out-of-place implementation |
| 39 | + |
| 40 | +Here we create a function to compute the production matrix with return type `Matrix{Float64}`. |
| 41 | + |
| 42 | +```@example NPZD |
| 43 | +using PositiveIntegrators # load ConservativePDSProblem |
| 44 | +
|
| 45 | +function prod(u, p, t) |
| 46 | + N, P, Z, D = u |
| 47 | +
|
| 48 | + p12 = 0.01 * P |
| 49 | + p13 = 0.01 * Z |
| 50 | + p14 = 0.003 * D |
| 51 | + p21 = N / (0.01 + N) * P |
| 52 | + p32 = 0.5 * (1.0 - exp(-1.21 * P^2)) * Z |
| 53 | + p42 = 0.05 * P |
| 54 | + p43 = 0.02 * Z |
| 55 | +
|
| 56 | + return [0.0 p12 p13 p14; |
| 57 | + p21 0.0 0.0 0.0; |
| 58 | + 0.0 p32 0.0 0.0; |
| 59 | + 0.0 p42 p43 0.0] |
| 60 | +end |
| 61 | +nothing #hide |
| 62 | +``` |
| 63 | +The solution of the NPZD model can now be computed as follows. |
| 64 | +```@example NPZD |
| 65 | +u0 = [8.0, 2.0, 1.0, 4.0] # initial values |
| 66 | +tspan = (0.0, 10.0) # time domain |
| 67 | +prob_oop = ConservativePDSProblem(prod, u0, tspan) # create the PDS |
| 68 | +
|
| 69 | +sol_oop = solve(prob_oop, MPRK43I(1.0, 0.5)) |
| 70 | +
|
| 71 | +nothing #hide |
| 72 | +``` |
| 73 | +Plotting the solution shows that the components ``N`` and ``P`` are in danger of becoming negative. |
| 74 | +```@example NPZD |
| 75 | +using Plots |
| 76 | +
|
| 77 | +plot(sol_oop; label = ["N" "P" "Z" "D"], xguide = "t") |
| 78 | +``` |
| 79 | +[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) provides the function [`isnonnegative`](@ref) (and also [`isnegative`](@ref)) to check if the solution is actually nonnegative, as expected from an MPRK scheme. |
| 80 | +```@example NPZD |
| 81 | +isnonnegative(sol_oop) |
| 82 | +``` |
| 83 | + |
| 84 | +### Standard in-place implementation |
| 85 | + |
| 86 | +Next we create an in-place function for the production matrix. |
| 87 | + |
| 88 | +```@example NPZD |
| 89 | +
|
| 90 | +function prod!(PMat, u, p, t) |
| 91 | + N, P, Z, D = u |
| 92 | +
|
| 93 | + p12 = 0.01 * P |
| 94 | + p13 = 0.01 * Z |
| 95 | + p14 = 0.003 * D |
| 96 | + p21 = N / (0.01 + N) * P |
| 97 | + p32 = 0.5 * (1.0 - exp(-1.21 * P^2)) * Z |
| 98 | + p42 = 0.05 * P |
| 99 | + p43 = 0.02 * Z |
| 100 | +
|
| 101 | + fill!(PMat, zero(eltype(PMat))) |
| 102 | +
|
| 103 | + PMat[1, 2] = p12 |
| 104 | + PMat[1, 3] = p13 |
| 105 | + PMat[1, 4] = p14 |
| 106 | + PMat[2, 1] = p21 |
| 107 | + PMat[3, 2] = p32 |
| 108 | + PMat[4, 2] = p42 |
| 109 | + PMat[4, 3] = p43 |
| 110 | +
|
| 111 | + return nothing |
| 112 | +end |
| 113 | +nothing #hide |
| 114 | +``` |
| 115 | + |
| 116 | +The solution of the in-place implementation of the NPZD model can now be computed as follows. |
| 117 | +```@example NPZD |
| 118 | +
|
| 119 | +prob_ip = ConservativePDSProblem(prod!, u0, tspan) |
| 120 | +sol_ip = solve(prob_ip, MPRK43I(1.0, 0.5)) |
| 121 | +nothing #hide |
| 122 | +``` |
| 123 | +```@example NPZD |
| 124 | +
|
| 125 | +plot(sol_ip; label = ["N" "P" "Z" "D"], xguide = "t") |
| 126 | +``` |
| 127 | + |
| 128 | +We also check that the in-place and out-of-place solutions are equivalent. |
| 129 | +```@example NPZD |
| 130 | +sol_oop.t ≈ sol_ip.t && sol_oop.u ≈ sol_ip.u |
| 131 | +``` |
| 132 | + |
| 133 | +### Using static arrays |
| 134 | +For PDS with a small number of differential equations like the NPZD model the use of static arrays will be more efficient. To create a function which computes the production matrix and returns a static matrix, we only need to add the `@SMatrix` macro. |
| 135 | + |
| 136 | +```@example NPZD |
| 137 | +using StaticArrays |
| 138 | +
|
| 139 | +function prod_static(u, p, t) |
| 140 | + N, P, Z, D = u |
| 141 | +
|
| 142 | + p12 = 0.01 * P |
| 143 | + p13 = 0.01 * Z |
| 144 | + p14 = 0.003 * D |
| 145 | + p21 = N / (0.01 + N) * P |
| 146 | + p32 = 0.5 * (1.0 - exp(-1.21 * P^2)) * Z |
| 147 | + p42 = 0.05 * P |
| 148 | + p43 = 0.02 * Z |
| 149 | +
|
| 150 | + return @SMatrix [0.0 p12 p13 p14; |
| 151 | + p21 0.0 0.0 0.0; |
| 152 | + 0.0 p32 0.0 0.0; |
| 153 | + 0.0 p42 p43 0.0] |
| 154 | +end |
| 155 | +nothing #hide |
| 156 | +``` |
| 157 | +In addition we also want to use a static vector to hold the initial conditions. |
| 158 | +```@example NPZD |
| 159 | +u0_static = @SVector [8.0, 2.0, 1.0, 4.0] # initial values |
| 160 | +prob_static = ConservativePDSProblem(prod_static, u0_static, tspan) # create the PDS |
| 161 | +
|
| 162 | +sol_static = solve(prob_static, MPRK43I(1.0, 0.5)) |
| 163 | +
|
| 164 | +nothing #hide |
| 165 | +``` |
| 166 | +```@example NPZD |
| 167 | +using Plots |
| 168 | +
|
| 169 | +plot(sol_static; label = ["N" "P" "Z" "D"], xguide = "t") |
| 170 | +``` |
| 171 | +This solution is also nonnegative. |
| 172 | +```@example NPZD |
| 173 | +isnonnegative(sol_static) |
| 174 | +``` |
| 175 | + |
| 176 | +The above implementation of the NPZD model using `StaticArrays` can also be found in the [Example Problems](https://skopecz.github.io/PositiveIntegrators.jl/dev/api_reference/#Example-problems) as [`prob_pds_npzd`](@ref). |
| 177 | + |
| 178 | +### Performance comparison |
| 179 | + |
| 180 | +Finally, we use [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) |
| 181 | +to show the benefit of using static arrays. |
| 182 | + |
| 183 | +```@example NPZD |
| 184 | +using BenchmarkTools |
| 185 | +@benchmark solve(prob_oop, MPRK43I(1.0, 0.5)) |
| 186 | +``` |
| 187 | + |
| 188 | +```@example NPZD |
| 189 | +using BenchmarkTools |
| 190 | +@benchmark solve(prob_ip, MPRK43I(1.0, 0.5)) |
| 191 | +``` |
| 192 | + |
| 193 | +```@example NPZD |
| 194 | +@benchmark solve(prob_static, MPRK43I(1.0, 0.5)) |
| 195 | +``` |
| 196 | + |
| 197 | +## Package versions |
| 198 | + |
| 199 | +These results were obtained using the following versions. |
| 200 | +```@example NPZD |
| 201 | +using InteractiveUtils |
| 202 | +versioninfo() |
| 203 | +println() |
| 204 | +
|
| 205 | +using Pkg |
| 206 | +Pkg.status(["PositiveIntegrators", "StaticArrays", "LinearSolve", "OrdinaryDiffEq"], |
| 207 | + mode=PKGMODE_MANIFEST) |
| 208 | +nothing # hide |
| 209 | +``` |
0 commit comments