Skip to content

Commit 10352b4

Browse files
SKopeczSKopecz
committed
We now test 'Different matrix types' with constant step size and adaptive time stepping. For adaptive time stepping isapprox is used with a larger relative tolerance.
1 parent 86a791b commit 10352b4

1 file changed

Lines changed: 225 additions & 12 deletions

File tree

‎test/runtests.jl‎

Lines changed: 225 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -552,12 +552,14 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [
552552
prob_sparse_op = ConservativePDSProblem(prod, u0, tspan;
553553
p_prototype = P_sparse)
554554

555-
sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt)
556-
sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt)
557-
sol_dense_ip = solve(prob_dense_ip, alg; dt)
558-
sol_dense_op = solve(prob_dense_op, alg; dt)
559-
sol_sparse_ip = solve(prob_sparse_ip, alg; dt)
560-
sol_sparse_op = solve(prob_sparse_op, alg; dt)
555+
sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt,
556+
adaptive = false)
557+
sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt,
558+
adaptive = false)
559+
sol_dense_ip = solve(prob_dense_ip, alg; dt, adaptive = false)
560+
sol_dense_op = solve(prob_dense_op, alg; dt, adaptive = false)
561+
sol_sparse_ip = solve(prob_sparse_ip, alg; dt, adaptive = false)
562+
sol_sparse_op = solve(prob_sparse_op, alg; dt, adaptive = false)
561563

562564
@test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t
563565
@test sol_dense_ip.t ≈ sol_dense_op.t
@@ -574,6 +576,89 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [
574576
end
575577
end
576578

579+
@testset "Different matrix types (conservative, adaptive)" begin
580+
prod_1! = (P, u, p, t) -> begin
581+
fill!(P, zero(eltype(P)))
582+
for i in 1:(length(u) - 1)
583+
P[i, i + 1] = i * u[i]
584+
end
585+
return nothing
586+
end
587+
588+
prod_2! = (P, u, p, t) -> begin
589+
fill!(P, zero(eltype(P)))
590+
for i in 1:(length(u) - 1)
591+
P[i + 1, i] = i * u[i + 1]
592+
end
593+
return nothing
594+
end
595+
596+
prod_3! = (P, u, p, t) -> begin
597+
fill!(P, zero(eltype(P)))
598+
for i in 1:(length(u) - 1)
599+
P[i, i + 1] = i * u[i]
600+
P[i + 1, i] = i * u[i + 1]
601+
end
602+
return nothing
603+
end
604+
605+
n = 4
606+
P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3],
607+
zeros(n),
608+
[0.4, 0.5, 0.6])
609+
P_dense = Matrix(P_tridiagonal)
610+
P_sparse = sparse(P_tridiagonal)
611+
u0 = [1.0, 1.5, 2.0, 2.5]
612+
tspan = (0.0, 1.0)
613+
dt = 0.25
614+
615+
rtol = sqrt(eps(Float32))
616+
617+
@testset "$alg" for alg in (MPE(),
618+
MPRK22(0.5), MPRK22(1.0),
619+
MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75),
620+
MPRK43II(2.0 / 3.0), MPRK43II(0.5), SSPMPRK43())
621+
for prod! in (prod_1!, prod_2!, prod_3!)
622+
prod = (u, p, t) -> begin
623+
P = similar(u, (length(u), length(u)))
624+
prod!(P, u, p, t)
625+
return P
626+
end
627+
prob_tridiagonal_ip = ConservativePDSProblem(prod!, u0, tspan;
628+
p_prototype = P_tridiagonal)
629+
prob_tridiagonal_op = ConservativePDSProblem(prod, u0, tspan;
630+
p_prototype = P_tridiagonal)
631+
prob_dense_ip = ConservativePDSProblem(prod!, u0, tspan;
632+
p_prototype = P_dense)
633+
prob_dense_op = ConservativePDSProblem(prod, u0, tspan;
634+
p_prototype = P_dense)
635+
prob_sparse_ip = ConservativePDSProblem(prod!, u0, tspan;
636+
p_prototype = P_sparse)
637+
prob_sparse_op = ConservativePDSProblem(prod, u0, tspan;
638+
p_prototype = P_sparse)
639+
640+
sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg; dt)
641+
sol_tridiagonal_op = solve(prob_tridiagonal_op, alg; dt)
642+
sol_dense_ip = solve(prob_dense_ip, alg; dt)
643+
sol_dense_op = solve(prob_dense_op, alg; dt)
644+
sol_sparse_ip = solve(prob_sparse_ip, alg; dt)
645+
sol_sparse_op = solve(prob_sparse_op, alg; dt)
646+
647+
@test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol)
648+
@test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol)
649+
@test isapprox(sol_sparse_ip.t, sol_sparse_op.t; rtol)
650+
@test isapprox(sol_tridiagonal_ip.t, sol_dense_ip.t; rtol)
651+
@test isapprox(sol_tridiagonal_ip.t, sol_sparse_ip.t; rtol)
652+
653+
@test isapprox(sol_tridiagonal_ip.u, sol_tridiagonal_op.u; rtol)
654+
@test isapprox(sol_dense_ip.u, sol_dense_op.u; rtol)
655+
@test isapprox(sol_sparse_ip.u, sol_sparse_op.u; rtol)
656+
@test isapprox(sol_tridiagonal_ip.u, sol_dense_ip.u; rtol)
657+
@test isapprox(sol_tridiagonal_ip.u, sol_sparse_ip.u; rtol)
658+
end
659+
end
660+
end
661+
577662
# Here we check that in-place and out-of-place implementations
578663
# deliver the same results
579664
@testset "Different matrix types (nonconservative)" begin
@@ -674,17 +759,17 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [
674759
p_prototype = P_sparse)
675760

676761
sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg;
677-
dt)
762+
dt, adaptive = false)
678763
sol_tridiagonal_op = solve(prob_tridiagonal_op, alg;
679-
dt)
764+
dt, adaptive = false)
680765
sol_dense_ip = solve(prob_dense_ip, alg;
681-
dt)
766+
dt, adaptive = false)
682767
sol_dense_op = solve(prob_dense_op, alg;
683-
dt)
768+
dt, adaptive = false)
684769
sol_sparse_ip = solve(prob_sparse_ip, alg;
685-
dt)
770+
dt, adaptive = false)
686771
sol_sparse_op = solve(prob_sparse_op, alg;
687-
dt)
772+
dt, adaptive = false)
688773

689774
@test sol_tridiagonal_ip.t ≈ sol_tridiagonal_op.t
690775
@test sol_dense_ip.t ≈ sol_dense_op.t
@@ -701,6 +786,134 @@ const prob_pds_linmod_nonconservative_inplace = PDSProblem(linmodP!, linmodD!, [
701786
end
702787
end
703788

789+
@testset "Different matrix types (nonconservative, adaptive)" begin
790+
prod_1! = (P, u, p, t) -> begin
791+
fill!(P, zero(eltype(P)))
792+
for i in 1:(length(u) - 1)
793+
P[i, i + 1] = i * u[i]
794+
end
795+
for i in 1:length(u)
796+
P[i, i] = i * u[i]
797+
end
798+
return nothing
799+
end
800+
dest_1! = (D, u, p, t) -> begin
801+
fill!(D, zero(eltype(D)))
802+
for i in 1:length(u)
803+
D[i] = (i + 1) * u[i]
804+
end
805+
return nothing
806+
end
807+
808+
prod_2! = (P, u, p, t) -> begin
809+
fill!(P, zero(eltype(P)))
810+
for i in 1:(length(u) - 1)
811+
P[i + 1, i] = i * u[i + 1]
812+
end
813+
for i in 1:length(u)
814+
P[i, i] = (i - 1) * u[i]
815+
end
816+
return nothing
817+
end
818+
dest_2! = (D, u, p, t) -> begin
819+
fill!(D, zero(eltype(D)))
820+
for i in 1:length(u)
821+
D[i] = i * u[i]
822+
end
823+
return nothing
824+
end
825+
826+
prod_3! = (P, u, p, t) -> begin
827+
fill!(P, zero(eltype(P)))
828+
for i in 1:(length(u) - 1)
829+
P[i, i + 1] = i * u[i]
830+
P[i + 1, i] = i * u[i + 1]
831+
end
832+
for i in 1:length(u)
833+
P[i, i] = (i + 1) * u[i]
834+
end
835+
return nothing
836+
end
837+
dest_3! = (D, u, p, t) -> begin
838+
fill!(D, zero(eltype(D)))
839+
for i in 1:length(u)
840+
D[i] = (i - 1) * u[i]
841+
end
842+
return nothing
843+
end
844+
845+
n = 4
846+
P_tridiagonal = Tridiagonal([0.1, 0.2, 0.3],
847+
zeros(n),
848+
[0.4, 0.5, 0.6])
849+
P_dense = Matrix(P_tridiagonal)
850+
P_sparse = sparse(P_tridiagonal)
851+
u0 = [1.0, 1.5, 2.0, 2.5]
852+
D = u0
853+
tspan = (0.0, 1.0)
854+
dt = 0.25
855+
856+
rtol = sqrt(eps(Float32))
857+
@testset "$alg" for alg in (MPE(),
858+
MPRK22(0.5), MPRK22(1.0),
859+
MPRK43I(1.0, 0.5), MPRK43I(0.5, 0.75),
860+
MPRK43II(2.0 / 3.0), MPRK43II(0.5),
861+
SSPMPRK22(0.5, 1.0), SSPMPRK43())
862+
for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!),
863+
(dest_1!, dest_2!, dest_3!))
864+
prod! = prod_3!
865+
dest! = dest_3!
866+
prod = (u, p, t) -> begin
867+
P = similar(u, (length(u), length(u)))
868+
prod!(P, u, p, t)
869+
return P
870+
end
871+
dest = (u, p, t) -> begin
872+
D = similar(u)
873+
dest!(D, u, p, t)
874+
return D
875+
end
876+
prob_tridiagonal_ip = PDSProblem(prod!, dest!, u0, tspan;
877+
p_prototype = P_tridiagonal)
878+
prob_tridiagonal_op = PDSProblem(prod, dest, u0, tspan;
879+
p_prototype = P_tridiagonal)
880+
prob_dense_ip = PDSProblem(prod!, dest!, u0, tspan;
881+
p_prototype = P_dense)
882+
prob_dense_op = PDSProblem(prod, dest, u0, tspan;
883+
p_prototype = P_dense)
884+
prob_sparse_ip = PDSProblem(prod!, dest!, u0, tspan;
885+
p_prototype = P_sparse)
886+
prob_sparse_op = PDSProblem(prod, dest, u0, tspan;
887+
p_prototype = P_sparse)
888+
889+
sol_tridiagonal_ip = solve(prob_tridiagonal_ip, alg;
890+
dt)
891+
sol_tridiagonal_op = solve(prob_tridiagonal_op, alg;
892+
dt)
893+
sol_dense_ip = solve(prob_dense_ip, alg;
894+
dt)
895+
sol_dense_op = solve(prob_dense_op, alg;
896+
dt)
897+
sol_sparse_ip = solve(prob_sparse_ip, alg;
898+
dt)
899+
sol_sparse_op = solve(prob_sparse_op, alg;
900+
dt)
901+
902+
@test isapprox(sol_tridiagonal_ip.t, sol_tridiagonal_op.t; rtol)
903+
@test isapprox(sol_dense_ip.t, sol_dense_op.t; rtol)
904+
@test isapprox(sol_sparse_ip.t, sol_sparse_op.t; rtol)
905+
@test isapprox(sol_tridiagonal_ip.t, sol_dense_ip.t; rtol)
906+
@test isapprox(sol_tridiagonal_ip.t, sol_sparse_ip.t; rtol)
907+
908+
@test isapprox(sol_tridiagonal_ip.u, sol_tridiagonal_op.u; rtol)
909+
@test isapprox(sol_dense_ip.u, sol_dense_op.u; rtol)
910+
@test isapprox(sol_sparse_ip.u, sol_sparse_op.u; rtol)
911+
@test isapprox(sol_tridiagonal_ip.u, sol_dense_ip.u; rtol)
912+
@test isapprox(sol_tridiagonal_ip.u, sol_sparse_ip.u; rtol)
913+
end
914+
end
915+
end
916+
704917
# Here we check the convergence order of pth-order schemes for which
705918
# also an interpolation of order p is available
706919
@testset "Convergence tests (conservative)" begin

0 commit comments

Comments
 (0)