From 258c798273cf4dcc445fad6dfe3c4d3b7a6e433f Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 18 Jul 2024 15:25:09 +0200 Subject: [PATCH] improve performance with sparse matrices (#101) * implement sparse optimization for out-of-place build_mprk_matrix! * fix broadcasting * format * format * add package versions to tutorial * add warning to docs * fix nonconservative systems * fix comment * mention KLUFactorization * add LinearSolve * skip broken tests * Update docs/src/faq.md Co-authored-by: Stefan Kopecz --------- Co-authored-by: Stefan Kopecz --- docs/Project.toml | 6 + docs/make.jl | 1 + docs/src/faq.md | 30 ++++ docs/src/linear_advection.md | 50 +++++-- src/PositiveIntegrators.jl | 3 +- src/mprk.jl | 265 +++++++++++++++++++++++++++++++---- src/proddest.jl | 6 +- src/sspmprk.jl | 144 +++++++++++++++++-- test/runtests.jl | 22 ++- 9 files changed, 468 insertions(+), 59 deletions(-) create mode 100644 docs/src/faq.md diff --git a/docs/Project.toml b/docs/Project.toml index 286d97ce..8b3b9da2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,13 +1,19 @@ [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] BenchmarkTools = "1" Documenter = "1" +InteractiveUtils = "1" +LinearSolve = "2.21" OrdinaryDiffEq = "6.59" +Pkg = "1" Plots = "1" SparseArrays = "1.7" diff --git a/docs/make.jl b/docs/make.jl index e4307ce1..fba3333e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -76,6 +76,7 @@ makedocs(modules = [PositiveIntegrators], "Tutorials" => [ "Linear Advection" => "linear_advection.md", ], + "Troubleshooting, FAQ" => "faq.md", "API reference" => "api_reference.md", "Contributing" => "contributing.md", "Code of conduct" => "code_of_conduct.md", diff --git a/docs/src/faq.md b/docs/src/faq.md new file mode 100644 index 00000000..a112682b --- /dev/null +++ b/docs/src/faq.md @@ -0,0 +1,30 @@ +# Troubleshooting and frequently asked questions + +## Sparse matrices + +You can use sparse matrices for the linear systems arising in +[PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl), +as described, e.g., in the [tutorial on linear advection](@ref tutorial-linear-advection). +However, you need to make sure that you do not change the sparsity pattern +of the production term matrix since we assume that the structural nonzeros +are kept fixed. This is a [known issue](https://github.com/JuliaSparse/SparseArrays.jl/issues/190). +For example, you should avoid something like + +```@repl +using SparseArrays +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. + +```@repl +using SparseArrays +p = spdiagm(0 => ones(4), 1 => zeros(3)) +for j in axes(p, 2) + for idx in nzrange(p, j) + i = rowvals(p)[idx] + nonzeros(p)[idx] = 10 * i + j # value p[i, j] + end +end; p +``` diff --git a/docs/src/linear_advection.md b/docs/src/linear_advection.md index 6694d441..982d6c46 100644 --- a/docs/src/linear_advection.md +++ b/docs/src/linear_advection.md @@ -1,7 +1,7 @@ -# Tutorial: Solution of the linear advection equation +# [Tutorial: Solution of the linear advection equation](@id tutorial-linear-advection) -This tutorial is about the efficient solution of production-destruction systems (PDS) with a large number of differential equations. -We will explore several ways to represent such large systems and assess their efficiency. +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 @@ -11,7 +11,7 @@ One example of the occurrence of a PDS with a large number of equations is the s \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 +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 ones(eltype(p), size(p, 1))) +end ##################################################################### # out-of-place for dense and static arrays @@ -26,7 +36,7 @@ end # in-place for dense arrays @muladd function build_mprk_matrix!(M, P, sigma, dt, d = nothing) - # M[i,i] = (sigma[i] + dt * sum_j P[j,i]) / sigma[i] + # M[i,i] = (sigma[i] + dt * d[i] + dt * sum_j≠i P[j,i]) / sigma[i] # M[i,j] = -dt * P[i,j] / sigma[j] # TODO: the performance of this can likely be improved Base.require_one_based_indexing(M, P, sigma) @@ -71,15 +81,14 @@ end return nothing end -# optimized versions for Tridiagonal matrices +# optimized version for Tridiagonal matrices @muladd function build_mprk_matrix!(M::Tridiagonal, P::Tridiagonal, σ, dt, d = nothing) - # M[i,i] = (sigma[i] + dt * sum_j P[j,i]) / sigma[i] + # M[i,i] = (sigma[i] + dt * d[i] + dt * sum_j≠i P[j,i]) / sigma[i] # M[i,j] = -dt * P[i,j] / sigma[j] Base.require_one_based_indexing(M.dl, M.d, M.du, P.dl, P.d, P.du, σ) @assert length(M.dl) + 1 == length(M.d) == length(M.du) + 1 == length(P.dl) + 1 == length(P.d) == length(P.du) + 1 == length(σ) - if !isnothing(d) Base.require_one_based_indexing(d) @assert length(σ) == length(d) @@ -114,6 +123,100 @@ end return nothing end +# optimized version for sparse matrices +@muladd function build_mprk_matrix!(M::AbstractSparseMatrix, P::AbstractSparseMatrix, + σ, dt, d = nothing) + # M[i,i] = (sigma[i] + dt * d[i] + dt * sum_j≠i P[j,i]) / sigma[i] + # M[i,j] = -dt * P[i,j] / sigma[j] + Base.require_one_based_indexing(M, P, σ) + @assert size(M, 1) == size(M, 2) == size(P, 1) == size(P, 2) == length(σ) + if !isnothing(d) + Base.require_one_based_indexing(d) + @assert length(σ) == length(d) + end + + M_rows = rowvals(M) + M_vals = nonzeros(M) + P_rows = rowvals(P) + P_vals = nonzeros(P) + n = size(M, 2) + + # Diagonal entries + for j in 1:n + for idx_M in nzrange(M, j) + # First, we look for the index corresponding to the + # diagonal entry of M + i = M_rows[idx_M] + if i != j + continue + end + + # Next, we sum over all production terms in the + # i-th column of P except for the diagonal entry + val = zero(eltype(P)) + for idx_P in nzrange(P, i) + if P_rows[idx_P] != i + val += P_vals[idx_P] + end + end + + # Finally, we set the diagonal entry of M + σi = σ[i] + if isnothing(d) + M_vals[idx_M] = (σi + dt * val) / σi + else + M_vals[idx_M] = (σi + dt * d[i] + dt * val) / σi + end + end + end + + # Fill remaining entries + for j in 1:n + iter_M = iterate(nzrange(M, j)) + iter_P = iterate(nzrange(P, j)) + + while iter_M !== nothing + idx_M = first(iter_M) + i = M_rows[idx_M] + + if iter_P === nothing + # If there are no more nonzeros in P, we set the + # entry in M to zero unless we are on the diagonal, + # where we do nothing since we have already set the + # diagonal entries. + if i != j + M_vals[idx_M] = zero(eltype(P)) + end + else + idx_P = first(iter_P) + i_P = P_rows[idx_P] + if i_P < i + throw(ArgumentError("The production matrix has more non-zero entries than the MPRK matrix.")) + elseif i_P == i + # P has a nonzero entry at this position + if i != j + M_vals[idx_M] = -dt * P_vals[idx_P] / σ[j] + end + iter_P = iterate(nzrange(P, j), last(iter_P)) + else + # P has a zero entry at this position + if i != j + M_vals[idx_M] = zero(eltype(P)) + end + end + end + + iter_M = iterate(nzrange(M, j), last(iter_M)) + + if (iter_M === nothing) && (iter_P !== nothing) + throw(ArgumentError("The production matrix has more non-zero entries than the MPRK matrix.")) + end + end + end + + return nothing +end + ### MPE ##################################################################################### """ MPE([linsolve = ..., small_constant = ...]) @@ -138,7 +241,7 @@ You can also choose the parameter `small_constant` which is added to all Patanka to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to `floatmin` of the floating point type used. -The current implementation only supports fixed time steps. +The current implementation only supports fixed time steps. ## References @@ -334,12 +437,12 @@ end A family of second-order modified Patankar-Runge-Kutta algorithms for production-destruction systems. Each member of this family is an adaptive, one-step, two-stage method which is second-order accurate, unconditionally positivity-preserving, and linearly -implicit. In this implementation the stage-values are conservative as well. +implicit. In this implementation the stage-values are conservative as well. The parameter `α` is described by Kopecz and Meister (2018) and studied by Izgin, Kopecz and Meister (2022) as well as -Torlo, Öffner and Ranocha (2022). +Torlo, Öffner and Ranocha (2022). -This method supports adaptive time stepping, using the Patankar-weight denominators +This method supports adaptive time stepping, using the Patankar-weight denominators ``σ_i``, see Kopecz and Meister (2018), as first order approximations to estimate the error. The scheme was introduced by Kopecz and Meister for conservative production-destruction systems. @@ -591,7 +694,16 @@ end f.p(P, uprev, p, t) # evaluate production terms f.d(D, uprev, p, t) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - @.. broadcast=false P2=a21 * P + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=a21 * nz_P + else + @.. broadcast=false P2=a21 * P + end @.. broadcast=false D2=a21 * D # avoid division by zero due to zero Patankar weights @@ -623,7 +735,16 @@ end f.d(D2, u, p, t + a21 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - @.. broadcast=false P2=b1 * P + b2 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=b1 * nz_P + b2 * nz_P2 + else + @.. broadcast=false P2=b1 * P + b2 * P2 + end @.. broadcast=false D2=b1 * D + b2 * D2 # tmp holds the right hand side of the linear system @@ -666,7 +787,16 @@ end # as well as to store the system matrix of the linear system f.p(P, uprev, p, t) # evaluate production terms integrator.stats.nf += 1 - @.. broadcast=false P2=a21 * P + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=a21 * nz_P + else + @.. broadcast=false P2=a21 * P + end # Avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant @@ -690,7 +820,16 @@ end f.p(P2, u, p, t + a21 * dt) # evaluate production terms integrator.stats.nf += 1 - @.. broadcast=false P2=b1 * P + b2 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=b1 * nz_P + b2 * nz_P2 + else + @.. broadcast=false P2=b1 * P + b2 * P2 + end build_mprk_matrix!(P2, P2, σ, dt) @@ -725,7 +864,7 @@ implicit. In this implementation the stage-values are conservative as well. The parameters `α` and `β` must be chosen such that the Runge--Kutta coefficients are nonnegative, see Kopecz and Meister (2018) for details. -These methods support adaptive time stepping, using the Patankar-weight denominators +These methods support adaptive time stepping, using the Patankar-weight denominators ``σ_i``, see Kopecz and Meister (2018), as second order approximations to estimate the error. The scheme was introduced by Kopecz and Meister for conservative production-destruction systems. @@ -828,7 +967,7 @@ implicit. In this implementation the stage-values are conservative as well. The `3/8 ≤ γ ≤ 3/4`. Further details are given in Kopecz and Meister (2018). -This method supports adaptive time stepping, using the Patankar-weight denominators +This method supports adaptive time stepping, using the Patankar-weight denominators ``σ_i``, see Kopecz and Meister (2018), as second order approximations to estimate the error. The scheme was introduced by Kopecz and Meister for conservative production-destruction systems. @@ -1134,7 +1273,16 @@ end f.p(P, uprev, p, t) # evaluate production terms f.d(D, uprev, p, t) # evaluate nonconservative destruction terms - @.. broadcast=false P3=a21 * P + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=a21 * nz_P + else + @.. broadcast=false P3=a21 * P + end @.. broadcast=false D3=a21 * D integrator.stats.nf += 1 @@ -1164,7 +1312,17 @@ end f.p(P2, u, p, t + c2 * dt) # evaluate production terms f.d(D2, u, p, t + c2 * dt) # evaluate nonconservative destruction terms - @.. broadcast=false P3=a31 * P + a32 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=a31 * nz_P + a32 * nz_P2 + else + @.. broadcast=false P3=a31 * P + a32 * P2 + end @.. broadcast=false D3=a31 * D + a32 * D2 integrator.stats.nf += 1 @@ -1188,7 +1346,17 @@ end @.. broadcast=false σ=σ + small_constant end - @.. broadcast=false P3=beta1 * P + beta2 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=beta1 * nz_P + beta2 * nz_P2 + else + @.. broadcast=false P3=beta1 * P + beta2 * P2 + end @.. broadcast=false D3=beta1 * D + beta2 * D2 # tmp holds the right hand side of the linear system @@ -1210,7 +1378,17 @@ end f.p(P3, u, p, t + c3 * dt) # evaluate production terms f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms - @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=b1 * nz_P + b2 * nz_P2 + b3 * nz_P3 + else + @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 + end @.. broadcast=false D3=b1 * D + b2 * D2 + b3 * D3 integrator.stats.nf += 1 @@ -1251,7 +1429,16 @@ 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 - @.. broadcast=false P3=a21 * P + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=a21 * nz_P + else + @.. broadcast=false P3=a21 * P + end integrator.stats.nf += 1 # avoid division by zero due to zero Patankar weights @@ -1273,7 +1460,17 @@ end @.. broadcast=false σ=σ + small_constant f.p(P2, u, p, t + c2 * dt) # evaluate production terms - @.. broadcast=false P3=a31 * P + a32 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=a31 * nz_P + a32 * nz_P2 + else + @.. broadcast=false P3=a31 * P + a32 * P2 + end integrator.stats.nf += 1 build_mprk_matrix!(P3, P3, σ, dt) @@ -1290,7 +1487,17 @@ end @.. broadcast=false σ=σ + small_constant end - @.. broadcast=false P3=beta1 * P + beta2 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=beta1 * nz_P + beta2 * nz_P2 + else + @.. broadcast=false P3=beta1 * P + beta2 * P2 + end build_mprk_matrix!(P3, P3, σ, dt) @@ -1304,7 +1511,17 @@ end @.. broadcast=false σ=σ + small_constant f.p(P3, u, p, t + c3 * dt) # evaluate production terms - @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=b1 * nz_P + b2 * nz_P2 + b3 * nz_P3 + else + @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 + end integrator.stats.nf += 1 build_mprk_matrix!(P3, P3, σ, dt) diff --git a/src/proddest.jl b/src/proddest.jl index 01bb76c3..58e32c3a 100644 --- a/src/proddest.jl +++ b/src/proddest.jl @@ -3,9 +3,9 @@ abstract type AbstractPDSProblem end """ PDSProblem(P, D, u0, tspan, p = NullParameters(); - p_prototype = nothing, - d_prototype = nothing, - analytic = nothing) + p_prototype = nothing, + d_prototype = nothing, + analytic = nothing) A structure describing a system of ordinary differential equations in form of a production-destruction system (PDS). `P` denotes the production matrix. diff --git a/src/sspmprk.jl b/src/sspmprk.jl index bd871d33..89dc362b 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -11,8 +11,8 @@ The difference to [`MPRK22`](@ref) is that this method is based on the SSP formu an explicit second-order Runge-Kutta method. This family of schemes contains the [`MPRK22`](@ref) family, where `MPRK22(α) = SSMPRK22(0, α)` applies. -This method supports adaptive time stepping, using the first order approximations -``(σ_i - u_i^n) / τ + u_i^n`` with ``τ=1+(α_{21}β_{10}^2)/(β_{20}+β_{21})``, +This method supports adaptive time stepping, using the first order approximations +``(σ_i - u_i^n) / τ + u_i^n`` with ``τ=1+(α_{21}β_{10}^2)/(β_{20}+β_{21})``, see (2.7) in Huang and Shu (2019), to estimate the error. The scheme was introduced by Huang and Shu for conservative production-destruction systems. @@ -276,7 +276,16 @@ end f.p(P, uprev, p, t) # evaluate production terms f.d(D, uprev, p, t) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - @.. broadcast=false P2=b10 * P + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=b10 * nz_P + else + @.. broadcast=false P2=b10 * P + end @.. broadcast=false D2=b10 * D # avoid division by zero due to zero Patankar weights @@ -308,7 +317,16 @@ end f.d(D2, u, p, t + b10 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - @.. broadcast=false P2=b20 * P + b21 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=b20 * nz_P + b21 * nz_P2 + else + @.. broadcast=false P2=b20 * P + b21 * P2 + end @.. broadcast=false D2=b20 * D + b21 * D2 # tmp holds the right hand side of the linear system @@ -352,7 +370,16 @@ end # as well as to store the system matrix of the linear system f.p(P, uprev, p, t) # evaluate production terms integrator.stats.nf += 1 - @.. broadcast=false P2=b10 * P + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=b10 * nz_P + else + @.. broadcast=false P2=b10 * P + end # Avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant @@ -379,7 +406,16 @@ end f.p(P2, u, p, t + b10 * dt) # evaluate production terms integrator.stats.nf += 1 - @.. broadcast=false P2=b20 * P + b21 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + @.. broadcast=false nz_P2=b20 * nz_P + b21 * nz_P2 + else + @.. broadcast=false P2=b20 * P + b21 * P2 + end @.. broadcast=false tmp=a20 * uprev + a21 * u @@ -433,7 +469,7 @@ You can also choose the parameter `small_constant` which is added to all Patanka to avoid divisions by zero. You can pass a value explicitly, otherwise `small_constant` is set to `floatmin` of the floating point type used. -The current implementation only supports fixed time steps. +The current implementation only supports fixed time steps. ## References @@ -751,7 +787,16 @@ end f.p(P, uprev, p, t) # evaluate production terms f.d(D, uprev, p, t) # evaluate nonconservative destruction terms - @.. broadcast=false P3=β10 * P + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=β10 * nz_P + else + @.. broadcast=false P3=β10 * P + end @.. broadcast=false D3=β10 * D integrator.stats.nf += 1 @@ -779,7 +824,17 @@ end f.p(P2, u, p, t + β10 * dt) # evaluate production terms f.d(D2, u, p, t + β10 * dt) # evaluate nonconservative destruction terms - @.. broadcast=false P3=β20 * P + β21 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=β20 * nz_P + β21 * nz_P2 + else + @.. broadcast=false P3=β20 * P + β21 * P2 + end @.. broadcast=false D3=β20 * D + β21 * D2 integrator.stats.nf += 1 @@ -801,7 +856,17 @@ end @.. broadcast=false σ=σ^(1 - s) * tmp2^s @.. broadcast=false σ=σ + small_constant - @.. broadcast=false P3=η3 * P + η4 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=η3 * nz_P + η4 * nz_P2 + else + @.. broadcast=false P3=η3 * P + η4 * P2 + end @.. broadcast=false D3=η3 * D + η4 * D2 # tmp holds the right hand side of the linear system @@ -826,7 +891,17 @@ end f.p(P3, u, p, t + c3 * dt) # evaluate production terms f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms - @.. broadcast=false P3=β30 * P + β31 * P2 + β32 * P3 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=β30 * nz_P + β31 * nz_P2 + β32 * nz_P3 + else + @.. broadcast=false P3=β30 * P + β31 * P2 + β32 * P3 + end @.. broadcast=false D3=β30 * D + β31 * D2 + β32 * D3 integrator.stats.nf += 1 @@ -868,7 +943,16 @@ 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 - @.. broadcast=false P3=β10 * P + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=β10 * nz_P + else + @.. broadcast=false P3=β10 * P + end integrator.stats.nf += 1 # avoid division by zero due to zero Patankar weights @@ -891,7 +975,17 @@ end @.. broadcast=false ρ=ρ + small_constant f.p(P2, u, p, t + β10 * dt) # evaluate production terms - @.. broadcast=false P3=β20 * P + β21 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=β20 * nz_P + β21 * nz_P2 + else + @.. broadcast=false P3=β20 * P + β21 * P2 + end integrator.stats.nf += 1 @.. broadcast=false tmp=α20 * uprev + α21 * tmp2 @@ -907,7 +1001,17 @@ end @.. broadcast=false σ=σ^(1 - s) * tmp2^s @.. broadcast=false σ=σ + small_constant - @.. broadcast=false P3=η3 * P + η4 * P2 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=η3 * nz_P + η4 * nz_P2 + else + @.. broadcast=false P3=η3 * P + η4 * P2 + end @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 build_mprk_matrix!(P3, P3, σ, dt) @@ -923,7 +1027,17 @@ end @.. broadcast=false σ=σ + small_constant f.p(P3, u, p, t + c3 * dt) # evaluate production terms - @.. broadcast=false P3=β30 * P + β31 * P2 + β32 * P3 + if issparse(P) + # We need to keep the structural nonzeros of the production terms. + # However, this is not guaranteed by broadcasting, see + # https://github.com/JuliaSparse/SparseArrays.jl/issues/190 + nz_P = nonzeros(P) + nz_P2 = nonzeros(P2) + nz_P3 = nonzeros(P3) + @.. broadcast=false nz_P3=β30 * nz_P + β31 * nz_P2 + β32 * nz_P3 + else + @.. broadcast=false P3=β30 * P + β31 * P2 + β32 * P3 + end integrator.stats.nf += 1 @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u diff --git a/test/runtests.jl b/test/runtests.jl index 3b14f06a..1e13ea52 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,7 +19,7 @@ using Aqua: Aqua ref_alg = TRBDF2(autodiff = false)) Solve `prob` with `alg` and fixed time steps taken from `dts`, and compute -the errors at `test_time`. If`test_time` is not specified the error is computed +the errors at `test_time`. If`test_time` is not specified the error is computed at the final time. Return the associated experimental orders of convergence. @@ -137,7 +137,7 @@ end # This is the usual conservative linear model problem, rewritten as # 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₂) +# 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]] @@ -345,9 +345,9 @@ end end end - # Here we check that solutions of equivalent ODEProblems, PDSProblems or - # ConservativePDS Problems are approximately equal. - # We also check that solvers from OrdinaryDiffEq can solve PDSProblems and + # Here we check that solutions of equivalent ODEProblems, PDSProblems or + # ConservativePDS Problems are approximately equal. + # We also check that solvers from OrdinaryDiffEq can solve PDSProblems and # ConservativePDSProblems. @testset "Check compatibility of PositiveIntegrators and OrdinaryDiffEq" begin @testset "Linear model" begin @@ -485,7 +485,7 @@ end N = 1000 # number of nodes u0 = sin.(π * LinRange(0.0, 1.0, N + 1))[2:end] # initial values tspan = (0.0, 1.0) - # Linear advection discretized with finite differences and upwind, periodic boundary conditions + # Linear advection discretized with finite differences and upwind, periodic boundary conditions linear_advection_fd_upwind_ODE = ODEProblem(linear_advection_fd_upwind_f!, u0, tspan) # problem with dense matrices @@ -897,6 +897,14 @@ 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()) + + # Something is going wrong on macOS with sparse matrices for MPRK43II, see + # https://github.com/SKopecz/PositiveIntegrators.jl/pull/101 + if Sys.isapple() && (alg isa MPRK43II) + @test_skip false + continue + end + for prod! in (prod_1!, prod_2!, prod_3!) prod = (u, p, t) -> begin P = similar(u, (length(u), length(u))) @@ -1503,7 +1511,7 @@ end @testset "$alg" for alg in algs @testset "$i" for (i, prob) in enumerate(probs) if prob == prob_pds_stratreac && alg == SSPMPRK22(0.5, 1.0) - #TODO: SSPMPRK22(0.5, 1.0) is unstable for prob_pds_stratreac. + #TODO: SSPMPRK22(0.5, 1.0) is unstable for prob_pds_stratreac. #Need to figure out if this is a problem of the algorithm or not. break elseif prob == prob_pds_stratreac && alg == MPRK43I(0.5, 0.75)