Skip to content

add J, W cache and linsolve interface for rosenbrock oop methods and add J cache for firk oop #873

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 8 commits into from
Aug 16, 2019
Merged

add J, W cache and linsolve interface for rosenbrock oop methods and add J cache for firk oop #873

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
c4b6c93
add J and W cache for rosenbrock and firk oop
huanglangwen Aug 14, 2019
c2d9f45
refactor init J and W
huanglangwen Aug 14, 2019
1330fa7
fix test error
huanglangwen Aug 14, 2019
9a46866
fix test error
huanglangwen Aug 14, 2019
1653cb9
Update src/integrators/integrator_utils.jl
huanglangwen Aug 15, 2019
a59f756
add linsolve interface to oop rosenbrock methods
huanglangwen Aug 15, 2019
108d091
split tab and oop cache for ros23/32
huanglangwen Aug 15, 2019
dee2e87
#TODO: Rename other tableau structs
huanglangwen Aug 15, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 8 additions & 2 deletions src/caches/firk_caches.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
mutable struct RadauIIA5ConstantCache{F,Tab,Tol,Dt,U} <: OrdinaryDiffEqConstantCache
mutable struct RadauIIA5ConstantCache{F,Tab,Tol,Dt,U,JType} <: OrdinaryDiffEqConstantCache
uf::F
tab::Tab
κ::Tol
Expand All @@ -10,6 +10,7 @@ mutable struct RadauIIA5ConstantCache{F,Tab,Tol,Dt,U} <: OrdinaryDiffEqConstantC
dtprev::Dt
W_dt::Dt
status::DiffEqBase.NLStatus
J::JType
end

function alg_cache(alg::RadauIIA5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
Expand All @@ -19,8 +20,13 @@ function alg_cache(alg::RadauIIA5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeN
tab = RadauIIA5Tableau(uToltype, real(tTypeNoUnits))

κ = alg.κ !== nothing ? convert(uToltype, alg.κ) : convert(uToltype, 1//100)
if rate_prototype isa Number
J = 0
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)'
end

RadauIIA5ConstantCache(uf, tab, κ, one(uToltype), 10000, u, u, u, dt, dt, Convergence)
RadauIIA5ConstantCache(uf, tab, κ, one(uToltype), 10000, u, u, u, dt, dt, Convergence, J)
end

mutable struct RadauIIA5Cache{uType,cuType,uNoUnitsType,rateType,JType,W1Type,W2Type,UF,JC,F1,F2,Tab,Tol,Dt,rTol,aTol} <: OrdinaryDiffEqMutableCache
Expand Down
148 changes: 48 additions & 100 deletions src/caches/rosenbrock_caches.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -64,19 +64,10 @@ function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,uBottomElty
fsalfirst = zero(rate_prototype)
fsallast = zero(rate_prototype)
dT = zero(rate_prototype)
if ArrayInterface.isstructured(f.jac_prototype) || f.jac_prototype isa SparseMatrixCSC
J = similar(f.jac_prototype)
W = similar(J)
elseif DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
W = WOperator(f, dt, true)
J = nothing # is J = W.J better?
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
W = similar(J)
end
J,W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
tmp = zero(rate_prototype)
atmp = similar(u, uEltypeNoUnits)
tab = Rosenbrock23ConstantCache(real(uBottomEltypeNoUnits),identity,identity)
tab = Rosenbrock23ConstantCache(real(uBottomEltypeNoUnits),identity,identity,false,false)
tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
linsolve_tmp = zero(rate_prototype)
Expand All @@ -101,19 +92,10 @@ function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,uBottomElty
fsalfirst = zero(rate_prototype)
fsallast = zero(rate_prototype)
dT = zero(rate_prototype)
if ArrayInterface.isstructured(f.jac_prototype) || f.jac_prototype isa SparseMatrixCSC
J = similar(f.jac_prototype)
W = similar(J)
elseif DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
W = WOperator(f, dt, true)
J = nothing # is J = W.J better?
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
W = similar(J)
end
J,W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
tmp = zero(rate_prototype)
atmp = similar(u, uEltypeNoUnits)
tab = Rosenbrock32ConstantCache(real(uBottomEltypeNoUnits),identity,identity)
tab = Rosenbrock32ConstantCache(real(uBottomEltypeNoUnits),identity,identity,false,false)

tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
Expand All @@ -124,52 +106,60 @@ function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,uBottomElty
Rosenbrock32Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,fsalfirst,fsallast,dT,J,W,tmp,atmp,tab,tf,uf,linsolve_tmp,linsolve,jac_config,grad_config)
end

struct Rosenbrock23ConstantCache{T,TF,UF} <: OrdinaryDiffEqConstantCache
struct Rosenbrock23ConstantCache{T,TF,UF,JType,WType} <: OrdinaryDiffEqConstantCache
c₃₂::T
d::T
tf::TF
uf::UF
J::JType
W::WType
end

function Rosenbrock23ConstantCache(T::Type,tf,uf)
function Rosenbrock23ConstantCache(T::Type,tf,uf,J,W)
c₃₂ = convert(T,6 + sqrt(2))
d = convert(T,1/(2+sqrt(2)))
Rosenbrock23ConstantCache(c₃₂,d,tf,uf)
Rosenbrock23ConstantCache(c₃₂,d,tf,uf,J,W)
end

function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
Rosenbrock23ConstantCache(real(uBottomEltypeNoUnits),tf,uf)
J,W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
Rosenbrock23ConstantCache(real(uBottomEltypeNoUnits),tf,uf,J,W)
end

struct Rosenbrock32ConstantCache{T,TF,UF} <: OrdinaryDiffEqConstantCache
struct Rosenbrock32ConstantCache{T,TF,UF,JType,WType} <: OrdinaryDiffEqConstantCache
c₃₂::T
d::T
tf::TF
uf::UF
J::JType
W::WType
end

function Rosenbrock32ConstantCache(T::Type,tf,uf)
function Rosenbrock32ConstantCache(T::Type,tf,uf,J,W)
c₃₂ = convert(T,6 + sqrt(2))
d = convert(T,1/(2+sqrt(2)))
Rosenbrock32ConstantCache(c₃₂,d,tf,uf)
Rosenbrock32ConstantCache(c₃₂,d,tf,uf,J,W)
end

function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
Rosenbrock32ConstantCache(real(uBottomEltypeNoUnits),tf,uf)
J,W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
Rosenbrock32ConstantCache(real(uBottomEltypeNoUnits),tf,uf,J,W)
end

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

### 3rd order specialized Rosenbrocks

struct Rosenbrock33ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
struct Rosenbrock33ConstantCache{TF,UF,Tab,JType,WType} <: OrdinaryDiffEqConstantCache
tf::TF
uf::UF
tab::Tab
J::JType
W::WType
end

@cache mutable struct Rosenbrock33Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
Expand Down Expand Up @@ -209,16 +199,7 @@ function alg_cache(alg::ROS3P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUni
fsalfirst = zero(rate_prototype)
fsallast = zero(rate_prototype)
dT = zero(rate_prototype)
if ArrayInterface.isstructured(f.jac_prototype) || f.jac_prototype isa SparseMatrixCSC
J = similar(f.jac_prototype)
W = similar(J)
elseif DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
W = WOperator(f, dt, true)
J = nothing # is J = W.J better?
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
W = similar(J)
end
J,W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
tmp = zero(rate_prototype)
atmp = similar(u, uEltypeNoUnits)
tab = ROS3PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
Expand All @@ -236,7 +217,8 @@ end
function alg_cache(alg::ROS3P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
Rosenbrock33ConstantCache(tf,uf,ROS3PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
J,W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
Rosenbrock33ConstantCache(tf,uf,ROS3PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W)
end

@cache mutable struct Rosenbrock34Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
Expand Down Expand Up @@ -276,16 +258,7 @@ function alg_cache(alg::Rodas3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUn
fsalfirst = zero(rate_prototype)
fsallast = zero(rate_prototype)
dT = zero(rate_prototype)
if ArrayInterface.isstructured(f.jac_prototype) || f.jac_prototype isa SparseMatrixCSC
J = similar(f.jac_prototype)
W = similar(J)
elseif DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
W = WOperator(f, dt, true)
J = nothing # is J = W.J better?
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
W = similar(J)
end
J,W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
tmp = zero(rate_prototype)
atmp = similar(u, uEltypeNoUnits)
tab = Rodas3ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
Expand All @@ -301,16 +274,19 @@ function alg_cache(alg::Rodas3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUn
linsolve,jac_config,grad_config)
end

struct Rosenbrock34ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
struct Rosenbrock34ConstantCache{TF,UF,Tab,JType,WType} <: OrdinaryDiffEqConstantCache
tf::TF
uf::UF
tab::Tab
J::JType
W::WType
end

function alg_cache(alg::Rodas3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
Rosenbrock34ConstantCache(tf,uf,Rodas3ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
J,W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
Rosenbrock34ConstantCache(tf,uf,Rodas3ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W)
end

################################################################################
Expand All @@ -330,10 +306,12 @@ jac_cache(c::Rosenbrock4Cache) = (c.J,c.W)

### Rodas methods

struct Rodas4ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
struct Rodas4ConstantCache{TF,UF,Tab,JType,WType} <: OrdinaryDiffEqConstantCache
tf::TF
uf::UF
tab::Tab
J::JType
W::WType
end

@cache mutable struct Rodas4Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
Expand Down Expand Up @@ -381,16 +359,7 @@ function alg_cache(alg::Rodas4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUn
fsalfirst = zero(rate_prototype)
fsallast = zero(rate_prototype)
dT = zero(rate_prototype)
if ArrayInterface.isstructured(f.jac_prototype) || f.jac_prototype isa SparseMatrixCSC
J = similar(f.jac_prototype)
W = similar(J)
elseif DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
W = WOperator(f, dt, true)
J = nothing # is J = W.J better?
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
W = similar(J)
end
J,W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
tmp = zero(rate_prototype)
atmp = similar(u, uEltypeNoUnits)
tab = Rodas4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
Expand All @@ -410,7 +379,8 @@ end
function alg_cache(alg::Rodas4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
Rodas4ConstantCache(tf,uf,Rodas4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
J,W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
Rodas4ConstantCache(tf,uf,Rodas4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W)
end

function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
Expand All @@ -428,16 +398,7 @@ function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoU
fsalfirst = zero(rate_prototype)
fsallast = zero(rate_prototype)
dT = zero(rate_prototype)
if ArrayInterface.isstructured(f.jac_prototype) || f.jac_prototype isa SparseMatrixCSC
J = similar(f.jac_prototype)
W = similar(J)
elseif DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
W = WOperator(f, dt, true)
J = nothing # is J = W.J better?
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
W = similar(J)
end
J,W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
tmp = zero(rate_prototype)
atmp = similar(u, uEltypeNoUnits)
tab = Rodas42ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
Expand All @@ -457,7 +418,8 @@ end
function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
Rodas4ConstantCache(tf,uf,Rodas42ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
J,W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
Rodas4ConstantCache(tf,uf,Rodas42ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W)
end

function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
Expand All @@ -475,16 +437,7 @@ function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoU
fsalfirst = zero(rate_prototype)
fsallast = zero(rate_prototype)
dT = zero(rate_prototype)
if ArrayInterface.isstructured(f.jac_prototype) || f.jac_prototype isa SparseMatrixCSC
J = similar(f.jac_prototype)
W = similar(J)
elseif DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
W = WOperator(f, dt, true)
J = nothing # is J = W.J better?
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
W = similar(J)
end
J,W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
tmp = zero(rate_prototype)
atmp = similar(u, uEltypeNoUnits)
tab = Rodas4PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
Expand All @@ -504,17 +457,20 @@ end
function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
Rodas4ConstantCache(tf,uf,Rodas4PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
J,W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
Rodas4ConstantCache(tf,uf,Rodas4PConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W)
end

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

### Rosenbrock5

struct Rosenbrock5ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
struct Rosenbrock5ConstantCache{TF,UF,Tab,JType,WType} <: OrdinaryDiffEqConstantCache
tf::TF
uf::UF
tab::Tab
J::JType
W::WType
end

@cache mutable struct Rosenbrock5Cache{uType,rateType,uNoUnitsType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
Expand Down Expand Up @@ -566,16 +522,7 @@ function alg_cache(alg::Rodas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUn
fsalfirst = zero(rate_prototype)
fsallast = zero(rate_prototype)
dT = zero(rate_prototype)
if ArrayInterface.isstructured(f.jac_prototype) || f.jac_prototype isa SparseMatrixCSC
J = similar(f.jac_prototype)
W = similar(J)
elseif DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
W = WOperator(f, dt, true)
J = nothing # is J = W.J better?
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
W = similar(J)
end
J,W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
tmp = zero(rate_prototype)
atmp = similar(u, uEltypeNoUnits)
tab = Rodas5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))
Expand All @@ -595,7 +542,8 @@ end
function alg_cache(alg::Rodas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
Rosenbrock5ConstantCache(tf,uf,Rodas5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
J,W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
Rosenbrock5ConstantCache(tf,uf,Rodas5ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W)
end

################################################################################
Expand Down
15 changes: 6 additions & 9 deletions src/generic_rosenbrock.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,10 +162,12 @@ because in the end of Rosenbrock's method, we have: `y_{n+1}=y_n+ki*bi`.
"""
function gen_cache_struct(tab::RosenbrockTableau,cachename::Symbol,constcachename::Symbol)
kstype=[:($(Symbol(:k,i))::rateType) for i in 1:length(tab.b)]
constcacheexpr=quote struct $constcachename{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
constcacheexpr=quote struct $constcachename{TF,UF,Tab,JType,WType} <: OrdinaryDiffEqConstantCache
tf::TF
uf::UF
tab::Tab
J::JType
W::WType
end
end
cacheexpr=quote
Expand Down Expand Up @@ -207,7 +209,8 @@ function gen_algcache(cacheexpr::Expr,cachename::Symbol,constcachename::Symbol,a
function alg_cache(alg::$algname,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
$constcachename(tf,uf,$tabname(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)))
J,W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
$constcachename(tf,uf,$tabname(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)),J,W)
end
function alg_cache(alg::$algname,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
du = zero(rate_prototype)
Expand All @@ -217,13 +220,7 @@ function gen_algcache(cacheexpr::Expr,cachename::Symbol,constcachename::Symbol,a
fsalfirst = zero(rate_prototype)
fsallast = zero(rate_prototype)
dT = zero(rate_prototype)
if DiffEqBase.has_jac(f) && !DiffEqBase.has_Wfact(f) && f.jac_prototype !== nothing
W = WOperator(f, dt, true)
J = nothing # is J = W.J better?
else
J = false .* vec(rate_prototype) .* vec(rate_prototype)' # uEltype?
W = similar(J)
end
J,W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
tmp = zero(rate_prototype)
tab = $tabname(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits))

Expand Down
Loading