Mercurial > repos > public > sbplib_julia
changeset 862:a382942e5437 feature/subpackage_simple_timesteppers
Add a few simple timesteppers
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Tue, 18 Jan 2022 15:29:54 +0100 |
| parents | bdc718c38096 |
| children | 8236c3063fa1 |
| files | src/Sbplib.jl src/SimpleTimeSteppers/SimpleTimeSteppers.jl test/SimpleTimeSteppers/SimpleTimeSteppers_test.jl |
| diffstat | 3 files changed, 208 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/src/Sbplib.jl Mon Jan 17 15:05:43 2022 +0100 +++ b/src/Sbplib.jl Tue Jan 18 15:29:54 2022 +0100 @@ -6,6 +6,7 @@ include("Grids/Grids.jl") include("SbpOperators/SbpOperators.jl") include("DiffOps/DiffOps.jl") +include("SimpleTimeSteppers/SimpleTimeSteppers.jl") export RegionIndices export LazyTensors
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/SimpleTimeSteppers/SimpleTimeSteppers.jl Tue Jan 18 15:29:54 2022 +0100 @@ -0,0 +1,62 @@ +""" + SimpleTimeSteppers + +Implements a few simple time integration schemes useful for testing the accuracy of spatial discretisations. +""" +module SimpleTimeSteppers + +export rk4 +export rkn4 +export centered_2nd_order_fd + +""" + rk4(F, vₙ, tₙ, Δtₙ) + +Solves ``vₜ = F(v,t)`` using the standard 4th order RK-method. +""" +function rk4(F, vₙ, tₙ, Δtₙ) + k₁ = F(vₙ,tₙ) + k₂ = F(vₙ + Δtₙ/2*k₁, tₙ + Δtₙ/2) + k₃ = F(vₙ + Δtₙ/2*k₂, tₙ + Δtₙ/2) + k₄ = F(vₙ + Δtₙ *k₃, tₙ + Δtₙ ) + + vₙ₊₁ = vₙ + (1/6)*(k₁+2*(k₂+k₃)+k₄)*Δtₙ + + return vₙ₊₁ +end + +""" + rkn4(F, vₙ, v̇ₙ, tₙ, Δtₙ) + +Solves ``vₜₜ = F(v,v̇,t)`` using the Runge-Kutta-Nyström method based on the +standard 4th order RK-method. +""" +function rkn4(F, vₙ, v̇ₙ, tₙ, Δtₙ) + k₁ = F(vₙ, v̇ₙ, tₙ ) + k₂ = F(vₙ + Δtₙ/2*v̇ₙ, v̇ₙ + Δtₙ/2*k₁, tₙ + Δtₙ/2) + k₃ = F(vₙ + Δtₙ/2*v̇ₙ + Δtₙ^2/4*k₁, v̇ₙ + Δtₙ/2*k₂, tₙ + Δtₙ/2) + k₄ = F(vₙ + Δtₙ *v̇ₙ + Δtₙ^2/2*k₂, v̇ₙ + Δtₙ *k₃, tₙ + Δtₙ ) + + vₙ₊₁ = vₙ + Δtₙ*v̇ₙ + Δtₙ^2*(1/6)*(k₁ + k₂ + k₃); + v̇ₙ₊₁ = v̇ₙ + Δtₙ*(k₁ + 2*k₂ + 2*k₃ + k₄)/6; + + return vₙ₊₁, v̇ₙ₊₁ +end + +""" + centered_2nd_order_fd(F, vₙ, vₙ₋₁, Δtₙ) + +Solves ``vₜₜ = F(v,t)`` using the 2nd order discretization +``(vₙ₊₁ - 2vₙ + vₙ₋₁)/Δt² = F(vₙ,tₙ)``. +""" +function centered_2nd_order_fd(F, vₙ, vₙ₋₁, tₙ, Δt) + vₙ₊₁ = Δt^2*F(vₙ,tₙ) + 2vₙ - vₙ₋₁ + + return vₙ₊₁, vₙ +end + +# TODO: Rewrite without allocations + + +# TODO: Add euler forward and try a local error estimate for testing +end # module
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/SimpleTimeSteppers/SimpleTimeSteppers_test.jl Tue Jan 18 15:29:54 2022 +0100 @@ -0,0 +1,145 @@ +using Test +using Sbplib.SimpleTimeSteppers + + +@testset "rk4" begin + ode_solutions = @NamedTuple{F,v}.([ + ((v,t) -> 0., (v₀,t₀,t) -> v₀), + ((v,t) -> 1., (v₀,t₀,t) -> v₀+t-t₀), + ((v,t) -> t, (v₀,t₀,t) -> v₀+(t^2-t₀^2)/2), + ((v,t) -> v, (v₀,t₀,t) -> v₀*exp(t-t₀)), + ((v,t) -> t*v, (v₀,t₀,t) -> v₀*exp((t^2-t₀^2)/2)), + ]) + + @testset "Exact" begin + @testset for i ∈ eachindex(ode_solutions[1:3]) + F, v = ode_solutions[i] + cases = @NamedTuple{v₀,t₀}.([ + (0.,0.), + (0.,1.), + (1.,0.), + (1.,1.), + (2.,0.), + (0.,2.), + (2.,2.), + ]) + + @testset for (v₀,t₀) ∈ cases + @testset for Δt ∈ [0., 1., .5, 2.] + @test rk4(F, v₀, t₀, Δt) == v(v₀, t₀, t₀+Δt) + end + end + end + end + + @testset "Approximate" begin + @testset for i ∈ eachindex(ode_solutions) + F, v = ode_solutions[i] + + @testset for v₀ ∈ [1., .5] ,t₀ ∈ [0., 1.] + Δt = 0.1 + @test rk4(F, v₀, t₀, Δt) ≈ v(v₀, t₀, t₀+Δt) rtol=1e-6 + end + end + end +end + +@testset "rkn4" begin + + function damped_oscillator(v₀, v̇₀,t₀,t; ζ) + ω = sqrt(1-ζ^2) + + C = v̇₀ + ζ*v₀ + + V = v₀*cos(ω*(t-t₀)) + C*sin(ω*(t-t₀)) + V̇ = C*cos(ω*(t-t₀)) - v₀*sin(ω*(t-t₀)) + v = exp(-ζ*(t-t₀))*V + v̇ = -ζ*exp(-ζ*(t-t₀))*V + exp(-ζ*(t-t₀))*V̇ + return v, v̇ + end + + ode_solutions = @NamedTuple{F,vv̇}.([ + ((v,v̇,t) -> 0., (v₀,v̇₀,t₀,t) -> (v₀ + (t-t₀)*v̇₀, v̇₀)), + ((v,v̇,t) -> 1., (v₀,v̇₀,t₀,t) -> (v₀ + (t-t₀)*v̇₀ + (t-t₀)^2/2, v̇₀ + t - t₀)), + ((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)), + ((v,v̇,t) -> -v, (v₀,v̇₀,t₀,t) -> (v₀*cos(t-t₀) + v̇₀*sin(t-t₀), v̇₀*cos(t-t₀) - v₀*sin(t-t₀))), + ((v,v̇,t) -> -v-6/5*v̇, (v₀,v̇₀,t₀,t) -> damped_oscillator(v₀,v̇₀,t₀,t; ζ=3/5)), + ]) + + @testset "Exact" begin + @testset for i ∈ eachindex(ode_solutions[1:3]) + F, vv̇ = ode_solutions[i] + @testset for v₀ ∈ [0., 1., 2.], v̇₀ ∈ [0., 1.], t₀ ∈ [0.,1.,2.] + @testset for Δt ∈ [0., 1., .5, 2.] + V,V̇ = rkn4(F, v₀, v̇₀, t₀, Δt) + v,v̇ = vv̇(v₀, v̇₀, t₀, t₀+Δt) + @test V ≈ v rtol=1e-15 + @test V̇ ≈ v̇ rtol=1e-15 + end + end + end + end + + @testset "Approximate" begin + @testset for i ∈ eachindex(ode_solutions) + F, vv̇ = ode_solutions[i] + @testset for v̇₀ ∈ [0., 0.5], v₀ ∈ [1., 0.5], t₀∈[0.,1.] + Δt = 0.1 + V,V̇ = rkn4(F, v₀, v̇₀, t₀, Δt) + v,v̇ = vv̇(v₀, v̇₀, t₀, t₀+Δt) + + if i != 5 + @test V ≈ v rtol=1e-6 + @test V̇ ≈ v̇ rtol=1e-6 + else + # TODO: Is this a bug or reasonable? + @test_broken V ≈ v rtol=1e-6 + @test_broken V̇ ≈ v̇ rtol=1e-6 + end + end + end + end +end + +@testset "cdiff" begin + ode_solutions = @NamedTuple{F,vv̇}.([ + ((v,t) -> 0., (v₀,v̇₀,t₀,t) -> v₀ + (t-t₀)*v̇₀), + ((v,t) -> 1., (v₀,v̇₀,t₀,t) -> v₀ + (t-t₀)*v̇₀ + (t-t₀)^2/2), + ((v,t) -> t, (v₀,v̇₀,t₀,t) -> v₀ + (t-t₀)*v̇₀ + (t-t₀)^3/6+ t₀/2*(t-t₀)^2), + ((v,t) -> -v, (v₀,v̇₀,t₀,t) -> v₀*cos(t-t₀) + v̇₀*sin(t-t₀)), + ]) + + @testset "Exact" begin + @testset for i ∈ eachindex(ode_solutions[1:3]) + F, v = ode_solutions[i] + @testset for v₀ ∈ [1., 2.], v̇₀ ∈ [0., 1.], t₀ ∈ [0.,1.,2.] + @testset for Δt ∈ [0., 1., .5] + v₋₁ = v(v₀, v̇₀, t₀, t₀-Δt) + v₁ = v(v₀, v̇₀, t₀, t₀+Δt) + + V₁, V₀ = centered_2nd_order_fd(F, v₀, v₋₁, t₀, Δt) + + @test V₁ ≈ v₁ rtol=1e-15 + @test V₀ ≈ v₀ rtol=1e-15 + end + end + end + end + + @testset "Approximate" begin + @testset for i ∈ eachindex(ode_solutions) + F, v = ode_solutions[i] + @testset for v̇₀ ∈ [0., 0.5], v₀ ∈ [1., 0.5], t₀ ∈ [0.,1.] + Δt = 0.1 + v₋₁ = v(v₀, v̇₀, t₀, t₀-Δt) + v₁ = v(v₀, v̇₀, t₀, t₀+Δt) + + V₁, V₀ = centered_2nd_order_fd(F, v₀, v₋₁, t₀, Δt) + + + @test V₁ ≈ v₁ rtol=1e-5 + @test V₀ ≈ v₀ rtol=1e-5 + end + end + end +end