Skip to content

Commit

Permalink
use LinearSolve.jl for in-place MPRK22 and set default alg to `LUFa…
Browse files Browse the repository at this point in the history 
…ctorization()` (#60)

* use LinearSolve.jl for in-place MPRK22 and set default alg to LUFactorization()

* format

* set version to v0.1.5
  • Loading branch information
@ranocha
ranocha authored Apr 6, 2024
1 parent 6af3b64 commit 95b1be4
Show file tree
Hide file tree
Showing 3 changed files with 87 additions and 36 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "PositiveIntegrators"
uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5"
authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"]
version = "0.1.4"
version = "0.1.5"

[deps]
FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
Expand Down
2 changes: 1 addition & 1 deletion src/PositiveIntegrators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ using OrdinaryDiffEq: OrdinaryDiffEq, OrdinaryDiffEqAlgorithm

using SymbolicIndexingInterface

using LinearSolve: LinearSolve, LinearProblem
using LinearSolve: LinearSolve, LinearProblem, LUFactorization

import SciMLBase: __has_mass_matrix, __has_analytic, __has_tgrad,
__has_jac, __has_jvp, __has_vjp, __has_jac_prototype,
Expand Down
119 changes: 85 additions & 34 deletions src/mprk.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,15 +178,20 @@ end

### MPE #####################################################################################
"""
MPE()
MPE([linsolve = ...])
The first-order modified Patankar-Euler algorithm for conservative production-destruction
systems. This one-step, one-stage method is first-order accurate, unconditionally
positivity-preserving, and linearly implicit.
The first-order modified Patankar-Euler algorithm for (conservative)
production-destruction systems. This one-step, one-stage method is
first-order accurate, unconditionally positivity-preserving, and
linearly implicit.
The modified Patankar-Euler method requires the special structure of a
[`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref).
You can optionally choose the linear solver to be used by passing an
algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
as keyword argument `linsolve`.
## References
- Hans Burchard, Eric Deleersnijder, and Andreas Meister.
Expand All @@ -199,22 +204,25 @@ struct MPE{F} <: OrdinaryDiffEqAlgorithm
linsolve::F
end

function MPE(; linsolve = nothing)
function MPE(; linsolve = LUFactorization())
MPE(linsolve)
end

# TODO: Consider switching to the interface of LinearSolve.jl directly,
# avoiding `dolinesolve` from OrdinaryDiffEq.jl.
# TODO: Think about adding preconditioners to the MPE algorithm
# This hack is currently required to make OrdinaryDiffEq.jl happy...
function Base.getproperty(mpe::MPE, f::Symbol)
function Base.getproperty(alg::MPE, f::Symbol)
# preconditioners
if f === :precs
return Returns((nothing, nothing))
else
return getfield(mpe, f)
return getfield(alg, f)
end
end

#@cache
OrdinaryDiffEq.alg_order(alg::MPE) = 1

struct MPECache{uType, rateType, PType, F, uNoUnitsType} <: OrdinaryDiffEqMutableCache
u::uType
uprev::uType
Expand All @@ -229,8 +237,8 @@ struct MPECache{uType, rateType, PType, F, uNoUnitsType} <: OrdinaryDiffEqMutabl
end

function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
dt, reltol, p, calck,
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
uprev, uprev2, f, t, dt, reltol, p, calck,
::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
tmp = zero(u)
P = p_prototype(u, f)
Expand All @@ -256,8 +264,8 @@ struct MPEConstantCache{T} <: OrdinaryDiffEqConstantCache
end

function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
dt, reltol, p, calck,
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
uprev, uprev2, f, t, dt, reltol, p, calck,
::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
MPEConstantCache(floatmin(uEltypeNoUnits))
end
Expand Down Expand Up @@ -325,9 +333,11 @@ function perform_step!(integrator, cache::MPECache, repeat_step = false)
f.p(P, uprev, p, t) # evaluate production terms
sum_destruction_terms!(D, P) # store destruction terms in D
build_mprk_matrix!(P, 1, P, D, uprev, dt)
#linres = P\uprev # TODO: needs to be implemented without allocations
linres = dolinsolve(integrator, cache.linsolve; A = P, b = _vec(uprev),
du = integrator.fsalfirst, u = u, p = p, t = t, weight = weight)
# linres = P\uprev # TODO: needs to be implemented without allocations
linres = dolinsolve(integrator, cache.linsolve;
A = P, b = _vec(uprev),
du = integrator.fsalfirst, u = u, p = p, t = t,
weight = weight)

u .= linres

Expand All @@ -337,17 +347,22 @@ end

### MPRK #####################################################################################
"""
MPRK22(α)
MPRK22(α; [linsolve = ...])
The second-order modified Patankar-Runge-Kutta algorithm for conservative production-destruction
systems. This one-step, two-stage method is second-order accurate, unconditionally
positivity-preserving, and linearly implicit. The parameter `α` is described by
Kopecz and Meister (2018) and studied by Izgin, Kopecz and Meister (2022) as well as
The second-order modified Patankar-Runge-Kutta algorithm for (conservative)
production-destruction systems. This one-step, two-stage method is
second-order accurate, unconditionally positivity-preserving, and linearly
implicit. The parameter `α` is described by Kopecz and Meister (2018) and
studied by Izgin, Kopecz and Meister (2022) as well as
Torlo, Öffner and Ranocha (2022).
This modified Patankar-Runge-Kutta method requires the special structure of a
[`PDSProblem`](@ref) or a [`ConservativePDSProblem`](@ref).
You can optionally choose the linear solver to be used by passing an
algorithm from [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl)
as keyword argument `linsolve`.
## References
- Hans Burchard, Eric Deleersnijder, and Andreas Meister.
Expand Down Expand Up @@ -375,14 +390,27 @@ struct MPRK22{T, Thread, F} <: OrdinaryDiffEqAdaptiveAlgorithm
linsolve::F
end

function MPRK22(alpha; thread = False(), linsolve = nothing)
function MPRK22(alpha; thread = False(), linsolve = LUFactorization())
MPRK22{typeof(alpha), typeof(thread), typeof(linsolve)}(alpha, thread, linsolve)
end

# TODO: Consider switching to the interface of LinearSolve.jl directly,
# avoiding `dolinesolve` from OrdinaryDiffEq.jl.
# TODO: Think about adding preconditioners to the MPRK22 algorithm
# This hack is currently required to make OrdinaryDiffEq.jl happy...
function Base.getproperty(alg::MPRK22, f::Symbol)
# preconditioners
if f === :precs
return Returns((nothing, nothing))
else
return getfield(alg, f)
end
end

OrdinaryDiffEq.alg_order(alg::MPRK22) = 2

#@cache
struct MPRK22Cache{uType, rateType, PType, tabType, Thread} <: OrdinaryDiffEqMutableCache
struct MPRK22Cache{uType, rateType, PType, tabType, Thread, F, uNoUnitsType} <:
OrdinaryDiffEqMutableCache
u::uType
uprev::uType
tmp::uType
Expand All @@ -397,6 +425,9 @@ struct MPRK22Cache{uType, rateType, PType, tabType, Thread} <: OrdinaryDiffEqMut
σ::uType
tab::tabType
thread::Thread
linsolve_tmp::uType # stores rhs of linear system
linsolve::F
weight::uNoUnitsType
end

function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
Expand All @@ -405,18 +436,30 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
tab = MPRK22ConstantCache(alg.alpha, 1 - 1 / (2 * alg.alpha), 1 / (2 * alg.alpha),
alg.alpha, floatmin(uEltypeNoUnits))
MPRK22Cache(u, uprev,
zero(u), # tmp

tmp = zero(u)

M = p_prototype(u, f)
linsolve_tmp = zero(u)
weight = similar(u, uEltypeNoUnits)
recursivefill!(weight, false)

linprob = LinearProblem(M, _vec(linsolve_tmp); u0 = _vec(tmp))
linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
assumptions = LinearSolve.OperatorAssumptions(true))

MPRK22Cache(u, uprev, tmp,
zero(u), # atmp
zero(rate_prototype), # k
zero(rate_prototype), #fsalfirst
p_prototype(u, f), # P
p_prototype(u, f), # P2
zero(u), # D
zero(u), # D2
p_prototype(u, f), # M
M,
zero(u), # σ
tab, alg.thread)
tab, alg.thread,
linsolve_tmp, linsolve, weight)
end

struct MPRK22ConstantCache{T} <: OrdinaryDiffEqConstantCache
Expand All @@ -428,8 +471,8 @@ struct MPRK22ConstantCache{T} <: OrdinaryDiffEqConstantCache
end

function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
dt, reltol, p, calck,
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
uprev, uprev2, f, t, dt, reltol, p, calck,
::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}

#TODO: Should assert alg.alpha >= 0.5
Expand Down Expand Up @@ -527,16 +570,20 @@ end

function perform_step!(integrator, cache::MPRK22Cache, repeat_step = false)
@unpack t, dt, uprev, u, f, p = integrator
@unpack tmp, atmp, P, P2, D, D2, M, σ, thread = cache
@unpack tmp, atmp, P, P2, D, D2, M, σ, thread, weight = cache
@unpack a21, b1, b2, c2, small_constant = cache.tab

uprev .= uprev .+ small_constant

f.p(P, uprev, p, t) # evaluate production terms
sum_destruction_terms!(D, P) # store destruction terms in D
build_mprk_matrix!(M, a21, P, D, uprev, dt)
tmp = M \ uprev #TODO: needs to be implemented without allocations.
u .= tmp
# linres = M \ uprev #TODO: needs to be implemented without allocations.
linres = dolinsolve(integrator, cache.linsolve;
A = M, b = _vec(uprev),
du = integrator.fsalfirst, u = u, p = p, t = t,
weight = weight)
u .= linres

u .= u .+ small_constant

Expand All @@ -545,8 +592,12 @@ function perform_step!(integrator, cache::MPRK22Cache, repeat_step = false)
f.p(P2, u, p, t + a21 * dt) # evaluate production terms
sum_destruction_terms!(D2, P2) # store destruction terms in D2
build_mprk_matrix!(M, b1, P, D, b2, P2, D2, σ, dt)
tmp = M \ uprev #TODO: needs to be implemented without allocations.
u .= tmp
# linres = M \ uprev #TODO: needs to be implemented without allocations.
linres = dolinsolve(integrator, cache.linsolve;
A = M, b = _vec(uprev),
du = integrator.fsalfirst, u = u, p = p, t = t,
weight = weight)
u .= linres

tmp .= u .- σ
calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol,
Expand Down

2 comments on commit 95b1be4

@ranocha
Copy link
Collaborator Author

@ranocha ranocha commented on 95b1be4 Apr 6, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/104380

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.1.5 -m "<description of version>" 95b1be42ddf833060e1e105e8ec0bf4fc3edca6f
git push origin v0.1.5

Please sign in to comment.