From b3ab1f23456ecc5913cbb6f80bf03fa08ef7d9ae Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 29 Dec 2025 13:51:14 +0100 Subject: [PATCH 01/29] created function basic_patankar_step for oop --- src/mprk.jl | 42 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 34 insertions(+), 8 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 1328d86d..fc1d3eb3 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -292,19 +292,34 @@ end function initialize!(integrator, cache::MPEConstantCache) end +function basic_patankar_step(uprev, P, σ, dt, linsolve, d = nothing) + if isnothing(d) + M = build_mprk_matrix(P, σ, dt) + linprob = LinearProblem(M, uprev) + else + rhs = uprev + dt * diag(P) + M = build_mprk_matrix(P, σ, dt, d) + linprob = LinearProblem(M, rhs) + end + + # solve linear system + sol = solve(linprob, linsolve) + return sol.u +end + @muladd function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false) @unpack alg, t, dt, uprev, f, p = integrator @unpack small_constant = cache - f = integrator.f - - # evaluate production matrix + # evaluate production matrix and destruction vector P = f.p(uprev, p, t) + d = f isa PDSFunction ? f.d(uprev, p, t) : nothing integrator.stats.nf += 1 # 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 @@ -320,6 +335,9 @@ end # solve linear system sol = solve(linprob, alg.linsolve) u = sol.u + =# + + u = basic_patankar_step(uprev, P, σ, dt, alg.linsolve, d) integrator.stats.nsolve += 1 integrator.u = u @@ -556,16 +574,17 @@ end @unpack alg, t, dt, uprev, f, p = integrator @unpack a21, b1, b2, small_constant = cache - f = integrator.f - # evaluate production matrix P = f.p(uprev, p, t) Ptmp = a21 * P + d = f isa PDSFunction ? f.d(uprev, p, t) : nothing + dtmp = f isa PDSFunction ? a21 * d : nothing integrator.stats.nf += 1 # 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 @@ -582,6 +601,9 @@ end # solve linear system sol = solve(linprob, alg.linsolve) u = sol.u + =# + + u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 # compute Patankar weight denominator @@ -594,13 +616,14 @@ 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 + Ptmp = b1 * P + b2 * f.p(u, p, t + a21 * dt) + dtmp = f isa PDSFunction ? b1 * d + b2 * f.d(u, p, t + a21 * dt) : nothing integrator.stats.nf += 1 + #= # build linear system matrix and rhs if f isa PDSFunction - d2 = f.d(u, p, t + a21 * dt) # evaluate nonconservative destruction terms + d2 = # evaluate nonconservative destruction terms dtmp = b1 * d + b2 * d2 rhs = uprev + dt * diag(Ptmp) M = build_mprk_matrix(Ptmp, σ, dt, dtmp) @@ -614,6 +637,9 @@ end # solve linear system sol = solve(linprob, alg.linsolve) u = sol.u + =# + + u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems. From 62292551366b8c8a3676b054fe4d69901ae74be8 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 29 Dec 2025 16:54:47 +0100 Subject: [PATCH 02/29] revised MPRK43 --- src/mprk.jl | 89 +++++++++++++++-------------------------------------- 1 file changed, 25 insertions(+), 64 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index fc1d3eb3..6de3c9c7 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -294,15 +294,15 @@ end function basic_patankar_step(uprev, P, σ, dt, linsolve, d = nothing) if isnothing(d) + rhs = uprev M = build_mprk_matrix(P, σ, dt) - linprob = LinearProblem(M, uprev) else rhs = uprev + dt * diag(P) M = build_mprk_matrix(P, σ, dt, d) - linprob = LinearProblem(M, rhs) end # solve linear system + linprob = LinearProblem(M, rhs) sol = solve(linprob, linsolve) return sol.u end @@ -319,24 +319,6 @@ end # 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 - rhs = uprev + dt * diag(P) - M = build_mprk_matrix(P, σ, dt, d) - linprob = LinearProblem(M, rhs) - else - # f isa ConservativePDSFunction - M = build_mprk_matrix(P, σ, dt) - linprob = LinearProblem(M, uprev) - end - - # solve linear system - sol = solve(linprob, alg.linsolve) - u = sol.u - =# - u = basic_patankar_step(uprev, P, σ, dt, alg.linsolve, d) integrator.stats.nsolve += 1 @@ -584,25 +566,6 @@ end # 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 - M = build_mprk_matrix(Ptmp, σ, dt) - linprob = LinearProblem(M, uprev) - end - - # solve linear system - sol = solve(linprob, alg.linsolve) - u = sol.u - =# - u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 @@ -620,25 +583,6 @@ end dtmp = f isa PDSFunction ? b1 * d + b2 * f.d(u, p, t + a21 * dt) : nothing integrator.stats.nf += 1 - #= - # build linear system matrix and rhs - if f isa PDSFunction - d2 = # 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 - M = build_mprk_matrix(Ptmp, σ, dt) - linprob = LinearProblem(M, uprev) - end - - # solve linear system - sol = solve(linprob, alg.linsolve) - u = sol.u - =# - u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 @@ -1126,11 +1070,11 @@ end @unpack alg, t, dt, uprev, f, p = integrator @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache - f = integrator.f - # evaluate production matrix P = f.p(uprev, p, t) Ptmp = a21 * P + d = f isa PDSFunction ? f.d(uprev, p, t) : nothing + dtmp = f isa PDSFunction ? a21 * d : nothing integrator.stats.nf += 1 # avoid division by zero due to zero Patankar weights @@ -1139,6 +1083,7 @@ end σ0 = σ end + #= # build linear system matrix and rhs if f isa PDSFunction d = f.d(uprev, p, t) @@ -1154,19 +1099,24 @@ end # solve linear system sol = solve(linprob, alg.linsolve) u2 = sol.u + =# + + u2 = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) u = u2 integrator.stats.nsolve += 1 # compute Patankar weight denominator σ = σ .^ (1 - q1) .* u .^ q1 - # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) P2 = f.p(u, p, t + c2 * dt) Ptmp = a31 * P + a32 * P2 + d2 = f isa PDSFunction ? f.d(u, p, t + c2 * dt) : nothing + dtmp = f isa PDSFunction ? a31 * d + a32 * d2 : nothing 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 @@ -1182,6 +1132,9 @@ end # solve linear system sol = solve(linprob, alg.linsolve) u = sol.u + =# + + u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 # compute Patankar weight denominator @@ -1193,7 +1146,8 @@ end end Ptmp = beta1 * P + beta2 * P2 - + dtmp = f isa PDSFunction ? beta1 * d + beta2 * d2 : nothing + #= # build linear system matrix and rhs if f isa PDSFunction dtmp = beta1 * d + beta2 * d2 @@ -1208,15 +1162,19 @@ end # solve linear system sol = solve(linprob, alg.linsolve) σ = sol.u + =# + + σ = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 # 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 + Ptmp = b1 * P + b2 * P2 + b3 * f.p(u, p, t + c3 * dt) + dtmp = f isa PDSFunction ? b1 * d + b2 * d2 + b3 * f.d(u, p, t + c3 * dt) : nothing integrator.stats.nf += 1 + #= # build linear system matrix if f isa PDSFunction d3 = f.d(u, p, t + c3 * dt) # evaluate nonconservative destruction terms @@ -1232,6 +1190,9 @@ end # solve linear system sol = solve(linprob, alg.linsolve) u = sol.u + =# + + u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 tmp = u - σ From 71f8d53e4cc9b5da7b79588b78b72314ac0c1516 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 29 Dec 2025 21:20:25 +0100 Subject: [PATCH 03/29] revised SSPMPRK oop implementations --- src/mprk.jl | 76 ++------------------------------------------------ src/sspmprk.jl | 48 ++++++++++++++++++++++++++----- 2 files changed, 44 insertions(+), 80 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 6de3c9c7..eeac65f3 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -292,12 +292,12 @@ end function initialize!(integrator, cache::MPEConstantCache) end -function basic_patankar_step(uprev, P, σ, dt, linsolve, d = nothing) +function basic_patankar_step(v, P, σ, dt, linsolve, d = nothing) if isnothing(d) - rhs = uprev + rhs = v M = build_mprk_matrix(P, σ, dt) else - rhs = uprev + dt * diag(P) + rhs = v + dt * diag(P) M = build_mprk_matrix(P, σ, dt, d) end @@ -1083,24 +1083,6 @@ end σ0 = σ 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) - linprob = LinearProblem(M, rhs) - else - M = build_mprk_matrix(Ptmp, σ, dt) - linprob = LinearProblem(M, uprev) - end - - # solve linear system - sol = solve(linprob, alg.linsolve) - u2 = sol.u - =# - u2 = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) u = u2 integrator.stats.nsolve += 1 @@ -1116,24 +1098,6 @@ end dtmp = f isa PDSFunction ? a31 * d + a32 * d2 : nothing 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) - linprob = LinearProblem(M, rhs) - else - M = build_mprk_matrix(Ptmp, σ, dt) - linprob = LinearProblem(M, uprev) - end - - # solve linear system - sol = solve(linprob, alg.linsolve) - u = sol.u - =# - u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 @@ -1147,22 +1111,6 @@ end Ptmp = beta1 * P + beta2 * P2 dtmp = f isa PDSFunction ? beta1 * d + beta2 * d2 : nothing - #= - # build linear system matrix and rhs - if f isa PDSFunction - dtmp = beta1 * d + beta2 * d2 - rhs = uprev + dt * diag(Ptmp) - M = build_mprk_matrix(Ptmp, σ, dt, dtmp) - linprob = LinearProblem(M, rhs) - else - M = build_mprk_matrix(Ptmp, σ, dt) - linprob = LinearProblem(M, uprev) - end - - # solve linear system - sol = solve(linprob, alg.linsolve) - σ = sol.u - =# σ = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 @@ -1174,24 +1122,6 @@ end dtmp = f isa PDSFunction ? b1 * d + b2 * d2 + b3 * f.d(u, p, t + c3 * dt) : nothing integrator.stats.nf += 1 - #= - # 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) - linprob = LinearProblem(M, rhs) - else - M = build_mprk_matrix(Ptmp, σ, dt) - linprob = LinearProblem(M, uprev) - end - - # solve linear system - sol = solve(linprob, alg.linsolve) - u = sol.u - =# - u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 diff --git a/src/sspmprk.jl b/src/sspmprk.jl index f3a325a5..4a62b376 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -124,16 +124,17 @@ end @unpack alg, t, dt, uprev, f, p = integrator @unpack a21, a10, a20, b10, b20, b21, s, τ, small_constant = cache - f = integrator.f - # evaluate production matrix P = f.p(uprev, p, t) Ptmp = b10 * P + d = f isa PDSFunction ? f.d(uprev, p, t) : nothing + dtmp = f isa PDSFunction ? b10 * d : nothing integrator.stats.nf += 1 # 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 @@ -150,6 +151,10 @@ end linprob = LinearProblem(M, rhs) sol = solve(linprob, alg.linsolve) u = sol.u + =# + + v = a10 * uprev + u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 # compute Patankar weight denominator @@ -161,10 +166,11 @@ 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 + Ptmp = b20 * P + b21 * f.p(u, p, t + b10 * dt) + dtmp = f isa PDSFunction ? b20 * d + b21 * f.d(u, p, t + b10 * dt) : nothing integrator.stats.nf += 1 + #= # build linear system matrix and rhs if f isa PDSFunction d2 = f.d(u, p, t + b10 * dt) # evaluate nonconservative destruction terms @@ -181,6 +187,10 @@ end linprob = LinearProblem(M, rhs) sol = solve(linprob, alg.linsolve) u = sol.u + =# + + v = a20 * uprev + a21 * u + u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 # Unless τ = 1, σ is not a first order approximation, since @@ -614,11 +624,14 @@ end # evaluate production matrix P = f.p(uprev, p, t) Ptmp = β10 * P + d = f isa PDSFunction ? f.d(uprev, p, t) : nothing + dtmp = f isa PDSFunction ? β10 * d : nothing integrator.stats.nf += 1 # 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) @@ -634,6 +647,10 @@ end linprob = LinearProblem(M, rhs) sol = solve(linprob, alg.linsolve) u2 = sol.u + =# + + v = α10 * uprev + u2 = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) u = u2 integrator.stats.nsolve += 1 @@ -644,8 +661,11 @@ end P2 = f.p(u, p, t + β10 * dt) Ptmp = β20 * P + β21 * P2 + d2 = f isa PDSFunction ? f.d(u, p, t + β10 * dt) : nothing + dtmp = f isa PDSFunction ? β20 * d + β21 * d2 : nothing integrator.stats.nf += 1 + #= # build linear system matrix and rhs if f isa PDSFunction d2 = f.d(u, p, t + β10 * dt) # evaluate nonconservative destruction terms @@ -662,6 +682,9 @@ end linprob = LinearProblem(M, rhs) sol = solve(linprob, alg.linsolve) u = sol.u + =# + v = α20 * uprev + α21 * u2 + u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 # compute Patankar weight denominator @@ -670,7 +693,9 @@ end σ = add_small_constant(σ, small_constant) Ptmp = η3 * P + η4 * P2 + dtmp = f isa PDSFunction ? η3 * d + η4 * d2 : nothing + #= # build linear system matrix and rhs if f isa PDSFunction dtmp = η3 * d + η4 * d2 @@ -688,6 +713,10 @@ end linprob = LinearProblem(M, rhs) sol = solve(linprob, alg.linsolve) σ = sol.u + =# + + v = η1 * uprev + η2 * u2 + σ = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 # compute Patankar weight denominator @@ -695,10 +724,11 @@ 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 + Ptmp = β30 * P + β31 * P2 + β32 * f.p(u, p, t + c3 * dt) + dtmp = f isa PDSFunction ? β30 * d + β31 * d2 + β32 * f.d(u, p, t + c3 * dt) : nothing integrator.stats.nf += 1 + #= # build linear system matrix if f isa PDSFunction d3 = f.d(u, p, t + c3 * dt) # evaluate nonconservative destruction terms @@ -714,6 +744,10 @@ end linprob = LinearProblem(M, rhs) sol = solve(linprob, alg.linsolve) u = sol.u + =# + + v = α30 * uprev + α31 * u2 + α32 * u + u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 #TODO: Figure out if a second order approximation of the solution @@ -724,7 +758,7 @@ end integrator.opts.reltol, integrator.opts.internalnorm, t) integrator.EEst = integrator.opts.internalnorm(atmp, t) =# - + integrator.u = u end From 75e3f639a5f4a4d19e963bc966b7f23271bb613f Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 29 Dec 2025 21:21:00 +0100 Subject: [PATCH 04/29] format --- src/sspmprk.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 4a62b376..57c5bfcb 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -154,7 +154,7 @@ end =# v = a10 * uprev - u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) + u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 # compute Patankar weight denominator @@ -693,7 +693,7 @@ end σ = add_small_constant(σ, small_constant) Ptmp = η3 * P + η4 * P2 - dtmp = f isa PDSFunction ? η3 * d + η4 * d2 : nothing + dtmp = f isa PDSFunction ? η3 * d + η4 * d2 : nothing #= # build linear system matrix and rhs @@ -758,7 +758,7 @@ end integrator.opts.reltol, integrator.opts.internalnorm, t) integrator.EEst = integrator.opts.internalnorm(atmp, t) =# - + integrator.u = u end From d14efde5dc3f5914e3264e78ac6ccb471887f4e1 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 30 Dec 2025 09:00:33 +0100 Subject: [PATCH 05/29] fix SSPMPRK43 --- src/mprk.jl | 33 ++++++++++++++++++--------------- src/sspmprk.jl | 2 +- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index eeac65f3..44c6333a 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -21,6 +21,24 @@ function _p_prototype(prototype::AbstractSparseMatrix) end ##################################################################### + +function basic_patankar_step(v, P, σ, dt, linsolve, d = nothing, P2 = P) + if isnothing(d) + rhs = v + M = build_mprk_matrix(P, σ, dt) + else + rhs = v + dt * diag(P2) + M = build_mprk_matrix(P, σ, dt, d) + end + + # solve linear system + linprob = LinearProblem(M, rhs) + sol = solve(linprob, linsolve) + return sol.u +end + +##################################################################### + # out-of-place for dense and static arrays function build_mprk_matrix(P, sigma, dt, d = nothing) # re-use the in-place version implemented below @@ -292,21 +310,6 @@ end function initialize!(integrator, cache::MPEConstantCache) end -function basic_patankar_step(v, P, σ, dt, linsolve, d = nothing) - if isnothing(d) - rhs = v - M = build_mprk_matrix(P, σ, dt) - else - rhs = v + dt * diag(P) - M = build_mprk_matrix(P, σ, dt, d) - end - - # solve linear system - linprob = LinearProblem(M, rhs) - sol = solve(linprob, linsolve) - return sol.u -end - @muladd function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false) @unpack alg, t, dt, uprev, f, p = integrator @unpack small_constant = cache diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 57c5bfcb..cc36a914 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -716,7 +716,7 @@ end =# v = η1 * uprev + η2 * u2 - σ = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) + σ = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp, η5 * P + η6 * P2) integrator.stats.nsolve += 1 # compute Patankar weight denominator From fae4106950646f6d38a247d9fe50240a99ffd7a0 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 30 Dec 2025 17:01:41 +0100 Subject: [PATCH 06/29] bugfix SSPMPRK43 --- src/sspmprk.jl | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index cc36a914..1276a0a8 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -631,6 +631,11 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(uprev, small_constant) + v = α10 * uprev + u2 = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) + u = u2 + integrator.stats.nsolve += 1 + #= # build linear system matrix and rhs if f isa PDSFunction @@ -646,13 +651,10 @@ end # solve linear system linprob = LinearProblem(M, rhs) sol = solve(linprob, alg.linsolve) - u2 = sol.u - =# + utmp1 = sol.u - v = α10 * uprev - u2 = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) - u = u2 - integrator.stats.nsolve += 1 + @assert utmp1 ≈ u2 + =# # compute Patankar weight denominator ρ = n1 * u + n2 * u .^ 2 ./ σ @@ -665,7 +667,11 @@ end dtmp = f isa PDSFunction ? β20 * d + β21 * d2 : nothing integrator.stats.nf += 1 - #= + v = α20 * uprev + α21 * u2 + u = basic_patankar_step(v, Ptmp, ρ, dt, alg.linsolve, dtmp) + integrator.stats.nsolve += 1 + + #= # build linear system matrix and rhs if f isa PDSFunction d2 = f.d(u, p, t + β10 * dt) # evaluate nonconservative destruction terms @@ -681,11 +687,10 @@ end # solve linear system linprob = LinearProblem(M, rhs) sol = solve(linprob, alg.linsolve) - u = sol.u + utmp2 = sol.u + + @assert utmp2 ≈ u =# - v = α20 * uprev + α21 * u2 - u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) - integrator.stats.nsolve += 1 # compute Patankar weight denominator σ = σ .^ (1 - s) .* u2 .^ s From 327b15ca9e27dba6bd6ffe349fb41b19988c6943 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 2 Jan 2026 12:42:30 +0100 Subject: [PATCH 07/29] more structure, more dispatch --- src/mprk.jl | 98 ++++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 70 insertions(+), 28 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 44c6333a..3f328acb 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -22,18 +22,63 @@ end ##################################################################### -function basic_patankar_step(v, P, σ, dt, linsolve, d = nothing, P2 = P) - if isnothing(d) - rhs = v - M = build_mprk_matrix(P, σ, dt) - else - rhs = v + dt * diag(P2) - M = build_mprk_matrix(P, σ, dt, d) +@inline function evaluate_pds(f, u, p, t) + # Always evaluate production rates + P = f.p(u, p, t) + + # Evaluate destruction rates only if the function type supports it + d = f isa PDSFunction ? f.d(u, p, t) : nothing + + return P, d +end + +# Linear combination of production matrices and destruction vectors for PDSFunctions +@inline function lincomb(c1::Number, P1, d1::AbstractArray, args...) + Ptmp = c1 .* P1 + dtmp = c1 .* d1 + + # Process additional stages provided in triplets (coeff, P, d) + for i in 1:3:length(args) + c, P, d = args[i], args[i + 1], args[i + 2] + Ptmp .+= c .* P + dtmp .+= c .* d end + return Ptmp, dtmp +end + +# Linear combination of production matrices and destruction vectors for ConservativePDSFunctions +@inline function lincomb(c1::Number, P1, d1::Nothing, args...) + Ptmp = c1 .* P1 + # For conservative systems, d remains nothing + for i in 1:3:length(args) + c, P = args[i], args[i + 1] + # args[i+2] is also nothing, so we skip it + Ptmp .+= c .* P + end + return Ptmp, nothing +end + +# Version for PDSFunctions +function basic_patankar_step(v, P, σ, dt, linsolve, d::AbstractVector, P2 = P) + rhs = v + dt * diag(P2) + M = build_mprk_matrix(P, σ, dt, d) # solve linear system linprob = LinearProblem(M, rhs) sol = solve(linprob, linsolve) + + return sol.u +end + +# Version for ConservativePDSFunctions +function basic_patankar_step(v, P, σ, dt, linsolve, d::Nothing, P2 = P) + rhs = v + M = build_mprk_matrix(P, σ, dt) + + # solve linear system + linprob = LinearProblem(M, rhs) + sol = solve(linprob, linsolve) + return sol.u end @@ -315,8 +360,7 @@ end @unpack small_constant = cache # evaluate production matrix and destruction vector - P = f.p(uprev, p, t) - d = f isa PDSFunction ? f.d(uprev, p, t) : nothing + P, d = evaluate_pds(f, uprev, p, t) integrator.stats.nf += 1 # avoid division by zero due to zero Patankar weights @@ -560,12 +604,11 @@ end @unpack a21, b1, b2, small_constant = cache # evaluate production matrix - P = f.p(uprev, p, t) - Ptmp = a21 * P - d = f isa PDSFunction ? f.d(uprev, p, t) : nothing - dtmp = f isa PDSFunction ? a21 * d : nothing + P, d = evaluate_pds(f, uprev, p, t) integrator.stats.nf += 1 + Ptmp, dtmp = lincomb(a21, P, d) + # avoid division by zero due to zero Patankar weights σ = add_small_constant(uprev, small_constant) @@ -582,10 +625,11 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - Ptmp = b1 * P + b2 * f.p(u, p, t + a21 * dt) - dtmp = f isa PDSFunction ? b1 * d + b2 * f.d(u, p, t + a21 * dt) : nothing + P2, d2 = evaluate_pds(f, u, p, t + a21 * dt) integrator.stats.nf += 1 + Ptmp, dtmp = lincomb(b1, P, d, b2, P2, d2) + u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 @@ -1073,13 +1117,12 @@ end @unpack alg, t, dt, uprev, f, p = integrator @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache - # evaluate production matrix - P = f.p(uprev, p, t) - Ptmp = a21 * P - d = f isa PDSFunction ? f.d(uprev, p, t) : nothing - dtmp = f isa PDSFunction ? a21 * d : nothing + # evaluate production and destruction terms + P, d = evaluate_pds(f, uprev, p, t) integrator.stats.nf += 1 + Ptmp, dtmp = lincomb(a21, P, d) + # avoid division by zero due to zero Patankar weights σ = add_small_constant(uprev, small_constant) if !(q1 ≈ q2) @@ -1095,12 +1138,11 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - P2 = f.p(u, p, t + c2 * dt) - Ptmp = a31 * P + a32 * P2 - d2 = f isa PDSFunction ? f.d(u, p, t + c2 * dt) : nothing - dtmp = f isa PDSFunction ? a31 * d + a32 * d2 : nothing + P2, d2 = evaluate_pds(f, u, p, t + c2 * dt) integrator.stats.nf += 1 + Ptmp, dtmp = lincomb(a31, P, d, a32, P2, d2) + u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 @@ -1112,8 +1154,7 @@ end σ = add_small_constant(σ, small_constant) end - Ptmp = beta1 * P + beta2 * P2 - dtmp = f isa PDSFunction ? beta1 * d + beta2 * d2 : nothing + Ptmp, dtmp = lincomb(beta1, P, d, beta2, P2, d2) σ = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 @@ -1121,10 +1162,11 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - Ptmp = b1 * P + b2 * P2 + b3 * f.p(u, p, t + c3 * dt) - dtmp = f isa PDSFunction ? b1 * d + b2 * d2 + b3 * f.d(u, p, t + c3 * dt) : nothing + P3, d3 = evaluate_pds(f, u, p, t + c3 * dt) integrator.stats.nf += 1 + Ptmp, dtmp = lincomb(b1, P, d, b2, P2, d2, b3, P3, d3) + u = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 From 915b7969e0bf4033d6dbce72a8076a250eeee09c Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 2 Jan 2026 13:42:32 +0100 Subject: [PATCH 08/29] bugfix --- src/mprk.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 3f328acb..0782b207 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -40,8 +40,8 @@ end # Process additional stages provided in triplets (coeff, P, d) for i in 1:3:length(args) c, P, d = args[i], args[i + 1], args[i + 2] - Ptmp .+= c .* P - dtmp .+= c .* d + Ptmp = Ptmp .+ c .* P + dtmp = dtmp .+ c .* d end return Ptmp, dtmp end @@ -53,7 +53,7 @@ end for i in 1:3:length(args) c, P = args[i], args[i + 1] # args[i+2] is also nothing, so we skip it - Ptmp .+= c .* P + Ptmp = Ptmp .+ c .* P end return Ptmp, nothing end From 6f0be7d1ca6f3bf8530287707cf5182efbc99757 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 2 Jan 2026 14:53:18 +0100 Subject: [PATCH 09/29] revised SSPMPRK schemes --- src/sspmprk.jl | 147 +++++-------------------------------------------- 1 file changed, 14 insertions(+), 133 deletions(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 1276a0a8..2b5a3f0b 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -125,34 +125,14 @@ end @unpack a21, a10, a20, b10, b20, b21, s, τ, small_constant = cache # evaluate production matrix - P = f.p(uprev, p, t) - Ptmp = b10 * P - d = f isa PDSFunction ? f.d(uprev, p, t) : nothing - dtmp = f isa PDSFunction ? b10 * d : nothing + P, d = evaluate_pds(f, uprev, p, t) integrator.stats.nf += 1 + Ptmp, dtmp = lincomb(b10, P, d) + # 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 - M = build_mprk_matrix(Ptmp, σ, dt) - rhs = a10 * uprev - end - - # solve linear system - linprob = LinearProblem(M, rhs) - sol = solve(linprob, alg.linsolve) - u = sol.u - =# - v = a10 * uprev u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 @@ -166,28 +146,10 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - Ptmp = b20 * P + b21 * f.p(u, p, t + b10 * dt) - dtmp = f isa PDSFunction ? b20 * d + b21 * f.d(u, p, t + b10 * dt) : nothing + P2, d2 = evaluate_pds(f, u, p, t + b10 * dt) integrator.stats.nf += 1 - #= - # 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 - M = build_mprk_matrix(Ptmp, σ, dt) - rhs = a20 * uprev + a21 * u - end - - # solve linear system - linprob = LinearProblem(M, rhs) - sol = solve(linprob, alg.linsolve) - u = sol.u - =# + Ptmp, tmp = lincomb(b20, P, d, b21, P2, d2) v = a20 * uprev + a21 * u u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) @@ -622,12 +584,11 @@ end f = integrator.f # evaluate production matrix - P = f.p(uprev, p, t) - Ptmp = β10 * P - d = f isa PDSFunction ? f.d(uprev, p, t) : nothing - dtmp = f isa PDSFunction ? β10 * d : nothing + P, d = evaluate_pds(f, uprev, p, t) integrator.stats.nf += 1 + Ptmp, dtmp = lincomb(β10, P, d) + # avoid division by zero due to zero Patankar weights σ = add_small_constant(uprev, small_constant) @@ -636,89 +597,26 @@ end u = u2 integrator.stats.nsolve += 1 - #= - # 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 - rhs = α10 * uprev - M = build_mprk_matrix(Ptmp, σ, dt) - end - - # solve linear system - linprob = LinearProblem(M, rhs) - sol = solve(linprob, alg.linsolve) - utmp1 = sol.u - - @assert utmp1 ≈ u2 - =# - # compute Patankar weight denominator ρ = n1 * u + n2 * u .^ 2 ./ σ # 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 - d2 = f isa PDSFunction ? f.d(u, p, t + β10 * dt) : nothing - dtmp = f isa PDSFunction ? β20 * d + β21 * d2 : nothing + P2, d2 = evaluate_pds(f, u, p, t + β10 * dt) integrator.stats.nf += 1 + Ptmp, dtmp = lincomb(β20, P, d, β21, P2, d2) + v = α20 * uprev + α21 * u2 u = basic_patankar_step(v, Ptmp, ρ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 - #= - # 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 - rhs = α20 * uprev + α21 * u2 - M = build_mprk_matrix(Ptmp, ρ, dt) - end - - # solve linear system - linprob = LinearProblem(M, rhs) - sol = solve(linprob, alg.linsolve) - utmp2 = sol.u - - @assert utmp2 ≈ u - =# - # compute Patankar weight denominator σ = σ .^ (1 - s) .* u2 .^ s # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - Ptmp = η3 * P + η4 * P2 - dtmp = f isa PDSFunction ? η3 * d + η4 * d2 : nothing - - #= - # build linear system matrix and rhs - if f isa PDSFunction - dtmp = η3 * d + η4 * d2 - - # see (3.25 f) in original paper - rhs = η1 * uprev + η2 * u2 + dt * (η5 * diag(P) + η6 * diag(P2)) - - M = build_mprk_matrix(Ptmp, σ, dt, dtmp) - else - rhs = η1 * uprev + η2 * u2 - M = build_mprk_matrix(Ptmp, σ, dt) - end - - # solve linear system - linprob = LinearProblem(M, rhs) - sol = solve(linprob, alg.linsolve) - σ = sol.u - =# + Ptmp, dtmp = lincomb(η3, P, d, η4, P2, d2) v = η1 * uprev + η2 * u2 σ = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp, η5 * P + η6 * P2) @@ -729,27 +627,10 @@ end # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) - Ptmp = β30 * P + β31 * P2 + β32 * f.p(u, p, t + c3 * dt) - dtmp = f isa PDSFunction ? β30 * d + β31 * d2 + β32 * f.d(u, p, t + c3 * dt) : nothing + P3, d3 = evaluate_pds(f, u, p, t + c3 * dt) integrator.stats.nf += 1 - #= - # 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 - rhs = α30 * uprev + α31 * u2 + α32 * u - M = build_mprk_matrix(Ptmp, σ, dt) - end - - # solve linear system - linprob = LinearProblem(M, rhs) - sol = solve(linprob, alg.linsolve) - u = sol.u - =# + Ptmp, dtmp = lincomb(β30, P, d, β31, P2, d2, β32, P3, d3) v = α30 * uprev + α31 * u2 + α32 * u u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) From facb30d7646a87c1666eed7cd9267e86321e8d6a Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 2 Jan 2026 16:10:27 +0100 Subject: [PATCH 10/29] bugfix --- src/sspmprk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 2b5a3f0b..1c645ebe 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -149,7 +149,7 @@ end P2, d2 = evaluate_pds(f, u, p, t + b10 * dt) integrator.stats.nf += 1 - Ptmp, tmp = lincomb(b20, P, d, b21, P2, d2) + Ptmp, dtmp = lincomb(b20, P, d, b21, P2, d2) v = a20 * uprev + a21 * u u = basic_patankar_step(v, Ptmp, σ, dt, alg.linsolve, dtmp) From 4e1bb6ab069091d0fc5968393d83dd1ae703597a Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Sun, 4 Jan 2026 11:58:08 +0100 Subject: [PATCH 11/29] revised lincomb function --- src/PositiveIntegrators.jl | 2 +- src/mprk.jl | 57 +++++++++++++++++++++++--------------- 2 files changed, 35 insertions(+), 24 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index f167f351..1653ec74 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -6,7 +6,7 @@ using Statistics: median using SparseArrays: SparseArrays, AbstractSparseMatrix, issparse, nonzeros, nzrange, rowvals, spdiagm -using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix, MMatrix +using StaticArrays: SVector, SMatrix, StaticArray, StaticMatrix, @SVector, @SMatrix, MMatrix using FastBroadcast: @.. using MuladdMacro: @muladd diff --git a/src/mprk.jl b/src/mprk.jl index 0782b207..e97b1d99 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -32,32 +32,30 @@ end return P, d end -# Linear combination of production matrices and destruction vectors for PDSFunctions -@inline function lincomb(c1::Number, P1, d1::AbstractArray, args...) - Ptmp = c1 .* P1 - dtmp = c1 .* d1 - - # Process additional stages provided in triplets (coeff, P, d) - for i in 1:3:length(args) - c, P, d = args[i], args[i + 1], args[i + 2] - Ptmp = Ptmp .+ c .* P - dtmp = dtmp .+ c .* d - end - return Ptmp, dtmp +### lincomb ##################################################################### +# --- Base cases (End of recursion) --- +# For PDSFunctions (with destruction vectors) +@inline lincomb(c1::Number, P1, d1::AbstractArray) = (c1 .* P1, c1 .* d1) + +# For ConservativePDSFunctions (without destruction vectors) +@inline lincomb(c1::Number, P1, d1::Nothing) = (c1 .* P1, nothing) + +# --- Recursive steps --- + +# For PDSFunctions: Processes triplets (coeff, P, d) +@inline function lincomb(c1::Number, P1, d1::AbstractArray, c2, P2, d2, args...) + P_tail, d_tail = lincomb(c2, P2, d2, args...) + return (c1 .* P1 .+ P_tail, c1 .* d1 .+ d_tail) end -# Linear combination of production matrices and destruction vectors for ConservativePDSFunctions -@inline function lincomb(c1::Number, P1, d1::Nothing, args...) - Ptmp = c1 .* P1 - # For conservative systems, d remains nothing - for i in 1:3:length(args) - c, P = args[i], args[i + 1] - # args[i+2] is also nothing, so we skip it - Ptmp = Ptmp .+ c .* P - end - return Ptmp, nothing +# For ConservativePDSFunctions: Processes triplets (coeff, P, nothing) +@inline function lincomb(c1::Number, P1, d1::Nothing, c2, P2, d2, args...) + P_tail, _ = lincomb(c2, P2, d2, args...) + return (c1 .* P1 .+ P_tail, nothing) end +###################################################################### + # Version for PDSFunctions function basic_patankar_step(v, P, σ, dt, linsolve, d::AbstractVector, P2 = P) rhs = v + dt * diag(P2) @@ -84,7 +82,7 @@ end ##################################################################### -# out-of-place for dense and static arrays +# out-of-place for dense arrays function build_mprk_matrix(P, sigma, dt, d = nothing) # re-use the in-place version implemented below M = similar(P) @@ -97,6 +95,19 @@ function build_mprk_matrix(P, sigma, dt, d = nothing) end end +# out-of-place for static arrays +@inline function build_mprk_matrix_new(P::StaticMatrix{N, N, T}, sigma, dt, + d = nothing) where {N, T} + return SMatrix{N, N, T}((i == j) ? + # diagonal + 1.0 + + (dt / sigma[i]) * + (sum(P[:, i]) - P[i, i] + (d === nothing ? zero(T) : d[i])) : + # off-diagonal + -(dt / sigma[j]) * P[i, j] + for i in 1:N, j in 1:N) +end + # in-place for dense arrays @muladd function build_mprk_matrix!(M, P, sigma, dt, d = nothing) # M[i,i] = (sigma[i] + dt * d[i] + dt * sum_j≠i P[j,i]) / sigma[i] From 1de95f57093eb417d3884d7df9ab7ed4b62b07ba Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Sun, 4 Jan 2026 12:11:35 +0100 Subject: [PATCH 12/29] added @muladd --- src/mprk.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index e97b1d99..b7406269 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -43,13 +43,13 @@ end # --- Recursive steps --- # For PDSFunctions: Processes triplets (coeff, P, d) -@inline function lincomb(c1::Number, P1, d1::AbstractArray, c2, P2, d2, args...) +@muladd @inline function lincomb(c1::Number, P1, d1::AbstractArray, c2, P2, d2, args...) P_tail, d_tail = lincomb(c2, P2, d2, args...) return (c1 .* P1 .+ P_tail, c1 .* d1 .+ d_tail) end # For ConservativePDSFunctions: Processes triplets (coeff, P, nothing) -@inline function lincomb(c1::Number, P1, d1::Nothing, c2, P2, d2, args...) +@muladd @inline function lincomb(c1::Number, P1, d1::Nothing, c2, P2, d2, args...) P_tail, _ = lincomb(c2, P2, d2, args...) return (c1 .* P1 .+ P_tail, nothing) end @@ -57,7 +57,7 @@ end ###################################################################### # Version for PDSFunctions -function basic_patankar_step(v, P, σ, dt, linsolve, d::AbstractVector, P2 = P) +@muladd function basic_patankar_step(v, P, σ, dt, linsolve, d::AbstractVector, P2 = P) rhs = v + dt * diag(P2) M = build_mprk_matrix(P, σ, dt, d) @@ -96,7 +96,7 @@ function build_mprk_matrix(P, sigma, dt, d = nothing) end # out-of-place for static arrays -@inline function build_mprk_matrix_new(P::StaticMatrix{N, N, T}, sigma, dt, +@muladd @inline function build_mprk_matrix_new(P::StaticMatrix{N, N, T}, sigma, dt, d = nothing) where {N, T} return SMatrix{N, N, T}((i == j) ? # diagonal From 5bae76e9c0cca56fd532154fcdea1509ae9dd309 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 9 Jan 2026 11:25:20 +0100 Subject: [PATCH 13/29] build_mprk_matrix for StaticMatrix --- src/mprk.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index b7406269..5b42a0b5 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -96,8 +96,8 @@ function build_mprk_matrix(P, sigma, dt, d = nothing) end # out-of-place for static arrays -@muladd @inline function build_mprk_matrix_new(P::StaticMatrix{N, N, T}, sigma, dt, - d = nothing) where {N, T} +@muladd @inline function build_mprk_matrix(P::StaticMatrix{N, N, T}, sigma, dt, + d = nothing) where {N, T} return SMatrix{N, N, T}((i == j) ? # diagonal 1.0 + From 505c2d011d816f43278b84b4b8ee7d03cd98055e Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 9 Jan 2026 16:19:58 +0100 Subject: [PATCH 14/29] bugfix build_mprk_matrix for StaticMatrix --- src/mprk.jl | 10 ++++++---- test/runtests.jl | 7 +++++++ 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 5b42a0b5..a5a069c9 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -89,6 +89,9 @@ function build_mprk_matrix(P, sigma, dt, d = nothing) build_mprk_matrix!(M, P, sigma, dt, d) if P isa StaticArray + #TODO: This is unnecessary once + # build_mprk_matrix(P::StaticMatrix{N, N, T}, ...) + # is tested successfully. return SMatrix(M) else return M @@ -100,9 +103,10 @@ end d = nothing) where {N, T} return SMatrix{N, N, T}((i == j) ? # diagonal - 1.0 + + one(T) + (dt / sigma[i]) * - (sum(P[:, i]) - P[i, i] + (d === nothing ? zero(T) : d[i])) : + (sum(P[1:(j - 1), i]) + sum(P[(j + 1):end, i]) + + (d === nothing ? zero(T) : d[i])) : # off-diagonal -(dt / sigma[j]) * P[i, j] for i in 1:N, j in 1:N) @@ -1160,7 +1164,6 @@ end # compute Patankar weight denominator if !(q1 ≈ q2) σ = σ0 .^ (1 - q2) .* u2 .^ q2 - # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) end @@ -1169,7 +1172,6 @@ end σ = basic_patankar_step(uprev, Ptmp, σ, dt, alg.linsolve, dtmp) integrator.stats.nsolve += 1 - # avoid division by zero due to zero Patankar weights σ = add_small_constant(σ, small_constant) diff --git a/test/runtests.jl b/test/runtests.jl index faf5862c..34fcfed2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2363,6 +2363,13 @@ end elseif prob == prob_pds_stratreac && alg == MPDeC(9; nodes = :lagrange) # unstable break + elseif prob == prob_pds_stratreac && alg == MPRK22(0.5) + # unstable + # MPRK22(0.5) is known to oscillate for stiff problems. + # See e.g. Izgin, Kopecz, Meister: On the stability of unconditionally + # positive and linear invariants preserving time integration schemes. + # SIAM J Numer Anal, 60, 2022. + break end # later versions of OrdinaryDiffEq.jl use dtmin = 0 by default, # see https://github.com/SciML/OrdinaryDiffEq.jl/pull/2098 From 852dcdab2c1501d33c5475c6d0ad4b8d88b5ff1a Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 9 Jan 2026 16:37:07 +0100 Subject: [PATCH 15/29] build_mprk_matrix for StaticMatrix without slicing allocations --- src/mprk.jl | 16 +++++++++++----- test/runtests.jl | 2 +- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index a5a069c9..5625108b 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -99,14 +99,20 @@ function build_mprk_matrix(P, sigma, dt, d = nothing) end # out-of-place for static arrays +@inline function sum_offdiagonal_col(P::StaticMatrix{N, N, T}, j) where {N, T} + res = zero(T) + for i in 1:N + if i != j + res += P[i, j] + end + end + return res +end @muladd @inline function build_mprk_matrix(P::StaticMatrix{N, N, T}, sigma, dt, d = nothing) where {N, T} return SMatrix{N, N, T}((i == j) ? - # diagonal - one(T) + - (dt / sigma[i]) * - (sum(P[1:(j - 1), i]) + sum(P[(j + 1):end, i]) + - (d === nothing ? zero(T) : d[i])) : + # diagonal: sum over P[k, i] where k != i + one(T) + (dt / sigma[i]) * (sum_offdiagonal_col(P, i) + (d === nothing ? zero(T) : d[i])) : # off-diagonal -(dt / sigma[j]) * P[i, j] for i in 1:N, j in 1:N) diff --git a/test/runtests.jl b/test/runtests.jl index 34fcfed2..cc59f785 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2368,7 +2368,7 @@ end # MPRK22(0.5) is known to oscillate for stiff problems. # See e.g. Izgin, Kopecz, Meister: On the stability of unconditionally # positive and linear invariants preserving time integration schemes. - # SIAM J Numer Anal, 60, 2022. + # SIAM Journal on Numerical Analysis, 60, 2022. break end # later versions of OrdinaryDiffEq.jl use dtmin = 0 by default, From 84ef5d3e2ffd6cffc49c73aa128008af6c7ee71f Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Fri, 9 Jan 2026 16:37:24 +0100 Subject: [PATCH 16/29] format --- src/mprk.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/mprk.jl b/src/mprk.jl index 5625108b..5140d947 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -112,7 +112,9 @@ end d = nothing) where {N, T} return SMatrix{N, N, T}((i == j) ? # diagonal: sum over P[k, i] where k != i - one(T) + (dt / sigma[i]) * (sum_offdiagonal_col(P, i) + (d === nothing ? zero(T) : d[i])) : + one(T) + + (dt / sigma[i]) * + (sum_offdiagonal_col(P, i) + (d === nothing ? zero(T) : d[i])) : # off-diagonal -(dt / sigma[j]) * P[i, j] for i in 1:N, j in 1:N) From 5e611f257223c8f62b30283f7fdaf9762b6e5144 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Sat, 10 Jan 2026 16:12:40 +0100 Subject: [PATCH 17/29] Revised implementation of build_mprk_matrix for StaticMatrix. Now tests pass as before. --- src/mprk.jl | 9 +++++---- test/runtests.jl | 7 ------- 2 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 5140d947..05caf264 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -112,11 +112,12 @@ end d = nothing) where {N, T} return SMatrix{N, N, T}((i == j) ? # diagonal: sum over P[k, i] where k != i - one(T) + - (dt / sigma[i]) * - (sum_offdiagonal_col(P, i) + (d === nothing ? zero(T) : d[i])) : + (sigma[i] + + dt * + (sum_offdiagonal_col(P, i) + (d === nothing ? zero(T) : d[i]))) / + sigma[i] : # off-diagonal - -(dt / sigma[j]) * P[i, j] + -(dt * P[i, j]) / sigma[j] for i in 1:N, j in 1:N) end diff --git a/test/runtests.jl b/test/runtests.jl index cc59f785..faf5862c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2363,13 +2363,6 @@ end elseif prob == prob_pds_stratreac && alg == MPDeC(9; nodes = :lagrange) # unstable break - elseif prob == prob_pds_stratreac && alg == MPRK22(0.5) - # unstable - # MPRK22(0.5) is known to oscillate for stiff problems. - # See e.g. Izgin, Kopecz, Meister: On the stability of unconditionally - # positive and linear invariants preserving time integration schemes. - # SIAM Journal on Numerical Analysis, 60, 2022. - break end # later versions of OrdinaryDiffEq.jl use dtmin = 0 by default, # see https://github.com/SciML/OrdinaryDiffEq.jl/pull/2098 From 49d43e3522c413ac6f84c10dbecefcd4edd26553 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 12 Jan 2026 11:05:59 +0100 Subject: [PATCH 18/29] lincomb! in MPRK schemes --- src/PositiveIntegrators.jl | 2 +- src/mprk.jl | 97 +++++++++++++++++++++++++++++++++++--- 2 files changed, 91 insertions(+), 8 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 1653ec74..e0a83f05 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -4,7 +4,7 @@ module PositiveIntegrators using LinearAlgebra: LinearAlgebra, Tridiagonal, I, diag, mul! using Statistics: median -using SparseArrays: SparseArrays, AbstractSparseMatrix, +using SparseArrays: SparseArrays, AbstractSparseMatrix, SparseMatrixCSC, issparse, nonzeros, nzrange, rowvals, spdiagm using StaticArrays: SVector, SMatrix, StaticArray, StaticMatrix, @SVector, @SMatrix, MMatrix diff --git a/src/mprk.jl b/src/mprk.jl index 05caf264..5c2b0513 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -32,7 +32,7 @@ end return P, d end -### lincomb ##################################################################### +### lincomb and lincomb! ############################################################## # --- Base cases (End of recursion) --- # For PDSFunctions (with destruction vectors) @inline lincomb(c1::Number, P1, d1::AbstractArray) = (c1 .* P1, c1 .* d1) @@ -54,6 +54,40 @@ end return (c1 .* P1 .+ P_tail, nothing) end +@muladd @generated function lincomb!(dest::AbstractArray, pairs...) + n = length(pairs) ÷ 2 + expr = :(pairs[1] * pairs[2]) + for i in 2:n + expr = :($expr + pairs[$(2 * i - 1)] * pairs[$(2 * i)]) + end + return quote + @.. broadcast=false dest=$expr + return dest + end +end + +@muladd @generated function lincomb!(dest::SparseMatrixCSC, pairs...) + n = length(pairs) ÷ 2 + nz_names = [Symbol("nz_", i) for i in 1:n] + setup_block = Expr(:block) + # 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 + for i in 1:n + push!(setup_block.args, :($(nz_names[i]) = nonzeros(pairs[$(2 * i)]))) + end + expr = :(pairs[1] * $(nz_names[1])) + for i in 2:n + expr = :($expr + pairs[$(2 * i - 1)] * $(nz_names[i])) + end + return quote + $setup_block + nz_dest = nonzeros(dest) + @.. broadcast=false nz_dest=$expr + return dest + end +end + ###################################################################### # Version for PDSFunctions @@ -749,6 +783,8 @@ end f.p(P, uprev, p, t) # evaluate production terms f.d(D, uprev, p, t) # evaluate nonconservative destruction terms integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -760,6 +796,9 @@ end @.. broadcast=false P2=a21 * P end @.. broadcast=false D2=a21 * D + =# + lincomb!(P2, a21, P2) + lincomb!(D2, a21, D2) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant @@ -790,6 +829,7 @@ end f.d(D2, u, p, t + a21 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -801,6 +841,9 @@ end @.. broadcast=false P2=b1 * P + b2 * P2 end @.. broadcast=false D2=b1 * D + b2 * D2 + =# + lincomb!(P2, b1, P, b2, P2) + lincomb!(D2, b1, D, b2, D2) # tmp holds the right hand side of the linear system tmp .= uprev @@ -842,6 +885,7 @@ 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 + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -852,6 +896,8 @@ end else @.. broadcast=false P2=a21 * P end + =# + lincomb!(P2, a21, P) # Avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant @@ -875,6 +921,7 @@ end f.p(P2, u, p, t + a21 * dt) # evaluate production terms integrator.stats.nf += 1 + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -885,6 +932,8 @@ end else @.. broadcast=false P2=b1 * P + b2 * P2 end + =# + lincomb!(P2, b1, P, b2, P2) build_mprk_matrix!(P2, P2, σ, dt) @@ -1289,6 +1338,9 @@ end f.p(P, uprev, p, t) # evaluate production terms f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1300,7 +1352,9 @@ end @.. broadcast=false P3=a21 * P end @.. broadcast=false D3=a21 * D - integrator.stats.nf += 1 + =# + lincomb!(P3, a21, P) + lincomb!(D3, a21, D) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant @@ -1328,6 +1382,9 @@ end f.p(P2, u, p, t + c2 * dt) # evaluate production terms f.d(D2, u, p, t + c2 * dt) # evaluate nonconservative destruction terms + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1340,7 +1397,9 @@ end @.. broadcast=false P3=a31 * P + a32 * P2 end @.. broadcast=false D3=a31 * D + a32 * D2 - integrator.stats.nf += 1 + =# + lincomb!(P3, a31, P, a32, P2) + lincomb!(D3, a31, D, a32, D2) # tmp holds the right hand side of the linear system tmp .= uprev @@ -1362,6 +1421,7 @@ end @.. broadcast=false σ=σ + small_constant end + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1374,6 +1434,9 @@ end @.. broadcast=false P3=beta1 * P + beta2 * P2 end @.. broadcast=false D3=beta1 * D + beta2 * D2 + =# + lincomb!(P3, beta1, P, beta2, P2) + lincomb!(D3, beta1, D, beta2, D2) # tmp holds the right hand side of the linear system tmp .= uprev @@ -1394,6 +1457,9 @@ end f.p(P3, u, p, t + c3 * dt) # evaluate production terms f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1406,7 +1472,9 @@ end @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 end @.. broadcast=false D3=b1 * D + b2 * D2 + b3 * D3 - integrator.stats.nf += 1 + =# + lincomb!(P3, b1, P, b2, P2, b3, P3) + lincomb!(D3, b1, D, b2, D2, b3, D3) # tmp holds the right hand side of the linear system tmp .= uprev @@ -1445,6 +1513,9 @@ 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 + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1455,7 +1526,8 @@ end else @.. broadcast=false P3=a21 * P end - integrator.stats.nf += 1 + =# + lincomb!(P3, a21, P) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant @@ -1476,6 +1548,9 @@ end @.. broadcast=false σ=σ + small_constant f.p(P2, u, p, t + c2 * dt) # evaluate production terms + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1487,7 +1562,8 @@ end else @.. broadcast=false P3=a31 * P + a32 * P2 end - integrator.stats.nf += 1 + =# + lincomb!(P3, a31, P, a32, P2) build_mprk_matrix!(P3, P3, σ, dt) @@ -1503,6 +1579,7 @@ end @.. broadcast=false σ=σ + small_constant end + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1514,6 +1591,8 @@ end else @.. broadcast=false P3=beta1 * P + beta2 * P2 end + =# + lincomb!(P3, beta1, P, beta2, P2) build_mprk_matrix!(P3, P3, σ, dt) @@ -1527,6 +1606,9 @@ end @.. broadcast=false σ=σ + small_constant f.p(P3, u, p, t + c3 * dt) # evaluate production terms + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -1538,7 +1620,8 @@ end else @.. broadcast=false P3=b1 * P + b2 * P2 + b3 * P3 end - integrator.stats.nf += 1 + =# + lincomb!(P3, b1, P, b2, P2, b3, P3) build_mprk_matrix!(P3, P3, σ, dt) From 23246713a83934fc43ba500a8de8c008cc1d03d4 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 12 Jan 2026 16:35:36 +0100 Subject: [PATCH 19/29] bugfix MPRK22 --- src/mprk.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 5c2b0513..25e99cc0 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -797,8 +797,8 @@ end end @.. broadcast=false D2=a21 * D =# - lincomb!(P2, a21, P2) - lincomb!(D2, a21, D2) + lincomb!(P2, a21, P) + lincomb!(D2, a21, D) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant From ed7e8dcd7a4cdb82fe6eb9c639331b31b074f780 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 12 Jan 2026 18:18:03 +0100 Subject: [PATCH 20/29] remove stale explicit imports --- src/PositiveIntegrators.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index e0a83f05..566e1a68 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -28,8 +28,7 @@ using LinearSolve: LinearSolve, LinearProblem, LUFactorization, solve! import SciMLBase: interp_summary -using OrdinaryDiffEqCore: @cache, - OrdinaryDiffEqAdaptiveAlgorithm, +using OrdinaryDiffEqCore: OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqConstantCache, OrdinaryDiffEqMutableCache, False, _vec From 4cc972123ffc63a7bddc4d46bfc3c172174fbe64 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Tue, 13 Jan 2026 11:25:00 +0100 Subject: [PATCH 21/29] introduced basic_patankar_step! and basic_patankar_step_conservative! --- src/mprk.jl | 247 ++++++++++++++++++---------------------------------- 1 file changed, 85 insertions(+), 162 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 25e99cc0..c9be2542 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -114,6 +114,39 @@ function basic_patankar_step(v, P, σ, dt, linsolve, d::Nothing, P2 = P) return sol.u end +# non-conservative PDS +@muladd @inline function basic_patankar_step!(u, uprev, P, d, σ, dt, linsolve) + b = linsolve.b + b .= uprev + + # TODO: improve performance for sparse matrices + @inbounds for i in eachindex(u) + b[i] += dt * P[i, i] + end + + build_mprk_matrix!(P, P, σ, dt, d) + + linsolve.A = P + sol = solve!(linsolve) + u .= sol.u + + return nothing +end + +# non-conservative PDS +@inline function basic_patankar_step_conservative!(u, uprev, P, σ, dt, linsolve) + b = linsolve.b + b .= uprev + + build_mprk_matrix!(P, P, σ, dt) + + linsolve.A = P + sol = solve!(linsolve) + u .= sol.u + + return nothing +end + ##################################################################### # out-of-place for dense arrays @@ -504,6 +537,7 @@ end # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant + #= linsolve_rhs .= uprev @inbounds for i in eachindex(linsolve_rhs) linsolve_rhs[i] += dt * P[i, i] @@ -516,6 +550,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, uprev, P, D, σ, dt, linsolve) integrator.stats.nsolve += 1 end @@ -535,6 +571,7 @@ end # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant + #= build_mprk_matrix!(P, P, σ, dt) # Same as linres = P \ uprev @@ -542,6 +579,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, uprev, P, σ, dt, linsolve) integrator.stats.nsolve += 1 end @@ -784,25 +823,13 @@ end f.d(D, uprev, p, t) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P2, a21, P) lincomb!(D2, a21, D) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant + #= # tmp holds the right hand side of the linear system tmp .= uprev @inbounds for i in eachindex(tmp) @@ -816,6 +843,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, uprev, P2, D2, σ, dt, linsolve) integrator.stats.nsolve += 1 if isone(a21) @@ -829,22 +858,10 @@ end f.d(D2, u, p, t + a21 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P2, b1, P, b2, P2) lincomb!(D2, b1, D, b2, D2) + #= # tmp holds the right hand side of the linear system tmp .= uprev @inbounds for i in eachindex(tmp) @@ -858,6 +875,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, uprev, P2, D2, σ, dt, linsolve) integrator.stats.nsolve += 1 # Now σ stores the error estimate @@ -878,30 +897,20 @@ end @unpack tmp, P, P2, σ, linsolve = cache @unpack a21, b1, b2, small_constant = cache.tab - # Set right hand side of linear system - tmp .= uprev - # 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 integrator.stats.nf += 1 - #= - 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 - =# + lincomb!(P2, a21, P) # Avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant + #= + # Set right hand side of linear system + tmp .= uprev + build_mprk_matrix!(P2, P2, σ, dt) # Same as linres = P2 \ tmp @@ -909,6 +918,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, uprev, P2, σ, dt, linsolve) integrator.stats.nsolve += 1 if isone(a21) @@ -921,20 +932,9 @@ end f.p(P2, u, p, t + a21 * dt) # evaluate production terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P2, b1, P, b2, P2) + #= build_mprk_matrix!(P2, P2, σ, dt) # Same as linres = P2 \ tmp @@ -942,6 +942,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, uprev, P2, σ, dt, linsolve) integrator.stats.nsolve += 1 # Now σ stores the error estimate @@ -1340,25 +1342,13 @@ end f.d(D, uprev, p, t) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, a21, P) lincomb!(D3, a21, D) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant + #= # tmp holds the right hand side of the linear system tmp .= uprev @inbounds for i in eachindex(tmp) @@ -1372,10 +1362,13 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, uprev, P3, D3, σ, dt, linsolve) + integrator.stats.nsolve += 1 + 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 @@ -1384,23 +1377,10 @@ end f.d(D2, u, p, t + c2 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, a31, P, a32, P2) lincomb!(D3, a31, D, a32, D2) + #= # tmp holds the right hand side of the linear system tmp .= uprev @inbounds for i in eachindex(tmp) @@ -1414,6 +1394,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, uprev, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 if !(q1 ≈ q2) @@ -1421,23 +1403,10 @@ end @.. broadcast=false σ=σ + small_constant end - #= - 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 - =# lincomb!(P3, beta1, P, beta2, P2) lincomb!(D3, beta1, D, beta2, D2) + #= # tmp holds the right hand side of the linear system tmp .= uprev @inbounds for i in eachindex(tmp) @@ -1449,9 +1418,11 @@ end # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) + σ .= linres + =# + basic_patankar_step!(σ, uprev, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 - σ .= linres # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=σ + small_constant @@ -1459,23 +1430,10 @@ end f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, b1, P, b2, P2, b3, P3) lincomb!(D3, b1, D, b2, D2, b3, D3) + #= # tmp holds the right hand side of the linear system tmp .= uprev @inbounds for i in eachindex(tmp) @@ -1489,6 +1447,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, uprev, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 # Now tmp stores the error estimate @@ -1515,23 +1475,12 @@ end f.p(P, uprev, p, t) # evaluate production terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, a21, P) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant + #= build_mprk_matrix!(P3, P3, σ, dt) # Same as linres = P3 \ tmp @@ -1539,6 +1488,10 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, uprev, P3, σ, dt, linsolve) + integrator.stats.nsolve += 1 + if !(q1 ≈ q2) tmp2 .= u #u2 in out-of-place version end @@ -1550,21 +1503,9 @@ end f.p(P2, u, p, t + c2 * dt) # evaluate production terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, a31, P, a32, P2) + #= build_mprk_matrix!(P3, P3, σ, dt) # Same as linres = P3 \ tmp @@ -1572,6 +1513,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, uprev, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 if !(q1 ≈ q2) @@ -1579,50 +1522,28 @@ end @.. broadcast=false σ=σ + small_constant end - #= - 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 - =# lincomb!(P3, beta1, P, beta2, P2) + #= build_mprk_matrix!(P3, P3, σ, dt) # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) + σ .= linres + =# + basic_patankar_step_conservative!(σ, uprev, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 - σ .= linres # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=σ + small_constant f.p(P3, u, p, t + c3 * dt) # evaluate production terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, b1, P, b2, P2, b3, P3) + #= build_mprk_matrix!(P3, P3, σ, dt) # Same as linres = P3 \ tmp @@ -1630,6 +1551,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, uprev, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 # Now tmp stores the error estimate From 302836f4c7a91b7342108682a0dd18b97829ee0b Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 14 Jan 2026 10:07:11 +0100 Subject: [PATCH 22/29] revised sspmprk --- src/mprk.jl | 197 +++++-------------------------------------------- src/sspmprk.jl | 110 ++++++++++++++++++++++++--- 2 files changed, 116 insertions(+), 191 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index c9be2542..86c85657 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -115,9 +115,9 @@ function basic_patankar_step(v, P, σ, dt, linsolve, d::Nothing, P2 = P) end # non-conservative PDS -@muladd @inline function basic_patankar_step!(u, uprev, P, d, σ, dt, linsolve) +@muladd @inline function basic_patankar_step!(u, v, P, d, σ, dt, linsolve) b = linsolve.b - b .= uprev + b .= v # TODO: improve performance for sparse matrices @inbounds for i in eachindex(u) @@ -134,9 +134,9 @@ end end # non-conservative PDS -@inline function basic_patankar_step_conservative!(u, uprev, P, σ, dt, linsolve) +@inline function basic_patankar_step_conservative!(u, v, P, σ, dt, linsolve) b = linsolve.b - b .= uprev + b .= v build_mprk_matrix!(P, P, σ, dt) @@ -537,20 +537,6 @@ end # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant - #= - linsolve_rhs .= uprev - @inbounds for i in eachindex(linsolve_rhs) - linsolve_rhs[i] += dt * P[i, i] - end - - build_mprk_matrix!(P, P, σ, dt, D) - - # Same as linres = P \ linsolve_rhs - linsolve.A = P - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step!(u, uprev, P, D, σ, dt, linsolve) integrator.stats.nsolve += 1 end @@ -571,15 +557,6 @@ end # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant - #= - build_mprk_matrix!(P, P, σ, dt) - - # Same as linres = P \ uprev - linsolve.A = P - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step_conservative!(u, uprev, P, σ, dt, linsolve) integrator.stats.nsolve += 1 end @@ -829,22 +806,8 @@ end # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant - #= - # tmp holds the right hand side of the linear system tmp .= uprev - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P2[i, i] - end - - build_mprk_matrix!(P2, P2, σ, dt, D2) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step!(u, uprev, P2, D2, σ, dt, linsolve) + basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve) integrator.stats.nsolve += 1 if isone(a21) @@ -861,22 +824,8 @@ end lincomb!(P2, b1, P, b2, P2) lincomb!(D2, b1, D, b2, D2) - #= - # tmp holds the right hand side of the linear system tmp .= uprev - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P2[i, i] - end - - build_mprk_matrix!(P2, P2, σ, dt, D2) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step!(u, uprev, P2, D2, σ, dt, linsolve) + basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve) integrator.stats.nsolve += 1 # Now σ stores the error estimate @@ -907,19 +856,8 @@ end # Avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant - #= - # Set right hand side of linear system tmp .= uprev - - build_mprk_matrix!(P2, P2, σ, dt) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step_conservative!(u, uprev, P2, σ, dt, linsolve) + basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve) integrator.stats.nsolve += 1 if isone(a21) @@ -934,16 +872,7 @@ end lincomb!(P2, b1, P, b2, P2) - #= - build_mprk_matrix!(P2, P2, σ, dt) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step_conservative!(u, uprev, P2, σ, dt, linsolve) + basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve) integrator.stats.nsolve += 1 # Now σ stores the error estimate @@ -1348,22 +1277,8 @@ end # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant - #= - # tmp holds the right hand side of the linear system tmp .= uprev - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P3[i, i] - end - - build_mprk_matrix!(P3, P3, σ, dt, D3) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step!(u, uprev, P3, D3, σ, dt, linsolve) + basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 if !(q1 ≈ q2) @@ -1380,22 +1295,8 @@ end lincomb!(P3, a31, P, a32, P2) lincomb!(D3, a31, D, a32, D2) - #= - # tmp holds the right hand side of the linear system tmp .= uprev - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P3[i, i] - end - - build_mprk_matrix!(P3, P3, σ, dt, D3) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step!(u, uprev, P3, D3, σ, dt, linsolve) + basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 if !(q1 ≈ q2) @@ -1406,21 +1307,8 @@ end lincomb!(P3, beta1, P, beta2, P2) lincomb!(D3, beta1, D, beta2, D2) - #= - # tmp holds the right hand side of the linear system tmp .= uprev - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P3[i, i] - end - - build_mprk_matrix!(P3, P3, σ, dt, D3) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - σ .= linres - =# - basic_patankar_step!(σ, uprev, P3, D3, σ, dt, linsolve) + basic_patankar_step!(σ, tmp, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 # avoid division by zero due to zero Patankar weights @@ -1433,22 +1321,8 @@ end lincomb!(P3, b1, P, b2, P2, b3, P3) lincomb!(D3, b1, D, b2, D2, b3, D3) - #= - # tmp holds the right hand side of the linear system tmp .= uprev - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P3[i, i] - end - - build_mprk_matrix!(P3, P3, σ, dt, D3) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step!(u, uprev, P3, D3, σ, dt, linsolve) + basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 # Now tmp stores the error estimate @@ -1467,9 +1341,6 @@ end @unpack tmp, tmp2, P, P2, P3, σ, linsolve = cache @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache.tab - # Set right hand side of linear system - tmp .= uprev - # 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 @@ -1480,16 +1351,8 @@ end # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant - #= - build_mprk_matrix!(P3, P3, σ, dt) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step_conservative!(u, uprev, P3, σ, dt, linsolve) + tmp .= uprev + basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 if !(q1 ≈ q2) @@ -1505,16 +1368,7 @@ end lincomb!(P3, a31, P, a32, P2) - #= - build_mprk_matrix!(P3, P3, σ, dt) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step_conservative!(u, uprev, P3, σ, dt, linsolve) + basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 if !(q1 ≈ q2) @@ -1524,15 +1378,7 @@ end lincomb!(P3, beta1, P, beta2, P2) - #= - build_mprk_matrix!(P3, P3, σ, dt) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - σ .= linres - =# - basic_patankar_step_conservative!(σ, uprev, P3, σ, dt, linsolve) + basic_patankar_step_conservative!(σ, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 # avoid division by zero due to zero Patankar weights @@ -1543,16 +1389,7 @@ end lincomb!(P3, b1, P, b2, P2, b3, P3) - #= - build_mprk_matrix!(P3, P3, σ, dt) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# - basic_patankar_step_conservative!(u, uprev, P3, σ, dt, linsolve) + basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 # Now tmp stores the error estimate diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 1c645ebe..4d50625b 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -252,6 +252,8 @@ end f.p(P, uprev, p, t) # evaluate production terms f.d(D, uprev, p, t) # evaluate nonconservative destruction terms integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -263,12 +265,17 @@ end @.. broadcast=false P2=b10 * P end @.. broadcast=false D2=b10 * D + =# + lincomb!(P2, b10, P) + lincomb!(D2, b10, D) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant # tmp holds the right hand side of the linear system @.. broadcast=false tmp=a10 * uprev + + #= @inbounds for i in eachindex(tmp) tmp[i] += dt * P2[i, i] end @@ -280,6 +287,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve) integrator.stats.nsolve += 1 if isone(s) @@ -293,6 +302,7 @@ end f.d(D2, u, p, t + b10 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -304,9 +314,13 @@ end @.. broadcast=false P2=b20 * P + b21 * P2 end @.. broadcast=false D2=b20 * D + b21 * D2 + =# + lincomb!(P2, b20, P, b21, P2) + lincomb!(D2, b20, D, b21, D2) # tmp holds the right hand side of the linear system @.. broadcast=false tmp=a20 * uprev + a21 * u + #= @inbounds for i in eachindex(tmp) tmp[i] += dt * P2[i, i] end @@ -318,6 +332,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve) integrator.stats.nsolve += 1 # Unless τ = 1, σ is not a first order approximation, since @@ -346,6 +362,8 @@ 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 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -356,6 +374,8 @@ end else @.. broadcast=false P2=b10 * P end + =# + lincomb!(P2, b10, P) # Avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant @@ -363,6 +383,7 @@ end # tmp holds the right hand side of the linear system @.. broadcast=false tmp=a10 * uprev + #= build_mprk_matrix!(P2, P2, σ, dt) # Same as linres = P2 \ tmp @@ -370,6 +391,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve) integrator.stats.nsolve += 1 if isone(s) @@ -382,6 +405,7 @@ end f.p(P2, u, p, t + b10 * dt) # evaluate production terms integrator.stats.nf += 1 + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -392,9 +416,11 @@ end else @.. broadcast=false P2=b20 * P + b21 * P2 end + =# + lincomb!(P2, b20, P, b21, P2) @.. broadcast=false tmp=a20 * uprev + a21 * u - + #= build_mprk_matrix!(P2, P2, σ, dt) # Same as linres = P2 \ tmp @@ -402,6 +428,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve) integrator.stats.nsolve += 1 # Unless τ = 1, σ is not a first order approximation, since @@ -734,6 +762,9 @@ end f.p(P, uprev, p, t) # evaluate production terms f.d(D, uprev, p, t) # evaluate nonconservative destruction terms + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -745,13 +776,16 @@ end @.. broadcast=false P3=β10 * P end @.. broadcast=false D3=β10 * D - integrator.stats.nf += 1 + =# + lincomb!(P3, β10, P) + lincomb!(D3, β10, D) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant # tmp holds the right hand side of the linear system @.. broadcast=false tmp=α10 * uprev + #= @inbounds for i in eachindex(tmp) tmp[i] += dt * P3[i, i] end @@ -763,14 +797,20 @@ end linres = solve!(linsolve) u .= linres - tmp2 .= u + =# + basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 + tmp2 .= u + @.. 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 + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -783,10 +823,14 @@ end @.. broadcast=false P3=β20 * P + β21 * P2 end @.. broadcast=false D3=β20 * D + β21 * D2 - integrator.stats.nf += 1 + =# + lincomb!(P3, β20, P, β21, P2) + lincomb!(D3, β20, D, β21, D2) # tmp holds the right hand side of the linear system @.. broadcast=false tmp=α20 * uprev + α21 * tmp2 + + #= @inbounds for i in eachindex(tmp) tmp[i] += dt * P3[i, i] end @@ -798,11 +842,14 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, tmp, P3, D3, ρ, dt, linsolve) integrator.stats.nsolve += 1 @.. broadcast=false σ=σ^(1 - s) * tmp2^s @.. broadcast=false σ=σ + small_constant + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -815,6 +862,9 @@ end @.. broadcast=false P3=η3 * P + η4 * P2 end @.. broadcast=false D3=η3 * D + η4 * D2 + =# + lincomb!(P3, η3, P, η4, P2) + lincomb!(D3, η3, D, η4, D2) # tmp holds the right hand side of the linear system @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 @@ -828,6 +878,7 @@ end # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) + integrator.stats.nsolve += 1 σ .= linres @@ -838,6 +889,9 @@ end f.p(P3, u, p, t + c3 * dt) # evaluate production terms f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -850,10 +904,13 @@ end @.. broadcast=false P3=β30 * P + β31 * P2 + β32 * P3 end @.. broadcast=false D3=β30 * D + β31 * D2 + β32 * D3 - integrator.stats.nf += 1 + =# + lincomb!(P3, β30, P, β31, P2, β32, P3) + lincomb!(D3, β30, D, β31, D2, β32, D3) # tmp holds the right hand side of the linear system @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u + #= @inbounds for i in eachindex(tmp) tmp[i] += dt * P3[i, i] end @@ -865,6 +922,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 #TODO: Figure out if a second order approximation of the solution @@ -890,6 +949,9 @@ 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 + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -900,7 +962,8 @@ end else @.. broadcast=false P3=β10 * P end - integrator.stats.nf += 1 + =# + lincomb!(P3, β10, P) # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=uprev + small_constant @@ -908,6 +971,7 @@ end # tmp holds the right hand side of the linear system @.. broadcast=false tmp=α10 * uprev + #= build_mprk_matrix!(P3, P3, σ, dt) # Same as linres = P3 \ tmp @@ -915,13 +979,19 @@ end linres = solve!(linsolve) u .= linres - tmp2 .= u + =# + basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 + tmp2 .= u + @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ @.. broadcast=false ρ=ρ + small_constant f.p(P2, u, p, t + β10 * dt) # evaluate production terms + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -933,9 +1003,11 @@ end else @.. broadcast=false P3=β20 * P + β21 * P2 end - integrator.stats.nf += 1 + =# + lincomb!(P3, β20, P, β21, P2) @.. broadcast=false tmp=α20 * uprev + α21 * tmp2 + #= build_mprk_matrix!(P3, P3, ρ, dt) # Same as linres = P3 \ tmp @@ -943,11 +1015,14 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, tmp, P3, ρ, dt, linsolve) integrator.stats.nsolve += 1 @.. broadcast=false σ=σ^(1 - s) * tmp2^s @.. broadcast=false σ=σ + small_constant + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -959,21 +1034,30 @@ end else @.. broadcast=false P3=η3 * P + η4 * P2 end - @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 + =# + lincomb!(P3, η3, P, η4, P2) + @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 + #= build_mprk_matrix!(P3, P3, σ, dt) # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) - integrator.stats.nsolve += 1 σ .= linres + =# + basic_patankar_step_conservative!(σ, tmp, P3, σ, dt, linsolve) + integrator.stats.nsolve += 1 + @.. broadcast=false σ=σ + z * uprev * u / ρ # avoid division by zero due to zero Patankar weights @.. broadcast=false σ=σ + small_constant f.p(P3, u, p, t + c3 * dt) # evaluate production terms + integrator.stats.nf += 1 + + #= if issparse(P) # We need to keep the structural nonzeros of the production terms. # However, this is not guaranteed by broadcasting, see @@ -985,9 +1069,11 @@ end else @.. broadcast=false P3=β30 * P + β31 * P2 + β32 * P3 end - integrator.stats.nf += 1 + =# + lincomb!(P3, β30, P, β31, P2, β32, P3) @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u + #= build_mprk_matrix!(P3, P3, σ, dt) # Same as linres = P3 \ tmp @@ -995,6 +1081,8 @@ end linres = solve!(linsolve) u .= linres + =# + basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 #TODO: Figure out if a second order approximation of the solution From f06680383c6757692dafa2f5bfce5e38806a13d2 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 14 Jan 2026 13:13:55 +0100 Subject: [PATCH 23/29] removed comments --- src/sspmprk.jl | 275 ------------------------------------------------- 1 file changed, 275 deletions(-) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 4d50625b..7daad55b 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -253,19 +253,6 @@ end f.d(D, uprev, p, t) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P2, b10, P) lincomb!(D2, b10, D) @@ -274,20 +261,6 @@ end # tmp holds the right hand side of the linear system @.. broadcast=false tmp=a10 * uprev - - #= - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P2[i, i] - end - - build_mprk_matrix!(P2, P2, σ, dt, D2) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve) integrator.stats.nsolve += 1 @@ -302,37 +275,11 @@ end f.d(D2, u, p, t + b10 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P2, b20, P, b21, P2) lincomb!(D2, b20, D, b21, D2) # tmp holds the right hand side of the linear system @.. broadcast=false tmp=a20 * uprev + a21 * u - #= - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P2[i, i] - end - - build_mprk_matrix!(P2, P2, σ, dt, D2) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve) integrator.stats.nsolve += 1 @@ -363,18 +310,6 @@ end f.p(P, uprev, p, t) # evaluate production terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P2, b10, P) # Avoid division by zero due to zero Patankar weights @@ -382,16 +317,6 @@ end # tmp holds the right hand side of the linear system @.. broadcast=false tmp=a10 * uprev - - #= - build_mprk_matrix!(P2, P2, σ, dt) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve) integrator.stats.nsolve += 1 @@ -405,30 +330,9 @@ end f.p(P2, u, p, t + b10 * dt) # evaluate production terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P2, b20, P, b21, P2) @.. broadcast=false tmp=a20 * uprev + a21 * u - #= - build_mprk_matrix!(P2, P2, σ, dt) - - # Same as linres = P2 \ tmp - linsolve.A = P2 - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve) integrator.stats.nsolve += 1 @@ -764,19 +668,6 @@ end f.d(D, uprev, p, t) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, β10, P) lincomb!(D3, β10, D) @@ -785,19 +676,6 @@ end # tmp holds the right hand side of the linear system @.. broadcast=false tmp=α10 * uprev - #= - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P3[i, i] - end - - build_mprk_matrix!(P3, P3, σ, dt, D3) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 @@ -810,59 +688,17 @@ end f.d(D2, u, p, t + β10 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, β20, P, β21, P2) lincomb!(D3, β20, D, β21, D2) # tmp holds the right hand side of the linear system @.. broadcast=false tmp=α20 * uprev + α21 * tmp2 - - #= - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P3[i, i] - end - - build_mprk_matrix!(P3, P3, ρ, dt, D3) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step!(u, tmp, P3, D3, ρ, dt, linsolve) integrator.stats.nsolve += 1 @.. broadcast=false σ=σ^(1 - s) * tmp2^s @.. broadcast=false σ=σ + small_constant - #= - 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 - =# lincomb!(P3, η3, P, η4, P2) lincomb!(D3, η3, D, η4, D2) @@ -891,38 +727,11 @@ end f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, β30, P, β31, P2, β32, P3) lincomb!(D3, β30, D, β31, D2, β32, D3) # tmp holds the right hand side of the linear system @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u - #= - @inbounds for i in eachindex(tmp) - tmp[i] += dt * P3[i, i] - end - - build_mprk_matrix!(P3, P3, σ, dt, D3) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve) integrator.stats.nsolve += 1 @@ -951,18 +760,6 @@ end f.p(P, uprev, p, t) # evaluate production terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, β10, P) # avoid division by zero due to zero Patankar weights @@ -971,15 +768,6 @@ end # tmp holds the right hand side of the linear system @.. broadcast=false tmp=α10 * uprev - #= - build_mprk_matrix!(P3, P3, σ, dt) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 @@ -991,62 +779,20 @@ end f.p(P2, u, p, t + β10 * dt) # evaluate production terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, β20, P, β21, P2) @.. broadcast=false tmp=α20 * uprev + α21 * tmp2 - #= - build_mprk_matrix!(P3, P3, ρ, dt) - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - u .= linres - =# basic_patankar_step_conservative!(u, tmp, P3, ρ, dt, linsolve) integrator.stats.nsolve += 1 @.. broadcast=false σ=σ^(1 - s) * tmp2^s @.. broadcast=false σ=σ + small_constant - #= - 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 - =# lincomb!(P3, η3, P, η4, P2) @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 - #= - build_mprk_matrix!(P3, P3, σ, dt) - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - - σ .= linres - =# basic_patankar_step_conservative!(σ, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 @@ -1057,31 +803,10 @@ end f.p(P3, u, p, t + c3 * dt) # evaluate production terms integrator.stats.nf += 1 - #= - 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 - =# lincomb!(P3, β30, P, β31, P2, β32, P3) @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u - #= - build_mprk_matrix!(P3, P3, σ, dt) - - # Same as linres = P3 \ tmp - linsolve.A = P3 - linres = solve!(linsolve) - u .= linres - =# basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 From 3249f15efc1d5627335ffdcce6ef0cecbb76ae0e Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 14 Jan 2026 13:31:15 +0100 Subject: [PATCH 24/29] add diagonal of sparse matrix efficiently --- src/mprk.jl | 22 +++++++++++++++++++++- src/sspmprk.jl | 16 ++++++++++------ 2 files changed, 31 insertions(+), 7 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 86c85657..efdb4402 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -115,14 +115,34 @@ function basic_patankar_step(v, P, σ, dt, linsolve, d::Nothing, P2 = P) end # non-conservative PDS +@inline function add_diagonal!(b, P, dt) + if P isa SparseMatrixCSC + for j in 1:size(P, 2) + for idx in nzrange(P, j) + if rowvals(P)[idx] == j + b[j] += dt * nonzeros(P)[idx] + break + end + end + end + else + @inbounds for i in eachindex(b) + b[i] += dt * P[i, i] + end + end + return nothing +end + @muladd @inline function basic_patankar_step!(u, v, P, d, σ, dt, linsolve) b = linsolve.b b .= v - # TODO: improve performance for sparse matrices + #= @inbounds for i in eachindex(u) b[i] += dt * P[i, i] end + =# + add_diagonal!(b, P, dt) build_mprk_matrix!(P, P, σ, dt, d) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 7daad55b..497cdaec 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -702,22 +702,27 @@ end lincomb!(P3, η3, P, η4, P2) lincomb!(D3, η3, D, η4, D2) - # tmp holds the right hand side of the linear system + # The next stage is the only stage that is not suited + # for a direct application of basic_patankar_step! + # We therefore build the system matrix explicitly. @.. broadcast=false tmp=η1 * uprev + η2 * tmp2 + # see (3.25 f) in original paper + #= @inbounds for i in eachindex(tmp) - # see (3.25 f) in original paper + tmp[i] += dt * (η5 * P[i, i] + η6 * P2[i, i]) end + =# + add_diagonal!(tmp, P, dt * η5) + add_diagonal!(tmp, P2, dt * η6) build_mprk_matrix!(P3, P3, σ, dt, D3) # Same as linres = P3 \ tmp linsolve.A = P3 linres = solve!(linsolve) - - integrator.stats.nsolve += 1 - σ .= linres + integrator.stats.nsolve += 1 @.. broadcast=false σ=σ + z * uprev * u / ρ # avoid division by zero due to zero Patankar weights @@ -806,7 +811,6 @@ end lincomb!(P3, β30, P, β31, P2, β32, P3) @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u - basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve) integrator.stats.nsolve += 1 From 5faece70a183b94a90a121ecf62ff426a85efafc Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Wed, 14 Jan 2026 16:19:53 +0100 Subject: [PATCH 25/29] use naive loop to add diagonal of sparse matrices --- src/mprk.jl | 22 +++------------------- 1 file changed, 3 insertions(+), 19 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index efdb4402..27f28708 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -115,20 +115,9 @@ function basic_patankar_step(v, P, σ, dt, linsolve, d::Nothing, P2 = P) end # non-conservative PDS -@inline function add_diagonal!(b, P, dt) - if P isa SparseMatrixCSC - for j in 1:size(P, 2) - for idx in nzrange(P, j) - if rowvals(P)[idx] == j - b[j] += dt * nonzeros(P)[idx] - break - end - end - end - else - @inbounds for i in eachindex(b) - b[i] += dt * P[i, i] - end +@muladd @inline function add_diagonal!(b, P, dt) + @inbounds for i in eachindex(b) + b[i] += dt * P[i, i] end return nothing end @@ -137,11 +126,6 @@ end b = linsolve.b b .= v - #= - @inbounds for i in eachindex(u) - b[i] += dt * P[i, i] - end - =# add_diagonal!(b, P, dt) build_mprk_matrix!(P, P, σ, dt, d) From 27e683fb9eb9cc71b2d49bdb98774da3c44f1d3f Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 15 Jan 2026 09:56:08 +0100 Subject: [PATCH 26/29] consolidate broadcasts --- src/mprk.jl | 22 ++++++++-------------- src/sspmprk.jl | 16 ++++++---------- 2 files changed, 14 insertions(+), 24 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 27f28708..b4393d3c 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -815,11 +815,10 @@ end integrator.stats.nsolve += 1 if isone(a21) - σ .= u + @.. broadcast=false σ=u + small_constant else - @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) + @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) + small_constant 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 @@ -865,11 +864,10 @@ end integrator.stats.nsolve += 1 if isone(a21) - σ .= u + @.. broadcast=false σ=u + small_constant else - @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) + @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) + small_constant end - @.. broadcast=false σ=σ + small_constant f.p(P2, u, p, t + a21 * dt) # evaluate production terms integrator.stats.nf += 1 @@ -1289,8 +1287,7 @@ end tmp2 .= u #u2 in out-of-place version end - @.. broadcast=false σ=σ^(1 - q1) * u^q1 - @.. broadcast=false σ=σ + small_constant + @.. broadcast=false σ=σ^(1 - q1) * u^q1 + 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 @@ -1304,8 +1301,7 @@ end integrator.stats.nsolve += 1 if !(q1 ≈ q2) - @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 - @.. broadcast=false σ=σ + small_constant + @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 + small_constant end lincomb!(P3, beta1, P, beta2, P2) @@ -1364,8 +1360,7 @@ end end integrator.stats.nsolve += 1 - @.. broadcast=false σ=σ^(1 - q1) * u^q1 - @.. broadcast=false σ=σ + small_constant + @.. broadcast=false σ=σ^(1 - q1) * u^q1 + small_constant f.p(P2, u, p, t + c2 * dt) # evaluate production terms integrator.stats.nf += 1 @@ -1376,8 +1371,7 @@ end integrator.stats.nsolve += 1 if !(q1 ≈ q2) - @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 - @.. broadcast=false σ=σ + small_constant + @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 + small_constant end lincomb!(P3, beta1, P, beta2, P2) diff --git a/src/sspmprk.jl b/src/sspmprk.jl index 497cdaec..a24135a6 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -265,11 +265,10 @@ end integrator.stats.nsolve += 1 if isone(s) - σ .= u + @.. broadcast=false σ=u + small_constant else - @.. broadcast=false σ=σ^(1 - s) * u^s + @.. broadcast=false σ=σ^(1 - s) * u^s + small_constant 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 @@ -321,11 +320,10 @@ end integrator.stats.nsolve += 1 if isone(s) - σ .= u + @.. broadcast=false σ=u + small_constant else - @.. broadcast=false σ=σ^(1 - s) * u^s + @.. broadcast=false σ=σ^(1 - s) * u^s + small_constant end - @.. broadcast=false σ=σ + small_constant f.p(P2, u, p, t + b10 * dt) # evaluate production terms integrator.stats.nf += 1 @@ -696,8 +694,7 @@ end basic_patankar_step!(u, tmp, P3, D3, ρ, dt, linsolve) integrator.stats.nsolve += 1 - @.. broadcast=false σ=σ^(1 - s) * tmp2^s - @.. broadcast=false σ=σ + small_constant + @.. broadcast=false σ=σ^(1 - s) * tmp2^s + small_constant lincomb!(P3, η3, P, η4, P2) lincomb!(D3, η3, D, η4, D2) @@ -791,8 +788,7 @@ end basic_patankar_step_conservative!(u, tmp, P3, ρ, dt, linsolve) integrator.stats.nsolve += 1 - @.. broadcast=false σ=σ^(1 - s) * tmp2^s - @.. broadcast=false σ=σ + small_constant + @.. broadcast=false σ=σ^(1 - s) * tmp2^s + small_constant lincomb!(P3, η3, P, η4, P2) From 0ab817c2a8f0ae6a914f7278e88c04c111010207 Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Thu, 15 Jan 2026 12:04:05 +0100 Subject: [PATCH 27/29] removed @unpack --- Project.toml | 2 -- src/PositiveIntegrators.jl | 1 - src/mpdec.jl | 16 +++++----- src/mprk.jl | 62 +++++++++++++++++--------------------- src/sspmprk.jl | 32 ++++++++++---------- 5 files changed, 52 insertions(+), 61 deletions(-) diff --git a/Project.toml b/Project.toml index 6ca792d4..cfb82651 100644 --- a/Project.toml +++ b/Project.toml @@ -12,7 +12,6 @@ OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" -SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -34,7 +33,6 @@ OrdinaryDiffEqCore = "1.21, 2" RecipesBase = "1.3.4" Reexport = "1.2.2" SciMLBase = "2.78" -SimpleUnPack = "1" SparseArrays = "1" StaticArrays = "1.9.7" Statistics = "1" diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl index 566e1a68..eba70eb3 100644 --- a/src/PositiveIntegrators.jl +++ b/src/PositiveIntegrators.jl @@ -10,7 +10,6 @@ using StaticArrays: SVector, SMatrix, StaticArray, StaticMatrix, @SVector, @SMat using FastBroadcast: @.. using MuladdMacro: @muladd -using SimpleUnPack: @unpack using Reexport: @reexport diff --git a/src/mpdec.jl b/src/mpdec.jl index 5fb83dbb..b6f46bbf 100644 --- a/src/mpdec.jl +++ b/src/mpdec.jl @@ -564,8 +564,8 @@ function _build_mpdec_matrix_and_rhs!(M::AbstractSparseMatrix, rhs, P::AbstractS end @muladd function perform_step!(integrator, cache::MPDeCConstantCache, repeat_step = false) - @unpack alg, t, dt, uprev, f, p = integrator - @unpack K, M, nodes, theta, small_constant, ValM = cache + (; alg, t, dt, uprev, f, p) = integrator + (; K, M, nodes, theta, small_constant, ValM) = cache C = cmatrix(uprev, ValM) C2 = cmatrix(uprev, ValM) @@ -702,9 +702,9 @@ function initialize!(integrator, cache::Union{MPDeCCache, MPDeCConservativeCache end @muladd function perform_step!(integrator, cache::MPDeCCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, P, P2, d, σ, C, C2, linsolve_rhs, linsolve = cache - @unpack K, M, nodes, theta, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, P, P2, d, σ, C, C2, linsolve_rhs, linsolve) = cache + (; K, M, nodes, theta, small_constant) = cache.tab # Initialize C matrices for i in 1:(M + 1) @@ -743,9 +743,9 @@ end @muladd function perform_step!(integrator, cache::MPDeCConservativeCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, P, P2, σ, C, C2, linsolve_rhs, linsolve = cache - @unpack K, M, nodes, theta, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, P, P2, σ, C, C2, linsolve_rhs, linsolve) = cache + (; K, M, nodes, theta, small_constant) = cache.tab # Initialize right hand side of linear system linsolve_rhs .= uprev diff --git a/src/mprk.jl b/src/mprk.jl index b4393d3c..3ae3bd8a 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -153,23 +153,17 @@ end ##################################################################### -# out-of-place for dense arrays +# out-of-place for dense matrices function build_mprk_matrix(P, sigma, dt, d = nothing) # re-use the in-place version implemented below M = similar(P) + build_mprk_matrix!(M, P, sigma, dt, d) - if P isa StaticArray - #TODO: This is unnecessary once - # build_mprk_matrix(P::StaticMatrix{N, N, T}, ...) - # is tested successfully. - return SMatrix(M) - else - return M - end + return M end -# out-of-place for static arrays +# out-of-place for StaticMatrix @inline function sum_offdiagonal_col(P::StaticMatrix{N, N, T}, j) where {N, T} res = zero(T) for i in 1:N @@ -451,8 +445,8 @@ function initialize!(integrator, cache::MPEConstantCache) end @muladd function perform_step!(integrator, cache::MPEConstantCache, repeat_step = false) - @unpack alg, t, dt, uprev, f, p = integrator - @unpack small_constant = cache + (; alg, t, dt, uprev, f, p) = integrator + (; small_constant) = cache # evaluate production matrix and destruction vector P, d = evaluate_pds(f, uprev, p, t) @@ -525,9 +519,9 @@ function initialize!(integrator, cache::Union{MPECache, MPEConservativeCache}) end @muladd function perform_step!(integrator, cache::MPECache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack P, D, σ, linsolve, linsolve_rhs = cache - @unpack small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; P, D, σ, linsolve, linsolve_rhs) = cache + (; small_constant) = cache.tab # We use P to store the evaluation of the PDS # as well as to store the system matrix of the linear system @@ -546,9 +540,9 @@ end end @muladd function perform_step!(integrator, cache::MPEConservativeCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack P, σ, linsolve = cache - @unpack small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; P, σ, linsolve) = cache + (; small_constant) = cache.tab # We use P to store the evaluation of the PDS # as well as to store the system matrix of the linear system @@ -678,8 +672,8 @@ function initialize!(integrator, cache::MPRK22ConstantCache) end @muladd function perform_step!(integrator, cache::MPRK22ConstantCache, repeat_step = false) - @unpack alg, t, dt, uprev, f, p = integrator - @unpack a21, b1, b2, small_constant = cache + (; alg, t, dt, uprev, f, p) = integrator + (; a21, b1, b2, small_constant) = cache # evaluate production matrix P, d = evaluate_pds(f, uprev, p, t) @@ -793,9 +787,9 @@ function initialize!(integrator, cache::Union{MPRK22Cache, MPRK22ConservativeCac end @muladd function perform_step!(integrator, cache::MPRK22Cache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, P, P2, D, D2, σ, linsolve = cache - @unpack a21, b1, b2, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, P, P2, D, D2, σ, linsolve) = cache + (; a21, b1, b2, small_constant) = cache.tab # We use P2 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system @@ -845,9 +839,9 @@ end @muladd function perform_step!(integrator, cache::MPRK22ConservativeCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, P, P2, σ, linsolve = cache - @unpack a21, b1, b2, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, P, P2, σ, linsolve) = cache + (; a21, b1, b2, small_constant) = cache.tab # We use P2 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system @@ -1120,8 +1114,8 @@ function initialize!(integrator, cache::MPRK43ConstantCache) end @muladd function perform_step!(integrator, cache::MPRK43ConstantCache, repeat_step = false) - @unpack alg, t, dt, uprev, f, p = integrator - @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache + (; alg, t, dt, uprev, f, p) = integrator + (; a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant) = cache # evaluate production and destruction terms P, d = evaluate_pds(f, uprev, p, t) @@ -1262,9 +1256,9 @@ function initialize!(integrator, cache::Union{MPRK43Cache, MPRK43ConservativeCac end @muladd function perform_step!(integrator, cache::MPRK43Cache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, tmp2, P, P2, P3, D, D2, D3, σ, linsolve = cache - @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, tmp2, P, P2, P3, D, D2, D3, σ, linsolve) = cache + (; a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant) = cache.tab # We use P3 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system @@ -1337,9 +1331,9 @@ end @muladd function perform_step!(integrator, cache::MPRK43ConservativeCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, tmp2, P, P2, P3, σ, linsolve = cache - @unpack a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, tmp2, P, P2, P3, σ, linsolve) = cache + (; a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2, small_constant) = cache.tab # We use P3 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system diff --git a/src/sspmprk.jl b/src/sspmprk.jl index a24135a6..5b615257 100644 --- a/src/sspmprk.jl +++ b/src/sspmprk.jl @@ -121,8 +121,8 @@ end @muladd function perform_step!(integrator, cache::SSPMPRK22ConstantCache, repeat_step = false) - @unpack alg, t, dt, uprev, f, p = integrator - @unpack a21, a10, a20, b10, b20, b21, s, τ, small_constant = cache + (; alg, t, dt, uprev, f, p) = integrator + (; a21, a10, a20, b10, b20, b21, s, τ, small_constant) = cache # evaluate production matrix P, d = evaluate_pds(f, uprev, p, t) @@ -242,9 +242,9 @@ function initialize!(integrator, cache::Union{SSPMPRK22Cache, SSPMPRK22Conservat end @muladd function perform_step!(integrator, cache::SSPMPRK22Cache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, P, P2, D, D2, σ, linsolve = cache - @unpack a21, a10, a20, b10, b20, b21, s, τ, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, P, P2, D, D2, σ, linsolve) = cache + (; a21, a10, a20, b10, b20, b21, s, τ, small_constant) = cache.tab # We use P2 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system @@ -300,9 +300,9 @@ end @muladd function perform_step!(integrator, cache::SSPMPRK22ConservativeCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, P, P2, σ, linsolve = cache - @unpack a21, a10, a20, b10, b20, b21, s, τ, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, P, P2, σ, linsolve) = cache + (; a21, a10, a20, b10, b20, b21, s, τ, small_constant) = cache.tab # We use P2 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system @@ -508,8 +508,8 @@ end @muladd function perform_step!(integrator, cache::SSPMPRK43ConstantCache, repeat_step = false) - @unpack alg, t, dt, uprev, f, p = integrator - @unpack n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache + (; alg, t, dt, uprev, f, p) = integrator + (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) = cache f = integrator.f @@ -655,9 +655,9 @@ function initialize!(integrator, cache::Union{SSPMPRK43ConservativeCache, SSPMPR end @muladd function perform_step!(integrator, cache::SSPMPRK43Cache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, linsolve = cache - @unpack n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, linsolve) = cache + (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) = cache.tab # We use P3 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system @@ -753,9 +753,9 @@ end @muladd function perform_step!(integrator, cache::SSPMPRK43ConservativeCache, repeat_step = false) - @unpack t, dt, uprev, u, f, p = integrator - @unpack tmp, tmp2, P, P2, P3, σ, ρ, linsolve = cache - @unpack n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant = cache.tab + (; t, dt, uprev, u, f, p) = integrator + (; tmp, tmp2, P, P2, P3, σ, ρ, linsolve) = cache + (; n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) = cache.tab # We use P3 to store the last evaluation of the PDS # as well as to store the system matrix of the linear system From 7427043656ce43082856e6f3ab63788fd02a384f Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 19 Jan 2026 16:25:25 +0100 Subject: [PATCH 28/29] additional comments in source file --- src/mprk.jl | 60 +++++++++++++++++++++++++++++++++++------------------ 1 file changed, 40 insertions(+), 20 deletions(-) diff --git a/src/mprk.jl b/src/mprk.jl index 3ae3bd8a..4ac12b8b 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -1,16 +1,25 @@ -# Helper functions +##################################################################### +### Helper functions +##################################################################### + +### add_small_constant ############################################## +# add_small_constant is used in out-of-place implementations to add +# the safety parameter small_constant to avoid divisions by zero. add_small_constant(v, small_constant) = v .+ small_constant function add_small_constant(v::SVector{N, T}, small_constant::T) where {N, T} v + SVector{N, T}(ntuple(i -> small_constant, N)) end -##################################################################### +### p_prototype ###################################################### +# p_prototype is used to allocate memory for production matrices +# in in-place implementations. p_prototype(u, f) = zeros(eltype(u), length(u), length(u)) + p_prototype(u, f::ConservativePDSFunction) = _p_prototype(f.p_prototype) p_prototype(u, f::PDSFunction) = _p_prototype(f.p_prototype) - _p_prototype(prototype) = zero(prototype) + function _p_prototype(prototype::AbstractSparseMatrix) # We need to ensure that we store all structural nonzeros that # are required for the linear system. In particular, we need to @@ -20,19 +29,25 @@ function _p_prototype(prototype::AbstractSparseMatrix) return p + spdiagm(0 => ones(eltype(p), size(p, 1))) end -##################################################################### - -@inline function evaluate_pds(f, u, p, t) - # Always evaluate production rates +### evaluate_pds ###################################################### +# evaluate_pds is used in out-of-place implementations to evaluate +# the production and destruction terms +@inline function evaluate_pds(f::PDSFunction, u, p, t) P = f.p(u, p, t) - - # Evaluate destruction rates only if the function type supports it d = f isa PDSFunction ? f.d(u, p, t) : nothing - return P, d end -### lincomb and lincomb! ############################################################## +@inline function evaluate_pds(f::ConservativePDSFunction, u, p, t) + P = f.p(u, p, t) + return P, nothing +end + +### lincomb and lincomb! ############################################### +# These functions are used to compute linear combinations of production +# matrices and/or destruction vectors. + +# out-of-place versions # --- Base cases (End of recursion) --- # For PDSFunctions (with destruction vectors) @inline lincomb(c1::Number, P1, d1::AbstractArray) = (c1 .* P1, c1 .* d1) @@ -41,7 +56,6 @@ end @inline lincomb(c1::Number, P1, d1::Nothing) = (c1 .* P1, nothing) # --- Recursive steps --- - # For PDSFunctions: Processes triplets (coeff, P, d) @muladd @inline function lincomb(c1::Number, P1, d1::AbstractArray, c2, P2, d2, args...) P_tail, d_tail = lincomb(c2, P2, d2, args...) @@ -54,6 +68,7 @@ end return (c1 .* P1 .+ P_tail, nothing) end +# in-place version for general arrays @muladd @generated function lincomb!(dest::AbstractArray, pairs...) n = length(pairs) ÷ 2 expr = :(pairs[1] * pairs[2]) @@ -66,6 +81,7 @@ end end end +# in-place version for sparse matrices @muladd @generated function lincomb!(dest::SparseMatrixCSC, pairs...) n = length(pairs) ÷ 2 nz_names = [Symbol("nz_", i) for i in 1:n] @@ -88,9 +104,11 @@ end end end -###################################################################### +### basic_patankar_step, basic_patankar_step!, basic_patankar_step_conservative! ########### +# These functions implement the core Patankar step by building and solving the linear system. -# Version for PDSFunctions +# out-of-place implementations +# Version for PDSFunctions (non-conservative PDS) @muladd function basic_patankar_step(v, P, σ, dt, linsolve, d::AbstractVector, P2 = P) rhs = v + dt * diag(P2) M = build_mprk_matrix(P, σ, dt, d) @@ -102,7 +120,7 @@ end return sol.u end -# Version for ConservativePDSFunctions +# Version for ConservativePDSFunctions (conservative PDS) function basic_patankar_step(v, P, σ, dt, linsolve, d::Nothing, P2 = P) rhs = v M = build_mprk_matrix(P, σ, dt) @@ -114,6 +132,7 @@ function basic_patankar_step(v, P, σ, dt, linsolve, d::Nothing, P2 = P) return sol.u end +# in-place implementations # non-conservative PDS @muladd @inline function add_diagonal!(b, P, dt) @inbounds for i in eachindex(b) @@ -121,7 +140,6 @@ end end return nothing end - @muladd @inline function basic_patankar_step!(u, v, P, d, σ, dt, linsolve) b = linsolve.b b .= v @@ -137,7 +155,7 @@ end return nothing end -# non-conservative PDS +# conservative PDS @inline function basic_patankar_step_conservative!(u, v, P, σ, dt, linsolve) b = linsolve.b b .= v @@ -151,7 +169,9 @@ end return nothing end -##################################################################### +### build_mprk_matrix ############################################################ +# These functions build the system matrix M that needs to be solved in each +# Patankar step. # out-of-place for dense matrices function build_mprk_matrix(P, sigma, dt, d = nothing) @@ -233,7 +253,7 @@ end return nothing end -# optimized version for Tridiagonal matrices +# in-place optimized version for Tridiagonal matrices @muladd function build_mprk_matrix!(M::Tridiagonal, P::Tridiagonal, σ, 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] @@ -275,7 +295,7 @@ end return nothing end -# optimized version for sparse matrices +# in-place 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] From f52a7c722283019968e40085eec41a8098f89daf Mon Sep 17 00:00:00 2001 From: Stefan Kopecz Date: Mon, 19 Jan 2026 16:35:57 +0100 Subject: [PATCH 29/29] Update src/mprk.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> --- src/mprk.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mprk.jl b/src/mprk.jl index 4ac12b8b..ed4a90ef 100644 --- a/src/mprk.jl +++ b/src/mprk.jl @@ -34,7 +34,7 @@ end # the production and destruction terms @inline function evaluate_pds(f::PDSFunction, u, p, t) P = f.p(u, p, t) - d = f isa PDSFunction ? f.d(u, p, t) : nothing + d = f.d(u, p, t) return P, d end