Skip to content

Commit 7d858a9

Browse files
SKopeczSKopecz
committed
basic_patankar_step! with linear system matrix different from P & in-place versions of MPLM22, MPLM33
1 parent 43e72d7 commit 7d858a9

2 files changed

Lines changed: 300 additions & 13 deletions

File tree

‎src/mplm.jl‎

Lines changed: 257 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -62,8 +62,6 @@ function alg_cache(alg::MPLM22, u, rate_prototype, ::Type{uEltypeNoUnits},
6262
::Val{true},
6363
verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
6464
uprevprev = zero(u)
65-
uprevprev[1] = 111.0
66-
uprevprev[2] = 222.0
6765

6866
step = 1
6967
small_constant = alg.small_constant_function(uEltypeNoUnits)
@@ -111,6 +109,7 @@ struct MPLM33{F, T} <: OrdinaryDiffEqAlgorithm
111109
end
112110

113111
alg_order(alg::MPLM33) = 3
112+
isfsal(::MPLM33) = false
114113
alg_extrapolates(alg::MPLM33) = true # TODO: Should probably be false
115114

116115
@cache mutable struct MPLM33oopCache{uType, PType, dType, T, T2} <:
@@ -157,13 +156,90 @@ function alg_cache(alg::MPLM33, u, rate_prototype, ::Type{uEltypeNoUnits},
157156
MPLM33oopCache(u, u, P, P, d, d, αβ, 1, alg.small_constant_function(uEltypeNoUnits))
158157
end
159158

159+
@cache mutable struct MPLM33Cache{uType, dType, T, PType, F, TabType} <: MPLMMutableCache
160+
uprevprev::uType
161+
uprev3::uType
162+
v::uType
163+
vprev::uType
164+
vprev2::uType
165+
step::Int
166+
small_constant::T
167+
b::uType # rhs of the linear system
168+
P::PType
169+
P2::PType
170+
P3::PType
171+
A::PType # system matrix of the linear system
172+
d::dType
173+
d2::dType
174+
d3::dType
175+
d_sys::dType
176+
σ::uType
177+
linsolve::F
178+
αβ::TabType
179+
end
180+
181+
function alg_cache(alg::MPLM33, u, rate_prototype, ::Type{uEltypeNoUnits},
182+
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
183+
uprev, uprev2, f, t, dt, reltol, p, calck,
184+
::Val{true},
185+
verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
186+
uprevprev = zero(u)
187+
uprev3 = zero(u)
188+
v = zero(u)
189+
vprev = zero(u)
190+
vprev2 = zero(u)
191+
step = 1
192+
small_constant = alg.small_constant_function(uEltypeNoUnits)
193+
b = zero(u)
194+
P = p_prototype(u, f)
195+
P2 = p_prototype(u, f)
196+
P3 = p_prototype(u, f)
197+
A = p_prototype(u, f)
198+
σ = zero(u)
199+
200+
# MPLM33 coefficients
201+
α1 = zero(uEltypeNoUnits)
202+
α2 = zero(uEltypeNoUnits)
203+
α3 = one(uEltypeNoUnits)
204+
β1 = 9 / 4 * one(uEltypeNoUnits)
205+
β2 = zero(uEltypeNoUnits)
206+
β3 = 3 / 4 * one(uEltypeNoUnits)
207+
αβ = (α1, α2, α3, β1, β2, β3)
208+
209+
if f isa ConservativePDSFunction
210+
linprob = LinearProblem(A, _vec(b))
211+
linsolve = init(linprob, alg.linsolve,
212+
alias = LinearSolve.LinearAliasSpecifier(; alias_A = true,
213+
alias_b = true),
214+
assumptions = LinearSolve.OperatorAssumptions(true))
215+
216+
MPLM33Cache(uprevprev, uprev3, v, vprev, vprev2, step, small_constant, b, P, P2, P3,
217+
A, nothing, nothing, nothing, nothing, σ,
218+
linsolve, αβ)
219+
elseif f isa PDSFunction
220+
linprob = LinearProblem(A, _vec(b))
221+
linsolve = init(linprob, alg.linsolve,
222+
alias = LinearSolve.LinearAliasSpecifier(; alias_A = true,
223+
alias_b = true),
224+
assumptions = LinearSolve.OperatorAssumptions(true))
225+
226+
MPLM33Cache(uprevprev, uprev3, v, vprev, vprev2, step, small_constant, b, P, P2, P3,
227+
A,
228+
similar(u), similar(u), similar(u), similar(u),
229+
σ, linsolve, αβ)
230+
else
231+
throw(ArgumentError("MPLM33 can only be applied to production-destruction systems"))
232+
end
233+
end
234+
160235
#### MPLM43 ############################################################################
161236
struct MPLM43{F, T} <: OrdinaryDiffEqAlgorithm
162237
linsolve::F
163238
small_constant_function::T
164239
end
165240

166241
alg_order(alg::MPLM43) = 3
242+
isfsal(::MPLM43) = false
167243
alg_extrapolates(alg::MPLM43) = true # TODO: Should probably be false
168244

169245
@cache mutable struct MPLM43oopCache{uType, PType, dType, T, T2} <:
@@ -222,6 +298,7 @@ struct MPLM54{F, T} <: OrdinaryDiffEqAlgorithm
222298
end
223299

224300
alg_order(alg::MPLM54) = 4
301+
isfsal(::MPLM54) = false
225302
alg_extrapolates(alg::MPLM54) = true # TODO: Should probably be false
226303

227304
@cache mutable struct MPLM54oopCache{uType, PType, dType, T, T2} <:
@@ -285,6 +362,7 @@ struct MPLM75{F, T} <: OrdinaryDiffEqAlgorithm
285362
end
286363

287364
alg_order(alg::MPLM75) = 5
365+
isfsal(::MPLM75) = false
288366
alg_extrapolates(alg::MPLM75) = true # TODO: Should probably be false
289367

290368
@cache mutable struct MPLM75oopCache{uType, PType, dType, T, T2} <:
@@ -357,6 +435,7 @@ struct MPLM106{F, T} <: OrdinaryDiffEqAlgorithm
357435
end
358436

359437
alg_order(alg::MPLM106) = 6
438+
isfsal(::MPLM106) = false
360439
alg_extrapolates(alg::MPLM106) = true # TODO: Should probably be false
361440

362441
@cache mutable struct MPLM106oopCache{uType, PType, dType, T, T2} <:
@@ -493,6 +572,7 @@ end
493572
#TODO Use αβ in MPLM22
494573
@muladd function perform_step_MPLM22_oop(P, d, dt, uprev, uprevprev, linsolve,
495574
small_constant)
575+
496576
# First σ approximation
497577
σ = add_small_constant(uprev, small_constant)
498578

@@ -502,7 +582,6 @@ end
502582
σ = add_small_constant(σ, small_constant)
503583

504584
u = basic_patankar_step(uprevprev, P, σ, 2 * dt, linsolve, d)
505-
506585
# statistics: 2 nsolve
507586

508587
return u
@@ -514,18 +593,16 @@ end
514593
@.. broadcast=false σ=uprev + small_constant
515594

516595
# use lincomb! to handle cases in which d2 is nothing
517-
P2 .= P
518-
lincomb!(d2, 1, d)
596+
#P2 .= P
597+
#lincomb!(d2, 1, d)
519598

520-
basic_patankar_step!(σ, uprev, P2, d2, σ, dt, linsolve)
599+
#basic_patankar_step!(σ, uprev, P2, d2, linsolve.A, σ, dt, linsolve)
600+
basic_patankar_step!(σ, uprev, P, d, linsolve.A, σ, dt, linsolve)
521601

522602
# Main step
523603
@.. broadcast=false σ=σ + small_constant
524604

525-
P2 .= P
526-
lincomb!(d2, 1, d)
527-
528-
basic_patankar_step!(u, uprevprev, P2, d2, σ, 2 * dt, linsolve)
605+
basic_patankar_step!(u, uprevprev, P, d, linsolve.A, σ, 2 * dt, linsolve)
529606

530607
# statistics: 2 nsolve
531608

@@ -555,6 +632,9 @@ end
555632
integrator.stats.nf += nf
556633
integrator.stats.nsolve += ns
557634
else
635+
# increase step counter
636+
cache.step += 1
637+
558638
# evaluate production matrix
559639
P, d = evaluate_pds(f, uprev, p, t)
560640
integrator.stats.nf += 1
@@ -593,6 +673,9 @@ end
593673
integrator.stats.nsolve += ns
594674
else
595675

676+
# increase step counter
677+
cache.step += 1
678+
596679
# evaluate production matrix
597680
evaluate_pds!(P, D, f, uprev, p, t)
598681
integrator.stats.nf += 1
@@ -647,6 +730,50 @@ end
647730

648731
return (v, u), t, nf, ns
649732
end
733+
734+
@muladd function start_MPLM33!(v, tmp, P, P2, d, d2, t, dt, vprev, vprev2, σ, f, p,
735+
small_constant, linsolve)
736+
737+
# 1 macro step consists of 4 substeps
738+
dts = dt / 4
739+
740+
### first macro step ###############################################################
741+
# substep 1
742+
nf, ns = start_MPLM22!(v, P, d, t, dts, vprev, σ, f, p, small_constant, linsolve)
743+
744+
# substeps 2 - 4
745+
for _ in 1:3
746+
vprev2 .= vprev
747+
vprev .= v
748+
749+
evaluate_pds!(P, d, f, vprev, p, t)
750+
nf += 1
751+
752+
perform_step_MPLM22!(v, P, P2, d, d2, dts, vprev, vprev2, σ, linsolve,
753+
small_constant)
754+
t += dts
755+
ns += 2
756+
end
757+
758+
tmp .= v
759+
760+
### second macro step ############################################################
761+
for _ in 1:4
762+
vprev2 .= vprev
763+
vprev .= v
764+
765+
evaluate_pds!(P, d, f, vprev, p, t)
766+
nf += 1
767+
768+
perform_step_MPLM22!(v, P, P2, d, d2, dts, vprev, vprev2, σ, linsolve,
769+
small_constant)
770+
t += dts
771+
ns += 2
772+
end
773+
774+
return nf, ns
775+
end
776+
650777
@muladd function perform_step_MPLM33_oop(P_tup, d_tup, dt, u_tup, linsolve, αβ,
651778
small_constant)
652779
P, P2, P3 = P_tup
@@ -668,14 +795,48 @@ end
668795
σ = add_small_constant(σ, small_constant)
669796

670797
Ptmp, dtmp = lincomb(β1, P, d, β2, P2, d2, β3, P3, d3)
798+
671799
v = α1 * uprev + α2 * uprevprev + α3 * uprev3
800+
672801
u = basic_patankar_step(v, Ptmp, σ, dt, linsolve, dtmp)
673802

674803
# statistics: 3 nsolve
675804

676805
return u
677806
end
678807

808+
@muladd function perform_step_MPLM33!(u, P_tup, d_tup, dt, u_tup, σ, linsolve, αβ,
809+
small_constant)
810+
P, P2, P3 = P_tup
811+
d, d2, d3 = d_tup
812+
uprev, uprevprev, uprev3 = u_tup
813+
α1, α2, α3, β1, β2, β3 = αβ
814+
815+
# First σ approximation
816+
@.. broadcast=false σ=uprev + small_constant
817+
818+
basic_patankar_step!(σ, uprev, P, d, linsolve.A, σ, dt, linsolve)
819+
820+
# Second σ approximation
821+
@.. broadcast=false σ=σ + small_constant
822+
823+
basic_patankar_step!(σ, uprevprev, P, d, linsolve.A, σ, 2 * dt, linsolve)
824+
825+
# Main step
826+
@.. broadcast=false σ=σ + small_constant
827+
828+
lincomb!(P3, β1, P, β2, P2, β3, P3)
829+
lincomb!(d3, β1, d, β2, d2, β3, d3)
830+
831+
@.. broadcast=false linsolve.b=α1 * uprev + α2 * uprevprev + α3 * uprev3
832+
833+
basic_patankar_step!(u, linsolve.b, P3, d3, linsolve.A, σ, dt, linsolve)
834+
835+
# statistics: 3 nsolve
836+
837+
return nothing
838+
end
839+
679840
@muladd function perform_step!(integrator, cache::MPLM33oopCache, repeat_step = false)
680841
(; alg, t, dt, uprev, uprev2, f, p) = integrator
681842
(; uprevprev, uprev3, P2, P3, d2, d3, αβ, small_constant) = cache
@@ -720,7 +881,10 @@ end
720881

721882
cache.uprev3 = uprevprev
722883
cache.uprevprev = uprev
884+
723885
else
886+
# increase step count
887+
cache.step += 1
724888

725889
# evaluate production matrix
726890
P, d = evaluate_pds(f, uprev, p, t)
@@ -746,6 +910,89 @@ end
746910
cache.d2 = d
747911
end
748912

913+
@muladd function perform_step!(integrator, cache::MPLM33Cache, repeat_step = false)
914+
(; alg, t, dt, uprev, uprev2, u, f, p) = integrator
915+
(; uprevprev, uprev3, v, vprev, vprev2, P, P2, P3, d, d2, d3, σ, αβ, small_constant, linsolve) = cache
916+
917+
#TODO: is this necessary?
918+
if integrator.u_modified
919+
cache.step = 1
920+
end
921+
922+
if cache.step == 1
923+
# increase step count
924+
cache.step += 1
925+
926+
# initilialze vprev
927+
vprev .= uprev
928+
929+
# evaluate production matrix at tspan[1]
930+
evaluate_pds!(P, d, f, uprev, p, t)
931+
integrator.stats.nf += 1
932+
933+
# save current P
934+
P3 .= P
935+
936+
# compute initial values for MPLM33
937+
938+
# we use uprev3 as temporary storage for the value of u needed in step 1.
939+
nf, ns = start_MPLM33!(v, uprev3, P, P2, d, d2, t, dt, vprev, vprev2, σ, f, p,
940+
small_constant,
941+
linsolve)
942+
integrator.stats.nf += nf
943+
integrator.stats.nsolve += ns
944+
945+
# reset P
946+
P .= P3
947+
948+
# u at time tspan[1] + dt
949+
u .= uprev3
950+
951+
# we use uprev3 as temporary storage for the value of u needed in step 2.
952+
uprev3 .= v
953+
954+
uprevprev .= uprev
955+
956+
elseif cache.step == 2
957+
# increase step count
958+
cache.step += 1
959+
960+
# evaluate production matrix at tspan[1] + dt
961+
evaluate_pds!(P, d, f, uprev, p, t)
962+
integrator.stats.nf += 1
963+
964+
# u at time tspan[1] + 2*dt (this was computed in step 1)
965+
u .= uprev3
966+
967+
uprev3 .= uprevprev
968+
uprevprev .= uprev
969+
970+
else
971+
# increase step count
972+
cache.step += 1
973+
974+
# evaluate production matrix
975+
evaluate_pds!(P, d, f, uprev, p, t)
976+
integrator.stats.nf += 1
977+
978+
P_tup = (P, P2, P3)
979+
d_tup = (d, d2, d3)
980+
u_tup = (uprev, uprevprev, uprev3)
981+
982+
perform_step_MPLM33!(u, P_tup, d_tup, dt, u_tup, σ, linsolve, αβ,
983+
small_constant)
984+
integrator.stats.nsolve += 3
985+
986+
uprev3 .= uprevprev
987+
uprevprev .= uprev
988+
end
989+
990+
P3 .= P2
991+
P2 .= P
992+
!isnothing(d2) && (d3 .= d2)
993+
!isnothing(d) && (d2 .= d)
994+
end
995+
749996
#### MPLM43 ############################################################################
750997
@muladd function start_MPLM43_oop(P, d, t, dt, uprev, f, p, small_constant, linsolve)
751998
αβ33 = (0, 0, 1, 9 / 4, 0, 3 / 4)

0 commit comments

Comments
 (0)