@@ -13,8 +13,8 @@ struct RungeKutta <: DiffSolver
1313 order:: Int
1414 supersample:: Int
1515 function RungeKutta (order:: Int , supersample:: Int )
16- if order ≠ 4
17- throw (ArgumentError (" only 4th order Runge-Kutta is supported." ))
16+ if order ≠ 4 && order ≠ 1
17+ throw (ArgumentError (" only 1st and 4th order Runge-Kutta is supported." ))
1818 end
1919 if order < 1
2020 throw (ArgumentError (" order must be greater than 0" ))
2929"""
3030 RungeKutta(order=4; supersample=1)
3131
32- Create a Runge-Kutta solver with optional super-sampling.
32+ Create an explicit Runge-Kutta solver with optional super-sampling.
3333
34- Only the 4th order Runge-Kutta is supported for now . The keyword argument `supersample`
35- provides the number of internal steps (default to 1 step). This solver is allocation-free if
36- the `f!` and `h!` functions do not allocate.
34+ Only the 1st and 4th order Runge-Kutta is supported. The keyword argument
35+ `supersample` provides the number of internal steps (default to 1 step). This solver is
36+ allocation-free if the `f!` and `h!` functions do not allocate.
3737"""
3838RungeKutta (order:: Int = 4 ; supersample:: Int = 1 ) = RungeKutta (order, supersample)
3939
4040" Get the `f!` and `h!` functions for Runge-Kutta solver."
4141function get_solver_functions (NT:: DataType , solver:: RungeKutta , fc!, hc!, Ts, _ , nx, _ , _ )
42+ order = solver. order
4243 Ts_inner = Ts/ solver. supersample
4344 Nc = nx + 2
4445 xcur_cache:: DiffCache{Vector{NT}, Vector{NT}} = DiffCache (zeros (NT, nx), Nc)
4546 k1_cache:: DiffCache{Vector{NT}, Vector{NT}} = DiffCache (zeros (NT, nx), Nc)
4647 k2_cache:: DiffCache{Vector{NT}, Vector{NT}} = DiffCache (zeros (NT, nx), Nc)
4748 k3_cache:: DiffCache{Vector{NT}, Vector{NT}} = DiffCache (zeros (NT, nx), Nc)
4849 k4_cache:: DiffCache{Vector{NT}, Vector{NT}} = DiffCache (zeros (NT, nx), Nc)
49- f! = function inner_solver_f! (xnext, x, u, d, p)
50- CT = promote_type (eltype (x), eltype (u), eltype (d))
51- xcur = get_tmp (xcur_cache, CT)
52- k1 = get_tmp (k1_cache, CT)
53- k2 = get_tmp (k2_cache, CT)
54- k3 = get_tmp (k3_cache, CT)
55- k4 = get_tmp (k4_cache, CT)
56- xterm = xnext
57- @. xcur = x
58- for i= 1 : solver. supersample
59- @. xterm = xcur
60- fc! (k1, xterm, u, d, p)
61- @. xterm = xcur + k1 * Ts_inner/ 2
62- fc! (k2, xterm, u, d, p)
63- @. xterm = xcur + k2 * Ts_inner/ 2
64- fc! (k3, xterm, u, d, p)
65- @. xterm = xcur + k3 * Ts_inner
66- fc! (k4, xterm, u, d, p)
67- @. xcur = xcur + (k1 + 2 k2 + 2 k3 + k4)* Ts_inner/ 6
50+ if order== 4
51+ f! = function rk4_solver! (xnext, x, u, d, p)
52+ CT = promote_type (eltype (x), eltype (u), eltype (d))
53+ xcur = get_tmp (xcur_cache, CT)
54+ k1 = get_tmp (k1_cache, CT)
55+ k2 = get_tmp (k2_cache, CT)
56+ k3 = get_tmp (k3_cache, CT)
57+ k4 = get_tmp (k4_cache, CT)
58+ xterm = xnext
59+ @. xcur = x
60+ for i= 1 : solver. supersample
61+ fc! (k1, xcur, u, d, p)
62+ @. xterm = xcur + k1 * Ts_inner/ 2
63+ fc! (k2, xterm, u, d, p)
64+ @. xterm = xcur + k2 * Ts_inner/ 2
65+ fc! (k3, xterm, u, d, p)
66+ @. xterm = xcur + k3 * Ts_inner
67+ fc! (k4, xterm, u, d, p)
68+ @. xcur = xcur + (k1 + 2 k2 + 2 k3 + k4)* Ts_inner/ 6
69+ end
70+ @. xnext = xcur
71+ return nothing
6872 end
69- @. xnext = xcur
70- return nothing
7173 end
7274 h! = hc!
7375 return f!, h!
0 commit comments