From 0dd83c06121077bab7158a6445c33af8d204d12d Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Sat, 12 Aug 2017 20:17:52 -0700 Subject: [PATCH 1/6] proper handling and testing of mass matrix in Rosenbrock23 --- src/OrdinaryDiffEq.jl | 3 - src/algorithms.jl | 69 +++++++++++++---------- src/initdt.jl | 14 ++++- src/integrators/rosenbrock_integrators.jl | 36 ++++++++++-- src/solve.jl | 4 +- test/mass_matrix_tests.jl | 33 +++++++++++ 6 files changed, 119 insertions(+), 40 deletions(-) create mode 100644 test/mass_matrix_tests.jl diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 212fad2fca..20f6721a44 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -120,9 +120,6 @@ module OrdinaryDiffEq # Reexport the Alg Types - export OrdinaryDiffEqAlgorithm, OrdinaryDiffEqAdaptiveAlgorithm, - OrdinaryDiffEqCompositeAlgorithm - export Discrete, FunctionMap, Euler, Heun, Ralston, Midpoint, SSPRK22, SSPRK33, SSPRK432, SSPRK104, RK4, ExplicitRK, OwrenZen3, OwrenZen4, OwrenZen5, BS3, BS5, diff --git a/src/algorithms.jl b/src/algorithms.jl index b787c69425..53afc0506b 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -1,8 +1,15 @@ abstract type OrdinaryDiffEqAlgorithm <: AbstractODEAlgorithm end abstract type OrdinaryDiffEqAdaptiveAlgorithm <: OrdinaryDiffEqAlgorithm end -abstract type OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} <: OrdinaryDiffEqAdaptiveAlgorithm end abstract type OrdinaryDiffEqCompositeAlgorithm <: OrdinaryDiffEqAlgorithm end +abstract type OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS,AD} <: OrdinaryDiffEqAdaptiveAlgorithm end +abstract type OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} <: OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS,AD} end +abstract type OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} <: OrdinaryDiffEqAdaptiveImplicitAlgorithm{CS,AD} end + +abstract type OrdinaryDiffEqImplicitAlgorithm{CS,AD} <: OrdinaryDiffEqAlgorithm end +abstract type OrdinaryDiffEqNewtonAlgorithm{CS,AD} <: OrdinaryDiffEqImplicitAlgorithm{CS,AD} end +abstract type OrdinaryDiffEqRosenbrockAlgorithm{CS,AD} <: OrdinaryDiffEqImplicitAlgorithm{CS,AD} end + struct Discrete{apply_map,scale_by_time} <: OrdinaryDiffEqAlgorithm end Base.@pure Discrete(;apply_map=false,scale_by_time=false) = Discrete{apply_map,scale_by_time}() @@ -109,7 +116,7 @@ struct StrangSplitting <: OrdinaryDiffEqAlgorithm end # SDIRK Methods -struct ImplicitEuler{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct ImplicitEuler{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -129,7 +136,7 @@ Base.@pure ImplicitEuler(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct Trapezoid{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct Trapezoid{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -149,7 +156,7 @@ Base.@pure Trapezoid(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct TRBDF2{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct TRBDF2{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -170,7 +177,7 @@ TRBDF2{chunk_size,autodiff,typeof(linsolve), linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct SDIRK2{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct SDIRK2{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -191,7 +198,7 @@ Base.@pure SDIRK2(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct SSPSDIRK2{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqAlgorithm # Not adaptive +struct SSPSDIRK2{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAlgorithm{CS,AD} # Not adaptive linsolve::F diff_type::Symbol κ::K @@ -212,7 +219,7 @@ Base.@pure SSPSDIRK2(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct Kvaerno3{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct Kvaerno3{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -233,7 +240,7 @@ Base.@pure Kvaerno3(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct KenCarp3{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct KenCarp3{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -255,7 +262,7 @@ Base.@pure KenCarp3(;chunk_size=0,autodiff=true,diff_type=:central, max_newton_iter,new_jac_conv_bound) -struct Cash4{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct Cash4{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -277,7 +284,7 @@ Base.@pure Cash4(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound,embedding) -struct Hairer4{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct Hairer4{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -298,7 +305,7 @@ Base.@pure Hairer4(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct Hairer42{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct Hairer42{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -319,7 +326,7 @@ Base.@pure Hairer42(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct Kvaerno4{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct Kvaerno4{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -340,7 +347,7 @@ Base.@pure Kvaerno4(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct Kvaerno5{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct Kvaerno5{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -361,7 +368,7 @@ Base.@pure Kvaerno5(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct KenCarp4{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct KenCarp4{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -382,7 +389,7 @@ Base.@pure KenCarp4(;chunk_size=0,autodiff=true,diff_type=:central, linsolve,diff_type,κ,tol,smooth_est,extrapolant,min_newton_iter, max_newton_iter,new_jac_conv_bound) -struct KenCarp5{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{Controller} +struct KenCarp5{CS,AD,F,K,T,T2,Controller} <: OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller} linsolve::F diff_type::Symbol κ::K @@ -407,91 +414,91 @@ Base.@pure KenCarp5(;chunk_size=0,autodiff=true,diff_type=:central, # Rosenbrock Methods -struct Rosenbrock23{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Rosenbrock23{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Rosenbrock23(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Rosenbrock23{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct Rosenbrock32{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Rosenbrock32{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Rosenbrock32(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Rosenbrock32{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct ROS3P{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct ROS3P{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure ROS3P(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = ROS3P{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct Rodas3{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Rodas3{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Rodas3(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Rodas3{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct RosShamp4{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct RosShamp4{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure RosShamp4(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = RosShamp4{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct Veldd4{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Veldd4{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Veldd4(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Veldd4{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct Velds4{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Velds4{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Velds4(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Velds4{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct GRK4T{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct GRK4T{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure GRK4T(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = GRK4T{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct GRK4A{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct GRK4A{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure GRK4A(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = GRK4A{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct Ros4LStab{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Ros4LStab{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Ros4LStab(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Ros4LStab{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct Rodas4{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Rodas4{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Rodas4(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Rodas4{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct Rodas42{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Rodas42{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Rodas42(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Rodas42{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct Rodas4P{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Rodas4P{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Rodas4P(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Rodas4P{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct Rodas5{CS,AD,F} <: OrdinaryDiffEqAdaptiveAlgorithm +struct Rodas5{CS,AD,F} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} linsolve::F diff_type::Symbol end Base.@pure Rodas5(;chunk_size=0,autodiff=true,diff_type=:central,linsolve=DEFAULT_LINSOLVE) = Rodas5{chunk_size,autodiff,typeof(linsolve)}(linsolve,diff_type) -struct GeneralRosenbrock{CS,AD,F,TabType} <: OrdinaryDiffEqAdaptiveAlgorithm +struct GeneralRosenbrock{CS,AD,F,TabType} <: OrdinaryDiffEqRosenbrockAdaptiveAlgorithm{CS,AD} tableau::TabType factorization::F end @@ -516,6 +523,8 @@ struct LawsonEuler <: OrdinaryDiffEqAlgorithm end struct NorsettEuler <: OrdinaryDiffEqAlgorithm end struct SplitEuler <: OrdinaryDiffEqAlgorithm end +######################################### + struct CompositeAlgorithm{T,F} <: OrdinaryDiffEqCompositeAlgorithm algs::T choice_function::F diff --git a/src/initdt.jl b/src/initdt.jl index c09f8ef826..5cf4109647 100644 --- a/src/initdt.jl +++ b/src/initdt.jl @@ -1,4 +1,4 @@ -@muladd function ode_determine_initdt{tType,uType}(u0,t::tType,tdir,dtmax,abstol,reltol,internalnorm,prob::AbstractODEProblem{uType,tType,true},order) +@muladd function ode_determine_initdt{tType,uType}(u0,t::tType,tdir,dtmax,abstol,reltol,internalnorm,prob::AbstractODEProblem{uType,tType,true},order,alg) f = prob.f oneunit_tType = oneunit(tType) dtmax_tdir = tdir*dtmax @@ -9,6 +9,13 @@ f₀ = zeros(u0./t) f(t,u0,f₀) + + if prob.mass_matrix != I + ftmp = similar(f₀) + alg.linsolve(ftmp, copy(prob.mass_matrix), f₀, true) + f₀ .= ftmp + end + if any((isnan(x) for x in f₀)) warn("First function call produced NaNs. Exiting.") end @@ -29,6 +36,11 @@ f₁ = similar(f₀) f(t+dt₀_tdir,u₁,f₁) + if prob.mass_matrix != I + alg.linsolve(ftmp, prob.mass_matrix, f₁, false) + f₁ .= ftmp + end + @. tmp = (f₁-f₀)/sk*oneunit_tType d₂ = internalnorm(tmp)/dt₀*oneunit_tType # Hairer has d₂ = sqrt(sum(abs2,tmp))/dt₀, note the lack of norm correction diff --git a/src/integrators/rosenbrock_integrators.jl b/src/integrators/rosenbrock_integrators.jl index 82619ad529..cbe615a0b5 100644 --- a/src/integrators/rosenbrock_integrators.jl +++ b/src/integrators/rosenbrock_integrators.jl @@ -71,7 +71,13 @@ end @. u = uprev + dto2*k₁ f(t+dto2, u, f₁) - @. linsolve_tmp = f₁ - k₁ + if mass_matrix == I + tmp .= k₁ + else + A_mul_B!(tmp,mass_matrix,k₁) # vectmp == k₁ + end + + @. linsolve_tmp = f₁ - tmp if has_invW(f) A_mul_B!(vectmp2, W, linsolve_tmp_vec) else @@ -85,7 +91,15 @@ end if integrator.opts.adaptive f(t+dt, u, fsallast) - @. linsolve_tmp = fsallast - c₃₂*(k₂-f₁) - 2(k₁-fsalfirst) + dt*dT + if mass_matrix == I + @. linsolve_tmp = fsallast - c₃₂*(k₂-f₁) - 2(k₁-fsalfirst) + dt*dT + else + @. du2 = c₃₂*k₂ + 2k₁ + A_mul_B!(du1,mass_matrix,du2) + @. linsolve_tmp = fsallast - du1 + c₃₂*f₁ + 2fsalfirst + dt*dT + end + + if has_invW(f) A_mul_B!(vectmp3, W, linsolve_tmp_vec) else @@ -173,7 +187,14 @@ end @. u = uprev + dto2*k₁ f(t+dto2, u, f₁) - @. linsolve_tmp = f₁ - k₁ + if mass_matrix == I + tmp .= k₁ + else + A_mul_B!(tmp,mass_matrix,k₁) # vectmp == k₁ + end + + @. linsolve_tmp = f₁ - tmp + if has_invW(f) A_mul_B!(vectmp2, W, linsolve_tmp_vec) else @@ -185,7 +206,14 @@ end @. tmp = uprev + dt*k₂ f(t+dt, tmp, fsallast) - @. linsolve_tmp = fsallast - c₃₂*(k₂-f₁) - 2(k₁-fsalfirst) + dt*dT + if mass_matrix == I + @. linsolve_tmp = fsallast - c₃₂*(k₂-f₁) - 2(k₁-fsalfirst) + dt*dT + else + @. du2 = c₃₂*k₂ + 2k₁ + A_mul_B!(du1,mass_matrix,du2) + @. linsolve_tmp = fsallast - du1 + c₃₂*f₁ + 2fsalfirst + dt*dT + end + if has_invW(f) A_mul_B!(vectmp3, W, linsolve_tmp_vec) else diff --git a/src/solve.jl b/src/solve.jl index 1ab168fb76..3473b56938 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -52,7 +52,7 @@ function init{algType<:OrdinaryDiffEqAlgorithm,recompile_flag}( error("This solver is not able to use mass matrices.") end elseif !(typeof(prob)<:DiscreteProblem) && - !(typeof(alg) <: Union{Rosenbrock23,Rosenbrock32}) && + !(typeof(alg) <: Union{OrdinaryDiffEqAdaptiveImplicitAlgorithm,OrdinaryDiffEqImplicitAlgorithm}) && prob.mass_matrix != I error("This solver is not able to use mass matrices.") end @@ -126,7 +126,7 @@ function init{algType<:OrdinaryDiffEqAlgorithm,recompile_flag}( dtmax > zero(dtmax) && tdir < 0 && (dtmax *= tdir) # Allow positive dtmax, but auto-convert # dtmin is all abs => does not care about sign already. if dt == zero(dt) && adaptive - dt = tType(ode_determine_initdt(u,t,tdir,dtmax,abstol_internal,reltol_internal,internalnorm,prob,order)) + dt = tType(ode_determine_initdt(u,t,tdir,dtmax,abstol_internal,reltol_internal,internalnorm,prob,order,alg)) if sign(dt)!=tdir && dt!=tType(0) && !isnan(dt) error("Automatic dt setting has the wrong sign. Exiting. Please report this error.") end diff --git a/test/mass_matrix_tests.jl b/test/mass_matrix_tests.jl new file mode 100644 index 0000000000..6e52e4c418 --- /dev/null +++ b/test/mass_matrix_tests.jl @@ -0,0 +1,33 @@ +using OrdinaryDiffEq, Base.Test + +A = [-2.0 1 4 + 4 -2 1 + 2 1 3] + +function f(t,u,du) + A_mul_B!(du,A,u) +end +function f(::Type{Val{:analytic}},t,u0) + u0.*exp(t) +end + +function g(t,u,du) + du.=u +end +function g(::Type{Val{:analytic}},t,u0) + u0.*exp(t) +end +prob2 = ODEProblem(g,ones(3),(0.0,1.0)) +prob = ODEProblem(f,ones(3),(0.0,1.0),mass_matrix=A) + +######################################### Test each method for exactness + +sol = solve(prob, Rosenbrock23()) +sol2 = solve(prob2,Rosenbrock23()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-11 + +sol = solve(prob, Rosenbrock32()) +sol2 = solve(prob2,Rosenbrock32()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-11 From 041f7ba01da1cbf1fd036ed2bba4517b91d5bdce Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Sat, 12 Aug 2017 22:16:35 -0700 Subject: [PATCH 2/6] all Rosenbrock MM --- src/algorithms.jl | 1 - src/integrators/rosenbrock_integrators.jl | 221 ++++++++++++++++++---- test/mass_matrix_tests.jl | 54 +++++- 3 files changed, 233 insertions(+), 43 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 53afc0506b..067565cd1c 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -55,7 +55,6 @@ struct Feagin14 <: OrdinaryDiffEqAdaptiveAlgorithm end # Symplectic methods -#struct Verlet <: OrdinaryDiffEqAlgorithm end struct SymplecticEuler <: OrdinaryDiffEqAlgorithm end struct VelocityVerlet <: OrdinaryDiffEqAlgorithm end struct VerletLeapfrog <: OrdinaryDiffEqAlgorithm end diff --git a/src/integrators/rosenbrock_integrators.jl b/src/integrators/rosenbrock_integrators.jl index cbe615a0b5..c2b9c3ff6d 100644 --- a/src/integrators/rosenbrock_integrators.jl +++ b/src/integrators/rosenbrock_integrators.jl @@ -74,7 +74,7 @@ end if mass_matrix == I tmp .= k₁ else - A_mul_B!(tmp,mass_matrix,k₁) # vectmp == k₁ + A_mul_B!(tmp,mass_matrix,k₁) end @. linsolve_tmp = f₁ - tmp @@ -190,11 +190,11 @@ end if mass_matrix == I tmp .= k₁ else - A_mul_B!(tmp,mass_matrix,k₁) # vectmp == k₁ + A_mul_B!(tmp,mass_matrix,k₁) end @. linsolve_tmp = f₁ - tmp - + if has_invW(f) A_mul_B!(vectmp2, W, linsolve_tmp_vec) else @@ -498,7 +498,14 @@ end @. u = uprev + a21*k1 f(t+c2*dt, u, du) - @. linsolve_tmp = du + dtd2*dT + dtC21*k1 + if mass_matrix == I + @. linsolve_tmp = du + dtd2*dT + dtC21*k1 + else + @. du1 = dtC21*k1 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd2*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp2, W, linsolve_tmp_vec) else @@ -509,7 +516,14 @@ end @. u = uprev + a31*k1 + a32*k2 f(t+c3*dt, u, du) - @. linsolve_tmp = du + dtd3*dT + dtC31*k1 + dtC32*k2 + if mass_matrix == I + @. linsolve_tmp = du + dtd3*dT + dtC31*k1 + dtC32*k2 + else + @. du1 = dtC31*k1 + dtC32*k2 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd3*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp3, W, linsolve_tmp_vec) else @@ -698,7 +712,14 @@ end f(t+c2*dt,u,du) =# - @. linsolve_tmp = fsalfirst + dtd2*dT + dtC21*k1 + if mass_matrix == I + @. linsolve_tmp = fsalfirst + dtd2*dT + dtC21*k1 + else + @. du1 = dtC21*k1 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = fsalfirst + dtd2*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp2, W, linsolve_tmp_vec) else @@ -709,7 +730,14 @@ end @. u = uprev + a31*k1 + a32*k2 f(t+c3*dt, u, du) - @. linsolve_tmp = du + dtd3*dT + dtC31*k1 + dtC32*k2 + if mass_matrix == I + @. linsolve_tmp = du + dtd3*dT + dtC31*k1 + dtC32*k2 + else + @. du1 = dtC31*k1 + dtC32*k2 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd3*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp3, W, linsolve_tmp_vec) else @@ -717,7 +745,14 @@ end end k3 = reshape(vectmp3, sizeu...) - @. linsolve_tmp = du + dtd4*dT + dtC41*k1 + dtC42*k2 + dtC43*k3 + + if mass_matrix == I + @. linsolve_tmp = du + dtd4*dT + dtC41*k1 + dtC42*k2 + dtC43*k3 + else + @. du1 = dtC41*k1 + dtC42*k2 + dtC43*k3 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd4*dT + du2 + end if has_invW(f) A_mul_B!(vectmp4, W, linsolve_tmp_vec) @@ -729,6 +764,8 @@ end @. u = uprev + b1*k1 + b2*k2 + b3*k3 + b4*k4 f(t, u, fsallast) + #@show k1,k2,k3,k4 + if integrator.opts.adaptive @. utilde = btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4 calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol, @@ -905,7 +942,14 @@ end @. u = uprev + a21*k1 f(t+c2*dt, u, du) - @. linsolve_tmp = du + dtd2*dT + dtC21*k1 + if mass_matrix == I + @. linsolve_tmp = du + dtd2*dT + dtC21*k1 + else + @. du1 = dtC21*k1 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd2*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp2, W, linsolve_tmp_vec) else @@ -916,7 +960,14 @@ end @. u = uprev + a31*k1 + a32*k2 f(t+c3*dt, u, du) - @. linsolve_tmp = du + dtd3*dT + dtC31*k1 + dtC32*k2 + if mass_matrix == I + @. linsolve_tmp = du + dtd3*dT + dtC31*k1 + dtC32*k2 + else + @. du1 = dtC31*k1 + dtC32*k2 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd3*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp3, W, linsolve_tmp_vec) else @@ -925,7 +976,14 @@ end k3 = reshape(vectmp3, sizeu...) - @. linsolve_tmp = du + dtd4*dT + dtC41*k1 + dtC42*k2 + dtC43*k3 + if mass_matrix == I + @. linsolve_tmp = du + dtd4*dT + dtC41*k1 + dtC42*k2 + dtC43*k3 + else + @. du1 = dtC41*k1 + dtC42*k2 + dtC43*k3 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd4*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp4, W, linsolve_tmp_vec) else @@ -1150,7 +1208,14 @@ end @. u = uprev + a21*k1 f(t+c2*dt, u, du) - @. linsolve_tmp = du + dtd2*dT + dtC21*k1 + if mass_matrix == I + @. linsolve_tmp = du + dtd2*dT + dtC21*k1 + else + @. du1 = dtC21*k1 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd2*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp2, W, linsolve_tmp_vec) else @@ -1161,7 +1226,14 @@ end @. u = uprev + a31*k1 + a32*k2 f(t+c3*dt, u, du) - @. linsolve_tmp = du + dtd3*dT + (dtC31*k1 + dtC32*k2) + if mass_matrix == I + @. linsolve_tmp = du + dtd3*dT + (dtC31*k1 + dtC32*k2) + else + @. du1 = dtC31*k1 + dtC32*k2 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd3*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp3, W, linsolve_tmp_vec) else @@ -1172,7 +1244,13 @@ end @. u = uprev + a41*k1 + a42*k2 + a43*k3 f(t+c4*dt, u, du) - @. linsolve_tmp = du + dtd4*dT + (dtC41*k1 + dtC42*k2 + dtC43*k3) + if mass_matrix == I + @. linsolve_tmp = du + dtd4*dT + (dtC41*k1 + dtC42*k2 + dtC43*k3) + else + @. du1 = dtC41*k1 + dtC42*k2 + dtC43*k3 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd4*dT + du2 + end if has_invW(f) A_mul_B!(vectmp4, W, linsolve_tmp_vec) @@ -1184,7 +1262,14 @@ end @. u = uprev + a51*k1 + a52*k2 + a53*k3 + a54*k4 f(t+dt, u, du) - @. linsolve_tmp = du + (dtC52*k2 + dtC54*k4 + dtC51*k1 + dtC53*k3) + if mass_matrix == I + @. linsolve_tmp = du + (dtC52*k2 + dtC54*k4 + dtC51*k1 + dtC53*k3) + else + @. du1 = dtC52*k2 + dtC54*k4 + dtC51*k1 + dtC53*k3 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + du2 + end + if has_invW(f) A_mul_B!(vectmp5, W, linsolve_tmp_vec) else @@ -1195,10 +1280,20 @@ end u .+= k5 f(t, u, du) - # @. linsolve_tmp = du + (dtC61*k1 + dtC62*k2 + dtC65*k5 + dtC64*k4 + dtC63*k3) - @tight_loop_macros for i in uidx - @inbounds linsolve_tmp[i] = du[i] + (dtC61*k1[i] + dtC62*k2[i] + dtC65*k5[i] + dtC64*k4[i] + dtC63*k3[i]) + if mass_matrix == I + # @. linsolve_tmp = du + (dtC61*k1 + dtC62*k2 + dtC65*k5 + dtC64*k4 + dtC63*k3) + @tight_loop_macros for i in uidx + @inbounds linsolve_tmp[i] = du[i] + (dtC61*k1[i] + dtC62*k2[i] + dtC65*k5[i] + dtC64*k4[i] + dtC63*k3[i]) + end + else + # @. linsolve_tmp = dtC61*k1 + dtC62*k2 + dtC65*k5 + dtC64*k4 + dtC63*k3 + @tight_loop_macros for i in uidx + @inbounds du1[i] = dtC61*k1[i] + dtC62*k2[i] + dtC65*k5[i] + dtC64*k4[i] + dtC63*k3[i] + end + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + du2 end + if has_invW(f) A_mul_B!(vectmp6, W, linsolve_tmp_vec) else @@ -1499,7 +1594,16 @@ end @. u = uprev + a21*k1 f(t+c2*dt, u, du) - @. linsolve_tmp = du + dtd2*dT + dtC21*k1 + + + if mass_matrix == I + @. linsolve_tmp = du + dtd2*dT + dtC21*k1 + else + @. du1 = dtC21*k1 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd2*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp2, W, linsolve_tmp_vec) else @@ -1510,7 +1614,14 @@ end @. u = uprev + a31*k1 + a32*k2 f(t+c3*dt, u, du) - @. linsolve_tmp = du + dtd3*dT + (dtC31*k1 + dtC32*k2) + if mass_matrix == I + @. linsolve_tmp = du + dtd3*dT + (dtC31*k1 + dtC32*k2) + else + @. du1 = dtC31*k1 + dtC32*k2 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd3*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp3, W, linsolve_tmp_vec) else @@ -1521,7 +1632,14 @@ end @. u = uprev + a41*k1 + a42*k2 + a43*k3 f(t+c4*dt, u, du) - @. linsolve_tmp = du + dtd4*dT + (dtC41*k1 + dtC42*k2 + dtC43*k3) + if mass_matrix == I + @. linsolve_tmp = du + dtd4*dT + (dtC41*k1 + dtC42*k2 + dtC43*k3) + else + @. du1 = dtC41*k1 + dtC42*k2 + dtC43*k3 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd4*dT + du2 + end + if has_invW(f) A_mul_B!(vectmp4, W, linsolve_tmp_vec) else @@ -1532,10 +1650,17 @@ end @. u = uprev + a51*k1 + a52*k2 + a53*k3 + a54*k4 f(t+c5*dt, u, du) -# @. linsolve_tmp = du + dtd5*dT + (dtC52*k2 + dtC54*k4 + dtC51*k1 + dtC53*k3) - @tight_loop_macros for i in uidx - @inbounds linsolve_tmp[i] = du[i] + dtd5*dT[i] + (dtC52*k2[i] + dtC54*k4[i] + dtC51*k1[i] + dtC53*k3[i]) + if mass_matrix == I + # @. linsolve_tmp = du + dtd5*dT + (dtC52*k2 + dtC54*k4 + dtC51*k1 + dtC53*k3) + @tight_loop_macros for i in uidx + @inbounds linsolve_tmp[i] = du[i] + dtd5*dT[i] + (dtC52*k2[i] + dtC54*k4[i] + dtC51*k1[i] + dtC53*k3[i]) + end + else + @. du1 = dtC52*k2 + dtC54*k4 + dtC51*k1 + dtC53*k3 + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + dtd5*dT + du2 end + if has_invW(f) A_mul_B!(vectmp5, W, linsolve_tmp_vec) else @@ -1549,10 +1674,20 @@ end end f(t+dt, u, du) - # @. linsolve_tmp = du + (dtC61*k1 + dtC62*k2 + dtC63*k3 + dtC64*k4 + dtC65*k5) - @tight_loop_macros for i in uidx - @inbounds linsolve_tmp[i] = du[i] + (dtC61*k1[i] + dtC62*k2[i] + dtC63*k3[i] + dtC64*k4[i] + dtC65*k5[i]) + if mass_matrix == I + # @. linsolve_tmp = du + (dtC61*k1 + dtC62*k2 + dtC63*k3 + dtC64*k4 + dtC65*k5) + @tight_loop_macros for i in uidx + @inbounds linsolve_tmp[i] = du[i] + (dtC61*k1[i] + dtC62*k2[i] + dtC63*k3[i] + dtC64*k4[i] + dtC65*k5[i]) + end + else + # @. du1 = dtC61*k1 + dtC62*k2 + dtC63*k3 + dtC64*k4 + dtC65*k5 + @tight_loop_macros for i in uidx + @inbounds du1[i] = dtC61*k1[i] + dtC62*k2[i] + dtC63*k3[i] + dtC64*k4[i] + dtC65*k5[i] + end + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + du2 end + if has_invW(f) A_mul_B!(vectmp6, W, linsolve_tmp_vec) else @@ -1563,10 +1698,20 @@ end u .+= k6 f(t+dt, u, du) - # @. linsolve_tmp = du + (dtC71*k1 + dtC72*k2 + dtC73*k3 + dtC74*k4 + dtC75*k5 + dtC76*k6) - @tight_loop_macros for i in uidx - @inbounds linsolve_tmp[i] = du[i] + (dtC71*k1[i] + dtC72*k2[i] + dtC73*k3[i] + dtC74*k4[i] + dtC75*k5[i] + dtC76*k6[i]) + if mass_matrix == I + # @. linsolve_tmp = du + (dtC71*k1 + dtC72*k2 + dtC73*k3 + dtC74*k4 + dtC75*k5 + dtC76*k6) + @tight_loop_macros for i in uidx + @inbounds linsolve_tmp[i] = du[i] + (dtC71*k1[i] + dtC72*k2[i] + dtC73*k3[i] + dtC74*k4[i] + dtC75*k5[i] + dtC76*k6[i]) + end + else + # @. du1 =dtC72*k2 + dtC73*k3 + dtC74*k4 + dtC75*k5 + dtC76*k6 + @tight_loop_macros for i in uidx + @inbounds du1[i] = dtC71*k1[i] + dtC72*k2[i] + dtC73*k3[i] + dtC74*k4[i] + dtC75*k5[i] + dtC76*k6[i] + end + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + du2 end + if has_invW(f) A_mul_B!(vectmp7, W, linsolve_tmp_vec) else @@ -1577,10 +1722,20 @@ end u .+= k7 f(t+dt, u, du) - # @. linsolve_tmp = du + (dtC81*k1 + dtC82*k2 + dtC83*k3 + dtC84*k4 + dtC85*k5 + dtC86*k6 + dtC87*k7) - @tight_loop_macros for i in uidx - @inbounds linsolve_tmp[i] = du[i] + (dtC81*k1[i] + dtC82*k2[i] + dtC83*k3[i] + dtC84*k4[i] + dtC85*k5[i] + dtC86*k6[i] + dtC87*k7[i]) + if mass_matrix == I + # @. linsolve_tmp = du + (dtC81*k1 + dtC82*k2 + dtC83*k3 + dtC84*k4 + dtC85*k5 + dtC86*k6 + dtC87*k7) + @tight_loop_macros for i in uidx + @inbounds linsolve_tmp[i] = du[i] + (dtC81*k1[i] + dtC82*k2[i] + dtC83*k3[i] + dtC84*k4[i] + dtC85*k5[i] + dtC86*k6[i] + dtC87*k7[i]) + end + else + # @. du1 = dtC81*k1 + dtC82*k2 + dtC83*k3 + dtC84*k4 + dtC85*k5 + dtC86*k6 + dtC87*k7 + @tight_loop_macros for i in uidx + @inbounds du1[i] = dtC81*k1[i] + dtC82*k2[i] + dtC83*k3[i] + dtC84*k4[i] + dtC85*k5[i] + dtC86*k6[i] + dtC87*k7[i] + end + A_mul_B!(du2,mass_matrix,du1) + @. linsolve_tmp = du + du2 end + if has_invW(f) A_mul_B!(vectmp8, W, linsolve_tmp_vec) else diff --git a/test/mass_matrix_tests.jl b/test/mass_matrix_tests.jl index 6e52e4c418..0079d83d51 100644 --- a/test/mass_matrix_tests.jl +++ b/test/mass_matrix_tests.jl @@ -1,28 +1,29 @@ using OrdinaryDiffEq, Base.Test -A = [-2.0 1 4 - 4 -2 1 - 2 1 3] - +const A = [-2.0 1 4 + 4 -2 1 + 2 1 3] +const b = A*ones(3) function f(t,u,du) A_mul_B!(du,A,u) + tmp = t*b + du .+= tmp end function f(::Type{Val{:analytic}},t,u0) - u0.*exp(t) + @. 2ones(3)*exp(t) - t - 1 end - function g(t,u,du) - du.=u + du .= u + t end function g(::Type{Val{:analytic}},t,u0) - u0.*exp(t) + @. 2ones(3)*exp(t) - t - 1 end prob2 = ODEProblem(g,ones(3),(0.0,1.0)) prob = ODEProblem(f,ones(3),(0.0,1.0),mass_matrix=A) ######################################### Test each method for exactness -sol = solve(prob, Rosenbrock23()) +sol = solve(prob ,Rosenbrock23()) sol2 = solve(prob2,Rosenbrock23()) @test norm(sol .- sol2) ≈ 0 atol=1e-11 @@ -31,3 +32,38 @@ sol = solve(prob, Rosenbrock32()) sol2 = solve(prob2,Rosenbrock32()) @test norm(sol .- sol2) ≈ 0 atol=1e-11 + +sol = solve(prob, ROS3P()) +sol2 = solve(prob2,ROS3P()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-11 + +sol = solve(prob, Rodas3()) +sol2 = solve(prob2,Rodas3()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-11 + +sol = solve(prob, RosShamp4()) +sol2 = solve(prob2,RosShamp4()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-10 + +sol = solve(prob, Rodas4()) +sol2 = solve(prob2,Rodas4()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-11 + +sol = solve(prob, Rodas5()) +sol2 = solve(prob2,Rodas5()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-9 + +#sol = solve(prob, SDIRK2()) +#sol2 = solve(prob2, SDIRK2()) + +#@test norm(sol .- sol2) ≈ 0 atol=1e-11 + + + +#sol = solve(prob, Rodas3(),adaptive=false,dt=1/10) +#sol2 = solve(prob2, Rodas3(),adaptive=false,dt=1/10) From a4e732c451f671ec188b7b6f052a60f3820008c3 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Sun, 13 Aug 2017 00:28:26 -0700 Subject: [PATCH 3/6] implicit euler mass matrix --- src/integrators/sdirk_integrators.jl | 15 ++++++++++++--- test/mass_matrix_tests.jl | 5 +++++ 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/src/integrators/sdirk_integrators.jl b/src/integrators/sdirk_integrators.jl index d9c81f40ce..f55467b748 100644 --- a/src/integrators/sdirk_integrators.jl +++ b/src/integrators/sdirk_integrators.jl @@ -161,7 +161,12 @@ end iter += 1 f(t+dt,u,k) scale!(k,dt) - k .-= z + if mass_matrix == I + k .-= z + else + A_mul_B!(du1,mass_matrix,z) + k .-= du1 + end if has_invW(f) A_mul_B!(vec(dz),W,vec(k)) # Here W is actually invW else @@ -183,7 +188,12 @@ end iter += 1 f(t+dt,u,k) scale!(k,dt) - k .-= z + if mass_matrix == I + k .-= z + else + A_mul_B!(du1,mass_matrix,z) + k .-= du1 + end if has_invW(f) A_mul_B!(dz,W,k) # Here W is actually invW else @@ -209,7 +219,6 @@ end cache.ηold = η cache.newton_iters = iter - @. u = uprev + z if integrator.opts.adaptive && integrator.success_iter > 0 # Use 2rd divided differences a la SPICE and Shampine diff --git a/test/mass_matrix_tests.jl b/test/mass_matrix_tests.jl index 0079d83d51..23fe9f27a3 100644 --- a/test/mass_matrix_tests.jl +++ b/test/mass_matrix_tests.jl @@ -58,6 +58,11 @@ sol2 = solve(prob2,Rodas5()) @test norm(sol .- sol2) ≈ 0 atol=1e-9 +sol = solve(prob, ImplicitEuler()) +sol2 = solve(prob2,ImplicitEuler()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-9 + #sol = solve(prob, SDIRK2()) #sol2 = solve(prob2, SDIRK2()) From ed6cfca6b11f2d3c71deb623bdff971775db072b Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Sun, 13 Aug 2017 01:10:14 -0700 Subject: [PATCH 4/6] Trapezoid mass matrix --- src/integrators/sdirk_integrators.jl | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/integrators/sdirk_integrators.jl b/src/integrators/sdirk_integrators.jl index f55467b748..05b4166aee 100644 --- a/src/integrators/sdirk_integrators.jl +++ b/src/integrators/sdirk_integrators.jl @@ -406,7 +406,12 @@ end iter += 1 f(t+dto2,u,k) scale!(k,dto2) - k .-= z + if mass_matrix == I + k .-= z + else + A_mul_B!(du1,mass_matrix,z) + k .-= du1 + end if has_invW(f) A_mul_B!(vec(dz),W,vec(k)) # Here W is actually invW else @@ -428,7 +433,12 @@ end iter += 1 f(t+dto2,u,k) scale!(k,dto2) - k .-= z + if mass_matrix == I + k .-= z + else + A_mul_B!(du1,mass_matrix,z) + k .-= du1 + end if has_invW(f) A_mul_B!(vec(dz),W,vec(k)) # Here W is actually invW else From eb0fd35ed4e8d75c3a437e6e73c034ca7dd3eb01 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Sun, 13 Aug 2017 03:16:15 -0700 Subject: [PATCH 5/6] fix mass matrix tests --- src/initdt.jl | 2 +- test/mass_matrix_tests.jl | 39 +++++++++++++++++++++++++------------ test/ode/ode_cache_tests.jl | 2 +- test/runtests.jl | 1 + 4 files changed, 30 insertions(+), 14 deletions(-) diff --git a/src/initdt.jl b/src/initdt.jl index 5cf4109647..aa8bf24869 100644 --- a/src/initdt.jl +++ b/src/initdt.jl @@ -54,7 +54,7 @@ dt = tdir*min(100dt₀,dt₁,dtmax_tdir) end -@muladd function ode_determine_initdt{uType,tType}(u0::uType,t,tdir,dtmax,abstol,reltol,internalnorm,prob::AbstractODEProblem{uType,tType,false},order) +@muladd function ode_determine_initdt{uType,tType}(u0::uType,t,tdir,dtmax,abstol,reltol,internalnorm,prob::AbstractODEProblem{uType,tType,false},order,alg) f = prob.f oneunit_tType = oneunit(tType) dtmax_tdir = tdir*dtmax diff --git a/test/mass_matrix_tests.jl b/test/mass_matrix_tests.jl index 23fe9f27a3..366bf052e6 100644 --- a/test/mass_matrix_tests.jl +++ b/test/mass_matrix_tests.jl @@ -1,25 +1,25 @@ using OrdinaryDiffEq, Base.Test -const A = [-2.0 1 4 +const mm_A = [-2.0 1 4 4 -2 1 2 1 3] -const b = A*ones(3) -function f(t,u,du) - A_mul_B!(du,A,u) - tmp = t*b +const mm_b = mm_A*ones(3) +function mm_f(t,u,du) + A_mul_B!(du,mm_A,u) + tmp = t*mm_b du .+= tmp end -function f(::Type{Val{:analytic}},t,u0) +function mm_f(::Type{Val{:analytic}},t,u0) @. 2ones(3)*exp(t) - t - 1 end -function g(t,u,du) +function mm_g(t,u,du) du .= u + t end -function g(::Type{Val{:analytic}},t,u0) +function mm_g(::Type{Val{:analytic}},t,u0) @. 2ones(3)*exp(t) - t - 1 end -prob2 = ODEProblem(g,ones(3),(0.0,1.0)) -prob = ODEProblem(f,ones(3),(0.0,1.0),mass_matrix=A) +prob2 = ODEProblem(mm_g,ones(3),(0.0,1.0)) +prob = ODEProblem(mm_f,ones(3),(0.0,1.0),mass_matrix=mm_A) ######################################### Test each method for exactness @@ -63,6 +63,19 @@ sol2 = solve(prob2,ImplicitEuler()) @test norm(sol .- sol2) ≈ 0 atol=1e-9 +sol = solve(prob, Trapezoid()) +sol2 = solve(prob2,Trapezoid()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-9 + + +#= + +sol = solve(prob, TRBDF2()) +sol2 = solve(prob2,TRBDF2()) + +@test norm(sol .- sol2) ≈ 0 atol=1e-9 + #sol = solve(prob, SDIRK2()) #sol2 = solve(prob2, SDIRK2()) @@ -70,5 +83,7 @@ sol2 = solve(prob2,ImplicitEuler()) -#sol = solve(prob, Rodas3(),adaptive=false,dt=1/10) -#sol2 = solve(prob2, Rodas3(),adaptive=false,dt=1/10) +sol = solve(prob, TRBDF2(),adaptive=false,dt=1/10) +sol2 = solve(prob2, TRBDF2(),adaptive=false,dt=1/10) + +=# diff --git a/test/ode/ode_cache_tests.jl b/test/ode/ode_cache_tests.jl index e98c14b8d1..cf8d23556d 100644 --- a/test/ode/ode_cache_tests.jl +++ b/test/ode/ode_cache_tests.jl @@ -1,6 +1,6 @@ using OrdinaryDiffEq, DiffEqBase -NON_IMPLICIT_ALGS = filter((x)->isleaftype(x) && !OrdinaryDiffEq.isimplicit(x()),union(subtypes(OrdinaryDiffEqAlgorithm),subtypes(OrdinaryDiffEqAdaptiveAlgorithm))) +NON_IMPLICIT_ALGS = filter((x)->isleaftype(x) && !OrdinaryDiffEq.isimplicit(x()),union(subtypes(OrdinaryDiffEq.OrdinaryDiffEqAlgorithm),subtypes(OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm))) f = function (t,u,du) for i in 1:length(u) diff --git a/test/runtests.jl b/test/runtests.jl index 3910bf4d33..b4d2fd7a88 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -20,6 +20,7 @@ tic() @time @testset "Initial Dt Tests" begin include("ode/ode_initdt_tests.jl") end @time @testset "OwrenZen Tests" begin include("owrenzen_tests.jl") end @time @testset "Rosenbrock Tests" begin include("ode/ode_rosenbrock_tests.jl") end +@time @testset "Mass Matrix Tests" begin include("mass_matrix_tests.jl") end @time @testset "Differentiation Trait Tests" begin include("differentiation_traits_tests.jl") end @time @testset "Partitioned Methods Tests" begin include("partitioned_methods_tests.jl") end @time @testset "Split Methods Tests" begin include("split_methods_tests.jl") end From 5fea73f99be3f57f1f903a11ce87f445fb33c799 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Sun, 13 Aug 2017 03:26:17 -0700 Subject: [PATCH 6/6] remove Rodas4 FSAL --- src/alg_utils.jl | 6 +++--- src/integrators/rosenbrock_integrators.jl | 16 ++++------------ 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 678151c4aa..cfe4368034 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -18,9 +18,9 @@ isfsal(alg::Velds4) = true isfsal(alg::GRK4T) = true isfsal(alg::GRK4A) = true isfsal(alg::Ros4LStab) = true -isfsal(alg::Rodas4) = true -isfsal(alg::Rodas42) = true -isfsal(alg::Rodas4P) = true +isfsal(alg::Rodas4) = false +isfsal(alg::Rodas42) = false +isfsal(alg::Rodas4P) = false isfsal(alg::Rodas5) = true isfsal(alg::LawsonEuler) = true isfsal(alg::NorsettEuler) = true diff --git a/src/integrators/rosenbrock_integrators.jl b/src/integrators/rosenbrock_integrators.jl index c2b9c3ff6d..bd1088be3f 100644 --- a/src/integrators/rosenbrock_integrators.jl +++ b/src/integrators/rosenbrock_integrators.jl @@ -1010,12 +1010,9 @@ function initialize!(integrator, cache::Rodas4ConstantCache) integrator.kshortsize = 2 k = eltype(integrator.sol.k)(2) integrator.k = k - integrator.fsalfirst = integrator.f(integrator.t, integrator.uprev) - # Avoid undefined entries if k is an array of arrays - integrator.fsallast = zero(integrator.fsalfirst) - integrator.k[1] = zero(integrator.fsalfirst) - integrator.k[2] = zero(integrator.fsalfirst) + integrator.k[1] = zero(integrator.u) + integrator.k[2] = zero(integrator.u) end @muladd function perform_step!(integrator, cache::Rodas4ConstantCache, repeat_step=false) @@ -1108,24 +1105,19 @@ end integrator.k[1] = @. h21*k1 + h22*k2 + h23*k3 + h24*k4 + h25*k5 integrator.k[2] = @. h31*k1 + h32*k2 + h33*k3 + h34*k4 + h35*k5 end - - integrator.fsallast = du integrator.u = u end function initialize!(integrator, cache::Rodas4Cache) integrator.kshortsize = 2 - @unpack fsalfirst,fsallast,dense1,dense2 = cache - integrator.fsalfirst = fsalfirst - integrator.fsallast = fsallast + @unpack dense1,dense2 = cache integrator.k = [dense1,dense2] - integrator.f(integrator.t, integrator.uprev, integrator.fsalfirst) end @muladd function perform_step!(integrator, cache::Rodas4Cache, repeat_step=false) @unpack t,dt,uprev,u,f = integrator - @unpack du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,vectmp5,vectmp6,fsalfirst,fsallast,dT,J,W,uf,tf,linsolve_tmp,linsolve_tmp_vec,jac_config = cache + @unpack du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,vectmp5,vectmp6,dT,J,W,uf,tf,linsolve_tmp,linsolve_tmp_vec,jac_config = cache @unpack 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,gamma,c2,c3,c4,d1,d2,d3,d4 = cache.tab # Assignments