@@ -41,7 +41,7 @@ mutable struct RosenbrockCache{uType, rateType, uNoUnitsType, JType, WType, TabT
4141 alg:: A
4242 step_limiter!:: StepLimiter
4343 stage_limiter!:: StageLimiter
44- order :: Int
44+ interp_order :: Int
4545end
4646function full_cache (c:: RosenbrockCache )
4747 return [c. u, c. uprev, c. dense... , c. du, c. du1, c. du2,
@@ -56,7 +56,7 @@ struct RosenbrockCombinedConstantCache{TF, UF, Tab, JType, WType, F, AD} <: Rose
5656 W:: WType
5757 linsolve:: F
5858 autodiff:: AD
59- order :: Int
59+ interp_order :: Int
6060end
6161
6262@cache mutable struct Rosenbrock23Cache{uType, rateType, uNoUnitsType, JType, WType,
@@ -718,8 +718,12 @@ tabtype(::Rodas4) = Rodas4Tableau
718718tabtype (:: Rodas42 ) = Rodas42Tableau
719719tabtype (:: Rodas4P ) = Rodas4PTableau
720720tabtype (:: Rodas4P2 ) = Rodas4P2Tableau
721+ tabtype (:: Rodas5 ) = Rodas5Tableau
722+ tabtype (:: Rodas5P ) = Rodas5PTableau
723+ tabtype (:: Rodas5Pr ) = Rodas5PTableau
724+ tabtype (:: Rodas5Pe ) = Rodas5PTableau
721725
722- function alg_cache (alg:: Union{Rodas4, Rodas42, Rodas4P, Rodas4P2} ,
726+ function alg_cache (alg:: Union{Rodas4, Rodas42, Rodas4P, Rodas4P2, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr } ,
723727 u, rate_prototype, :: Type{uEltypeNoUnits} ,
724728 :: Type{uBottomEltypeNoUnits} , :: Type{tTypeNoUnits} , uprev, uprev2, f, t,
725729 dt, reltol, p, calck,
@@ -729,21 +733,22 @@ function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2},
729733 J, W = build_J_W (alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val (false ))
730734 linprob = nothing # LinearProblem(W,copy(u); u0=copy(u))
731735 linsolve = nothing # init(linprob,alg.linsolve,alias_A=true,alias_b=true)
736+ tab = tabtype (alg)(constvalue (uBottomEltypeNoUnits), constvalue (tTypeNoUnits))
732737 RosenbrockCombinedConstantCache (tf, uf,
733- tabtype (alg)(constvalue (uBottomEltypeNoUnits),
734- constvalue (tTypeNoUnits)), J, W, linsolve,
735- alg_autodiff (alg), 4 )
738+ tab, J, W, linsolve,
739+ alg_autodiff (alg), size (tab. H, 1 ))
736740end
737741
738- function alg_cache (alg:: Union{Rodas4, Rodas42, Rodas4P, Rodas4P2} ,
742+ function alg_cache (alg:: Union{Rodas4, Rodas42, Rodas4P, Rodas4P2, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr } ,
739743 u, rate_prototype, :: Type{uEltypeNoUnits} ,
740744 :: Type{uBottomEltypeNoUnits} , :: Type{tTypeNoUnits} , uprev, uprev2, f, t,
741745 dt, reltol, p, calck,
742746 :: Val{true} ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
743747
748+ tab = tabtype (alg)(constvalue (uBottomEltypeNoUnits), constvalue (tTypeNoUnits))
744749 # Initialize vectors
745- dense = [zero (rate_prototype) for _ in 1 : 2 ]
746- ks = [zero (rate_prototype) for _ in 1 : 6 ]
750+ dense = [zero (rate_prototype) for _ in 1 : size (tab . H, 1 ) ]
751+ ks = [zero (rate_prototype) for _ in 1 : size (tab . A, 1 ) ]
747752 du = zero (rate_prototype)
748753 du1 = zero (rate_prototype)
749754 du2 = zero (rate_prototype)
@@ -762,7 +767,6 @@ function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2},
762767 recursivefill! (atmp, false )
763768 weight = similar (u, uEltypeNoUnits)
764769 recursivefill! (weight, false )
765- tab = tabtype (alg)(constvalue (uBottomEltypeNoUnits), constvalue (tTypeNoUnits))
766770
767771 tf = TimeGradientWrapper (f, uprev, p)
768772 uf = UJacobianWrapper (f, t, p)
@@ -785,120 +789,9 @@ function alg_cache(alg::Union{Rodas4, Rodas42, Rodas4P, Rodas4P2},
785789 u, uprev, dense, du, du1, du2, ks, fsalfirst, fsallast,
786790 dT, J, W, tmp, atmp, weight, tab, tf, uf, linsolve_tmp,
787791 linsolve, jac_config, grad_config, reltol, alg,
788- alg. step_limiter!, alg. stage_limiter!, 4 )
789- end
790-
791- # ###############################################################################
792-
793- # ## Rosenbrock5
794-
795- function alg_cache (alg:: Rodas5 , u, rate_prototype, :: Type{uEltypeNoUnits} ,
796- :: Type{uBottomEltypeNoUnits} , :: Type{tTypeNoUnits} , uprev, uprev2, f, t,
797- dt, reltol, p, calck,
798- :: Val{true} ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
799- dense = [zero (rate_prototype) for _ in 1 : 3 ]
800- du = zero (rate_prototype)
801- du1 = zero (rate_prototype)
802- du2 = zero (rate_prototype)
803- ks = [zero (rate_prototype) for _ in 1 : 7 ]
804- fsalfirst = zero (rate_prototype)
805- fsallast = zero (rate_prototype)
806- dT = zero (rate_prototype)
807- J, W = build_J_W (alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val (true ))
808- tmp = zero (rate_prototype)
809- atmp = similar (u, uEltypeNoUnits)
810- recursivefill! (atmp, false )
811- weight = similar (u, uEltypeNoUnits)
812- recursivefill! (weight, false )
813- tab = Rodas5Tableau (constvalue (uBottomEltypeNoUnits), constvalue (tTypeNoUnits))
814-
815- tf = TimeGradientWrapper (f, uprev, p)
816- uf = UJacobianWrapper (f, t, p)
817- linsolve_tmp = zero (rate_prototype)
818- linprob = LinearProblem (W, _vec (linsolve_tmp); u0 = _vec (tmp))
819- Pl, Pr = wrapprecs (
820- alg. precs (W, nothing , u, p, t, nothing , nothing , nothing ,
821- nothing )... , weight, tmp)
822- linsolve = init (linprob, alg. linsolve, alias_A = true , alias_b = true ,
823- Pl = Pl, Pr = Pr,
824- assumptions = LinearSolve. OperatorAssumptions (true ))
825- grad_config = build_grad_config (alg, f, tf, du1, t)
826- jac_config = build_jac_config (alg, f, uf, du1, uprev, u, tmp, du2)
827- RosenbrockCache (u, uprev, dense, du, du1, du2, ks,
828- fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf,
829- linsolve_tmp,
830- linsolve, jac_config, grad_config, reltol, alg, alg. step_limiter!,
831- alg. stage_limiter!, 5 )
832- end
833-
834- function alg_cache (alg:: Rodas5 , u, rate_prototype, :: Type{uEltypeNoUnits} ,
835- :: Type{uBottomEltypeNoUnits} , :: Type{tTypeNoUnits} , uprev, uprev2, f, t,
836- dt, reltol, p, calck,
837- :: Val{false} ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
838- tf = TimeDerivativeWrapper (f, u, p)
839- uf = UDerivativeWrapper (f, t, p)
840- J, W = build_J_W (alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val (false ))
841- linprob = nothing # LinearProblem(W,copy(u); u0=copy(u))
842- linsolve = nothing # init(linprob,alg.linsolve,alias_A=true,alias_b=true)
843- RosenbrockCombinedConstantCache (tf, uf,
844- Rodas5Tableau (constvalue (uBottomEltypeNoUnits),
845- constvalue (tTypeNoUnits)), J, W, linsolve, alg_autodiff (alg), 5 )
792+ alg. step_limiter!, alg. stage_limiter!, size (tab. H, 1 ))
846793end
847794
848- function alg_cache (
849- alg:: Union{Rodas5P, Rodas5Pe, Rodas5Pr} , u, rate_prototype, :: Type{uEltypeNoUnits} ,
850- :: Type{uBottomEltypeNoUnits} , :: Type{tTypeNoUnits} , uprev, uprev2, f, t,
851- dt, reltol, p, calck,
852- :: Val{true} ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
853- dense = [zero (rate_prototype) for _ in 1 : 3 ]
854- du = zero (rate_prototype)
855- du1 = zero (rate_prototype)
856- du2 = zero (rate_prototype)
857- ks = [zero (rate_prototype) for _ in 1 : 8 ]
858- fsalfirst = zero (rate_prototype)
859- fsallast = zero (rate_prototype)
860- dT = zero (rate_prototype)
861- J, W = build_J_W (alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val (true ))
862- tmp = zero (rate_prototype)
863- atmp = similar (u, uEltypeNoUnits)
864- recursivefill! (atmp, false )
865- weight = similar (u, uEltypeNoUnits)
866- recursivefill! (weight, false )
867- tab = Rodas5PTableau (constvalue (uBottomEltypeNoUnits), constvalue (tTypeNoUnits))
868-
869- tf = TimeGradientWrapper (f, uprev, p)
870- uf = UJacobianWrapper (f, t, p)
871- linsolve_tmp = zero (rate_prototype)
872- linprob = LinearProblem (W, _vec (linsolve_tmp); u0 = _vec (tmp))
873- Pl, Pr = wrapprecs (
874- alg. precs (W, nothing , u, p, t, nothing , nothing , nothing ,
875- nothing )... , weight, tmp)
876- linsolve = init (linprob, alg. linsolve, alias_A = true , alias_b = true ,
877- Pl = Pl, Pr = Pr,
878- assumptions = LinearSolve. OperatorAssumptions (true ))
879- grad_config = build_grad_config (alg, f, tf, du1, t)
880- jac_config = build_jac_config (alg, f, uf, du1, uprev, u, tmp, du2)
881- RosenbrockCache (u, uprev, dense, du, du1, du2, ks,
882- fsalfirst, fsallast, dT, J, W, tmp, atmp, weight, tab, tf, uf,
883- linsolve_tmp,
884- linsolve, jac_config, grad_config, reltol, alg, alg. step_limiter!,
885- alg. stage_limiter!, 5 )
886- end
887-
888- function alg_cache (
889- alg:: Union{Rodas5P, Rodas5Pe, Rodas5Pr} , u, rate_prototype, :: Type{uEltypeNoUnits} ,
890- :: Type{uBottomEltypeNoUnits} , :: Type{tTypeNoUnits} , uprev, uprev2, f, t,
891- dt, reltol, p, calck,
892- :: Val{false} ) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
893- tf = TimeDerivativeWrapper (f, u, p)
894- uf = UDerivativeWrapper (f, t, p)
895- J, W = build_J_W (alg, u, uprev, p, t, dt, f, uEltypeNoUnits, Val (false ))
896- linprob = nothing # LinearProblem(W,copy(u); u0=copy(u))
897- linsolve = nothing # init(linprob,alg.linsolve,alias_A=true,alias_b=true)
898- RosenbrockCombinedConstantCache (tf, uf,
899- Rodas5PTableau (constvalue (uBottomEltypeNoUnits),
900- constvalue (tTypeNoUnits)), J, W, linsolve, alg_autodiff (alg), 5 )
901- end
902795
903796function get_fsalfirstlast (
904797 cache:: Union {Rosenbrock23Cache, Rosenbrock32Cache, Rosenbrock33Cache,
0 commit comments