diff --git a/Project.toml b/Project.toml index 48f7eec0..3ca51dc4 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.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" diff --git a/docs/src/faq.md b/docs/src/faq.md index a112682b..19f472e9 100644 --- a/docs/src/faq.md +++ b/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 diff --git a/docs/src/heat_equation_dirichlet.md b/docs/src/heat_equation_dirichlet.md index e8015e2f..1b48b3cb 100644 --- a/docs/src/heat_equation_dirichlet.md +++ b/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) diff --git a/docs/src/index.md b/docs/src/index.md index 5e0a8aed..39211935 100644 --- a/docs/src/index.md +++ b/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)) diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md index 73e8ea74..c167efa9 100644 --- a/docs/src/linear_advection.md +++ b/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_00``, ``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, 1e4 .* y), 0, 2), (0, 3)], label = ["u₁" "10⁴u₂" "u₃"]) ``` -[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 robertson isnonnegative(sol) ``` ### Using callbacks to solve the Robertson problem with non-adatpive schemes -The `SSPMPRK43()` scheme is only available with fixed time stepping. With a scheme like this, it would take a huge amount of time to solve the Robertson problem, since the time step must be chosen very small to accurately solve the problem in its initial phase. However, the use of a `callback` allows us to modify the time step size after each step, which makes a solution with a fixed step method possible. +The `SSPMPRK43()` scheme is only available with fixed time stepping. With a scheme +like this, it would take a huge amount of time to solve the Robertson problem, +since the time step must be chosen very small to accurately solve the problem +in its initial phase. However, the use of a `callback` allows us to modify the +time step size after each step, which makes a solution with a fixed step method possible. -In the following example the `callback` increases the time step size by a factor of 1.5 after each time step. +In the following example the `callback` increases the time step size by +a factor of 1.5 after each time step. ```@example robertson using OrdinaryDiffEq diff --git a/docs/src/stratospheric_reaction.md b/docs/src/stratospheric_reaction.md index ef18a7b0..92a08548 100644 --- a/docs/src/stratospheric_reaction.md +++ b/docs/src/stratospheric_reaction.md @@ -1,11 +1,18 @@ # [Tutorial: Solution of a stratospheric reaction problem](@id tutorial-stratos) -This tutorial is about the efficient solution of a stiff non-autonomous and non-conservative 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 a stiff non-autonomous +and non-conservative 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 -This stratospheric reaction problem was described by Adrian Sandu in [Positive Numerical Integration Methods for Chemical Kinetic Systems](https://doi.org/10.1006/jcph.2001.6750), see also the paper [Positivity-preserving adaptive Runge–Kutta methods](https://doi.org/10.2140/camcos.2021.16.155) by Stefan Nüßlein, Hendrik Ranocha and David I. Ketcheson. The governing equations are +This stratospheric reaction problem was described by Adrian Sandu in +[Positive Numerical Integration Methods for Chemical Kinetic Systems](https://doi.org/10.1006/jcph.2001.6750), +see also the paper [Positivity-preserving adaptive Runge–Kutta methods](https://doi.org/10.2140/camcos.2021.16.155) +by Stefan Nüßlein, Hendrik Ranocha and David I. Ketcheson. The governing equations are ```math \begin{aligned} \frac{dO^{1D}}{dt} &= r_5 - r_6 - r_7,\\ @@ -32,10 +39,16 @@ T &= t/3600 \mod 24,\quad T_r=4.5,\quad T_s = 19.5,\\ σ(T) &= \begin{cases}1, & T_r≤ T≤ T_s,\\0, & \text{otherwise}.\end{cases} \end{aligned} ``` -Setting ``\mathbf u = (O^{1D}, O, O_3, O_2, NO, NO_2)`` the initial value is ``\mathbf{u}_0 = (9.906⋅10^1, 6.624⋅10^8, 5.326⋅10^{11}, 1.697⋅10^{16}, 4⋅10^6, 1.093⋅10^9)^T``. The time domain in seconds is ``[4.32⋅10^{4}, 3.024⋅10^5]``, which corresponds to ``[12.0, 84.0]`` in hours. -There are two independent linear invariants, e.g. ``u_1+u_2+3u_3+2u_4+u_5+2u_6=(1,1,3,2,1,2)\cdot\mathbf{u}_0`` and ``u_5+u_6 = 1.097⋅10^9``. - -The stratospheric reaction problem can be represented as a (non-conservative) PDS with production terms +Setting ``\mathbf u = (O^{1D}, O, O_3, O_2, NO, NO_2)`` the initial value +is ``\mathbf{u}_0 = (9.906⋅10^1, 6.624⋅10^8, 5.326⋅10^{11}, 1.697⋅10^{16}, 4⋅10^6, 1.093⋅10^9)^T``. +The time domain in seconds is ``[4.32⋅10^{4}, 3.024⋅10^5]``, which +corresponds to ``[12.0, 84.0]`` in hours. +There are two independent linear invariants, e.g. +``u_1+u_2+3u_3+2u_4+u_5+2u_6=(1,1,3,2,1,2)\cdot\mathbf{u}_0`` +and ``u_5+u_6 = 1.097⋅10^9``. + +The stratospheric reaction problem can be represented as a +(non-conservative) PDS with production terms ```math \begin{aligned} p_{13} &= r_5, & p_{21} &= r_6, & p_{22} &= r_1+r_{10},\\ @@ -54,21 +67,26 @@ In addition, all production and destruction terms not listed have the value zero ## Solution of the 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/). +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/). -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 an out-of-place function to compute the production matrix with return type `Matrix{Float64}` and a second out-of-place function for the destruction vector with return type `Vector{Float64}`. +Here we create an out-of-place function to compute the production matrix +with return type `Matrix{Float64}` and the destruction vector with return +type `Vector{Float64}`. ```@example stratreac using PositiveIntegrators # load PDSProblem -function prod(u, p, t) +function prod_dest(u, p, t) O1D, O, O3, O2, NO, NO2 = u Tr = 4.5 @@ -107,38 +125,30 @@ function prod(u, p, t) r10 = k10 * NO2 r11 = k11 * NO * O - return [0.0 0.0 r5 0.0 0.0 0.0; - r6 r1+r10 r3 r1 0.0 0.0; - 0.0 r2 0.0 0.0 0.0 0.0; - r7 r4+r9 r4+r7+r8 r3+r5 0.0 0.0; - 0.0 0.0 0.0 0.0 0.0 r9+r10; - 0.0 0.0 0.0 0.0 r8+r11 0.0] + P = [0.0 0.0 r5 0.0 0.0 0.0; + r6 r1+r10 r3 r1 0.0 0.0; + 0.0 r2 0.0 0.0 0.0 0.0; + r7 r4+r9 r4+r7+r8 r3+r5 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 r9+r10; + 0.0 0.0 0.0 0.0 r8+r11 0.0] + D = [0.0, r11, 0.0, r2, 0.0, 0.0] + return P, D end -function dest(u, p, t) - O1D, O, O3, O2, NO, NO2 = u - - k2 = 8.018e-17 - k11 = 1.0e-8 - - r2 = k2 * O * O2 - r11 = k11 * NO * O - - return [0.0, r11, 0.0, r2, 0.0, 0.0] -end nothing #hide ``` The solution of the stratospheric reaction problem can now be computed as follows. ```@example stratreac u0 = [9.906e1, 6.624e8, 5.326e11, 1.697e16, 4e6, 1.093e9] # initial values tspan = (4.32e4, 3.024e5) # time domain -prob_oop = PDSProblem(prod, dest, u0, tspan) # create the PDS +prob_oop = PDSProblem(prod_dest, u0, tspan) # create the PDS sol_oop = solve(prob_oop, MPRK43I(1.0, 0.5)) nothing #hide ``` -Plotting the solution shows that the components ``O¹ᴰ``, ``O`` and ``NO`` are in danger of becoming negative. +Plotting the solution shows that the components ``O¹ᴰ``, ``O`` and +``NO`` are in danger of becoming negative. ```@example stratreac using Plots @@ -149,11 +159,13 @@ plot(sol_oop, xticks = (range(first(tspan), last(tspan), 4), range(12.0, 84.0, 4)), yguide=["O¹ᴰ" "O" "O₃" "O₂" "NO" "NO₂"], tickfontsize = 7, - legend = :none, + legend = :none, widen = true ) ``` -[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 stratreac isnonnegative(sol_oop) ``` @@ -164,7 +176,7 @@ Next we create in-place functions for the production matrix and the destruction ```@example stratreac -function prod!(P, u, p, t) +function prod_dest!(P, D, u, p, t) O1D, O, O3, O2, NO, NO2 = u Tr = 4.5 @@ -216,31 +228,20 @@ function prod!(P, u, p, t) P[4, 4] = r3 + r5 P[5, 6] = r9 + r10 P[6, 5] = r8 + r11 - return nothing -end - -function dest!(D, u, p, t) - O1D, O, O3, O2, NO, NO2 = u - - k2 = 8.018e-17 - k11 = 1.0e-8 - - r2 = k2 * O * O2 - r11 = k11 * NO * O fill!(D, zero(eltype(D))) - D[2] = r11 D[4] = r2 - return nothing end + nothing #hide ``` -The solution of the in-place implementation of the stratospheric reaction problem can now be computed as follows. +The solution of the in-place implementation of the stratospheric +reaction problem can now be computed as follows. ```@example stratreac -prob_ip = PDSProblem(prod!, dest!, u0, tspan) # create the PDS +prob_ip = PDSProblem(prod_dest!, u0, tspan) # create the PDS sol_ip = solve(prob_ip, MPRK43I(1.0, 0.5)) nothing #hide ``` @@ -253,7 +254,7 @@ plot(sol_ip, xticks = (range(first(tspan), last(tspan), 4), range(12.0, 84.0, 4)), yguide=["O¹ᴰ" "O" "O₃" "O₂" "NO" "NO₂"], tickfontsize = 7, - legend = :none, + legend = :none, widen = true ) ``` @@ -264,12 +265,16 @@ 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 stratospheric reaction 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. Accordingly, we use the `@SVector` macro for the destruction vector. +For PDS with a small number of differential equations like the stratospheric +reaction 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. Accordingly, we use the `@SVector` +macro for the destruction vector. ```@example stratreac using StaticArrays -function prod_static(u, p, t) +function prod_dest_static(u, p, t) O1D, O, O3, O2, NO, NO2 = u Tr = 4.5 @@ -308,31 +313,22 @@ function prod_static(u, p, t) r10 = k10 * NO2 r11 = k11 * NO * O - return @SMatrix [0.0 0.0 r5 0.0 0.0 0.0; - r6 r1+r10 r3 r1 0.0 0.0; - 0.0 r2 0.0 0.0 0.0 0.0; - r7 r4+r9 r4+r7+r8 r3+r5 0.0 0.0; - 0.0 0.0 0.0 0.0 0.0 r9+r10; - 0.0 0.0 0.0 0.0 r8+r11 0.0] + P = @SMatrix [0.0 0.0 r5 0.0 0.0 0.0; + r6 r1+r10 r3 r1 0.0 0.0; + 0.0 r2 0.0 0.0 0.0 0.0; + r7 r4+r9 r4+r7+r8 r3+r5 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 r9+r10; + 0.0 0.0 0.0 0.0 r8+r11 0.0] + D = @SVector [0.0, r11, 0.0, r2, 0.0, 0.0] + return P, D end -function dest_static(u, p, t) - O1D, O, O3, O2, NO, NO2 = u - - k2 = 8.018e-17 - k11 = 1.0e-8 - - r2 = k2 * O * O2 - r11 = k11 * NO * O - - return @SVector [0.0, r11, 0.0, r2, 0.0, 0.0] -end nothing #hide ``` In addition we also want to use a static vector to hold the initial conditions. ```@example stratreac u0_static = @SVector [9.906e1, 6.624e8, 5.326e11, 1.697e16, 4e6, 1.093e9] # initial values -prob_static = PDSProblem(prod_static, dest_static, u0_static, tspan) # create the PDS +prob_static = PDSProblem(prod_dest_static, u0_static, tspan) # create the PDS sol_static = solve(prob_static, MPRK43I(1.0, 0.5)) @@ -353,15 +349,20 @@ plot(sol_static, xticks = (range(first(tspan), last(tspan), 4), range(12.0, 84.0, 4)), yguide=["O¹ᴰ" "O" "O₃" "O₂" "NO" "NO₂"], tickfontsize = 7, - legend = :none, + legend = :none, widen = true ) ``` -The above implementation of the stratospheric reaction problem using `StaticArrays` can also be found in the [Example Problems](https://skopecz.github.io/PositiveIntegrators.jl/dev/api_reference/#Example-problems) as [`prob_pds_stratreac`](@ref). +The above implementation of the stratospheric reaction problem using +`StaticArrays` can also be found in the +[Example Problems](https://skopecz.github.io/PositiveIntegrators.jl/dev/api_reference/#Example-problems) +as [`prob_pds_stratreac`](@ref). ### Preservation of linear invariants -As MPRK schemes do not preserve general linear invariants, especially when applied to non-conservative PDS, we compute and plot the relative errors with respect to both linear invariants to see how well these are preserved. +As MPRK schemes do not preserve general linear invariants, especially when +applied to non-conservative PDS, we compute and plot the relative errors +with respect to both linear invariants to see how well these are preserved. ```@example stratreac linear_invariant(a, u) = sum(a .* u) @@ -376,13 +377,16 @@ a2 = [0; 0; 0; 0; 1; 1] # second linear invariant p1 = plot(sol_oop.t, relerr_lininv(a1, u0, sol_oop)) p2 = plot(sol_oop.t, relerr_lininv(a2, u0, sol_oop)) -plot(p1, p2, +plot(p1, p2, xticks = (range(first(tspan), last(tspan), 4), range(12.0, 84.0, 4)), legend = :none) ``` -In contrast to MPRK schemes, Runge-Kutta and Rosenbrock methods preserve all linear invariants, but are not guaranteed to generate nonnegative solutions. -One way to enforce nonnegative solutions of such schemes is passing [`isnegative`](@ref) to the solver option [`isoutofdomain`](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). We show this using the Rosenbrock scheme `Rosenbrock23` as an example. +In contrast to MPRK schemes, Runge-Kutta and Rosenbrock methods preserve +all linear invariants, but are not guaranteed to generate nonnegative solutions. +One way to enforce nonnegative solutions of such schemes is passing [`isnegative`](@ref) +to the solver option [`isoutofdomain`](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/). +We show this using the Rosenbrock scheme `Rosenbrock23` as an example. ```@example stratreac using OrdinaryDiffEq diff --git a/src/PDSProblemLibrary.jl b/src/PDSProblemLibrary.jl index fb05d2fd..938c8bb3 100644 --- a/src/PDSProblemLibrary.jl +++ b/src/PDSProblemLibrary.jl @@ -75,7 +75,7 @@ u_3' &= 0.3 u_2, \\end{aligned} ``` with initial value ``\\mathbf{u}_0 = (9.98, 0.01, 0.01)^T`` and time domain ``(0.0, 30.0)``. -There is one independent linear invariant, e.g. ``u_1+u_2+u_3 = 10.0``. +There is one independent linear invariant, e.g. ``u_1 + u_2 + u_3 = 10.0``. ## References @@ -291,7 +291,7 @@ There is one independent linear invariant, e.g. ``u_1+u_2+u_3+u_4 = 15.0``. prob_pds_npzd = ConservativePDSProblem(P_npzd, u0_npzd, (0.0, 10.0), std_rhs = f_npzd) # stratospheric reaction problem -function P_stratreac(u, p, t) +function PD_stratreac(u, p, t) O1D, O, O3, O2, NO, NO2 = u Tr = 4.5 @@ -330,23 +330,14 @@ function P_stratreac(u, p, t) r10 = k10 * NO2 r11 = k11 * NO * O - return @SMatrix [0.0 0.0 r5 0.0 0.0 0.0; - r6 r1+r10 r3 r1 0.0 0.0; - 0.0 r2 0.0 0.0 0.0 0.0; - r7 r4+r9 r4+r7+r8 r3+r5 0.0 0.0; - 0.0 0.0 0.0 0.0 0.0 r9+r10; - 0.0 0.0 0.0 0.0 r8+r11 0.0] -end -function d_stratreac(u, p, t) - O1D, O, O3, O2, NO, NO2 = u - - k2 = 8.018e-17 - k11 = 1.0e-8 - - r2 = k2 * O * O2 - r11 = k11 * NO * O - - return @SVector [0.0, r11, 0.0, r2, 0.0, 0.0] + P = @SMatrix [0.0 0.0 r5 0.0 0.0 0.0; + r6 r1+r10 r3 r1 0.0 0.0; + 0.0 r2 0.0 0.0 0.0 0.0; + r7 r4+r9 r4+r7+r8 r3+r5 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0 r9+r10; + 0.0 0.0 0.0 0.0 r8+r11 0.0] + d = @SVector [0.0, r11, 0.0, r2, 0.0, 0.0] + return P, d end function f_stratreac(u, p, t) O1D, O, O3, O2, NO, NO2 = u @@ -436,7 +427,7 @@ There are two independent linear invariants, e.g. ``u_1+u_2+3u_3+2u_4+u_5+2u_6=( Communications in Applied Mathematics and Computer Science 16 (2021): 155-179. [DOI: 10.2140/camcos.2021.16.155](https://doi.org/10.2140/camcos.2021.16.155) """ -prob_pds_stratreac = PDSProblem(P_stratreac, d_stratreac, u0_stratreac, (4.32e4, 3.024e5), +prob_pds_stratreac = PDSProblem(PD_stratreac, u0_stratreac, (4.32e4, 3.024e5), std_rhs = f_stratreac) # mapk problem @@ -456,7 +447,7 @@ function f_minmapk(u, p, t) k3 * u[1] * u[3] - k4 * u[5]; k7 * u[1] - k6 * u[6]] end -function P_minmapk(u, p, t) +function PD_minmapk(u, p, t) k1 = 100 / 3 k2 = 1 / 3 k3 = 50 @@ -465,21 +456,17 @@ function P_minmapk(u, p, t) k6 = 0.1 k7 = 0.7 - return @SMatrix [0 0 0 k2*u[4] 0 k6*u[6] - 0 0 k5*u[3] 0 0 0; - 0 0 k2*u[4] 0 k4*u[5] 0; - k1*u[1]*u[2] 0 0 0 0 0; - 0 0 k3*u[1]*u[3] 0 0 0; - k7*u[1] 0 0 0 0 0] -end -function D_minmapk(u, p, t) - k1 = 100 / 3 - - return @SVector [0; k1 * u[1] * u[2]; 0; 0; 0; 0] + P = @SMatrix [0 0 0 k2*u[4] 0 k6*u[6] + 0 0 k5*u[3] 0 0 0; + 0 0 k2*u[4] 0 k4*u[5] 0; + k1*u[1]*u[2] 0 0 0 0 0; + 0 0 k3*u[1]*u[3] 0 0 0; + k7*u[1] 0 0 0 0 0] + d = @SVector [0; k1 * u[1] * u[2]; 0; 0; 0; 0] + return P, d end -u0 = @SVector [0.1; 0.175; 0.15; 1.15; 0.81; 0.5] -tspan = (0.0, 200.0) +u0_minmapk = @SVector [0.1; 0.175; 0.15; 1.15; 0.81; 0.5] """ prob_pds_minmapk @@ -518,4 +505,4 @@ There are two independent linear invariants, e.g. ``u_1+u_4+u_6=1.75`` and ``u_2 PLoS ONE 12 (2017): e0178457. [DOI: 10.1371/journal.pone.0178457](https://doi.org/10.1371/journal.pone.0178457) """ -prob_pds_minmapk = PDSProblem(P_minmapk, D_minmapk, u0, tspan; std_rhs = f_minmapk) +prob_pds_minmapk = PDSProblem(PD_minmapk, u0_minmapk, (0.0, 200.0); std_rhs = f_minmapk) diff --git a/src/mprk.jl b/src/mprk.jl index cd369eea..0f03edf4 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -298,8 +298,13 @@ end f = integrator.f - # evaluate production matrix - P = f.p(uprev, p, t) + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P, d = f.pd(uprev, p, t) + else # f isa ConservativePDSFunction + P = f.p(uprev, p, t) + end integrator.stats.nf += 1 # avoid division by zero due to zero Patankar weights @@ -307,12 +312,10 @@ end # build linear system matrix and rhs if f isa PDSFunction - d = f.d(uprev, p, t) # evaluate nonconservative destruction terms rhs = uprev + dt * diag(P) M = build_mprk_matrix(P, σ, dt, d) linprob = LinearProblem(M, rhs) - else - # f isa ConservativePDSFunction + else # f isa ConservativePDSFunction M = build_mprk_matrix(P, σ, dt) linprob = LinearProblem(M, uprev) end @@ -388,8 +391,7 @@ end # We require the users to set unused entries to zero! - f.p(P, uprev, p, t) # evaluate production terms - f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + f.pd(P, D, uprev, p, t) # evaluate production and destruction terms integrator.stats.nf += 1 # avoid division by zero due to zero Patankar weights @@ -405,9 +407,9 @@ end # Same as linres = P \ linsolve_rhs linsolve.A = P linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 end @muladd function perform_step!(integrator, cache::MPEConservativeCache, repeat_step = false) @@ -431,9 +433,9 @@ end # Same as linres = P \ uprev linsolve.A = P linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 end ### MPRK22 ##################################################################################### @@ -554,23 +556,26 @@ end f = integrator.f - # evaluate production matrix - P = f.p(uprev, p, t) - Ptmp = a21 * P + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P, d = f.pd(uprev, p, t) + else # f isa ConservativePDSFunction + P = f.p(uprev, p, t) + end integrator.stats.nf += 1 + Ptmp = a21 * P # avoid division by zero due to zero Patankar weights σ = add_small_constant(uprev, small_constant) # build linear system matrix and rhs if f isa PDSFunction - d = f.d(uprev, p, t) # evaluate nonconservative destruction terms dtmp = a21 * d rhs = uprev + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) linprob = LinearProblem(M, rhs) - else - # f isa ConservativePDSFunction + else # f isa ConservativePDSFunction M = build_mprk_matrix(Ptmp, σ, dt) linprob = LinearProblem(M, uprev) end @@ -590,19 +595,23 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - P2 = f.p(u, p, t + a21 * dt) - Ptmp = b1 * P + b2 * P2 + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P2, d2 = f.pd(u, p, t + a21 * dt) + else # f isa ConservativePDSFunction + P2 = f.p(u, p, t + a21 * dt) + end integrator.stats.nf += 1 + Ptmp = b1 * P + b2 * P2 # build linear system matrix and rhs if f isa PDSFunction - d2 = f.d(u, p, t + a21 * dt) # evaluate nonconservative destruction terms dtmp = b1 * d + b2 * d2 rhs = uprev + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) linprob = LinearProblem(M, rhs) - else - # f isa ConservativePDSFunction + else # f isa ConservativePDSFunction M = build_mprk_matrix(Ptmp, σ, dt) linprob = LinearProblem(M, uprev) end @@ -697,8 +706,7 @@ end # We use P2 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system - f.p(P, uprev, p, t) # evaluate production terms - f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + f.pd(P, D, uprev, p, t) # evaluate production and destruction terms integrator.stats.nf += 1 if issparse(P) # We need to keep the structural nonzeros of the production terms. @@ -726,9 +734,9 @@ end # Same as linres = P2 \ tmp linsolve.A = P2 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 if isone(a21) σ .= u @@ -737,8 +745,7 @@ end end @.. broadcast=false σ=σ + small_constant - f.p(P2, u, p, t + a21 * dt) # evaluate production terms - f.d(D2, u, p, t + a21 * dt) # evaluate nonconservative destruction terms + f.pd(P2, D2, u, p, t + a21 * dt) # evaluate production and destruction terms integrator.stats.nf += 1 if issparse(P) @@ -764,9 +771,9 @@ end # Same as linres = P2 \ tmp linsolve.A = P2 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 # Now σ stores the error estimate # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems. @@ -812,9 +819,9 @@ end # Same as linres = P2 \ tmp linsolve.A = P2 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 if isone(a21) σ .= u @@ -842,9 +849,9 @@ end # Same as linres = P2 \ tmp linsolve.A = P2 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 # Now σ stores the error estimate # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems. @@ -987,7 +994,7 @@ You can optionally choose the linear solver to be used by passing an algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) as keyword argument `linsolve`. You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators -to avoid divisions by zero. To display the default value for data type `type` evaluate +to avoid divisions by zero. To display the default value for data type `type` evaluate `MPRK43II(gamma).small_constant_function(type)`, where `type` can be, e.g., `Float64`. @@ -1094,10 +1101,15 @@ end f = integrator.f - # evaluate production matrix - P = f.p(uprev, p, t) - Ptmp = a21 * P + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P, d = f.pd(uprev, p, t) + else # f isa ConservativePDSFunction + P = f.p(uprev, p, t) + end integrator.stats.nf += 1 + Ptmp = a21 * P # avoid division by zero due to zero Patankar weights σ = add_small_constant(uprev, small_constant) @@ -1107,7 +1119,6 @@ end # build linear system matrix and rhs if f isa PDSFunction - d = f.d(uprev, p, t) dtmp = a21 * d rhs = uprev + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) @@ -1129,13 +1140,19 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - P2 = f.p(u, p, t + c2 * dt) + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P2, d2 = f.pd(u, p, t + c2 * dt) + else # f isa ConservativePDSFunction + P2 = f.p(u, p, t + c2 * dt) + end + integrator.stats.nf += 1 Ptmp = a31 * P + a32 * P2 integrator.stats.nf += 1 # build linear system matrix and rhs if f isa PDSFunction - d2 = f.d(u, p, t + c2 * dt) # evaluate nonconservative destruction terms dtmp = a31 * d + a32 * d2 rhs = uprev + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) @@ -1179,13 +1196,18 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - P3 = f.p(u, p, t + c3 * dt) - Ptmp = b1 * P + b2 * P2 + b3 * P3 + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P3, d3 = f.pd(u, p, t + c3 * dt) + else # f isa ConservativePDSFunction + P3 = f.p(u, p, t + c3 * dt) + end integrator.stats.nf += 1 + Ptmp = b1 * P + b2 * P2 + b3 * P3 # build linear system matrix if f isa PDSFunction - d3 = f.d(u, p, t + c3 * dt) # evaluate nonconservative destruction terms dtmp = b1 * d + b2 * d2 + b3 * d3 rhs = uprev + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) @@ -1291,8 +1313,7 @@ end # We use P3 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system - f.p(P, uprev, p, t) # evaluate production terms - f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + f.pd(P, D, uprev, p, t) # evaluate production and destruction terms if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1320,18 +1341,17 @@ end # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres if !(q1 ≈ q2) tmp2 .= u #u2 in out-of-place version end - integrator.stats.nsolve += 1 @.. broadcast=false σ=σ^(1 - q1) * u^q1 @.. broadcast=false σ=σ + small_constant - f.p(P2, u, p, t + c2 * dt) # evaluate production terms - f.d(D2, u, p, t + c2 * dt) # evaluate nonconservative destruction terms + f.pd(P2, D2, u, p, t + c2 * dt) # evaluate production and destruction terms if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1357,9 +1377,9 @@ end # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 if !(q1 ≈ q2) @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 @@ -1396,8 +1416,7 @@ end # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=σ + small_constant - f.p(P3, u, p, t + c3 * dt) # evaluate production terms - f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms + f.pd(P3, D3, u, p, t + c3 * dt) # evaluate production and destruction terms if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1423,9 +1442,9 @@ end # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 # Now tmp stores the error estimate @.. broadcast=false tmp=u - σ @@ -1469,12 +1488,12 @@ end # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres if !(q1 ≈ q2) tmp2 .= u #u2 in out-of-place version end - integrator.stats.nsolve += 1 @.. broadcast=false σ=σ^(1 - q1) * u^q1 @.. broadcast=false σ=σ + small_constant @@ -1498,9 +1517,9 @@ end # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 if !(q1 ≈ q2) @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 @@ -1549,9 +1568,9 @@ end # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) + integrator.stats.nsolve += 1 u .= linres - integrator.stats.nsolve += 1 # Now tmp stores the error estimate @.. broadcast=false tmp=u - σ diff --git a/src/proddest.jl b/src/proddest.jl index 646f466f..6ad6bc23 100644 --- a/src/proddest.jl +++ b/src/proddest.jl @@ -2,35 +2,42 @@ abstract type AbstractPDSProblem end """ - PDSProblem(P, D, u0, tspan, p = NullParameters(); + PDSProblem(PD, u0, tspan, p = NullParameters(); p_prototype = nothing, analytic = nothing, std_rhs = nothing) A structure describing a system of ordinary differential equations in form of a production-destruction system (PDS). -`P` denotes the function defining the production matrix ``P``. +`PD` denotes the function defining +- the production matrix ``P`` +- the destruction vector ``D`` + +The off-diagonal elements of ``P`` contain production terms with destruction counterparts +(that do not need to be stored explicitly but are extracted from ``P``). The diagonal of ``P`` contains production terms without destruction counterparts. -`D` is the function defining the vector of destruction terms ``D`` without production counterparts. +The vector of destruction terms ``D`` contains destruction terms without production counterparts. `u0` is the vector of initial conditions and `tspan` the time span `(t_initial, t_final)` of the problem. The optional argument `p` can be used -to pass additional parameters to the functions `P` and `D`. +to pass additional parameters to the function `PD`. -The functions `P` and `D` can be used either in the out-of-place form with signature -`production_terms = P(u, p, t)` or the in-place form `P(production_terms, u, p, t)`. +The function `PD` can be used either in the out-of-place form with signature +`production_matrix, destruction_vector = PD(u, p, t)` or the in-place form +`PD(production_matrix, destruction_vector, u, p, t)`. ### Keyword arguments: ### -- `p_prototype`: If `P` is given in in-place form, `p_prototype` or copies thereof are used to store evaluations of `P`. - If `p_prototype` is not specified explicitly and `P` is in-place, then `p_prototype` will be internally - set to `zeros(eltype(u0), (length(u0), length(u0)))`. +- `p_prototype`: If `PD` is given in in-place form, `p_prototype` or copies thereof are used + to store evaluations of `P`. + If `p_prototype` is not specified explicitly and `PD` is in-place, then `p_prototype` + will be internally set to `zeros(eltype(u0), (length(u0), length(u0)))`. - `analytic`: The analytic solution of a PDS must be given in the form `f(u0,p,t)`. Specifying the analytic solution can be useful for plotting and convergence tests. - `std_rhs`: The standard ODE right-hand side evaluation function callable as `du = std_rhs(u, p, t)` for the out-of-place form and as `std_rhs(du, u, p, t)` for the in-place form. Solvers that do not rely on - the production-destruction representation of the ODE, will use this function - instead to compute the solution. If not specified, - a default implementation calling `P` and `D` is used. + the production-destruction representation of the ODE will use this function + instead to compute the solution (e.g., standard ODE methods implemented in + OrdinaryDiffEq.jl). If not specified, a default implementation calling `PD` is used. ## References @@ -43,10 +50,9 @@ The functions `P` and `D` can be used either in the out-of-place form with signa struct PDSProblem{iip} <: AbstractPDSProblem end # New ODE function PDSFunction -struct PDSFunction{iip, specialize, P, D, PrototypeP, PrototypeD, StdRHS, Ta} <: +struct PDSFunction{iip, specialize, PD, PrototypeP, PrototypeD, StdRHS, Ta} <: AbstractODEFunction{iip} - p::P - d::D + pd::PD p_prototype::PrototypeP d_prototype::PrototypeD std_rhs::StdRHS @@ -73,21 +79,15 @@ function Base.getproperty(obj::PDSFunction, sym::Symbol) end # Most general constructor for PDSProblems -function PDSProblem(P, D, u0, tspan, p = NullParameters(); +function PDSProblem(PD, u0, tspan, p = NullParameters(); kwargs...) - Piip = isinplace(P, 4) - Diip = isinplace(D, 4) - if Piip == Diip - iip = Piip - else - error("Conflict due to the joint use of in-place and out-of-place functions.") - end - return PDSProblem{iip}(P, D, u0, tspan, p; kwargs...) + iip = isinplace(PD, 5, "PD"; outofplace_param_number = 3) + return PDSProblem{iip}(PD, u0, tspan, p; kwargs...) end # Specialized constructor for PDSProblems setting `iip` manually # (arbitrary functions) -function PDSProblem{iip}(P, D, u0, tspan, p = NullParameters(); +function PDSProblem{iip}(PD, u0, tspan, p = NullParameters(); p_prototype = nothing, analytic = nothing, std_rhs = nothing, @@ -101,98 +101,95 @@ function PDSProblem{iip}(P, D, u0, tspan, p = NullParameters(); # evaluations of D. d_prototype = similar(u0 ./ oneunit(first(tspan))) - PD = PDSFunction{iip}(P, D; p_prototype, d_prototype, - analytic, std_rhs) - PDSProblem{iip}(PD, u0, tspan, p; kwargs...) + pds = PDSFunction{iip}(PD; p_prototype, d_prototype, analytic, std_rhs) + PDSProblem{iip}(pds, u0, tspan, p; kwargs...) end # Specialized constructor for PDSProblems setting `iip` manually # (PDSFunction) -function PDSProblem{iip}(PD::PDSFunction, u0, tspan, p = NullParameters(); +function PDSProblem{iip}(pds::PDSFunction, u0, tspan, p = NullParameters(); kwargs...) where {iip} - ODEProblem(PD, u0, tspan, p, PDSProblem{iip}(); kwargs...) + ODEProblem(pds, u0, tspan, p, PDSProblem{iip}(); kwargs...) end # Specialized constructor for PDSFunction setting `iip` manually -function PDSFunction{iip}(P, D; kwargs...) where {iip} - PDSFunction{iip, FullSpecialize}(P, D; kwargs...) +function PDSFunction{iip}(PD; kwargs...) where {iip} + PDSFunction{iip, FullSpecialize}(PD; kwargs...) end # Most specific constructor for PDSFunction -function PDSFunction{iip, FullSpecialize}(P, D; +function PDSFunction{iip, FullSpecialize}(PD; p_prototype = nothing, d_prototype = nothing, analytic = nothing, std_rhs = nothing) where {iip} if std_rhs === nothing - std_rhs = PDSStdRHS(P, D, p_prototype, d_prototype) + std_rhs = PDSStdRHS(PD, p_prototype, d_prototype) end - PDSFunction{iip, FullSpecialize, typeof(P), typeof(D), typeof(p_prototype), + PDSFunction{iip, FullSpecialize, typeof(PD), typeof(p_prototype), typeof(d_prototype), - typeof(std_rhs), typeof(analytic)}(P, D, p_prototype, d_prototype, std_rhs, + typeof(std_rhs), typeof(analytic)}(PD, p_prototype, d_prototype, std_rhs, analytic) end # Evaluation of a PDSFunction -function (PD::PDSFunction)(u, p, t) - return PD.std_rhs(u, p, t) +function (pds::PDSFunction)(u, p, t) + return pds.std_rhs(u, p, t) end -function (PD::PDSFunction)(du, u, p, t) - return PD.std_rhs(du, u, p, t) +function (pds::PDSFunction)(du, u, p, t) + return pds.std_rhs(du, u, p, t) end # Default implementation of the standard right-hand side evaluation function -struct PDSStdRHS{P, D, PrototypeP, PrototypeD, TMP} <: Function - p::P - d::D +struct PDSStdRHS{PD, PrototypeP, PrototypeD, TMP} <: Function + pd::PD p_prototype::PrototypeP d_prototype::PrototypeD tmp::TMP end -function PDSStdRHS(P, D, p_prototype, d_prototype) +function PDSStdRHS(PD, p_prototype, d_prototype) if p_prototype isa AbstractSparseMatrix tmp = zeros(eltype(p_prototype), (size(p_prototype, 1),)) / oneunit(first(p_prototype)) # drop units else tmp = nothing end - PDSStdRHS(P, D, p_prototype, d_prototype, tmp) + PDSStdRHS(PD, p_prototype, d_prototype, tmp) end # Evaluation of a PDSStdRHS (out-of-place) function (PD::PDSStdRHS)(u, p, t) - P = PD.p(u, p, t) - D = PD.d(u, p, t) - diag(P) + vec(sum(P, dims = 2)) - - vec(sum(P, dims = 1)) - vec(D) + P, D = PD.pd(u, p, t) + rhs = diag(P) + vec(sum(P, dims = 2)) - vec(sum(P, dims = 1)) - vec(D) + return rhs end # Evaluation of a PDSStdRHS (in-place) function (PD::PDSStdRHS)(du, u, p, t) - PD.p(PD.p_prototype, u, p, t) + PD.pd(PD.p_prototype, PD.d_prototype, u, p, t) if PD.p_prototype isa AbstractSparseMatrix - # row sum coded as matrix-vector product + # row sum coded as matrix-vector product fill!(PD.tmp, one(eltype(PD.tmp))) mul!(vec(du), PD.p_prototype, PD.tmp) - for i in 1:length(u) #vec(du) .+= diag(PD.p_prototype) + vec(du) .-= PD.d_prototype + + for i in eachindex(u) #vec(du) .+= diag(PD.p_prototype) du[i] += PD.p_prototype[i, i] end sum!(PD.d_prototype', PD.p_prototype) vec(du) .-= PD.d_prototype - PD.d(PD.d_prototype, u, p, t) - vec(du) .-= PD.d_prototype else - PD.d(PD.d_prototype, u, p, t) # This implementation does not need any auxiliary vectors - for i in 1:length(u) - du[i] = PD.p_prototype[i, i] - PD.d_prototype[i] - for j in 1:length(u) - du[i] += PD.p_prototype[i, j] - PD.p_prototype[j, i] + for i in eachindex(u) + du_i = PD.p_prototype[i, i] - PD.d_prototype[i] + for j in eachindex(u) + du_i += PD.p_prototype[i, j] - PD.p_prototype[j, i] end + du[i] = du_i end end return nothing @@ -205,29 +202,33 @@ end analytic = nothing, std_rhs = nothing) -A structure describing a conservative system of ordinary differential equation in form of a production-destruction system (PDS). +A structure describing a conservative system of ordinary differential equation in form of +a production-destruction system (PDS). `P` denotes the function defining the production matrix ``P``. +The off-diagonal elements of ``P`` contain production terms with destruction counterparts +(that do not need to be stored explicitly but are extracted from ``P``). The diagonal of ``P`` contains production terms without destruction counterparts. `u0` is the vector of initial conditions and `tspan` the time span `(t_initial, t_final)` of the problem. The optional argument `p` can be used to pass additional parameters to the function `P`. The function `P` can be given either in the out-of-place form with signature -`production_terms = P(u, p, t)` or the in-place form `P(production_terms, u, p, t)`. +`production_matrix = P(u, p, t)` or the in-place form `P(production_matrix, u, p, t)`. ### Keyword arguments: ### -- `p_prototype`: If `P` is given in in-place form, `p_prototype` or copies thereof are used to store evaluations of `P`. - If `p_prototype` is not specified explicitly and `P` is in-place, then `p_prototype` will be internally - set to `zeros(eltype(u0), (length(u0), length(u0)))`. +- `p_prototype`: If `P` is given in in-place form, `p_prototype` or copies thereof + are used to store evaluations of `P`. + If `p_prototype` is not specified explicitly and `P` is in-place, then `p_prototype` + will be internally set to `zeros(eltype(u0), (length(u0), length(u0)))`. - `analytic`: The analytic solution of a PDS must be given in the form `f(u0,p,t)`. Specifying the analytic solution can be useful for plotting and convergence tests. - `std_rhs`: The standard ODE right-hand side evaluation function callable as `du = std_rhs(u, p, t)` for the out-of-place form and as `std_rhs(du, u, p, t)` for the in-place form. Solvers that do not rely on the production-destruction representation of the ODE, will use this function - instead to compute the solution. If not specified, - a default implementation calling `P` is used. + instead to compute the solution (e.g., standard ODE methods implemented in + OrdinaryDiffEq.jl). If not specified, a default implementation calling `P` is used. ## References diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 447eb0d4..4b84c1b3 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -126,22 +126,25 @@ end f = integrator.f - # evaluate production matrix - P = f.p(uprev, p, t) - Ptmp = b10 * P + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P, d = f.pd(uprev, p, t) + else # f isa ConservativePDSFunction + P = f.p(uprev, p, t) + end integrator.stats.nf += 1 + Ptmp = b10 * P # avoid division by zero due to zero Patankar weights σ = add_small_constant(uprev, small_constant) # build linear system matrix and rhs if f isa PDSFunction - d = f.d(uprev, p, t) # evaluate nonconservative destruction terms dtmp = b10 * d rhs = a10 * uprev + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) - else - # f isa ConservativePDSFunction + else # f isa ConservativePDSFunction M = build_mprk_matrix(Ptmp, σ, dt) rhs = a10 * uprev end @@ -161,18 +164,22 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - P2 = f.p(u, p, t + b10 * dt) - Ptmp = b20 * P + b21 * P2 + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P2, d2 = f.pd(u, p, t + b10 * dt) + else # f isa ConservativePDSFunction + P2 = f.p(u, p, t + b10 * dt) + end integrator.stats.nf += 1 + Ptmp = b20 * P + b21 * P2 # build linear system matrix and rhs if f isa PDSFunction - d2 = f.d(u, p, t + b10 * dt) # evaluate nonconservative destruction terms dtmp = b20 * d + b21 * d2 rhs = a20 * uprev + a21 * u + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) - else - # f isa ConservativePDSFunction + else # f isa ConservativePDSFunction M = build_mprk_matrix(Ptmp, σ, dt) rhs = a20 * uprev + a21 * u end @@ -273,8 +280,7 @@ end # We use P2 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system - f.p(P, uprev, p, t) # evaluate production terms - f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + f.pd(P, D, uprev, p, t) # evaluate production and destruction terms integrator.stats.nf += 1 if issparse(P) # We need to keep the structural nonzeros of the production terms. @@ -313,8 +319,7 @@ end end @.. broadcast=false σ=σ + small_constant - f.p(P2, u, p, t + b10 * dt) # evaluate production terms - f.d(D2, u, p, t + b10 * dt) # evaluate nonconservative destruction terms + f.pd(P2, D2, u, p, t + b10 * dt) # evaluate production and destruction terms integrator.stats.nf += 1 if issparse(P) @@ -466,7 +471,7 @@ You can optionally choose the linear solver to be used by passing an algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) as keyword argument `linsolve`. You can also choose the parameter `small_constant` which is added to all Patankar-weight denominators -to avoid divisions by zero. To display the default value for data type `type` evaluate +to avoid divisions by zero. To display the default value for data type `type` evaluate `SSPMPRK43. small_constant_function(type)`, where `type` can be, e.g., `Float64`. @@ -607,21 +612,25 @@ end f = integrator.f - # evaluate production matrix - P = f.p(uprev, p, t) - Ptmp = β10 * P + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P, d = f.pd(uprev, p, t) + else # f isa ConservativePDSFunction + P = f.p(uprev, p, t) + end integrator.stats.nf += 1 + Ptmp = β10 * P # avoid division by zero due to zero Patankar weights σ = add_small_constant(uprev, small_constant) # build linear system matrix and rhs if f isa PDSFunction - d = f.d(uprev, p, t) dtmp = β10 * d rhs = α10 * uprev + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) - else + else # f isa ConservativePDSFunction rhs = α10 * uprev M = build_mprk_matrix(Ptmp, σ, dt) end @@ -638,18 +647,22 @@ end # avoid division by zero due to zero Patankar weights ρ = add_small_constant(ρ, small_constant) - P2 = f.p(u, p, t + β10 * dt) - Ptmp = β20 * P + β21 * P2 + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P2, d2 = f.pd(u, p, t + β10 * dt) + else # f isa ConservativePDSFunction + P2 = f.p(u, p, t + β10 * dt) + end integrator.stats.nf += 1 + Ptmp = β20 * P + β21 * P2 # build linear system matrix and rhs if f isa PDSFunction - d2 = f.d(u, p, t + β10 * dt) # evaluate nonconservative destruction terms dtmp = β20 * d + β21 * d2 rhs = α20 * uprev + α21 * u2 + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, ρ, dt, dtmp) - - else + else # f isa ConservativePDSFunction rhs = α20 * uprev + α21 * u2 M = build_mprk_matrix(Ptmp, ρ, dt) end @@ -691,17 +704,22 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - P3 = f.p(u, p, t + c3 * dt) - Ptmp = β30 * P + β31 * P2 + β32 * P3 + # evaluate production matrix and nonconservative destruction terms + # (if present) + if f isa PDSFunction + P3, d3 = f.pd(u, p, t + c3 * dt) + else # f isa ConservativePDSFunction + P3 = f.p(u, p, t + c3 * dt) + end integrator.stats.nf += 1 + Ptmp = β30 * P + β31 * P2 + β32 * P3 # build linear system matrix if f isa PDSFunction - d3 = f.d(u, p, t + c3 * dt) # evaluate nonconservative destruction terms dtmp = β30 * d + β31 * d2 + β32 * d3 rhs = α30 * uprev + α31 * u2 + α32 * u + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) - else + else # f isa ConservativePDSFunction rhs = α30 * uprev + α31 * u2 + α32 * u M = build_mprk_matrix(Ptmp, σ, dt) end @@ -804,8 +822,7 @@ end # We use P3 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system - f.p(P, uprev, p, t) # evaluate production terms - f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + f.pd(P, D, uprev, p, t) # evaluate production and destruction terms if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -841,8 +858,7 @@ end @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ @.. broadcast=false ρ=ρ + small_constant - f.p(P2, u, p, t + β10 * dt) # evaluate production terms - f.d(D2, u, p, t + β10 * dt) # evaluate nonconservative destruction terms + f.pd(P2, D2, u, p, t + β10 * dt) # evaluate production and destruction terms if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -908,8 +924,7 @@ end # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=σ + small_constant - f.p(P3, u, p, t + c3 * dt) # evaluate production terms - f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms + f.pd(P3, D3, u, p, t + c3 * dt) # evaluate production and destruction terms if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see diff --git a/test/runtests.jl b/test/runtests.jl index 3c0b1829..92c9b798 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -142,26 +142,28 @@ end # u₁' = -3 u₁ + 0.5 u₂ - 2 u₁ + 0.5 u₂ (= -5 u₁ + u₂) # u₂' = 3 u₁ - 0.5 u₂ - 0.5 u₂ + 2 u₁ (= 5 u₁ - u₂) # linear model problem - nonconservative - out-of-place -linmodP(u, p, t) = [0.5*u[2] 0.5*u[2]; 3*u[1] 2*u[1]] -linmodD(u, p, t) = [2 * u[1]; 0.5 * u[2]] -const prob_pds_linmod_nonconservative = PDSProblem(linmodP, linmodD, [0.9, 0.1], (0.0, 2.0), +function linmodPD(u, p, t) + P = [0.5*u[2] 0.5*u[2]; 3*u[1] 2*u[1]] + d = [2 * u[1]; 0.5 * u[2]] + return P, d +end +const prob_pds_linmod_nonconservative = PDSProblem(linmodPD, [0.9, 0.1], (0.0, 2.0), [5.0, 1.0]; analytic = f_analytic) # linear model problem - nonconservative - in-place -function linmodP!(P, u, p, t) +function linmodPD!(P, D, u, p, t) P[1, 1] = 0.5 * u[2] P[1, 2] = 0.5 * u[2] P[2, 1] = 3 * u[1] P[2, 2] = 2 * u[1] - return nothing -end -function linmodD!(D, u, p, t) + D[1] = 2 * u[1] D[2] = 0.5 * u[2] + return nothing end -const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [0.9, 0.1], +const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodPD!, [0.9, 0.1], (0.0, 2.0), [5.0, 1.0]; analytic = f_analytic) @@ -200,8 +202,9 @@ function linear_advection_fd_upwind_P!(P::SparseMatrixCSC, u, p, t) end return nothing end -function linear_advection_fd_upwind_D!(D, u, p, t) - D .= 0.0 +function linear_advection_fd_upwind_PD!(P, D, u, p, t) + linear_advection_fd_upwind_P!(P, u, p, t) + fill!(D, 0.0) return nothing end @@ -219,15 +222,11 @@ end end @testset "ODE RHS" begin - let counter_p = Ref(1), counter_d = Ref(1), counter_rhs = Ref(1) + let counter_pd = Ref(1), counter_rhs = Ref(1) # out-of-place - prod1 = (u, p, t) -> begin - counter_p[] += 1 - return [0 u[2]; u[1] 0] - end - dest1 = (u, p, t) -> begin - counter_d[] += 1 - return zero(u) + prod_dest_1 = (u, p, t) -> begin + counter_pd[] += 1 + return [0 u[2]; u[1] 0], zero(u) end rhs1 = (u, p, t) -> begin counter_rhs[] += 1 @@ -235,36 +234,28 @@ end end u0 = [1.0, 0.0] tspan = (0.0, 1.0) - prob_default = PDSProblem(prod1, dest1, u0, tspan) - prob_special = PDSProblem(prod1, dest1, u0, tspan; std_rhs = rhs1) + prob_default = PDSProblem(prod_dest_1, u0, tspan) + prob_special = PDSProblem(prod_dest_1, u0, tspan; std_rhs = rhs1) - counter_p[] = 0 - counter_d[] = 0 + counter_pd[] = 0 counter_rhs[] = 0 @inferred prob_default.f(u0, nothing, 0.0) - @test counter_p[] == 1 - @test counter_d[] == 1 + @test counter_pd[] == 1 @test counter_rhs[] == 0 - counter_p[] = 0 - counter_d[] = 0 + counter_pd[] = 0 counter_rhs[] = 0 @inferred prob_special.f(u0, nothing, 0.0) - @test counter_p[] == 0 - @test counter_d[] == 0 + @test counter_pd[] == 0 @test counter_rhs[] == 1 # in-place - prod1! = (P, u, p, t) -> begin - counter_p[] += 1 + prod_dest_1! = (P, D, u, p, t) -> begin + counter_pd[] += 1 P[1, 1] = 0 P[1, 2] = u[2] P[2, 1] = u[1] P[2, 2] = 0 - return nothing - end - dest1! = (D, u, p, t) -> begin - counter_d[] += 1 fill!(D, 0) return nothing end @@ -276,50 +267,38 @@ end end u0 = [1.0, 0.0] tspan = (0.0, 1.0) - prob_default = PDSProblem(prod1!, dest1!, u0, tspan) - prob_special = PDSProblem(prod1!, dest1!, u0, tspan; std_rhs = rhs1!) + prob_default = PDSProblem(prod_dest_1!, u0, tspan) + prob_special = PDSProblem(prod_dest_1!, u0, tspan; std_rhs = rhs1!) du = similar(u0) - counter_p[] = 0 - counter_d[] = 0 + counter_pd[] = 0 counter_rhs[] = 0 @inferred prob_default.f(du, u0, nothing, 0.0) - @test counter_p[] == 1 - @test counter_d[] == 1 + @test counter_pd[] == 1 @test counter_rhs[] == 0 - counter_p[] = 0 - counter_d[] = 0 + counter_pd[] = 0 counter_rhs[] = 0 @inferred prob_special.f(du, u0, nothing, 0.0) - @test counter_p[] == 0 - @test counter_d[] == 0 + @test counter_pd[] == 0 @test counter_rhs[] == 1 - counter_p[] = 0 - counter_d[] = 0 + counter_pd[] = 0 counter_rhs[] = 0 @inferred solve(prob_default, MPE(); dt = 0.1) - @test 10 <= counter_p[] <= 11 - @test 10 <= counter_d[] <= 11 - @test counter_d[] == counter_p[] + @test 10 <= counter_pd[] <= 11 @test counter_rhs[] == 0 - counter_p[] = 0 - counter_d[] = 0 + counter_pd[] = 0 counter_rhs[] = 0 @inferred solve(prob_default, Euler(); dt = 0.1) - @test 10 <= counter_p[] <= 11 - @test 10 <= counter_d[] <= 11 - @test counter_d[] == counter_p[] + @test 10 <= counter_pd[] <= 11 @test counter_rhs[] == 0 - counter_p[] = 0 - counter_d[] = 0 + counter_pd[] = 0 counter_rhs[] = 0 @inferred solve(prob_special, Euler(); dt = 0.1) - @test counter_p[] == 0 - @test counter_d[] == 0 + @test counter_pd[] == 0 @test 10 <= counter_rhs[] <= 11 end end @@ -412,22 +391,18 @@ end linmod_ODE_ip = ODEProblem(linmod!, u0, tspan) # out-of-place syntax for PDS - linmodP(u, p, t) = [0 u[2]; 5*u[1] 0] - linmodD(u, p, t) = [0.0; 0.0] - linmod_PDS_op = PDSProblem(linmodP, linmodD, u0, tspan) + linmodPD(u, p, t) = [0 u[2]; 5*u[1] 0], [0.0; 0.0] + linmod_PDS_op = PDSProblem(linmodPD, u0, tspan) # in-place sytanx for PDS - function linmodP!(P, u, p, t) + function linmodPD!(P, u, p, t) fill!(P, zero(eltype(P))) P[1, 2] = u[2] P[2, 1] = 5 * u[1] - return nothing - end - function linmodD!(D, u, p, t) fill!(D, zero(eltype(D))) return nothing end - linmod_PDS_ip = PDSProblem(linmodP!, linmodD!, u0, tspan) + linmod_PDS_ip = PDSProblem(linmodPD!, u0, tspan) # solutions sol_linmod_ODE_op = solve(linmod_ODE_op, Tsit5()) @@ -449,14 +424,6 @@ end alloc2 = @allocated(solve(linmod_PDS_ip, Tsit5())) @test 0.95 < alloc1 / alloc2 < 1.05 end - - @testset "PDSProblem error handling" begin - P(u, p, t) = 0.0 - D(du, u, p, t) = 0.0 - if VERSION >= v"1.8" - @test_throws "in-place and out-of-place" PDSProblem(P, D, 0.0, (0.0, 1.0)) - end - end end # Here we check that solutions of equivalent ODEProblems, PDSProblems or @@ -480,9 +447,9 @@ end linmod_f_ip = ODEProblem(linmod!, u0, tspan) # out-of-place syntax for PDS linmodP(u, p, t) = [0.0 u[2]; 5.0*u[1] 0.0] - linmodD(u, p, t) = [0.0; 0.0] - linmod_PDS_op = PDSProblem(linmodP, linmodD, u0, tspan) - linmod_PDS_op_2 = PDSProblem{false}(linmodP, linmodD, u0, tspan) + linmodPD(u, p, t) = linmodP(u, p, t), [0.0; 0.0] + linmod_PDS_op = PDSProblem(linmodPD, u0, tspan) + linmod_PDS_op_2 = PDSProblem{false}(linmodPD, u0, tspan) linmod_ConsPDS_op = ConservativePDSProblem(linmodP, u0, tspan) linmod_ConsPDS_op_2 = ConservativePDSProblem{false}(linmodP, u0, tspan) # in-place sytanx for PDS @@ -492,12 +459,13 @@ end P[2, 1] = 5.0 * u[1] return nothing end - function linmodD!(D, u, p, t) + function linmodPD!(P, D, u, p, t) + linmodP!(P, u, p, t) D .= 0.0 return nothing end - linmod_PDS_ip = PDSProblem(linmodP!, linmodD!, u0, tspan) - linmod_PDS_ip_2 = PDSProblem{true}(linmodP!, linmodD!, u0, tspan) + linmod_PDS_ip = PDSProblem(linmodPD!, u0, tspan) + linmod_PDS_ip_2 = PDSProblem{true}(linmodPD!, u0, tspan) linmod_ConsPDS_ip = ConservativePDSProblem(linmodP!, u0, tspan) linmod_ConsPDS_ip_2 = ConservativePDSProblem{true}(linmodP!, u0, tspan) @@ -552,24 +520,20 @@ end end lotvol_f_ip = ODEProblem(lotvol!, u0, tspan) # out-of-place syntax for PDS - lotvolP(u, p, t) = [u[1] 0.0; u[1]*u[2] 0.0] - lotvolD(u, p, t) = [0.0; u[2]] - lotvol_PDS_op = PDSProblem(lotvolP, lotvolD, u0, tspan) - lotvol_PDS_op_2 = PDSProblem{false}(lotvolP, lotvolD, u0, tspan) + lotvolPD(u, p, t) = [u[1] 0.0; u[1]*u[2] 0.0], [0.0; u[2]] + lotvol_PDS_op = PDSProblem(lotvolPD, u0, tspan) + lotvol_PDS_op_2 = PDSProblem{false}(lotvolPD, u0, tspan) # in-place sytanx for PDS - function lotvolP!(P, u, p, t) + function lotvolPD!(P, D, u, p, t) P .= 0.0 P[1, 1] = u[1] P[2, 1] = u[2] * u[1] - return nothing - end - function lotvolD!(D, u, p, t) D .= 0.0 D[2] = u[2] return nothing end - lotvol_PDS_ip = PDSProblem(lotvolP!, lotvolD!, u0, tspan) - lotvol_PDS_ip_2 = PDSProblem{true}(lotvolP!, lotvolD!, u0, tspan) + lotvol_PDS_ip = PDSProblem(lotvolPD!, u0, tspan) + lotvol_PDS_ip_2 = PDSProblem{true}(lotvolPD!, u0, tspan) # solutions sol_lotvol_f_op = solve(lotvol_f_op, Tsit5()) @@ -603,18 +567,15 @@ end linear_advection_fd_upwind_ODE = ODEProblem(linear_advection_fd_upwind_f!, u0, tspan) # problem with dense matrices - linear_advection_fd_upwind_PDS_dense = PDSProblem(linear_advection_fd_upwind_P!, - linear_advection_fd_upwind_D!, + linear_advection_fd_upwind_PDS_dense = PDSProblem(linear_advection_fd_upwind_PD!, u0, tspan) # problem with sparse matrices p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1), N - 1 => ones(eltype(u0), 1)) - linear_advection_fd_upwind_PDS_sparse = PDSProblem(linear_advection_fd_upwind_P!, - linear_advection_fd_upwind_D!, + linear_advection_fd_upwind_PDS_sparse = PDSProblem(linear_advection_fd_upwind_PD!, u0, tspan; p_prototype = p_prototype) - linear_advection_fd_upwind_PDS_sparse_2 = PDSProblem{true}(linear_advection_fd_upwind_P!, - linear_advection_fd_upwind_D!, + linear_advection_fd_upwind_PDS_sparse_2 = PDSProblem{true}(linear_advection_fd_upwind_PD!, u0, tspan; p_prototype = p_prototype) linear_advection_fd_upwind_ConsPDS_sparse = ConservativePDSProblem(linear_advection_fd_upwind_P!, @@ -681,14 +642,16 @@ end f(u, p, t) = [u[2] - 5 * u[1]; -u[2] + 5 * u[1]] / u"s" P(u, p, t) = [0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"] D(u, p, t) = [0.0; 0.0]u"N/s" + PD(u, p, t) = P(u, p, t), D(u, p, t) g = PositiveIntegrators.ConservativePDSFunction{false}(P) - h = PositiveIntegrators.PDSFunction{false}(P, D) + h = PositiveIntegrators.PDSFunction{false}(PD) f_static(u, p, t) = SA[(u[2] - 5 * u[1]) / u"s"; (-u[2] + 5 * u[1]) / u"s"] P_static(u, p, t) = SA[0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"] D_static(u, p, t) = SA[0.0u"N/s"; 0.0u"N/s"] + PD_static(u, p, t) = P_static(u, p, t), D_static(u, p, t) g_static = PositiveIntegrators.ConservativePDSFunction{false}(P_static) - h_static = PositiveIntegrators.PDSFunction{false}(P_static, D_static) + h_static = PositiveIntegrators.PDSFunction{false}(PD_static) f1 = f(u0, nothing, t0) g1 = g(u0, nothing, t0) @@ -721,7 +684,8 @@ end end return nothing end - function D!(D, u, p, t) + function PD!(P, D, u, p, t) + P!(P, u, p, t) fill!(D, zero(eltype(D))) return nothing end @@ -738,13 +702,13 @@ end p_prototype = P_tridiagonal) f_sparse! = PositiveIntegrators.ConservativePDSFunction{true}(P!, p_prototype = P_sparse) - g_dense! = PositiveIntegrators.PDSFunction{true}(P!, D!, + g_dense! = PositiveIntegrators.PDSFunction{true}(PD!, p_prototype = P_dense, d_prototype = D_dense) - g_tridiagonal! = PositiveIntegrators.PDSFunction{true}(P!, D!, + g_tridiagonal! = PositiveIntegrators.PDSFunction{true}(PD!, p_prototype = P_tridiagonal, d_prototype = D_dense) - g_sparse! = PositiveIntegrators.PDSFunction{true}(P!, D!, + g_sparse! = PositiveIntegrators.PDSFunction{true}(PD!, p_prototype = P_sparse, d_prototype = D_dense) @@ -771,21 +735,23 @@ end f(u, p, t) = [u[2] - 5 * u[1]; -u[2] + 5 * u[1]] / u"s" P(u, p, t) = [0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"] D(u, p, t) = [0.0; 0.0]u"N/s" + PD(u, p, t) = P(u, p, t), D(u, p, t) f_static(u, p, t) = SA[(u[2] - 5 * u[1]) / u"s"; (-u[2] + 5 * u[1]) / u"s"] P_static(u, p, t) = SA[0u"N/s" u[2]/u"s"; 5 * u[1]/u"s" 0u"N/s"] D_static(u, p, t) = SA[0.0u"N/s"; 0.0u"N/s"] + PD_static(u, p, t) = P_static(u, p, t), D_static(u, p, t) prob_ode = ODEProblem(f, u0, tspan) prob_cpds = ConservativePDSProblem(P, u0, tspan) - prob_pds = PDSProblem(P, D, u0, tspan) + prob_pds = PDSProblem(PD, u0, tspan) prob_cpds_rhs = ConservativePDSProblem(P, u0, tspan; std_rhs = f) - prob_pds_rhs = PDSProblem(P, D, u0, tspan; std_rhs = f) + prob_pds_rhs = PDSProblem(PD, u0, tspan; std_rhs = f) prob_ode_static = ODEProblem(f, u0_static, tspan) prob_cpds_static = ConservativePDSProblem(P, u0_static, tspan) - prob_pds_static = PDSProblem(P, D, u0_static, tspan) + prob_pds_static = PDSProblem(PD, u0_static, tspan) prob_cpds_static_rhs = ConservativePDSProblem(P, u0_static, tspan; std_rhs = f_static) - prob_pds_static_rhs = PDSProblem(P, D, u0_static, tspan; std_rhs = f_static) + prob_pds_static_rhs = PDSProblem(PD, u0_static, tspan; std_rhs = f_static) # implicit solvers and units don't work # algs = (Euler(), ImplicitEuler(), Tsit5(), Rosenbrock23(), SDIRK2(),TRBDF2()) @@ -852,7 +818,8 @@ end end return nothing end - function D!(D, u, p, t) + function PD!(P, D, u, p, t) + P!(P, u, p, t) fill!(D, zero(eltype(D))) return nothing end @@ -881,15 +848,15 @@ end probs[6] = ConservativePDSProblem(P!, u0, tspan; p_prototype = P_sparse) probs[7] = ConservativePDSProblem(P!, u0, tspan; p_prototype = P_sparse, std_rhs = f!) - probs[8] = PDSProblem(P!, D!, u0, tspan) - probs[9] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_dense) - probs[10] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_dense, + probs[8] = PDSProblem(PD!, u0, tspan) + probs[9] = PDSProblem(PD!, u0, tspan; p_prototype = P_dense) + probs[10] = PDSProblem(PD!, u0, tspan; p_prototype = P_dense, std_rhs = f!) - probs[11] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_tridiagonal) - probs[12] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_tridiagonal, + probs[11] = PDSProblem(PD!, u0, tspan; p_prototype = P_tridiagonal) + probs[12] = PDSProblem(PD!, u0, tspan; p_prototype = P_tridiagonal, std_rhs = f!) - probs[13] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_sparse) - probs[14] = PDSProblem(P!, D!, u0, tspan; p_prototype = P_sparse, + probs[13] = PDSProblem(PD!, u0, tspan; p_prototype = P_sparse) + probs[14] = PDSProblem(PD!, u0, tspan; p_prototype = P_sparse, std_rhs = f!) for alg in algs @@ -951,7 +918,7 @@ end @testset "$alg" for prob in probs, alg in algs dt = (last(prob.tspan) - first(prob.tspan)) / 1e4 sol = solve(prob, alg; dt, isoutofdomain = isnegative) # use explicit f - sol2 = solve(PDSProblem(prob.f.p, prob.f.d, prob.u0, prob.tspan), alg; dt, + sol2 = solve(PDSProblem(prob.f.pd, prob.u0, prob.tspan), alg; dt, isoutofdomain = isnegative) # use p and d to compute f sol3 = solve(ODEProblem(prob.f.std_rhs, prob.u0, prob.tspan), alg; dt, isoutofdomain = isnegative) # use f to create ODEProblem @@ -974,7 +941,7 @@ end # Bertolazzi problem # Did not find any solver configuration to compute a reasonable solution and pass tests. # - constant time stepping requires very small dt - # - adaptive time stepping generates solutions with different number of time steps + # - adaptive time stepping generates solutions with different number of time steps # # Nevertheless, the following code shows that the same problem is solved in each case # prob = prob_pds_bertolazzi @@ -991,7 +958,7 @@ end @testset "$alg" for alg in algs dt = 1.0 sol = solve(prob, alg; dt) # use explicit f - sol2 = solve(PDSProblem(prob.f.p, prob.f.d, prob.u0, prob.tspan), alg; dt) # use p and d to compute f + sol2 = solve(PDSProblem(prob.f.pd, prob.u0, prob.tspan), alg; dt) # use p and d to compute f sol3 = solve(ODEProblem(prob.f.std_rhs, prob.u0, prob.tspan), alg; dt) # use f to create ODEProblem @test sol.t ≈ sol2.t ≈ sol3.t @test sol.u ≈ sol2.u ≈ sol3.u @@ -1078,7 +1045,7 @@ end dt = 0.1) end - # Here we check that algorithms which accept input parameters return constants + # Here we check that algorithms which accept input parameters return constants # of the same type as the inputs @testset "Constant types" begin algs = (MPRK22(0.5f0), MPRK22(1.0f0), MPRK22(2.0f0), MPRK43I(1.0f0, 0.5f0), @@ -1131,8 +1098,8 @@ end # linear model problem - out-of-place linmodP(u, p, t) = [0 p[2]*u[2]; p[1]*u[1] 0] - linmodD(u, p, t) = [0.0; 0.0] - prob_op = PDSProblem(linmodP, linmodD, u0, tspan, p; analytic = f_analytic) + linmodPD(u, p, t) = linmodP(u, p, t), [0.0; 0.0] + prob_op = PDSProblem(linmodPD, u0, tspan, p; analytic = f_analytic) prob_op_2 = ConservativePDSProblem(linmodP, u0, tspan, p; analytic = f_analytic) dt = 0.25 @@ -1150,11 +1117,12 @@ end P[2, 1] = p[1] * u[1] return nothing end - function linmodD!(D, u, p, t) + function linmodPD!(P, D, u, p, t) + linmodP!(P, u, p, t) fill!(D, zero(eltype(D))) return nothing end - prob_ip = PDSProblem(linmodP!, linmodD!, u0, tspan, p; analytic = f_analytic) + prob_ip = PDSProblem(linmodPD!, u0, tspan, p; analytic = f_analytic) prob_ip_2 = ConservativePDSProblem(linmodP!, u0, tspan, p; analytic = f_analytic) @@ -1209,11 +1177,12 @@ end P[2, 1] = p[1] * u[1] return nothing end - function linmodD!(D, u, p, t) + function linmodPD!(P, D, u, p, t) + linmodP!(P, u, p, t) fill!(D, zero(eltype(D))) return nothing end - prob_ip = PDSProblem(linmodP!, linmodD!, u0, tspan, p; analytic = f_analytic) + prob_ip = PDSProblem(linmodPD!, u0, tspan, p; analytic = f_analytic) prob_ip_2 = ConservativePDSProblem(linmodP!, u0, tspan, p; analytic = f_analytic) @@ -1418,7 +1387,7 @@ end # Here we check that in-place and out-of-place implementations # deliver the same results @testset "Different matrix types (nonconservative)" begin - prod_1! = (P, u, p, t) -> begin + prod_dest_1! = (P, D, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) P[i, i + 1] = i * u[i] @@ -1426,17 +1395,16 @@ end for i in 1:length(u) P[i, i] = i * u[i] end - return nothing - end - dest_1! = (D, u, p, t) -> begin + fill!(D, zero(eltype(D))) for i in 1:length(u) D[i] = (i + 1) * u[i] end + return nothing end - prod_2! = (P, u, p, t) -> begin + prod_dest_2! = (P, D, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) P[i + 1, i] = i * u[i + 1] @@ -1444,17 +1412,16 @@ end for i in 1:length(u) P[i, i] = (i - 1) * u[i] end - return nothing - end - dest_2! = (D, u, p, t) -> begin + fill!(D, zero(eltype(D))) for i in 1:length(u) D[i] = i * u[i] end + return nothing end - prod_3! = (P, u, p, t) -> begin + prod_dest_3! = (P, D, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) P[i, i + 1] = i * u[i] @@ -1463,13 +1430,12 @@ end for i in 1:length(u) P[i, i] = (i + 1) * u[i] end - return nothing - end - dest_3! = (D, u, p, t) -> begin + fill!(D, zero(eltype(D))) for i in 1:length(u) D[i] = (i - 1) * u[i] end + return nothing end @@ -1489,29 +1455,24 @@ end MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), SSPMPRK43()) - for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), - (dest_1!, dest_2!, dest_3!)) - prod = (u, p, t) -> begin + for (prod_dest!) in (prod_dest_1!, prod_dest_2!, prod_dest_3!) + prod_dest = (u, p, t) -> begin P = similar(u, (length(u), length(u))) - prod!(P, u, p, t) - return P - end - dest = (u, p, t) -> begin D = similar(u) - dest!(D, u, p, t) - return D + prod_dest!(P, D, u, p, t) + return P, D end - prob_tridiagonal_ip = PDSProblem(prod!, dest!, u0, tspan; + prob_tridiagonal_ip = PDSProblem(prod_dest!, u0, tspan; p_prototype = P_tridiagonal) - prob_tridiagonal_op = PDSProblem(prod, dest, u0, tspan; + prob_tridiagonal_op = PDSProblem(prod_dest, u0, tspan; p_prototype = P_tridiagonal) - prob_dense_ip = PDSProblem(prod!, dest!, u0, tspan; + prob_dense_ip = PDSProblem(prod_dest!, u0, tspan; p_prototype = P_dense) - prob_dense_op = PDSProblem(prod, dest, u0, tspan; + prob_dense_op = PDSProblem(prod_dest, u0, tspan; p_prototype = P_dense) - prob_sparse_ip = PDSProblem(prod!, dest!, u0, tspan; + prob_sparse_ip = PDSProblem(prod_dest!, u0, tspan; p_prototype = P_sparse) - prob_sparse_op = PDSProblem(prod, dest, u0, tspan; + prob_sparse_op = PDSProblem(prod_dest, u0, tspan; p_prototype = P_sparse) sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; @@ -1543,7 +1504,7 @@ end end @testset "Different matrix types (nonconservative, adaptive)" begin - prod_1! = (P, u, p, t) -> begin + prod_dest_1! = (P, D, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) P[i, i + 1] = i * u[i] @@ -1551,17 +1512,16 @@ end for i in 1:length(u) P[i, i] = i * u[i] end - return nothing - end - dest_1! = (D, u, p, t) -> begin + fill!(D, zero(eltype(D))) for i in 1:length(u) D[i] = (i + 1) * u[i] end + return nothing end - prod_2! = (P, u, p, t) -> begin + prod_dest_2! = (P, D, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) P[i + 1, i] = i * u[i + 1] @@ -1569,17 +1529,16 @@ end for i in 1:length(u) P[i, i] = (i - 1) * u[i] end - return nothing - end - dest_2! = (D, u, p, t) -> begin + fill!(D, zero(eltype(D))) for i in 1:length(u) D[i] = i * u[i] end + return nothing end - prod_3! = (P, u, p, t) -> begin + prod_dest_3! = (P, D, u, p, t) -> begin fill!(P, zero(eltype(P))) for i in 1:(length(u) - 1) P[i, i + 1] = i * u[i] @@ -1588,13 +1547,12 @@ end for i in 1:length(u) P[i, i] = (i + 1) * u[i] end - return nothing - end - dest_3! = (D, u, p, t) -> begin + fill!(D, zero(eltype(D))) for i in 1:length(u) D[i] = (i - 1) * u[i] end + return nothing end @@ -1615,31 +1573,24 @@ end MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0), SSPMPRK43()) - for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!), - (dest_1!, dest_2!, dest_3!)) - prod! = prod_3! - dest! = dest_3! - prod = (u, p, t) -> begin + for prod_dest! in (prod_dest_1!, prod_dest_2!, prod_dest_3!) + prod_dest = (u, p, t) -> begin P = similar(u, (length(u), length(u))) - prod!(P, u, p, t) - return P - end - dest = (u, p, t) -> begin D = similar(u) - dest!(D, u, p, t) - return D + prod_dest!(P, D, u, p, t) + return P, D end - prob_tridiagonal_ip = PDSProblem(prod!, dest!, u0, tspan; + prob_tridiagonal_ip = PDSProblem(prod_dest!, u0, tspan; p_prototype = P_tridiagonal) - prob_tridiagonal_op = PDSProblem(prod, dest, u0, tspan; + prob_tridiagonal_op = PDSProblem(prod_dest, u0, tspan; p_prototype = P_tridiagonal) - prob_dense_ip = PDSProblem(prod!, dest!, u0, tspan; + prob_dense_ip = PDSProblem(prod_dest!, u0, tspan; p_prototype = P_dense) - prob_dense_op = PDSProblem(prod, dest, u0, tspan; + prob_dense_op = PDSProblem(prod_dest, u0, tspan; p_prototype = P_dense) - prob_sparse_ip = PDSProblem(prod!, dest!, u0, tspan; + prob_sparse_ip = PDSProblem(prod_dest!, u0, tspan; p_prototype = P_sparse) - prob_sparse_op = PDSProblem(prod, dest, u0, tspan; + prob_sparse_op = PDSProblem(prod_dest, u0, tspan; p_prototype = P_sparse) sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; @@ -1687,18 +1638,33 @@ end prod_inner!(P, u, p, t) return nothing end + prod_dest_sparse! = (P, D, u, p, t) -> begin + @test P isa SparseMatrixCSC + prod_inner!(P, u, p, t) + fill!(D, zero(eltype(D))) + return nothing + end prod_tridiagonal! = (P, u, p, t) -> begin @test P isa Tridiagonal prod_inner!(P, u, p, t) return nothing end + prod_dest_tridiagonal! = (P, D, u, p, t) -> begin + @test P isa Tridiagonal + prod_inner!(P, u, p, t) + fill!(D, zero(eltype(D))) + return nothing + end prod_dense! = (P, u, p, t) -> begin @test P isa Matrix prod_inner!(P, u, p, t) return nothing end - dest! = (D, u, p, t) -> begin + prod_dest_dense! = (P, D, u, p, t) -> begin + @test P isa Matrix + prod_inner!(P, u, p, t) fill!(D, zero(eltype(D))) + return nothing end #prototypes P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3], @@ -1719,12 +1685,12 @@ end prob_sparse = ConservativePDSProblem(prod_sparse!, u0, tspan; p_prototype = P_sparse) ## nonconservative PDS - prob_default2 = PDSProblem(prod_dense!, dest!, u0, tspan) - prob_tridiagonal2 = PDSProblem(prod_tridiagonal!, dest!, u0, tspan; + prob_default2 = PDSProblem(prod_dest_dense!, u0, tspan) + prob_tridiagonal2 = PDSProblem(prod_dest_tridiagonal!, u0, tspan; p_prototype = P_tridiagonal) - prob_dense2 = PDSProblem(prod_dense!, dest!, u0, tspan; + prob_dense2 = PDSProblem(prod_dest_dense!, u0, tspan; p_prototype = P_dense) - prob_sparse2 = PDSProblem(prod_sparse!, dest!, u0, tspan; + prob_sparse2 = PDSProblem(prod_dest_sparse!, u0, tspan; p_prototype = P_sparse) #solve and test for alg in (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), @@ -1896,33 +1862,28 @@ end end @testset "Check convergence order (nonautonomous nonconservative PDS)" begin - prod! = (P, u, p, t) -> begin + prod_dest! = (P, d, u, p, t) -> begin fill!(P, zero(eltype(P))) P[1, 2] = sin(t)^2 * u[2] P[2, 1] = cos(2 * t)^2 * u[1] P[2, 2] = cos(t)^2 * u[2] - return nothing - end - dest! = (d, u, p, t) -> begin + fill!(d, zero(eltype(d))) d[1] = sin(2 * t)^2 * u[1] d[2] = sin(0.5 * t)^2 * u[2] + return nothing end - prod = (u, p, t) -> begin + prod_dest = (u, p, t) -> begin P = similar(u, (length(u), length(u))) - prod!(P, u, p, t) - return P - end - dest = (u, p, t) -> begin d = similar(u, (length(u),)) - dest!(d, u, p, t) - return d + prod_dest!(P, d, u, p, t) + return P, d end u0 = [1.0; 0.0] tspan = (0.0, 1.0) - prob_oop = PDSProblem(prod, dest, u0, tspan) #out-of-place - prob_ip = PDSProblem(prod!, dest!, u0, tspan) #in-place + prob_oop = PDSProblem(prod_dest, u0, tspan) #out-of-place + prob_ip = PDSProblem(prod_dest!, u0, tspan) #in-place dts = 0.5 .^ (4:15) algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), @@ -1964,24 +1925,27 @@ end fill!(P, zero(eltype(P))) P[2, 1] = λ * u[1] end - function dest!(D, u, p, t) + function prod_dest!(P, D, u, p, t) + prod!(P, u, p, t) fill!(D, zero(eltype(D))) + return nothing end function prod(u, p, t) P = similar(u, (length(u), length(u))) prod!(P, u, p, t) return P end - function dest(u, p, t) + function prod_dest(u, p, t) + P = similar(u, (length(u), length(u))) d = similar(u) - dest!(d, u, p, t) - return d + prod_dest!(P, d, u, p, t) + return P, d end prob_ip = ConservativePDSProblem(prod!, u0, tspan, p) - prob_ip_2 = PDSProblem(prod!, dest!, u0, tspan, p) + prob_ip_2 = PDSProblem(prod_dest!, u0, tspan, p) prob_oop = ConservativePDSProblem(prod, u0, tspan, p) - prob_oop_2 = PDSProblem(prod, dest, u0, tspan, p) + prob_oop_2 = PDSProblem(prod_dest, u0, tspan, p) algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), @@ -2012,7 +1976,8 @@ end fill!(P, zero(eltype(P))) P[2, 1] = λ * u[1] end - function dest!(D, u, p, t) + function prod_dest!(P, D, u, p, t) + prod!(P, u, p, t) fill!(D, zero(eltype(D))) end function prod(u, p, t) @@ -2020,16 +1985,17 @@ end prod!(P, u, p, t) return P end - function dest(u, p, t) + function prod_dest(u, p, t) + P = similar(u, (length(u), length(u))) d = similar(u) - dest!(d, u, p, t) - return d + prod_dest!(P, d, u, p, t) + return P, d end prob_ip = ConservativePDSProblem(prod!, u0, tspan, p) - prob_ip_2 = PDSProblem(prod!, dest!, u0, tspan, p) + prob_ip_2 = PDSProblem(prod_dest!, u0, tspan, p) prob_oop = ConservativePDSProblem(prod, u0, tspan, p) - prob_oop_2 = PDSProblem(prod, dest, u0, tspan, p) + prob_oop_2 = PDSProblem(prod_dest, u0, tspan, p) algs = (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK22(2.0), MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75), MPRK43II(0.5), @@ -2104,7 +2070,8 @@ end fill!(P, zero(eltype(P))) P[2, 1] = λ * u[1] end - function dest!(D, u, p, t) + function prod_dest!(P, D, u, p, t) + prod!(P, u, p, t) fill!(D, zero(eltype(D))) end function prod(u, p, t) @@ -2112,15 +2079,16 @@ end prod!(P, u, p, t) return P end - function dest(u, p, t) + function prod_dest(u, p, t) + P = similar(u, (length(u), length(u))) d = similar(u) - dest!(d, u, p, t) + prod_dest!(P, d, u, p, t) return d end prob_ip = ConservativePDSProblem(prod!, u0, tspan, p) - prob_ip_2 = PDSProblem(prod!, dest!, u0, tspan, p) + prob_ip_2 = PDSProblem(prod_dest!, u0, tspan, p) prob_oop = ConservativePDSProblem(prod, u0, tspan, p) - prob_oop_2 = PDSProblem(prod, dest, u0, tspan, p) + prob_oop_2 = PDSProblem(prod_dest, u0, tspan, p) probs = (prob_ip, prob_ip_2, prob_oop, prob_oop_2, prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod, @@ -2153,28 +2121,22 @@ end MPRK43II(2.0 / 3.0), SSPMPRK22(0.5, 1.0), SSPMPRK43()) @testset "$alg, $q" for alg in algs, q in 1:PositiveIntegrators.alg_order(alg) f(t) = q * t^(q - 1) - function prod!(P, u, p, t) + function prod_dest!(P, D, u, p, t) fill!(P, zero(eltype(P))) P[1, 1] = f(t) - end - function dest!(D, u, p, t) fill!(D, zero(eltype(D))) end - function prod(u, p, t) + function prod_dest(u, p, t) P = similar(u, (length(u), length(u))) - prod!(P, u, p, t) - return P - end - function dest(u, p, t) d = similar(u) - dest!(d, u, p, t) - return d + prod_dest!(P, d, u, p, t) + return P, d end u0 = [0.0; 2.0] - prob_oop = PDSProblem(prod, dest, u0, (0.0, 1.0)) + prob_oop = PDSProblem(prod_dest, u0, (0.0, 1.0)) sol_oop = solve(prob_oop, alg, dt = 0.1; adaptive = false) @test first(last(sol_oop.u)) ≈ 1.0 - prob_ip = PDSProblem(prod!, dest!, u0, (0.0, 1.0)) + prob_ip = PDSProblem(prod_dest!, u0, (0.0, 1.0)) sol_ip = solve(prob_ip, alg, dt = 0.1; adaptive = false) @test first(last(sol_ip.u)) ≈ 1.0 end @@ -2186,28 +2148,22 @@ end alg = MPE() u_exact(t) = 1 / (2 + 3 * t) f(t) = -3 / (2 + 3 * t)^2 - function prod!(P, u, p, t) + function prod_dest!(P, D, u, p, t) fill!(P, zero(eltype(P))) - end - function dest!(D, u, p, t) fill!(D, zero(eltype(D))) D[1, 1] = -f(t) end - function prod(u, p, t) + function prod_dest(u, p, t) P = similar(u, (length(u), length(u))) - prod!(P, u, p, t) - return P - end - function dest(u, p, t) d = similar(u) - dest!(d, u, p, t) - return d + prod_dest!(P, d, u, p, t) + return P end u0 = [0.5; 0.0] - prob_oop = PDSProblem(prod, dest, u0, (0.0, 1.0)) + prob_oop = PDSProblem(prod_dest, u0, (0.0, 1.0)) sol_oop = solve(prob_oop, alg, dt = 0.5; adaptive = false) @test first(last(sol_oop.u)) ≈ u_exact(last(prob_oop.tspan)) - prob_ip = PDSProblem(prod!, dest!, u0, (0.0, 1.0)) + prob_ip = PDSProblem(prod_dest!, u0, (0.0, 1.0)) sol_ip = solve(prob_ip, alg, dt = 0.1; adaptive = false) @test first(last(sol_ip.u)) ≈ u_exact(last(prob_ip.tspan)) end @@ -2236,4 +2192,4 @@ end @test isnegative(last(sol_bruss.u)) end end -end; +end