Skip to content

Commit f86c263

Browse files
SKopeczranocha
SKopecz
and
ranocha
authored
Remove d_prototype & create type tests for p_prototype (#104)
* Removed d_prototype * bug fix: removed d_prototype keywords in tests * Update test/runtests.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Set version to v0.2.0 * Updated News --------- Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
1 parent 69db9ab commit f86c263

File tree

6 files changed

+97
-31
lines changed

6 files changed

+97
-31
lines changed

‎NEWS.md

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,11 @@ PositiveIntegrators.jl.jl follows the interpretation of
55
used in the Julia ecosystem. Notable changes will be documented in this file
66
for human readability.
77

8+
## Changes when updating to v0.2 from v0.1.x
9+
10+
#### Removed
11+
12+
- The optional keyword argument `d_prototype` has been removed from `PDSProblem`
813

914
## Changes in the v0.1 lifecycle
1015

@@ -21,3 +26,5 @@ for human readability.
2126
`prob_pds_linmod`, `prob_pds_nonlinmod`, `prob_pds_npzd`, `prob_pds_robertson`,
2227
`prob_pds_sir`, `prob_pds_stratreac`
2328
- Modified Patankar methods `MPE` and `MPRK22`
29+
30+

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "PositiveIntegrators"
22
uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5"
33
authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"]
4-
version = "0.1.16"
4+
version = "0.2.0"
55

66
[deps]
77
FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"

‎src/mprk.jl

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -363,7 +363,7 @@ function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
363363
linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
364364
assumptions = LinearSolve.OperatorAssumptions(true))
365365

366-
MPECache(P, zero(u), σ, tab, linsolve_rhs, linsolve)
366+
MPECache(P, similar(u), σ, tab, linsolve_rhs, linsolve)
367367
else
368368
throw(ArgumentError("MPE can only be applied to production-destruction systems"))
369369
end
@@ -670,8 +670,8 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
670670
assumptions = LinearSolve.OperatorAssumptions(true))
671671

672672
MPRK22Cache(tmp, P, P2,
673-
zero(u), # D
674-
zero(u), # D2
673+
similar(u), # D
674+
similar(u), # D2
675675
σ,
676676
tab, #MPRK22ConstantCache
677677
linsolve)
@@ -1246,9 +1246,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt
12461246
assumptions = LinearSolve.OperatorAssumptions(true))
12471247
MPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, tab, linsolve)
12481248
elseif f isa PDSFunction
1249-
D = zero(u)
1250-
D2 = zero(u)
1251-
D3 = zero(u)
1249+
D = similar(u)
1250+
D2 = similar(u)
1251+
D3 = similar(u)
12521252

12531253
linprob = LinearProblem(P3, _vec(tmp))
12541254
linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,

‎src/proddest.jl

Lines changed: 7 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -3,9 +3,8 @@ abstract type AbstractPDSProblem end
33

44
"""
55
PDSProblem(P, D, u0, tspan, p = NullParameters();
6-
p_prototype = nothing,
7-
d_prototype = nothing,
8-
analytic = nothing)
6+
p_prototype = nothing,
7+
analytic = nothing)
98
109
A structure describing a system of ordinary differential equations in form of a production-destruction system (PDS).
1110
`P` denotes the production matrix.
@@ -20,13 +19,9 @@ The functions `P` and `D` can be used either in the out-of-place form with signa
2019
2120
### Keyword arguments: ###
2221
23-
- `p_prototype`: If `P` is given in in-place form, `p_prototype` is used to store evaluations of `P`.
22+
- `p_prototype`: If `P` is given in in-place form, `p_prototype` or copies thereof are used to store evaluations of `P`.
2423
If `p_prototype` is not specified explicitly and `P` is in-place, then `p_prototype` will be internally
2524
set to `zeros(eltype(u0), (length(u0), length(u0)))`.
26-
- `d_prototype`: If `D` is given in in-place form, `d_prototype` is used to store evaluations of `D`.
27-
If `d_prototype` is not specified explicitly and `D` is in-place, then `d_prototype` will be internally
28-
set to `zeros(eltype(u0), (length(u0),))`.
29-
3025
- `analytic`: The analytic solution of a PDS must be given in the form `f(u0,p,t)`.
3126
Specifying the analytic solution can be useful for plotting and convergence tests.
3227
@@ -86,18 +81,16 @@ end
8681
# (arbitrary functions)
8782
function PDSProblem{iip}(P, D, u0, tspan, p = NullParameters();
8883
p_prototype = nothing,
89-
d_prototype = nothing,
9084
analytic = nothing,
9185
kwargs...) where {iip}
9286

9387
# p_prototype is used to store evaluations of P, if P is in-place.
9488
if isnothing(p_prototype) && iip
9589
p_prototype = zeros(eltype(u0), (length(u0), length(u0)))
9690
end
97-
# d_prototype is used to store evaluations of D, if D is in-place.
98-
if isnothing(d_prototype) && iip
99-
d_prototype = zeros(eltype(u0), (length(u0),))
100-
end
91+
# If a PDSFunction is to be evaluated and D is in-place, then d_prototype is used to store
92+
# evaluations of D.
93+
d_prototype = similar(u0)
10194

10295
PD = PDSFunction{iip}(P, D; p_prototype = p_prototype, d_prototype = d_prototype,
10396
analytic = analytic)
@@ -177,7 +170,7 @@ The function `P` can be given either in the out-of-place form with signature
177170
178171
### Keyword arguments: ###
179172
180-
- `p_prototype`: If `P` is given in in-place form, `p_prototype` is used to store evaluations of `P`.
173+
- `p_prototype`: If `P` is given in in-place form, `p_prototype` or copies thereof are used to store evaluations of `P`.
181174
If `p_prototype` is not specified explicitly and `P` is in-place, then `p_prototype` will be internally
182175
set to `zeros(eltype(u0), (length(u0), length(u0)))`.
183176
- `analytic`: The analytic solution of a PDS must be given in the form `f(u0,p,t)`.

‎src/sspmprk.jl

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -252,8 +252,8 @@ function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
252252
assumptions = LinearSolve.OperatorAssumptions(true))
253253

254254
SSPMPRK22Cache(tmp, P, P2,
255-
zero(u), # D
256-
zero(u), # D2
255+
similar(u), # D
256+
similar(u), # D2
257257
σ,
258258
tab, #MPRK22ConstantCache
259259
linsolve)
@@ -760,9 +760,9 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
760760
assumptions = LinearSolve.OperatorAssumptions(true))
761761
SSPMPRK43ConservativeCache(tmp, tmp2, P, P2, P3, σ, ρ, tab, linsolve)
762762
elseif f isa PDSFunction
763-
D = zero(u)
764-
D2 = zero(u)
765-
D3 = zero(u)
763+
D = similar(u)
764+
D2 = similar(u)
765+
D3 = similar(u)
766766

767767
linprob = LinearProblem(P3, _vec(tmp))
768768
linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,

‎test/runtests.jl

Lines changed: 71 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -495,17 +495,14 @@ end
495495
# problem with sparse matrices
496496
p_prototype = spdiagm(-1 => ones(eltype(u0), N - 1),
497497
N - 1 => ones(eltype(u0), 1))
498-
d_prototype = zero(u0)
499498
linear_advection_fd_upwind_PDS_sparse = PDSProblem(linear_advection_fd_upwind_P!,
500499
linear_advection_fd_upwind_D!,
501500
u0, tspan;
502-
p_prototype = p_prototype,
503-
d_prototype = d_prototype)
501+
p_prototype = p_prototype)
504502
linear_advection_fd_upwind_PDS_sparse_2 = PDSProblem{true}(linear_advection_fd_upwind_P!,
505503
linear_advection_fd_upwind_D!,
506504
u0, tspan;
507-
p_prototype = p_prototype,
508-
d_prototype = d_prototype)
505+
p_prototype = p_prototype)
509506
linear_advection_fd_upwind_ConsPDS_sparse = ConservativePDSProblem(linear_advection_fd_upwind_P!,
510507
u0, tspan;
511508
p_prototype = p_prototype)
@@ -1201,6 +1198,75 @@ end
12011198
end
12021199
end
12031200

1201+
# Here we check that the type of p_prototype actually
1202+
# defines the types of the Ps inside the algorithm caches.
1203+
# We test sparse, tridiagonal, and dense matrices.
1204+
@testset "Prototype type check" begin
1205+
#prod and dest functions
1206+
prod_inner! = (P, u, p, t) -> begin
1207+
fill!(P, zero(eltype(P)))
1208+
for i in 1:(length(u) - 1)
1209+
P[i, i + 1] = i * u[i]
1210+
end
1211+
return nothing
1212+
end
1213+
prod_sparse! = (P, u, p, t) -> begin
1214+
@test P isa SparseMatrixCSC
1215+
prod_inner!(P, u, p, t)
1216+
return nothing
1217+
end
1218+
prod_tridiagonal! = (P, u, p, t) -> begin
1219+
@test P isa Tridiagonal
1220+
prod_inner!(P, u, p, t)
1221+
return nothing
1222+
end
1223+
prod_dense! = (P, u, p, t) -> begin
1224+
@test P isa Matrix
1225+
prod_inner!(P, u, p, t)
1226+
return nothing
1227+
end
1228+
dest! = (D, u, p, t) -> begin
1229+
fill!(D, zero(eltype(D)))
1230+
end
1231+
#prototypes
1232+
P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3],
1233+
[0.0, 0.0, 0.0, 0.0],
1234+
[0.4, 0.5, 0.6])
1235+
P_dense = Matrix(P_tridiagonal)
1236+
P_sparse = sparse(P_tridiagonal)
1237+
# problem definition
1238+
u0 = [1.0, 1.5, 2.0, 2.5]
1239+
tspan = (0.0, 1.0)
1240+
dt = 0.5
1241+
## conservative PDS
1242+
prob_default = ConservativePDSProblem(prod_dense!, u0, tspan)
1243+
prob_tridiagonal = ConservativePDSProblem(prod_tridiagonal!, u0, tspan;
1244+
p_prototype = P_tridiagonal)
1245+
prob_dense = ConservativePDSProblem(prod_dense!, u0, tspan;
1246+
p_prototype = P_dense)
1247+
prob_sparse = ConservativePDSProblem(prod_sparse!, u0, tspan;
1248+
p_prototype = P_sparse)
1249+
## nonconservative PDS
1250+
prob_default2 = PDSProblem(prod_dense!, dest!, u0, tspan)
1251+
prob_tridiagonal2 = PDSProblem(prod_tridiagonal!, dest!, u0, tspan;
1252+
p_prototype = P_tridiagonal)
1253+
prob_dense2 = PDSProblem(prod_dense!, dest!, u0, tspan;
1254+
p_prototype = P_dense)
1255+
prob_sparse2 = PDSProblem(prod_sparse!, dest!, u0, tspan;
1256+
p_prototype = P_sparse)
1257+
#solve and test
1258+
for alg in (MPE(), MPRK22(0.5), MPRK22(1.0), MPRK43I(1.0, 0.5),
1259+
MPRK43I(0.5, 0.75),
1260+
MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK22(0.5, 1.0),
1261+
SSPMPRK43())
1262+
for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse,
1263+
prob_default2,
1264+
prob_tridiagonal2, prob_dense2, prob_sparse2)
1265+
solve(prob, alg; dt, adaptive = false)
1266+
end
1267+
end
1268+
end
1269+
12041270
# Here we check the convergence order of pth-order schemes for which
12051271
# also an interpolation of order p is available
12061272
@testset "Convergence tests (conservative)" begin

0 commit comments

Comments
 (0)