production and destruction terms in one function

Project.toml

@@ -1,7 +1,7 @@
 name = "PositiveIntegrators"
 uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5"
 authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"]
-version = "0.2.6"
+version = "0.3.0-DEV"
 
 [deps]
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
@@ -23,7 +23,7 @@ LinearSolve = "2.21"
 MuladdMacro = "0.2.1"
 OrdinaryDiffEq = "6.59"
 Reexport = "1"
-SciMLBase = "2"
+SciMLBase = "2.52"
 SimpleUnPack = "1"
 SparseArrays = "1.7"
 StaticArrays = "1.5"

docs/src/faq.md

@@ -16,7 +16,8 @@ p = spdiagm(0 => ones(4), 1 => zeros(3))
 p .= 2 * p
 ```
 
-Instead, you should be able to use a pattern like the following, where the function `nonzeros` is used to modify the values of a sparse matrix.
+Instead, you should be able to use a pattern like the following, where the function
+`nonzeros` is used to modify the values of a sparse matrix.
 
 ```@repl
 using SparseArrays

docs/src/heat_equation_dirichlet.md

@@ -89,12 +89,13 @@ the resulting linear systems:
 ```@example HeatEquationDirichlet
 using PositiveIntegrators # load ConservativePDSProblem
 
-function heat_eq_P!(P, u, μ, t)
-    fill!(P, 0)
+function heat_eq_PD!(P, D, u, μ, t)
     N = length(u)
     Δx = 1 / N
     μ_Δx2 = μ / Δx^2
 
+    # production terms
+    fill!(P, 0)
     let i = 1
         # Dirichlet boundary condition
         P[i, i + 1] = u[i + 1] * μ_Δx2
@@ -111,15 +112,8 @@ function heat_eq_P!(P, u, μ, t)
         P[i, i - 1] = u[i - 1] * μ_Δx2
     end
 
-    return nothing
-end
-
-function heat_eq_D!(D, u, μ, t)
+    # nonconservative destruction terms
     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
@@ -128,7 +122,7 @@ function heat_eq_D!(D, u, μ, t)
 end
 
 μ = 1.0e-2
-prob = PDSProblem(heat_eq_P!, heat_eq_D!, u0, tspan, μ) # create the PDS
+prob = PDSProblem(heat_eq_PD!, u0, tspan, μ) # create the PDS
 
 sol = solve(prob, MPRK22(1.0); save_everystep = false)
 
@@ -153,7 +147,7 @@ you can use the keyword argument `p_prototype` of
 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, μ;
+prob_sparse = PDSProblem(heat_eq_PD!, u0, tspan, μ;
                          p_prototype = p_prototype)
 
 sol_sparse = solve(prob_sparse, MPRK22(1.0); save_everystep = false)
@@ -179,7 +173,7 @@ 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, μ;
+prob_tridiagonal = PDSProblem(heat_eq_PD!, u0, tspan, μ;
                               p_prototype = p_prototype)
 
 sol_tridiagonal = solve(prob_tridiagonal, MPRK22(1.0); save_everystep = false)

docs/src/index.md

@@ -29,15 +29,22 @@ julia> using Pkg; Pkg.update()
 
 ### Modified Patankar-Runge-Kutta schemes
 
-Modified Patankar-Runge-Kutta (MPRK) schemes are unconditionally positive and conservative time integration schemes for the solution of positive and conservative ODE systems. The application of these methods is based on the representation of the ODE system as a so-called production-destruction system (PDS).
+Modified Patankar-Runge-Kutta (MPRK) schemes are unconditionally positive and
+conservative time integration schemes for the solution of positive and
+conservative ODE systems. The application of these methods is based on the
+representation of the ODE system as a so-called production-destruction system (PDS).
 
 #### Production-destruction systems (PDS)
 
-The application of MPRK schemes requires the ODE system to be represented as a production-destruction system (PDS). A PDS takes the general form
+The application of MPRK schemes requires the ODE system to be represented as
+a production-destruction system (PDS). A PDS takes the general form
 ```math
     u_i'(t) = \sum_{j=1}^N \bigl(p_{ij}(t,\boldsymbol u) - d_{ij}(t,\boldsymbol u)\bigr),\quad i=1,\dots,N,
 ```
-where ``\boldsymbol u=(u_1,\dots,u_n)^T`` is the vector of unknowns and both production terms ``p_{ij}(t,\boldsymbol u)`` and destruction terms ``d_{ij}(t,\boldsymbol u)`` must be nonnegative for all ``i,j=1,\dots,N``. The meaning behind ``p_{ij}`` and ``d_{ij}`` is as follows:
+where ``\boldsymbol u=(u_1,\dots,u_n)^T`` is the vector of unknowns and
+both production terms ``p_{ij}(t,\boldsymbol u)`` and destruction terms
+``d_{ij}(t,\boldsymbol u)`` must be nonnegative for all ``i,j=1,\dots,N``.
+The meaning behind ``p_{ij}`` and ``d_{ij}`` is as follows:
 * ``p_{ij}`` with ``i\ne j`` represents the sum of all nonnegative terms which
   appear in equation ``i`` with a positive sign and in equation ``j`` with a negative sign.
 * ``d_{ij}`` with ``i\ne j`` represents the sum of all nonnegative terms which
@@ -47,7 +54,9 @@ where ``\boldsymbol u=(u_1,\dots,u_n)^T`` is the vector of unknowns and both pro
 * ``d_{ii}`` represents the sum of all negative terms which appear in
   equation ``i`` and don't have a positive counterpart in one of the other equations.
 
-This naming convention leads to ``p_{ij} = d_{ji}`` for ``i≠ j`` and therefore a PDS is completely defined by the production matrix ``\mathbf{P}=(p_{ij})_{i,j=1,\dots,N}`` and the destruction vector ``\mathbf{d}=(d_{ii})_{i=1,\dots,N}``.
+This naming convention leads to ``p_{ij} = d_{ji}`` for ``i≠ j`` and therefore
+a PDS is completely defined by the production matrix ``\mathbf{P}=(p_{ij})_{i,j=1,\dots,N}``
+and the destruction vector ``\mathbf{d}=(d_{ii})_{i=1,\dots,N}``.
 
 As an example we consider the Lotka-Volterra model
 ```math
@@ -68,27 +77,36 @@ d_{22}(u_1,u_2) &= u_2,
 where all remaining production and destruction terms are zero.
 Consequently the production matrix ``\mathbf P`` and destruction vector ``\mathbf d`` are
 ```math
-\mathbf P(u_1,u_2) = \begin{pmatrix}2u_1 & 0\\ u_1u_2 & 0\end{pmatrix},\quad \mathbf d(u_1,u_2) = \begin{pmatrix}0\\ u_2\end{pmatrix}.
+\mathbf P(u_1,u_2) = \begin{pmatrix}2u_1 & 0\\ u_1u_2 & 0\end{pmatrix},\quad
+\mathbf d(u_1,u_2) = \begin{pmatrix}0\\ u_2\end{pmatrix}.
 ```
 
 ```@setup LotkaVolterra
 import Pkg; Pkg.add("OrdinaryDiffEq");  Pkg.add("Plots")
 ```
-To solve this PDS together with initial values ``u_1(0)=u_2(0)=2`` on the time domain ``(0,10)``, we first need to create a `PDSProblem`.
+To solve this PDS together with initial values ``u_1(0)=u_2(0)=2``
+on the time domain ``(0,10)``, we first need to create a `PDSProblem`.
 ```@example LotkaVolterra
 using PositiveIntegrators # load PDSProblem
 
-P(u, p, t) = [2*u[1]  0; u[1]*u[2]  0] # Production matrix
-d(u, p, t) = [0; u[2]] # Destruction vector
+function Pd(u, p, t)
+  P = [2*u[1]  0; u[1]*u[2]  0] # Production matrix
+  d = [0; u[2]] # Destruction vector
+  return P, d
+end
 
 u0 = [2.0; 2.0] # initial values
 tspan = (0.0, 10.0) # time span
 
 # Create PDS
-prob = PDSProblem(P, d, u0, tspan)
+prob = PDSProblem(Pd, u0, tspan)
 nothing #hide
 ```
-Now that the problem has been created, we can solve it with any method of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl). In the following, we use the method `MPRK22(1.0)`. In addition, we could also use any method provided by [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/), but these might possibly generate negative approximations.
+Now that the problem has been created, we can solve it with any method
+of [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl).
+In the following, we use the method `MPRK22(1.0)`. In addition, we could
+also use any method provided by [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/),
+but these might possibly generate negative approximations.
 
 ```@example LotkaVolterra
 sol = solve(prob, MPRK22(1.0))
@@ -162,8 +180,11 @@ tspan = (0.0, 100.0); # time span
 prob = ConservativePDSProblem(P, u0, tspan)
 nothing # hide
 ```
-Since the SIR model is not only conservative but also positive, we can use any scheme from [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl) to solve it. Here we use `MPRK22(1.0)`.
-Please note that any method from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) can be used as well, but might possibly generate negative approximations.
+Since the SIR model is not only conservative but also positive, we can use
+any scheme from [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl)
+to solve it. Here we use `MPRK22(1.0)`.
+Please note that any method from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/)
+can be used as well, but might possibly generate negative approximations.
 
 ```@example SIR
 sol = solve(prob, MPRK22(1.0))

docs/src/linear_advection.md

@@ -1,18 +1,24 @@
 # [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.
+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.
 
 ## Definition of the production-destruction system
 
-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
+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
 
 ```math
 \partial_t u(t,x)=-a\partial_x u(t,x),\quad u(0,x)=u_0(x)
 ```
 
-with ``a>0``, ``t≥ 0``, ``x\in[0,1]`` and periodic boundary conditions. To keep things as simple as possible, we
-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
+with ``a>0``, ``t≥ 0``, ``x\in[0,1]`` and periodic boundary conditions.
+To keep things as simple as possible, we
+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
 
 ```math
 \begin{aligned}
@@ -22,33 +28,44 @@ discretize the space domain as ``0=x_0<x_1\dots <x_{N-1}<x_N=1`` with ``x_i = i
 ```
 
 where ``u_i(t)`` is an approximation of ``u(t,x_i)`` for ``i=1,\dots, N``.
-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
+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
 
 ```math
 \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}.
 ```
 
-In particular the matrix ``\mathbf A`` shows that there is a single production term and a single destruction term per equation.
+In particular the matrix ``\mathbf A`` shows that there is a single production
+term and a single destruction term per equation.
 Furthermore, the system is conservative as ``\mathbf A`` has column sum zero.
-To be precise, the production matrix ``\mathbf P = (p_{i,j})`` of this conservative PDS is given by
+To be precise, the production matrix ``\mathbf P = (p_{i,j})`` of this
+conservative PDS is given by
 
 ```math
 \begin{aligned}
 &p_{1,N}(t,\mathbf u(t)) = \frac{a}{Δ x}u_N(t),\\
 &p_{i,i-1}(t,\mathbf u(t)) = \frac{a}{Δ x}u_{i-1}(t),\quad i=2,\dots,N.
 \end{aligned}
 ```
-In addition, all production and destruction 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 ``\mathbf P``.
+In addition, all production and destruction 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 ``\mathbf P``.
 
 ## Solution of the production-destruction system
 
-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/). In the following we use ``a=1``, ``N=1000`` and the time domain ``t\in[0,1]``. Moreover, we choose the step function
+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/).
+In the following we use ``a=1``, ``N=1000`` and the time domain ``t\in[0,1]``.
+Moreover, we choose the step function
 
 ```math
 u_0(x)=\begin{cases}1, & 0.4 ≤ x ≤ 0.6,\\ 0,& \text{elsewhere}\end{cases}
 ```
 
-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)``.
+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)``.
 
 ```@example LinearAdvection
 N = 1000 # number of subintervals
@@ -60,7 +77,8 @@ tspan = (0.0, 1.0) # time domain
 nothing #hide
 ```
 
-As mentioned above, we will try different approaches to solve this PDS and compare their efficiency. These are
+As mentioned above, we will try different approaches to solve this PDS
+and compare their efficiency. These are
 1. an in-place implementation with a dense matrix,
 2. an in-place implementation with a sparse matrix.
 
@@ -98,7 +116,8 @@ plot(x, u0; label = "u0", xguide = "x", yguide = "u")
 plot!(x, last(sol.u); label = "u")
 ```
 
-We can use [`isnonnegative`](@ref) to check that the computed solution is nonnegative,  as expected from an MPRK scheme.
+We can use [`isnonnegative`](@ref) to check that the computed solution
+is nonnegative,  as expected from an MPRK scheme.
 ```@example LinearAdvection
 isnonnegative(sol)
 ```

docs/src/npzd_model.md

@@ -1,11 +1,16 @@
 # [Tutorial: Solution of an NPZD model](@id tutorial-npzd)
 
-This tutorial is about the efficient solution of production-destruction systems (PDS) with a small number of differential equations.
-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.
+This tutorial is about the efficient solution of production-destruction systems (PDS)
+with a small number of differential equations.
+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.
 
 ## Definition of the production-destruction system
 
-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
+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
 ```math
 \begin{aligned}
 N' &= 0.01P + 0.01Z + 0.003D - \frac{NP}{0.01 + N},\\
@@ -14,7 +19,8 @@ Z' &= 0.5(1 - e^{-1.21P^2})Z - 0.01Z - 0.02Z,\\
 D' &= 0.05P + 0.02Z - 0.003D,
 \end{aligned}
 ```
-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]``. 
+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]``.
 
 The model can be represented as a conservative PDS with production terms
 ```math
@@ -24,20 +30,27 @@ p_{21} &= \frac{NP}{0.01 + N}, & p_{32} &= 0.5  (1.0 - e^{-1.21  P^2})  Z,& p_{4
 p_{43} &= 0.02  Z,
 \end{aligned}
 ```
-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``.
+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``.
 
 ## Solution of the production-destruction system
 
-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/). 
+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/).
 
-As mentioned above, we will try different approaches to solve this PDS and compare their efficiency. These are
+As mentioned above, we will try different approaches to solve this PDS and
+compare their efficiency. These are
 1. an out-of-place implementation with standard (dynamic) matrices and vectors,
 2. an in-place implementation with standard (dynamic) matrices and vectors,
 3. an out-of-place implementation with static matrices and vectors from [StaticArrays.jl](https://juliaarrays.github.io/StaticArrays.jl/stable/).
 
 ### Standard out-of-place implementation
 
-Here we create a function to compute the production matrix with return type `Matrix{Float64}`.
+Here we create a function to compute the production matrix with
+return type `Matrix{Float64}`.
 
 ```@example NPZD
 using PositiveIntegrators # load ConservativePDSProblem
@@ -70,13 +83,16 @@ sol_oop = solve(prob_oop, MPRK43I(1.0, 0.5))
 
 nothing #hide
 ```
-Plotting the solution shows that the components ``N`` and ``P`` are in danger of becoming negative. 
+Plotting the solution shows that the components ``N`` and ``P`` are
+in danger of becoming negative.
 ```@example NPZD
 using Plots
 
 plot(sol_oop; label = ["N" "P" "Z" "D"], xguide = "t")
 ```
-[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.
+[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.
 ```@example NPZD
 isnonnegative(sol_oop)
 ```
@@ -131,7 +147,10 @@ sol_oop.t ≈ sol_ip.t && sol_oop.u ≈ sol_ip.u
 ```
 
 ### Using static arrays
-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.
+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.
 
 ```@example NPZD
 using StaticArrays
@@ -173,7 +192,9 @@ This solution is also nonnegative.
 isnonnegative(sol_static)
 ```
 
-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).
+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).
 
 ### Performance comparison