Skip to content

Commit 52fe8ec

Browse files
SKopeczranocha
SKopecz
and
ranocha
authored
additional error handling and tests (#84)
* assert alpha greater or equal 0.5 for MPRK22(alpha) * check that MPRK22 and MPRK43 are applied to a PDS * Fixed bug: MPRK43I instead of MPRK43 * additional check that MPRK43II can only be applied to PDS * Update src/mprk.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Update src/mprk.jl Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * Added get_constant_paraemters for MPRK22. The parameter check is now inside this function and thus valid for in-place and out-of-place implementations. * Replaced @asserts dealing with correct parameters choices for MPRK43 by throw(ArgumentError( )) * format --------- Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
1 parent 81f1bb0 commit 52fe8ec

2 files changed

Lines changed: 71 additions & 15 deletions

File tree

‎src/mprk.jl‎

Lines changed: 56 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -477,6 +477,23 @@ end
477477
alg_order(::MPRK22) = 2
478478
isfsal(::MPRK22) = false
479479

480+
function get_constant_parameters(alg::MPRK22)
481+
if !(alg.alpha ≥ 1 / 2)
482+
throw(ArgumentError("MPRK22 requires α ≥ 1/2."))
483+
end
484+
485+
a21 = alg.alpha
486+
b2 = 1 / (2 * a21)
487+
b1 = 1 - b2
488+
c2 = a21
489+
490+
# This should never happen
491+
if !all((a21, b1, b2, c2) .≥ 0)
492+
throw(ArgumentError("MPRK22 requires nonnegative RK coefficients."))
493+
end
494+
return a21, b1, b2, c2
495+
end
496+
480497
struct MPRK22Cache{uType, rateType, PType, tabType, Thread, F, uNoUnitsType} <:
481498
OrdinaryDiffEqMutableCache
482499
u::uType
@@ -501,8 +518,11 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
501518
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
502519
uprev, uprev2, f, t, dt, reltol, p, calck,
503520
::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
504-
tab = MPRK22ConstantCache(alg.alpha, 1 - 1 / (2 * alg.alpha), 1 / (2 * alg.alpha),
505-
alg.alpha, floatmin(uEltypeNoUnits))
521+
if !(f isa PDSFunction || f isa ConservativePDSFunction)
522+
throw(ArgumentError("MPRK22 can only be applied to production-destruction systems"))
523+
end
524+
a21, b1, b2, c2 = get_constant_parameters(alg)
525+
tab = MPRK22ConstantCache(a21, b1, b2, c2, floatmin(uEltypeNoUnits))
506526

507527
tmp = zero(u)
508528

@@ -542,11 +562,12 @@ function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
542562
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
543563
uprev, uprev2, f, t, dt, reltol, p, calck,
544564
::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
565+
if !(f isa PDSFunction || f isa ConservativePDSFunction)
566+
throw(ArgumentError("MPRK22 can only be applied to production-destruction systems"))
567+
end
545568

546-
#TODO: Should assert alg.alpha >= 0.5
547-
548-
MPRK22ConstantCache(alg.alpha, 1 - 1 / (2 * alg.alpha), 1 / (2 * alg.alpha), alg.alpha,
549-
floatmin(uEltypeNoUnits))
569+
a21, b1, b2, c2 = get_constant_parameters(alg)
570+
MPRK22ConstantCache(a21, b1, b2, c2, floatmin(uEltypeNoUnits))
550571
end
551572

552573
function initialize!(integrator, cache::MPRK22ConstantCache)
@@ -811,14 +832,22 @@ struct MPRK43Cache{uType, rateType, PType, tabType, Thread, F, uNoUnitsType} <:
811832
end
812833

813834
function get_constant_parameters(alg::MPRK43I)
814-
@assert alg.alpha ≥ 1 / 3&&alg.alpha ≠ 2 / 3 "MPRK43I requires α ≥ 1/3 and α ≠ 2/3."
835+
if !(alg.alpha ≥ 1 / 3 && alg.alpha ≠ 2 / 3)
836+
throw(ArgumentError("MPRK43I requires α ≥ 1/3 and α ≠ 2/3."))
837+
end
815838
α0 = 1 / 6 * (3 + (3 - 2 * sqrt(2))^(1 / 3) + (3 + 2 * sqrt(2))^(1 / 3))
816839
if 1 / 3 ≤ alg.alpha < 2 / 3
817-
@assert 2/3≤alg.beta≤3*alg.alpha*(1-alg.alpha) "For this choice of α MPRK43I requires 2/3 ≤ β ≤ 3α(1-α)."
840+
if !(2 / 3 ≤ alg.beta ≤ 3 * alg.alpha * (1 - alg.alpha))
841+
throw(ArgumentError("For this choice of α MPRK43I requires 2/3 ≤ β ≤ 3α(1-α)."))
842+
end
818843
elseif 2 / 3 < alg.alpha ≤ α0
819-
@assert 3*alg.alpha*(1-alg.alpha)≤alg.beta≤2/3 "For this choice of α MPRK43I requires 3α(1-α) ≤ β ≤ 2/3."
844+
if !(3 * alg.alpha * (1 - alg.alpha) ≤ alg.beta ≤ 2 / 3)
845+
throw(ArgumentError("For this choice of α MPRK43I requires 3α(1-α) ≤ β ≤ 2/3."))
846+
end
820847
else
821-
@assert (3 * alg.alpha - 2)/(6 * alg.alpha - 3)≤alg.beta≤2/3 "For this choice of α MPRK43I requires (3α-2)/(6α-3) ≤ β ≤ 2/3."
848+
if !((3 * alg.alpha - 2) / (6 * alg.alpha - 3) ≤ alg.beta ≤ 2 / 3)
849+
throw(ArgumentError("For this choice of α MPRK43I requires (3α-2)/(6α-3) ≤ β ≤ 2/3."))
850+
end
822851
end
823852

824853
a21 = alg.alpha
@@ -837,12 +866,17 @@ function get_constant_parameters(alg::MPRK43I)
837866
q1 = 1 / (3 * a21 * (a31 + a32) * b3)
838867
q2 = 1 / a21
839868

840-
@assert all((a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2) .≥ 0) "MPRK43I requires nonnegative RK coefficients."
869+
#This should never happen
870+
if !all((a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2) .≥ 0)
871+
throw(ArgumentError("MPRK43I requires nonnegative RK coefficients."))
872+
end
841873
return a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2
842874
end
843875

844876
function get_constant_parameters(alg::MPRK43II)
845-
@assert 3/8≤alg.gamma≤3/4 "MPRK43II requires 3/8 ≤ γ ≤ 3/4."
877+
if !(3 / 8 ≤ alg.gamma ≤ 3 / 4)
878+
throw(ArgumentError("MPRK43II requires 3/8 ≤ γ ≤ 3/4."))
879+
end
846880

847881
a21 = 2 * one(alg.gamma) / 3
848882
a31 = a21 - 1 / (4 * alg.gamma)
@@ -859,7 +893,10 @@ function get_constant_parameters(alg::MPRK43II)
859893
q1 = 1 / (3 * a21 * (a31 + a32) * b3)
860894
q2 = 1 / a21
861895

862-
@assert all((a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2) .≥ 0) "MPRK43II requires nonnegative RK coefficients."
896+
#This should never happen
897+
if !all((a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2) .≥ 0)
898+
throw(ArgumentError("MPRK43II requires nonnegative RK coefficients."))
899+
end
863900
return a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2
864901
end
865902

@@ -883,6 +920,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt
883920
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
884921
uprev, uprev2, f, t, dt, reltol, p, calck,
885922
::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
923+
if !(f isa PDSFunction || f isa ConservativePDSFunction)
924+
throw(ArgumentError("MPRK43 can only be applied to production-destruction systems"))
925+
end
886926
a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg)
887927
tab = MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3,
888928
beta1, beta2, q1, q2, floatmin(uEltypeNoUnits))
@@ -1010,6 +1050,9 @@ function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uElt
10101050
::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
10111051
uprev, uprev2, f, t, dt, reltol, p, calck,
10121052
::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
1053+
if !(f isa PDSFunction || f isa ConservativePDSFunction)
1054+
throw(ArgumentError("MPRK43 can only be applied to production-destruction systems"))
1055+
end
10131056
a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg)
10141057
tab = MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3,
10151058
beta1, beta2, q1, q2, floatmin(uEltypeNoUnits))

‎test/runtests.jl‎

Lines changed: 15 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -286,8 +286,21 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [
286286
@test_throws "MPE can only be applied to production-destruction systems" solve(prob_ip,
287287
MPE(),
288288
dt = 0.25)
289-
#TODO: Test that MPRK22 requires α ≥ 1/2 (not yet implemented)
290-
#TODO: Test that MPRK22, MPRK43I, MPRK43II can only be applied to PDS (not yet implemented)
289+
@test_throws "MPRK22 can only be applied to production-destruction systems" solve(prob_oop,
290+
MPRK22(1.0))
291+
@test_throws "MPRK22 can only be applied to production-destruction systems" solve(prob_ip,
292+
MPRK22(1.0))
293+
@test_throws "MPRK22 requires α ≥ 1/2." solve(prob_pds_linmod, MPRK22(0.25))
294+
@test_throws "MPRK43 can only be applied to production-destruction systems" solve(prob_oop,
295+
MPRK43I(1.0,
296+
0.5))
297+
@test_throws "MPRK43 can only be applied to production-destruction systems" solve(prob_ip,
298+
MPRK43I(1.0,
299+
0.5))
300+
@test_throws "MPRK43 can only be applied to production-destruction systems" solve(prob_oop,
301+
MPRK43II(0.5))
302+
@test_throws "MPRK43 can only be applied to production-destruction systems" solve(prob_ip,
303+
MPRK43II(0.5))
291304
@test_throws "MPRK43I requires α ≥ 1/3 and α ≠ 2/3." solve(prob_pds_linmod,
292305
MPRK43I(0.0, 0.5))
293306
@test_throws "MPRK43I requires α ≥ 1/3 and α ≠ 2/3." solve(prob_pds_linmod,

0 commit comments

Comments
 (0)