|
1 | 1 | using BoundaryValueDiffEqFIRK |
2 | 2 | using Test |
3 | 3 |
|
4 | | -nested = false |
5 | | - |
6 | | -for stage in (2, 3, 4, 5) |
7 | | - s = Symbol("LobattoIIIa$(stage)") |
8 | | - @eval lobattoIIIa_solver(::Val{$stage}, args...; kwargs...) = $(s)(args...; kwargs...) |
9 | | -end |
10 | | - |
11 | | -for stage in (2, 3, 4, 5) |
12 | | - s = Symbol("LobattoIIIb$(stage)") |
13 | | - @eval lobattoIIIb_solver(::Val{$stage}, args...; kwargs...) = $(s)(args...; kwargs...) |
14 | | -end |
15 | | - |
16 | | -for stage in (2, 3, 4, 5) |
17 | | - s = Symbol("LobattoIIIc$(stage)") |
18 | | - @eval lobattoIIIc_solver(::Val{$stage}, args...; kwargs...) = $(s)(args...; kwargs...) |
19 | | -end |
20 | | - |
21 | | -for stage in (2, 3, 5, 7) |
22 | | - s = Symbol("RadauIIa$(stage)") |
23 | | - @eval radau_solver(::Val{$stage}, args...; kwargs...) = $(s)(args...; kwargs...) |
24 | | -end |
25 | | - |
26 | | -# First order test |
27 | | -function f1!(du, u, p, t) |
28 | | - du[1] = u[2] |
29 | | - return du[2] = 0 |
30 | | -end |
31 | | -f1(u, p, t) = [u[2], 0] |
32 | | - |
33 | | -# Second order linear test |
34 | | -function f2!(du, u, p, t) |
35 | | - du[1] = u[2] |
36 | | - return du[2] = -u[1] |
37 | | -end |
38 | | -f2(u, p, t) = [u[2], -u[1]] |
39 | | - |
40 | | -function boundary!(residual, u, p, t) |
41 | | - residual[1] = u(0.0)[1] - 5 |
42 | | - return residual[2] = u(5.0)[1] |
43 | | -end |
44 | | -boundary(u, p, t) = [u(0.0)[1] - 5, u(5.0)[1]] |
45 | | - |
46 | | -# Array indexing for boundary conditions |
47 | | -function boundary_indexing!(residual, u, p, t) |
48 | | - residual[1] = u[:, 1][1] - 5 |
49 | | - return residual[2] = u[:, end][1] |
50 | | -end |
51 | | -boundary_indexing(u, p, t) = [u[:, 1][1] - 5, u[:, end][1]] |
52 | | - |
53 | | -function boundary_two_point_a!(resida, ua, p) |
54 | | - return resida[1] = ua[1] - 5 |
55 | | -end |
56 | | -function boundary_two_point_b!(residb, ub, p) |
57 | | - return residb[1] = ub[1] |
58 | | -end |
59 | | - |
60 | | -boundary_two_point_a(ua, p) = [ua[1] - 5] |
61 | | -boundary_two_point_b(ub, p) = [ub[1]] |
62 | | - |
63 | | -# Not able to change the initial condition. |
64 | | -# Hard coded solution. |
65 | | -odef1! = ODEFunction(f1!, analytic = (u0, p, t) -> [5 - t, -1]) |
66 | | -odef1 = ODEFunction(f1, analytic = (u0, p, t) -> [5 - t, -1]) |
67 | | - |
68 | | -odef2! = ODEFunction( |
69 | | - f2!, analytic = ( |
70 | | - u0, p, t, |
71 | | - ) -> [5 * (cos(t) - cot(5) * sin(t)), 5 * (-cos(t) * cot(5) - sin(t))] |
72 | | -) |
73 | | -odef2 = ODEFunction( |
74 | | - f2, analytic = ( |
75 | | - u0, p, t, |
76 | | - ) -> [5 * (cos(t) - cot(5) * sin(t)), 5 * (-cos(t) * cot(5) - sin(t))] |
77 | | -) |
78 | | - |
79 | | -bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1)) |
80 | | - |
81 | | -tspan = (0.0, 5.0) |
82 | | -u0 = [5.0, -3.5] |
83 | | - |
84 | | -probArr = [ |
85 | | - BVProblem(odef1!, boundary!, u0, tspan, nlls = Val(false)), |
86 | | - BVProblem(odef1, boundary, u0, tspan, nlls = Val(false)), |
87 | | - BVProblem(odef2!, boundary!, u0, tspan, nlls = Val(false)), |
88 | | - BVProblem(odef2, boundary, u0, tspan, nlls = Val(false)), |
89 | | - BVProblem(odef2!, boundary_indexing!, u0, tspan, nlls = Val(false)), |
90 | | - BVProblem(odef2, boundary_indexing, u0, tspan, nlls = Val(false)), |
91 | | - TwoPointBVProblem( |
92 | | - odef1!, (boundary_two_point_a!, boundary_two_point_b!), |
93 | | - u0, tspan; bcresid_prototype, nlls = Val(false) |
94 | | - ), |
95 | | - TwoPointBVProblem( |
96 | | - odef1, (boundary_two_point_a, boundary_two_point_b), |
97 | | - u0, tspan; bcresid_prototype, nlls = Val(false) |
98 | | - ), |
99 | | - TwoPointBVProblem( |
100 | | - odef2!, (boundary_two_point_a!, boundary_two_point_b!), |
101 | | - u0, tspan; bcresid_prototype, nlls = Val(false) |
102 | | - ), |
103 | | - TwoPointBVProblem( |
104 | | - odef2, (boundary_two_point_a, boundary_two_point_b), |
105 | | - u0, tspan; bcresid_prototype, nlls = Val(false) |
106 | | - ), |
107 | | -] |
108 | | - |
109 | | -testTol = 0.3 |
110 | | -affineTol = 1.0e-2 |
111 | | -dts = 1 .// 2 .^ (5:-1:3) |
112 | | -# Stage-4 Lobatto methods reach ~1e-14 error (the double-precision floor) at the |
113 | | -# finest `dts` step, which corrupts the Richardson order estimate. A coarser step |
114 | | -# range keeps the finest error above the round-off floor (~1e-12) so the measured |
115 | | -# order reflects the true order 6 rather than floating-point noise. |
116 | | -dts_stage4 = 1 .// 2 .^ (4:-1:2) |
117 | | - |
118 | | -@testset "Affineness" begin |
119 | | - using LinearAlgebra |
120 | | - |
121 | | - @testset "Problem: $i" for i in (1, 2, 7, 8) |
122 | | - prob = probArr[i] |
123 | | - |
124 | | - @testset "LobattoIIIa$stage" for stage in (2, 3, 4, 5) |
125 | | - @time sol = solve(prob, lobattoIIIa_solver(Val(stage)); dt = 0.2, adaptive = false) |
126 | | - @test norm(diff(first.(sol.u)) .+ 0.2, Inf) + abs(sol.u[1][1] - 5) < affineTol |
127 | | - end |
128 | | - @testset "LobattoIIIb$stage" for stage in (2, 3, 4, 5) |
129 | | - @time sol = solve(prob, lobattoIIIb_solver(Val(stage)); dt = 0.2, adaptive = false) |
130 | | - @test norm(diff(first.(sol.u)) .+ 0.2, Inf) + abs(sol.u[1][1] - 5) < affineTol |
131 | | - end |
132 | | - @testset "LobattoIIIc$stage" for stage in (2, 3, 4, 5) |
133 | | - @time sol = solve(prob, lobattoIIIc_solver(Val(stage)); dt = 0.2, adaptive = false) |
134 | | - @test norm(diff(first.(sol.u)) .+ 0.2, Inf) + abs(sol.u[1][1] - 5) < affineTol |
135 | | - end |
136 | | - |
137 | | - @testset "RadauIIa$stage" for stage in (2, 3, 5, 7) |
138 | | - @time sol = solve(prob, radau_solver(Val(stage)); dt = 0.2, adaptive = false) |
139 | | - @test norm(diff(first.(sol.u)) .+ 0.2, Inf) + abs(sol.u[1][1] - 5) < affineTol |
140 | | - end |
141 | | - end |
142 | | -end |
143 | | - |
144 | | -# JET tests have been moved to the separate QA test group (test/qa/) |
145 | | - |
146 | | -@testset "Convergence on Linear" begin |
147 | | - using LinearAlgebra, DiffEqDevTools |
148 | | - |
149 | | - @testset "Problem: $i" for i in (3, 4, 9, 10) |
150 | | - prob = probArr[i] |
151 | | - |
152 | | - @testset "LobattoIIIa$stage" for stage in (2, 3, 4, 5) |
153 | | - stepsizes = stage == 4 ? dts_stage4 : dts |
154 | | - @time sim = test_convergence(stepsizes, prob, lobattoIIIa_solver(Val(stage)); abstol = 1.0e-8) |
155 | | - if stage == 5 |
156 | | - @test_broken sim.𝒪est[:final] ≈ 2 * stage - 2 atol = testTol |
157 | | - else |
158 | | - @test sim.𝒪est[:final] ≈ 2 * stage - 2 atol = testTol |
159 | | - end |
160 | | - end |
161 | | - |
162 | | - @testset "LobattoIIIb$stage" for stage in (2, 3, 4, 5) |
163 | | - stepsizes = stage == 4 ? dts_stage4 : dts |
164 | | - @time sim = test_convergence( |
165 | | - stepsizes, prob, lobattoIIIb_solver(Val(stage)); abstol = 1.0e-8, reltol = 1.0e-8 |
166 | | - ) |
167 | | - if stage == 5 |
168 | | - @test_broken sim.𝒪est[:final] ≈ 2 * stage - 2 atol = testTol |
169 | | - else |
170 | | - @test sim.𝒪est[:final] ≈ 2 * stage - 2 atol = testTol |
171 | | - end |
172 | | - end |
173 | | - |
174 | | - @testset "LobattoIIIc$stage" for stage in (2, 3, 4, 5) |
175 | | - stepsizes = stage == 4 ? dts_stage4 : dts |
176 | | - @time sim = test_convergence( |
177 | | - stepsizes, prob, lobattoIIIc_solver(Val(stage)); abstol = 1.0e-8, reltol = 1.0e-8 |
178 | | - ) |
179 | | - if stage != 4 && first(sim.errors[:final]) < 1.0e-12 |
180 | | - @test_broken sim.𝒪est[:final] ≈ 2 * stage - 2 atol = testTol |
181 | | - else |
182 | | - @test sim.𝒪est[:final] ≈ 2 * stage - 2 atol = testTol |
183 | | - end |
184 | | - end |
185 | | - |
186 | | - @testset "RadauIIa$stage" for stage in (2, 3, 5, 7) |
187 | | - @time sim = test_convergence( |
188 | | - dts, prob, radau_solver(Val(stage)); abstol = 1.0e-8, reltol = 1.0e-8 |
189 | | - ) |
190 | | - if first(sim.errors[:final]) < 1.0e-12 |
191 | | - @test_broken sim.𝒪est[:final] ≈ 2 * stage - 1 atol = testTol |
192 | | - else |
193 | | - @test sim.𝒪est[:final] ≈ 2 * stage - 1 atol = testTol |
194 | | - end |
195 | | - end |
196 | | - end |
197 | | -end |
| 4 | +include("firk_test_setup.jl") |
198 | 5 |
|
199 | 6 | # FIXME: This is a really bad test. Needs interpolation |
200 | 7 | @testset "Simple Pendulum" begin |
|
0 commit comments