diff --git a/docs/src/manual/ensemblegpukernel.md b/docs/src/manual/ensemblegpukernel.md index 8950d30a..d828ff56 100644 --- a/docs/src/manual/ensemblegpukernel.md +++ b/docs/src/manual/ensemblegpukernel.md @@ -15,7 +15,9 @@ GPUVern9 GPUEM GPUSIEA GPURosenbrock23 +GPURodas3 GPURodas4 +GPURodas42 GPURodas5P GPUKvaerno3 GPUKvaerno5 diff --git a/src/DiffEqGPU.jl b/src/DiffEqGPU.jl index 3160b4e5..3539d3c3 100644 --- a/src/DiffEqGPU.jl +++ b/src/DiffEqGPU.jl @@ -62,6 +62,7 @@ include("ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl") include("ensemblegpukernel/perform_step/gpu_em_perform_step.jl") include("ensemblegpukernel/perform_step/gpu_siea_perform_step.jl") include("ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl") +include("ensemblegpukernel/perform_step/gpu_rodas3_perform_step.jl") include("ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl") include("ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl") include("ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl") @@ -78,7 +79,8 @@ export EnsembleCPUArray, EnsembleGPUArray, EnsembleGPUKernel, LinSolveGPUSplitFa export GPUTsit5, GPUVern7, GPUVern9, GPUEM, GPUSIEA ## Stiff ODE solvers -export GPURosenbrock23, GPURodas4, GPURodas5P, GPUKvaerno3, GPUKvaerno5 +export GPURosenbrock23, GPURodas3, GPURodas4, GPURodas42, GPURodas5P, GPUKvaerno3, + GPUKvaerno5 export terminate! diff --git a/src/ensemblegpukernel/alg_utils.jl b/src/ensemblegpukernel/alg_utils.jl index b9701aaa..be527050 100644 --- a/src/ensemblegpukernel/alg_utils.jl +++ b/src/ensemblegpukernel/alg_utils.jl @@ -6,7 +6,9 @@ alg_order(alg::GPUTsit5) = 5 alg_order(alg::GPUVern7) = 7 alg_order(alg::GPUVern9) = 9 alg_order(alg::GPURosenbrock23) = 2 +alg_order(alg::GPURodas3) = 3 alg_order(alg::GPURodas4) = 4 +alg_order(alg::GPURodas42) = 4 alg_order(alg::GPURodas5P) = 5 alg_order(alg::GPUKvaerno3) = 3 alg_order(alg::GPUKvaerno5) = 5 diff --git a/src/ensemblegpukernel/gpukernel_algorithms.jl b/src/ensemblegpukernel/gpukernel_algorithms.jl index ce6e4d3f..60302ae3 100644 --- a/src/ensemblegpukernel/gpukernel_algorithms.jl +++ b/src/ensemblegpukernel/gpukernel_algorithms.jl @@ -31,6 +31,14 @@ generation with EnsembleGPUKernel. """ struct GPURosenbrock23{AD} <: GPUODEImplicitAlgorithm{AD} end +""" +GPURodas3() + +A specialized implementation of the `Rodas3` method specifically for kernel +generation with EnsembleGPUKernel. +""" +struct GPURodas3{AD} <: GPUODEImplicitAlgorithm{AD} end + """ GPURodas4() @@ -39,6 +47,14 @@ generation with EnsembleGPUKernel. """ struct GPURodas4{AD} <: GPUODEImplicitAlgorithm{AD} end +""" +GPURodas42() + +A specialized implementation of the `Rodas42` method specifically for kernel +generation with EnsembleGPUKernel. +""" +struct GPURodas42{AD} <: GPUODEImplicitAlgorithm{AD} end + """ GPURodas5P() @@ -63,7 +79,7 @@ generation with EnsembleGPUKernel. """ struct GPUKvaerno5{AD} <: GPUODEImplicitAlgorithm{AD} end -for Alg in [:GPURosenbrock23, :GPURodas4, :GPURodas5P, :GPUKvaerno3, :GPUKvaerno5] +for Alg in [:GPURosenbrock23, :GPURodas3, :GPURodas4, :GPURodas42, :GPURodas5P, :GPUKvaerno3, :GPUKvaerno5] @eval begin function $Alg(; autodiff = Val{true}()) return $Alg{SciMLBase._unwrap_val(autodiff)}() diff --git a/src/ensemblegpukernel/integrators/stiff/types.jl b/src/ensemblegpukernel/integrators/stiff/types.jl index 96bf1772..2e291186 100644 --- a/src/ensemblegpukernel/integrators/stiff/types.jl +++ b/src/ensemblegpukernel/integrators/stiff/types.jl @@ -365,6 +365,295 @@ const GPUARodas4I = GPUARodas4Integrator ) end +@inline function init( + alg::GPURodas42, f::F, IIP::Bool, u0::S, t0::T, dt::T, + p::P, tstops::TS, + callback::CB, + save_everystep::Bool, + saveat::ST + ) where { + F, P, T, + S, + TS, CB, ST, + } + !IIP && @assert S <: SArray + event_last_time = 1 + vector_event_last_time = 0 + last_event_error = zero(T) + + tab = Rodas42Tableau(T, T) + + return integ = GPURodas4I{IIP, S, T, ST, P, F, TS, CB, typeof(tab), typeof(alg)}( + alg, f, + copy(u0), + copy(u0), + copy(u0), t0, + t0, + t0, + dt, + sign(dt), p, + true, + tstops, 1, + callback, + save_everystep, + saveat, 1, 1, + event_last_time, + vector_event_last_time, + last_event_error, + copy(u0), + copy(u0), tab, + DiffEqBase.ReturnCode.Default + ) +end + +@inline function init( + alg::GPURodas42, f::F, IIP::Bool, u0::S, t0::T, tf::T, + dt::T, p::P, + abstol::TOL, reltol::TOL, + internalnorm::N, tstops::TS, + callback::CB, + saveat::ST + ) where { + F, P, S, T, N, TOL, TS, + CB, ST, + } + !IIP && @assert S <: SArray + qoldinit = T(1.0e-4) + event_last_time = 1 + vector_event_last_time = 0 + last_event_error = zero(T) + + tab = Rodas42Tableau(T, T) + + return integ = GPUARodas4I{ + IIP, S, T, ST, P, F, N, TOL, typeof(qoldinit), TS, CB, typeof(tab), + typeof(alg), + }( + alg, + f, + copy(u0), + copy(u0), + copy(u0), + t0, + t0, + t0, + tf, + dt, + dt, + sign( + tf - + t0 + ), + p, + true, + tstops, + 1, + callback, + false, + saveat, + 1, 1, + event_last_time, + vector_event_last_time, + last_event_error, + copy(u0), + copy(u0), + tab, + qoldinit, + abstol, + reltol, + internalnorm, + DiffEqBase.ReturnCode.Default + ) +end + +########################## +# Rodas 3 +########################## + +mutable struct GPURodas3Integrator{IIP, S, T, ST, P, F, TS, CB, TabType, AlgType} <: + DiffEqBase.AbstractODEIntegrator{AlgType, IIP, S, T} + alg::AlgType + f::F + uprev::S + u::S + tmp::S + tprev::T + t::T + t0::T + dt::T + tdir::T + p::P + u_modified::Bool + tstops::TS + tstops_idx::Int + callback::CB + save_everystep::Bool + saveat::ST + cur_t::Int + step_idx::Int + event_last_time::Int + vector_event_last_time::Int + last_event_error::T + k1::S + k2::S + tab::TabType + retcode::DiffEqBase.ReturnCode.T +end +const GPURodas3I = GPURodas3Integrator + +@inline function init( + alg::GPURodas3, f::F, IIP::Bool, u0::S, t0::T, dt::T, + p::P, tstops::TS, + callback::CB, + save_everystep::Bool, + saveat::ST + ) where { + F, P, T, + S, + TS, CB, ST, + } + !IIP && @assert S <: SArray + event_last_time = 1 + vector_event_last_time = 0 + last_event_error = zero(T) + + tab = Rodas3Tableau(T, T) + + return integ = GPURodas3I{IIP, S, T, ST, P, F, TS, CB, typeof(tab), typeof(alg)}( + alg, f, + copy(u0), + copy(u0), + copy(u0), t0, + t0, + t0, + dt, + sign(dt), p, + true, + tstops, 1, + callback, + save_everystep, + saveat, 1, 1, + event_last_time, + vector_event_last_time, + last_event_error, + copy(u0), + copy(u0), tab, + DiffEqBase.ReturnCode.Default + ) +end + +mutable struct GPUARodas3Integrator{ + IIP, + S, + T, + ST, + P, + F, + N, + TOL, + Q, + TS, + CB, + TabType, + AlgType, + } <: + DiffEqBase.AbstractODEIntegrator{AlgType, IIP, S, T} + alg::AlgType + f::F + uprev::S + u::S + tmp::S + tprev::T + t::T + t0::T + tf::T + dt::T + dtnew::T + tdir::T + p::P + u_modified::Bool + tstops::TS + tstops_idx::Int + callback::CB + save_everystep::Bool + saveat::ST + cur_t::Int + step_idx::Int + event_last_time::Int + vector_event_last_time::Int + last_event_error::T + k1::S + k2::S + tab::TabType + qold::Q + abstol::TOL + reltol::TOL + internalnorm::N + retcode::DiffEqBase.ReturnCode.T +end + +const GPUARodas3I = GPUARodas3Integrator + +@inline function init( + alg::GPURodas3, f::F, IIP::Bool, u0::S, t0::T, tf::T, + dt::T, p::P, + abstol::TOL, reltol::TOL, + internalnorm::N, tstops::TS, + callback::CB, + saveat::ST + ) where { + F, P, S, T, N, TOL, TS, + CB, ST, + } + !IIP && @assert S <: SArray + qoldinit = T(1.0e-4) + event_last_time = 1 + vector_event_last_time = 0 + last_event_error = zero(T) + + tab = Rodas3Tableau(T, T) + + return integ = GPUARodas3I{ + IIP, S, T, ST, P, F, N, TOL, typeof(qoldinit), TS, CB, typeof(tab), + typeof(alg), + }( + alg, + f, + copy(u0), + copy(u0), + copy(u0), + t0, + t0, + t0, + tf, + dt, + dt, + sign( + tf - + t0 + ), + p, + true, + tstops, + 1, + callback, + false, + saveat, + 1, 1, + event_last_time, + vector_event_last_time, + last_event_error, + copy(u0), + copy(u0), + tab, + qoldinit, + abstol, + reltol, + internalnorm, + DiffEqBase.ReturnCode.Default + ) +end + ########################## # Rodas 5P ########################## @@ -937,6 +1226,26 @@ end return integrator.u_modified = bool end +# GPURodas3Integrator +@inline function (integrator::GPURodas3I)(t) + Θ = (t - integrator.tprev) / integrator.dt + return _ode_interpolant(Θ, integrator.dt, integrator.uprev, integrator) +end + +@inline function DiffEqBase.u_modified!(integrator::GPURodas3I, bool::Bool) + return integrator.u_modified = bool +end + +# GPUARodas3Integrator +@inline function (integrator::GPUARodas3I)(t) + Θ = (t - integrator.tprev) / integrator.dt + return _ode_interpolant(Θ, integrator.dt, integrator.uprev, integrator) +end + +@inline function DiffEqBase.u_modified!(integrator::GPUARodas3I, bool::Bool) + return integrator.u_modified = bool +end + # GPURodas4Integrator @inline function (integrator::GPURodas4I)(t) Θ = (t - integrator.tprev) / integrator.dt diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas3_perform_step.jl new file mode 100644 index 00000000..3e0dfee4 --- /dev/null +++ b/src/ensemblegpukernel/perform_step/gpu_rodas3_perform_step.jl @@ -0,0 +1,193 @@ +@inline function step!(integ::GPURodas3I{false, S, T}, ts, us) where {T, S} + dt = integ.dt + t = integ.t + p = integ.p + f = integ.f + integ.uprev = integ.u + uprev = integ.u + @unpack a21, a31, a32, a41, a42, a43, C21, C31, C32, C41, C42, C43, + γ, c2, c3, d1, d2, d3, d4 = integ.tab + + integ.tprev = t + saved_in_cb = false + if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && + (integ.tstops[integ.tstops_idx] - integ.t - integ.dt - 100 * eps(T) < 0) + integ.t = integ.tstops[integ.tstops_idx] + dt = integ.t - integ.tprev + integ.tstops_idx += 1 + else + integ.t += dt + end + + Jf, _ = build_J_W(integ.alg, f, γ, dt) + J = Jf(uprev, p, t) + + Tgrad = build_tgrad(integ.alg, f) + dT = Tgrad(uprev, p, t) + + dtC21 = C21 / dt + dtC31 = C31 / dt + dtC32 = C32 / dt + dtC41 = C41 / dt + dtC42 = C42 / dt + dtC43 = C43 / dt + + dtd1 = dt * d1 + dtd2 = dt * d2 + dtd3 = dt * d3 + dtd4 = dt * d4 + dtgamma = dt * γ + + mass_matrix = f.mass_matrix + W = J - mass_matrix * inv(dtgamma) + du0 = f(uprev, p, t) + + linsolve_tmp = du0 + dtd1 * dT + k1 = linear_solve(W, -linsolve_tmp) + u = uprev + a21 * k1 + du = f(u, p, t + c2 * dt) + + linsolve_tmp = du + dtd2 * dT + mass_matrix * (dtC21 * k1) + k2 = linear_solve(W, -linsolve_tmp) + u = uprev + a31 * k1 + a32 * k2 + du = f(u, p, t + c3 * dt) + + linsolve_tmp = du + dtd3 * dT + mass_matrix * (dtC31 * k1 + dtC32 * k2) + k3 = linear_solve(W, -linsolve_tmp) + u = uprev + a41 * k1 + a42 * k2 + a43 * k3 + du = f(u, p, t + dt) + + linsolve_tmp = du + dtd4 * dT + mass_matrix * (dtC41 * k1 + dtC42 * k2 + dtC43 * k3) + k4 = linear_solve(W, -linsolve_tmp) + integ.u = u + k4 + + @inbounds begin + integ.k1 = du0 + integ.k2 = f(integ.u, p, t + dt) + end + + _, saved_in_cb = handle_callbacks!(integ, ts, us) + + return saved_in_cb +end + +@inline function step!(integ::GPUARodas3I{false, S, T}, ts, us) where {T, S} + beta1, beta2, qmax, qmin, gamma, qoldinit, + _ = build_adaptive_controller_cache( + integ.alg, + T + ) + + dt = integ.dtnew + t = integ.t + p = integ.p + tf = integ.tf + + f = integ.f + integ.uprev = integ.u + uprev = integ.u + + qold = integ.qold + abstol = integ.abstol + reltol = integ.reltol + + @unpack a21, a31, a32, a41, a42, a43, C21, C31, C32, C41, C42, C43, + γ, c2, c3, d1, d2, d3, d4 = integ.tab + + Jf, _ = build_J_W(integ.alg, f, γ, dt) + J = Jf(uprev, p, t) + + Tgrad = build_tgrad(integ.alg, f) + dT = Tgrad(uprev, p, t) + + EEst = convert(T, Inf) + + while EEst > convert(T, 1.0) + dt < convert(T, 1.0f-14) && error("dt 1 + dt = dt / min(inv(qmin), q11 / gamma) + else + q = max(inv(qmax), min(inv(qmin), q / gamma)) + qold = max(EEst, qoldinit) + dtnew = dt / q + dtnew = min(abs(dtnew), abs(tf - t - dt)) + + @inbounds begin + integ.k1 = du0 + integ.k2 = f(u, p, t + dt) + end + + integ.dt = dt + integ.dtnew = dtnew + integ.qold = qold + integ.tprev = t + integ.u = u + + if (tf - t - dt) < convert(T, 1.0f-14) + integ.t = tf + else + if integ.tstops !== nothing && integ.tstops_idx <= length(integ.tstops) && + integ.tstops[integ.tstops_idx] - integ.t - integ.dt - + 100 * eps(T) < 0 + integ.t = integ.tstops[integ.tstops_idx] + integ.u = integ(integ.t) + dt = integ.t - integ.tprev + integ.tstops_idx += 1 + else + integ.t += dt + end + end + end + end + + _, saved_in_cb = handle_callbacks!(integ, ts, us) + + return saved_in_cb +end diff --git a/src/ensemblegpukernel/tableaus/rodas_tableaus.jl b/src/ensemblegpukernel/tableaus/rodas_tableaus.jl index 2d3c278b..b66841fb 100644 --- a/src/ensemblegpukernel/tableaus/rodas_tableaus.jl +++ b/src/ensemblegpukernel/tableaus/rodas_tableaus.jl @@ -103,6 +103,117 @@ function Rodas4Tableau(T::Type{T1}, T2::Type{T1}) where {T1} ) end +function Rodas42Tableau(T::Type{T1}, T2::Type{T1}) where {T1} + γ = convert(T, 0.25) + + a21 = convert(T, 1.4028884) + a31 = convert(T, 0.6581212688557198) + a32 = -convert(T, 1.320936088384301) + a41 = convert(T, 7.131197445744498) + a42 = convert(T, 16.02964143958207) + a43 = -convert(T, 5.561572550509766) + a51 = convert(T, 22.73885722420363) + a52 = convert(T, 67.38147284535289) + a53 = -convert(T, 31.2187749303856) + a54 = convert(T, 0.7285641833203814) + + C21 = -convert(T, 5.1043536) + C31 = -convert(T, 2.899967805418783) + C32 = convert(T, 4.040399359702244) + C41 = -convert(T, 32.64449927841361) + C42 = -convert(T, 99.35311008728094) + C43 = convert(T, 49.99119122405989) + C51 = -convert(T, 76.46023087151691) + C52 = -convert(T, 278.5942120829058) + C53 = convert(T, 153.9294840910643) + C54 = convert(T, 10.97101866258358) + C61 = -convert(T, 76.29701586804983) + C62 = -convert(T, 294.2795630511232) + C63 = convert(T, 162.0029695867566) + C64 = convert(T, 23.6516690309527) + C65 = -convert(T, 7.652977706771382) + + c2 = convert(T2, 0.3507221) + c3 = convert(T2, 0.2557041) + c4 = convert(T2, 0.681779) + + d1 = convert(T, 0.25) + d2 = -convert(T, 0.0690221) + d3 = -convert(T, 0.0009672) + d4 = -convert(T, 0.087979) + + h21 = -convert(T, 38.71940424117216) + h22 = -convert(T, 135.8025833007622) + h23 = convert(T, 64.51068857505875) + h24 = -convert(T, 4.192663174613162) + h25 = -convert(T, 2.53193205033506) + h31 = -convert(T, 14.99268484949843) + h32 = -convert(T, 76.30242396627033) + h33 = convert(T, 58.65928432851416) + h34 = convert(T, 16.61359034616402) + h35 = -convert(T, 0.6758691794084156) + + return Rodas4Tableau( + a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, + C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, C61, C62, C63, C64, C65, + γ, c2, c3, c4, d1, d2, d3, d4, + h21, h22, h23, h24, h25, h31, h32, h33, h34, h35 + ) +end + +struct Rodas3Tableau{T, T2} + a21::T + a31::T + a32::T + a41::T + a42::T + a43::T + C21::T + C31::T + C32::T + C41::T + C42::T + C43::T + γ::T2 + c2::T2 + c3::T2 + d1::T + d2::T + d3::T + d4::T +end + +function Rodas3Tableau(T, T2) + a21 = zero(T) + a31 = convert(T, 2) + a32 = zero(T) + a41 = convert(T, 2) + a42 = zero(T) + a43 = one(T) + + C21 = convert(T, 4) + C31 = one(T) + C32 = -one(T) + C41 = one(T) + C42 = -one(T) + C43 = -convert(T, 8 / 3) + + γ = convert(T2, 1 / 2) + c2 = zero(T2) + c3 = one(T2) + + d1 = convert(T, 1 / 2) + d2 = convert(T, 3 / 2) + d3 = zero(T) + d4 = zero(T) + + return Rodas3Tableau( + a21, a31, a32, a41, a42, a43, + C21, C31, C32, C41, C42, C43, + γ, c2, c3, d1, d2, d3, d4 + ) +end + struct Rodas5PTableau{T, T2} a21::T a31::T diff --git a/test/gpu_kernel_de/finite_diff.jl b/test/gpu_kernel_de/finite_diff.jl index 8f3527ed..ec94b152 100644 --- a/test/gpu_kernel_de/finite_diff.jl +++ b/test/gpu_kernel_de/finite_diff.jl @@ -20,7 +20,8 @@ monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy = false) osol = solve(prob, Rodas5P(), dt = 0.01f0, save_everystep = false) for alg in ( - GPURosenbrock23(autodiff = false), GPURodas4(autodiff = false), + GPURosenbrock23(autodiff = false), GPURodas3(autodiff = false), + GPURodas4(autodiff = false), GPURodas42(autodiff = false), GPURodas5P(autodiff = false), GPUKvaerno3(autodiff = false), GPUKvaerno5(autodiff = false), ) diff --git a/test/gpu_kernel_de/forward_diff.jl b/test/gpu_kernel_de/forward_diff.jl index 5e3f5d7f..7944e818 100644 --- a/test/gpu_kernel_de/forward_diff.jl +++ b/test/gpu_kernel_de/forward_diff.jl @@ -34,7 +34,8 @@ monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy = false) for alg in ( GPUTsit5(), GPUVern7(), GPUVern9(), GPURosenbrock23(autodiff = false), - GPURodas4(autodiff = false), GPURodas5P(autodiff = false), + GPURodas3(autodiff = false), GPURodas4(autodiff = false), + GPURodas42(autodiff = false), GPURodas5P(autodiff = false), GPUKvaerno3(autodiff = false), GPUKvaerno5(autodiff = false), ) @info alg diff --git a/test/gpu_kernel_de/stiff_ode/gpu_ode_modelingtoolkit_dae.jl b/test/gpu_kernel_de/stiff_ode/gpu_ode_modelingtoolkit_dae.jl index 596e2631..80cb5128 100644 --- a/test/gpu_kernel_de/stiff_ode/gpu_ode_modelingtoolkit_dae.jl +++ b/test/gpu_kernel_de/stiff_ode/gpu_ode_modelingtoolkit_dae.jl @@ -57,6 +57,18 @@ monteprob = SciMLBase.EnsembleProblem(prob, safetycopy = false) @test abs(sol.u[1][end][1] + sol.u[1][end][2] - 1.0f0) < 0.01f0 end +@testset "GPURodas3 DAE" begin + sol = solve( + monteprob, GPURodas3(), EnsembleGPUKernel(backend), + trajectories = 2, + dt = 0.001f0, + adaptive = false + ) + @test length(sol.u) == 2 + @test !any(isnan, sol.u[1][end]) + @test abs(sol.u[1][end][1] + sol.u[1][end][2] - 1.0f0) < 0.01f0 +end + @testset "GPURodas4 DAE" begin sol = solve( monteprob, GPURodas4(), EnsembleGPUKernel(backend), @@ -69,6 +81,18 @@ end @test abs(sol.u[1][end][1] + sol.u[1][end][2] - 1.0f0) < 0.01f0 end +@testset "GPURodas42 DAE" begin + sol = solve( + monteprob, GPURodas42(), EnsembleGPUKernel(backend), + trajectories = 2, + dt = 0.001f0, + adaptive = false + ) + @test length(sol.u) == 2 + @test !any(isnan, sol.u[1][end]) + @test abs(sol.u[1][end][1] + sol.u[1][end][2] - 1.0f0) < 0.01f0 +end + @testset "GPURodas5P DAE" begin sol = solve( monteprob, GPURodas5P(), EnsembleGPUKernel(backend), diff --git a/test/jet_tests.jl b/test/jet_tests.jl index dbe991cb..55949453 100644 --- a/test/jet_tests.jl +++ b/test/jet_tests.jl @@ -19,7 +19,9 @@ using ForwardDiff # Implicit algorithms (with default AD) @test_opt GPURosenbrock23() + @test_opt GPURodas3() @test_opt GPURodas4() + @test_opt GPURodas42() @test_opt GPURodas5P() @test_opt GPUKvaerno3() @test_opt GPUKvaerno5() @@ -27,7 +29,9 @@ using ForwardDiff # Implicit algorithms with explicit autodiff setting @test_opt GPURosenbrock23(; autodiff = Val{true}()) @test_opt GPURosenbrock23(; autodiff = Val{false}()) + @test_opt GPURodas3(; autodiff = Val{true}()) @test_opt GPURodas4(; autodiff = Val{true}()) + @test_opt GPURodas42(; autodiff = Val{false}()) @test_opt GPUKvaerno5(; autodiff = Val{false}()) end @@ -43,7 +47,9 @@ using ForwardDiff @test_opt DiffEqGPU.alg_order(GPUSIEA()) # Implicit algorithm orders @test_opt DiffEqGPU.alg_order(GPURosenbrock23()) + @test_opt DiffEqGPU.alg_order(GPURodas3()) @test_opt DiffEqGPU.alg_order(GPURodas4()) + @test_opt DiffEqGPU.alg_order(GPURodas42()) @test_opt DiffEqGPU.alg_order(GPURodas5P()) @test_opt DiffEqGPU.alg_order(GPUKvaerno3()) @test_opt DiffEqGPU.alg_order(GPUKvaerno5()) @@ -51,7 +57,9 @@ using ForwardDiff @testset "alg_autodiff type stability" begin @test_opt DiffEqGPU.alg_autodiff(GPURosenbrock23()) + @test_opt DiffEqGPU.alg_autodiff(GPURodas3()) @test_opt DiffEqGPU.alg_autodiff(GPURodas4()) + @test_opt DiffEqGPU.alg_autodiff(GPURodas42()) @test_opt DiffEqGPU.alg_autodiff(GPURodas5P()) @test_opt DiffEqGPU.alg_autodiff(GPUKvaerno3()) @test_opt DiffEqGPU.alg_autodiff(GPUKvaerno5()) @@ -128,6 +136,14 @@ using ForwardDiff GPURodas4(), f_test, false, u0, t0, dt, p, tstops, callback, save_everystep, saveat ) + @test_opt DiffEqGPU.init( + GPURodas3(), f_test, false, u0, t0, dt, p, tstops, + callback, save_everystep, saveat + ) + @test_opt DiffEqGPU.init( + GPURodas42(), f_test, false, u0, t0, dt, p, tstops, + callback, save_everystep, saveat + ) @test_opt DiffEqGPU.init( GPURodas5P(), f_test, false, u0, t0, dt, p, tstops, callback, save_everystep, saveat @@ -168,6 +184,14 @@ using ForwardDiff GPURodas4(), f_test, false, u0, t0, tf, dt, p, abstol, reltol, internalnorm, tstops, callback, saveat_adaptive ) + @test_opt DiffEqGPU.init( + GPURodas3(), f_test, false, u0, t0, tf, dt, p, + abstol, reltol, internalnorm, tstops, callback, saveat_adaptive + ) + @test_opt DiffEqGPU.init( + GPURodas42(), f_test, false, u0, t0, tf, dt, p, + abstol, reltol, internalnorm, tstops, callback, saveat_adaptive + ) @test_opt DiffEqGPU.init( GPURodas5P(), f_test, false, u0, t0, tf, dt, p, abstol, reltol, internalnorm, tstops, callback, saveat_adaptive @@ -192,8 +216,12 @@ using ForwardDiff @test_opt DiffEqGPU.Vern9Tableau(Float64, Float64) # Rodas tableaus + @test_opt DiffEqGPU.Rodas3Tableau(Float32, Float32) + @test_opt DiffEqGPU.Rodas3Tableau(Float64, Float64) @test_opt DiffEqGPU.Rodas4Tableau(Float32, Float32) @test_opt DiffEqGPU.Rodas4Tableau(Float64, Float64) + @test_opt DiffEqGPU.Rodas42Tableau(Float32, Float32) + @test_opt DiffEqGPU.Rodas42Tableau(Float64, Float64) @test_opt DiffEqGPU.Rodas5PTableau(Float32, Float32) @test_opt DiffEqGPU.Rodas5PTableau(Float64, Float64) @@ -210,6 +238,8 @@ using ForwardDiff @test_opt DiffEqGPU.build_adaptive_controller_cache(GPUVern7(), Float32) @test_opt DiffEqGPU.build_adaptive_controller_cache(GPUVern9(), Float32) @test_opt DiffEqGPU.build_adaptive_controller_cache(GPURosenbrock23(), Float32) + @test_opt DiffEqGPU.build_adaptive_controller_cache(GPURodas3(), Float32) @test_opt DiffEqGPU.build_adaptive_controller_cache(GPURodas4(), Float32) + @test_opt DiffEqGPU.build_adaptive_controller_cache(GPURodas42(), Float32) end end diff --git a/test/lower_level_api.jl b/test/lower_level_api.jl index b5e87f73..f4015037 100644 --- a/test/lower_level_api.jl +++ b/test/lower_level_api.jl @@ -56,7 +56,7 @@ gpu_probs = adapt(backend, probs) ## Finally use the lower API for faster solves! (Fixed time-stepping) -algs = (GPUTsit5(), GPUVern7(), GPUVern9(), GPURosenbrock23(), GPURodas4()) +algs = (GPUTsit5(), GPUVern7(), GPUVern9(), GPURosenbrock23(), GPURodas3(), GPURodas4(), GPURodas42()) for alg in algs @info alg