Skip to content

Commit

Permalink
Rodas4 tests pass
Browse files Browse the repository at this point in the history
  • Loading branch information
@ChrisRackauckas
ChrisRackauckas committed Jul 6, 2017
1 parent 61d2c99 commit 0698191
Show file tree
Hide file tree
Showing 4 changed files with 172 additions and 7 deletions.
8 changes: 6 additions & 2 deletions src/alg_utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,9 @@ get_chunksize{CS,AD}(alg::Velds4{CS,AD}) = CS
get_chunksize{CS,AD}(alg::GRK4T{CS,AD}) = CS
get_chunksize{CS,AD}(alg::GRK4A{CS,AD}) = CS
get_chunksize{CS,AD}(alg::Ros4LStab{CS,AD}) = CS

get_chunksize{CS,AD}(alg::Rodas4{CS,AD}) = CS
get_chunksize{CS,AD}(alg::Rodas42{CS,AD}) = CS
get_chunksize{CS,AD}(alg::Rodas4P{CS,AD}) = CS

alg_extrapolates(alg::OrdinaryDiffEqAlgorithm) = false
alg_extrapolates(alg::ImplicitEuler) = true
Expand All @@ -115,7 +117,9 @@ alg_autodiff{CS,AD}(alg::Velds4{CS,AD}) = AD
alg_autodiff{CS,AD}(alg::GRK4T{CS,AD}) = AD
alg_autodiff{CS,AD}(alg::GRK4A{CS,AD}) = AD
alg_autodiff{CS,AD}(alg::Ros4LStab{CS,AD}) = AD

alg_autodiff{CS,AD}(alg::Rodas4{CS,AD}) = AD
alg_autodiff{CS,AD}(alg::Rodas42{CS,AD}) = AD
alg_autodiff{CS,AD}(alg::Rodas4P{CS,AD}) = AD

alg_order(alg::OrdinaryDiffEqAlgorithm) = error("Order is not defined for this algorithm")
alg_adaptive_order(alg::OrdinaryDiffEqAdaptiveAlgorithm) = error("Algorithm is adaptive with no order")
Expand Down
80 changes: 80 additions & 0 deletions src/caches/rosenbrock_caches.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -594,6 +594,10 @@ end
type Rodas4Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,TabType,TFType,UFType,F,JCType} <: OrdinaryDiffEqMutableCache
u::uType
uprev::uType
k1::rateType
k2::rateType
k3::rateType
k4::rateType
du::rateType
du1::rateType
du2::du2Type
Expand Down Expand Up @@ -667,12 +671,88 @@ function alg_cache(alg::Rodas4,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,upre
Rodas4ConstantCache(tf,uf,Rodas4ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits)))
end

function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
k1 = zeros(rate_prototype)
k2 = zeros(rate_prototype)
k3 = zeros(rate_prototype)
k4 = zeros(rate_prototype)
du = zeros(rate_prototype)
du1 = zeros(rate_prototype)
du2 = zeros(rate_prototype)
vectmp = vec(similar(u,indices(u)))
vectmp2 = vec(similar(u,indices(u)))
vectmp3 = vec(similar(u,indices(u)))
vectmp4 = vec(similar(u,indices(u)))
vectmp5 = vec(similar(u,indices(u)))
vectmp6 = vec(similar(u,indices(u)))
fsalfirst = zeros(rate_prototype)
fsallast = zeros(rate_prototype)
dT = similar(u,indices(u))
J = zeros(uEltypeNoUnits,length(u),length(u)) # uEltype?
W = similar(J);
tmp = similar(u,indices(u))
tab = Rodas42ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
vf = VectorF(f,size(u))
vfr = VectorFReturn(f,size(u))
tf = TimeGradientWrapper(vf,uprev)
uf = UJacobianWrapper(vfr,t)
linsolve_tmp = similar(u,indices(u))
linsolve_tmp_vec = vec(linsolve_tmp)
if alg_autodiff(alg)
jac_config = ForwardDiff.JacobianConfig(uf,vec(du1),vec(uprev),ForwardDiff.Chunk{determine_chunksize(u,alg)}())
else
jac_config = nothing
end
Rodas4Cache(u,uprev,k1,k2,k3,k4,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
vectmp5,vectmp6,
fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
linsolve_tmp_vec,alg.linsolve,jac_config)
end

function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}})
tf = TimeDerivativeWrapper(f,u)
uf = UDerivativeWrapper(f,t)
Rodas4ConstantCache(tf,uf,Rodas42ConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits)))
end

function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{true}})
k1 = zeros(rate_prototype)
k2 = zeros(rate_prototype)
k3 = zeros(rate_prototype)
k4 = zeros(rate_prototype)
du = zeros(rate_prototype)
du1 = zeros(rate_prototype)
du2 = zeros(rate_prototype)
vectmp = vec(similar(u,indices(u)))
vectmp2 = vec(similar(u,indices(u)))
vectmp3 = vec(similar(u,indices(u)))
vectmp4 = vec(similar(u,indices(u)))
vectmp5 = vec(similar(u,indices(u)))
vectmp6 = vec(similar(u,indices(u)))
fsalfirst = zeros(rate_prototype)
fsallast = zeros(rate_prototype)
dT = similar(u,indices(u))
J = zeros(uEltypeNoUnits,length(u),length(u)) # uEltype?
W = similar(J);
tmp = similar(u,indices(u))
tab = Rodas4PConstantCache(realtype(uEltypeNoUnits),realtype(tTypeNoUnits))
vf = VectorF(f,size(u))
vfr = VectorFReturn(f,size(u))
tf = TimeGradientWrapper(vf,uprev)
uf = UJacobianWrapper(vfr,t)
linsolve_tmp = similar(u,indices(u))
linsolve_tmp_vec = vec(linsolve_tmp)
if alg_autodiff(alg)
jac_config = ForwardDiff.JacobianConfig(uf,vec(du1),vec(uprev),ForwardDiff.Chunk{determine_chunksize(u,alg)}())
else
jac_config = nothing
end
Rodas4Cache(u,uprev,k1,k2,k3,k4,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
vectmp5,vectmp6,
fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
linsolve_tmp_vec,alg.linsolve,jac_config)
end

function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,::Type{Val{false}})
tf = TimeDerivativeWrapper(f,u)
uf = UDerivativeWrapper(f,t)
Expand Down
48 changes: 43 additions & 5 deletions src/integrators/rosenbrock_integrators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -878,7 +878,7 @@ end

################################################################################

#### ROS4 type method
#### Rodas4 type method

@inline function initialize!(integrator,cache::Rodas4ConstantCache,f=integrator.f)
integrator.kshortsize = 2
Expand Down Expand Up @@ -977,7 +977,6 @@ end
@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
mass_matrix = integrator.sol.prob.mass_matrix

utilde = du
atmp = du

# Setup Jacobian Calc
Expand Down Expand Up @@ -1061,6 +1060,12 @@ end

k3 = reshape(vectmp3,sizeu...)

@tight_loop_macros for i in uidx
@inbounds u[i] = uprev[i] + a41*k1[i] + a42*k2[i] + a43*k3[i]
end

f(t+c4*dt,u,du)

@tight_loop_macros for i in uidx
@inbounds linsolve_tmp[i] = du[i] + dt*d4*dT[i] + C41*k1[i]/dt + C42*k2[i]/dt + C43*k3[i]/dt
end
Expand All @@ -1074,15 +1079,48 @@ end
k4 = reshape(vectmp4,sizeu...)

@tight_loop_macros for i in uidx
@inbounds u[i] = uprev[i] + b1*k1[i] + b2*k2[i] + b3*k3[i] + b4*k4[i]
@inbounds u[i] = uprev[i] + a51*k1[i] + a52*k2[i] + a53*k3[i] + a54*k4[i]
end

f(t+dt,u,du)

@tight_loop_macros for i in uidx
@inbounds linsolve_tmp[i] = du[i] + C52*k2[i]/dt + C54*k4[i]/dt + C51*k1[i]/dt + C53*k3[i]/dt
end

if has_invW(f)
A_mul_B!(vectmp5,W,linsolve_tmp_vec)
else
integrator.alg.linsolve(vectmp5,W,linsolve_tmp_vec)
end

k5 = reshape(vectmp5,sizeu...)

@tight_loop_macros for i in uidx
@inbounds u[i] = u[i] + k5[i]
end

f(t,u,fsallast)

@tight_loop_macros for i in uidx
@inbounds linsolve_tmp[i] = fsallast[i] + C61*k1[i]/dt + C62*k2[i]/dt + C63*k3[i]/dt + C64*k4[i]/dt + C65*k5[i]/dt
end

if has_invW(f)
A_mul_B!(vectmp6,W,linsolve_tmp_vec)
else
integrator.alg.linsolve(vectmp6,W,linsolve_tmp_vec)
end

k6 = reshape(vectmp6,sizeu...)

@tight_loop_macros for i in uidx
@inbounds u[i] = u[i] + k6[i]
end

if integrator.opts.adaptive
@tight_loop_macros for (i,atol,rtol) in zip(uidx,Iterators.cycle(integrator.opts.abstol),Iterators.cycle(integrator.opts.reltol))
@inbounds utilde[i] = @muladd btilde1*k1[i] + btilde2*k2[i] + btilde3*k3[i] + btilde4*k4[i]
@inbounds atmp[i] = ((utilde[i])./@muladd(atol+max(abs(uprev[i]),abs(u[i])).*rtol))
@inbounds atmp[i] = ((k6[i])./@muladd(atol+max(abs(uprev[i]),abs(u[i])).*rtol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down
43 changes: 43 additions & 0 deletions test/ode/ode_rosenbrock_tests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -284,6 +284,49 @@ sim = test_convergence(dts,prob,Ros4LStab(linsolve=LinSolveFactorize(qrfact!)))
sol = solve(prob,Ros4LStab(linsolve=LinSolveFactorize(qrfact!)))
@test length(sol) < 20

### Rodas4 Algorithms

prob = prob_ode_linear

sim = test_convergence(dts,prob,Rodas4())
@test abs(sim.𝒪est[:final]-4) < testTol

sol = solve(prob,Rodas4())
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas42())
@test abs(sim.𝒪est[:final]-4) < testTol

sol = solve(prob,Rodas42())
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas4P())
@test abs(sim.𝒪est[:final]-4) < testTol

sol = solve(prob,Rodas4P())
@test length(sol) < 20

prob = prob_ode_2Dlinear

sim = test_convergence(dts,prob,Rodas4())
@test abs(sim.𝒪est[:final]-4) < testTol

sol = solve(prob,Rodas4())
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas42())
@test abs(sim.𝒪est[:final]-4) < testTol

sol = solve(prob,Rodas42())
@test length(sol) < 20

sim = test_convergence(dts,prob,Rodas4P())
@test abs(sim.𝒪est[:final]-4) < testTol

sol = solve(prob,Rodas4P())
@test length(sol) < 20


### Test on Stiff

prob = deepcopy(prob_ode_rober)
Expand Down

0 comments on commit 0698191

Please sign in to comment.