Mercurial > repos > public > sbplib_julia

comparison test/SimpleTimeSteppers/SimpleTimeSteppers_test.jl @ 862:a382942e5437 feature/subpackage_simple_timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add a few simple timesteppers
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 18 Jan 2022 15:29:54 +0100
parents
children
comparison
equal deleted inserted replaced
859:bdc718c38096 862:a382942e5437
1 using Test
2 using Sbplib.SimpleTimeSteppers
3
4
5 @testset "rk4" begin
6 ode_solutions = @NamedTuple{F,v}.([
7 ((v,t) -> 0., (v₀,t₀,t) -> v₀),
8 ((v,t) -> 1., (v₀,t₀,t) -> v₀+t-t₀),
9 ((v,t) -> t, (v₀,t₀,t) -> v₀+(t^2-t₀^2)/2),
10 ((v,t) -> v, (v₀,t₀,t) -> v₀*exp(t-t₀)),
11 ((v,t) -> t*v, (v₀,t₀,t) -> v₀*exp((t^2-t₀^2)/2)),
12 ])
13
14 @testset "Exact" begin
15 @testset for i ∈ eachindex(ode_solutions[1:3])
16 F, v = ode_solutions[i]
17 cases = @NamedTuple{v₀,t₀}.([
18 (0.,0.),
19 (0.,1.),
20 (1.,0.),
21 (1.,1.),
22 (2.,0.),
23 (0.,2.),
24 (2.,2.),
25 ])
26
27 @testset for (v₀,t₀) ∈ cases
28 @testset for Δt ∈ [0., 1., .5, 2.]
29 @test rk4(F, v₀, t₀, Δt) == v(v₀, t₀, t₀+Δt)
30 end
31 end
32 end
33 end
34
35 @testset "Approximate" begin
36 @testset for i ∈ eachindex(ode_solutions)
37 F, v = ode_solutions[i]
38
39 @testset for v₀ ∈ [1., .5] ,t₀ ∈ [0., 1.]
40 Δt = 0.1
41 @test rk4(F, v₀, t₀, Δt) ≈ v(v₀, t₀, t₀+Δt) rtol=1e-6
42 end
43 end
44 end
45 end
46
47 @testset "rkn4" begin
48
49 function damped_oscillator(v₀, v̇₀,t₀,t; ζ)
50 ω = sqrt(1-ζ^2)
51
52 C = v̇₀ + ζ*v₀
53
54 V = v₀*cos(ω*(t-t₀)) + C*sin(ω*(t-t₀))
55 V̇ = C*cos(ω*(t-t₀)) - v₀*sin(ω*(t-t₀))
56 v = exp(-ζ*(t-t₀))*V
57 v̇ = -ζ*exp(-ζ*(t-t₀))*V + exp(-ζ*(t-t₀))*V̇
58 return v, v̇
59 end
60
61 ode_solutions = @NamedTuple{F,vv̇}.([
62 ((v,v̇,t) -> 0., (v₀,v̇₀,t₀,t) -> (v₀ + (t-t₀)*v̇₀, v̇₀)),
63 ((v,v̇,t) -> 1., (v₀,v̇₀,t₀,t) -> (v₀ + (t-t₀)*v̇₀ + (t-t₀)^2/2, v̇₀ + t - t₀)),
64 ((v,v̇,t) -> t, (v₀,v̇₀,t₀,t) -> (v₀ + (t-t₀)*v̇₀ + (t-t₀)^3/6+ t₀/2*(t-t₀)^2, v̇₀ + (t^2-t₀^2)/2)),
65 ((v,v̇,t) -> -v, (v₀,v̇₀,t₀,t) -> (v₀*cos(t-t₀) + v̇₀*sin(t-t₀), v̇₀*cos(t-t₀) - v₀*sin(t-t₀))),
66 ((v,v̇,t) -> -v-6/5*v̇, (v₀,v̇₀,t₀,t) -> damped_oscillator(v₀,v̇₀,t₀,t; ζ=3/5)),
67 ])
68
69 @testset "Exact" begin
70 @testset for i ∈ eachindex(ode_solutions[1:3])
71 F, vv̇ = ode_solutions[i]
72 @testset for v₀ ∈ [0., 1., 2.], v̇₀ ∈ [0., 1.], t₀ ∈ [0.,1.,2.]
73 @testset for Δt ∈ [0., 1., .5, 2.]
74 V,V̇ = rkn4(F, v₀, v̇₀, t₀, Δt)
75 v,v̇ = vv̇(v₀, v̇₀, t₀, t₀+Δt)
76 @test V ≈ v rtol=1e-15
77 @test V̇ ≈ v̇ rtol=1e-15
78 end
79 end
80 end
81 end
82
83 @testset "Approximate" begin
84 @testset for i ∈ eachindex(ode_solutions)
85 F, vv̇ = ode_solutions[i]
86 @testset for v̇₀ ∈ [0., 0.5], v₀ ∈ [1., 0.5], t₀∈[0.,1.]
87 Δt = 0.1
88 V,V̇ = rkn4(F, v₀, v̇₀, t₀, Δt)
89 v,v̇ = vv̇(v₀, v̇₀, t₀, t₀+Δt)
90
91 if i != 5
92 @test V ≈ v rtol=1e-6
93 @test V̇ ≈ v̇ rtol=1e-6
94 else
95 # TODO: Is this a bug or reasonable?
96 @test_broken V ≈ v rtol=1e-6
97 @test_broken V̇ ≈ v̇ rtol=1e-6
98 end
99 end
100 end
101 end
102 end
103
104 @testset "cdiff" begin
105 ode_solutions = @NamedTuple{F,vv̇}.([
106 ((v,t) -> 0., (v₀,v̇₀,t₀,t) -> v₀ + (t-t₀)*v̇₀),
107 ((v,t) -> 1., (v₀,v̇₀,t₀,t) -> v₀ + (t-t₀)*v̇₀ + (t-t₀)^2/2),
108 ((v,t) -> t, (v₀,v̇₀,t₀,t) -> v₀ + (t-t₀)*v̇₀ + (t-t₀)^3/6+ t₀/2*(t-t₀)^2),
109 ((v,t) -> -v, (v₀,v̇₀,t₀,t) -> v₀*cos(t-t₀) + v̇₀*sin(t-t₀)),
110 ])
111
112 @testset "Exact" begin
113 @testset for i ∈ eachindex(ode_solutions[1:3])
114 F, v = ode_solutions[i]
115 @testset for v₀ ∈ [1., 2.], v̇₀ ∈ [0., 1.], t₀ ∈ [0.,1.,2.]
116 @testset for Δt ∈ [0., 1., .5]
117 v₋₁ = v(v₀, v̇₀, t₀, t₀-Δt)
118 v₁ = v(v₀, v̇₀, t₀, t₀+Δt)
119
120 V₁, V₀ = centered_2nd_order_fd(F, v₀, v₋₁, t₀, Δt)
121
122 @test V₁ ≈ v₁ rtol=1e-15
123 @test V₀ ≈ v₀ rtol=1e-15
124 end
125 end
126 end
127 end
128
129 @testset "Approximate" begin
130 @testset for i ∈ eachindex(ode_solutions)
131 F, v = ode_solutions[i]
132 @testset for v̇₀ ∈ [0., 0.5], v₀ ∈ [1., 0.5], t₀ ∈ [0.,1.]
133 Δt = 0.1
134 v₋₁ = v(v₀, v̇₀, t₀, t₀-Δt)
135 v₁ = v(v₀, v̇₀, t₀, t₀+Δt)
136
137 V₁, V₀ = centered_2nd_order_fd(F, v₀, v₋₁, t₀, Δt)
138
139
140 @test V₁ ≈ v₁ rtol=1e-5
141 @test V₀ ≈ v₀ rtol=1e-5
142 end
143 end
144 end
145 end