rkc shift

SciML · Jul 15, 2019 · bbcdf53 · bbcdf53
1 parent 25cdfa9
commit bbcdf53
Show file tree

Hide file tree

Showing 3 changed files with 14 additions and 24 deletions.
diff --git a/src/caches/rkc_caches.jl b/src/caches/rkc_caches.jl
@@ -111,35 +111,24 @@ function alg_cache(alg::RKC,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits
   RKCConstantCache(u)
 end
 
-@cache mutable struct IRKCConstantCache{uType,rateType,F,N} <: OrdinaryDiffEqConstantCache
+@cache mutable struct IRKCConstantCache{uType,rateType,N} <: OrdinaryDiffEqConstantCache
   minm::Int64
   zprev::uType
-  uf::F
   nlsolver::N
   du₁::rateType
   du₂::rateType
 end
 
-@cache mutable struct IRKCCache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,N,F} <: OrdinaryDiffEqMutableCache
+@cache mutable struct IRKCCache{uType,rateType,uNoUnitsType,N} <: OrdinaryDiffEqMutableCache
   u::uType
   uprev::uType
   gprev::uType
   gprev2::uType
   fsalfirst::rateType
-  k::rateType
-  du1::rateType
   f1ⱼ₋₁::rateType
   f1ⱼ₋₂::rateType
   f2ⱼ₋₁::rateType
-  z::uType
-  dz::uType
-  tmp::uType
   atmp::uNoUnitsType
-  J::JType
-  W::WType
-  uf::UF
-  jac_config::JC
-  linsolve::F
   nlsolver::N
   du₁::rateType
   du₂::rateType
@@ -150,32 +139,28 @@ function alg_cache(alg::IRKC,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnit
   γ, c = 1.0, 1.0
   J, W = oop_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
   nlsolver = oopnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
-  @getoopnlsolvefields
   zprev = u
   du₁ = rate_prototype; du₂ = rate_prototype
-  IRKCConstantCache(50,zprev,uf,nlsolver,du₁,du₂)
+  IRKCConstantCache(50,zprev,nlsolver,du₁,du₂)
 end
 
 function alg_cache(alg::IRKC,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
   γ, c = 1.0, 1.0
   J, W = iip_generate_W(alg,u,uprev,p,t,dt,f,uEltypeNoUnits)
   nlsolver = iipnlsolve(alg,u,uprev,p,t,dt,f,W,J,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,γ,c)
-  @getiipnlsolvefields
 
   gprev = similar(u)
   gprev2 = similar(u)
-  tmp = similar(u)
   atmp = similar(u,uEltypeNoUnits)
   fsalfirst = zero(rate_prototype)
-  k  = zero(rate_prototype)
   zprev = similar(u)
   f1ⱼ₋₁ = zero(rate_prototype)
   f1ⱼ₋₂ = zero(rate_prototype)
   f2ⱼ₋₁ = zero(rate_prototype)
   du₁ = zero(rate_prototype)
   du₂ = zero(rate_prototype)
-  constantcache = IRKCConstantCache(50,zprev,uf,nlsolver,du₁,du₂)
-  IRKCCache(u,uprev,gprev,gprev2,fsalfirst,k,du1,f1ⱼ₋₁,f1ⱼ₋₂,f2ⱼ₋₁,z,dz,tmp,atmp,J,W,uf,jac_config,linsolve,nlsolver,du₁,du₂,constantcache)
+  constantcache = IRKCConstantCache(50,zprev,nlsolver,du₁,du₂)
+  IRKCCache(u,uprev,gprev,gprev2,fsalfirst,f1ⱼ₋₁,f1ⱼ₋₂,f2ⱼ₋₁,atmp,nlsolver,du₁,du₂,constantcache)
 end
 
 mutable struct ESERK4ConstantCache{T, zType} <: OrdinaryDiffEqConstantCache

diff --git a/src/perform_step/rkc_perform_step.jl b/src/perform_step/rkc_perform_step.jl
@@ -653,7 +653,7 @@ function initialize!(integrator, cache::IRKCCache)
   @unpack f1, f2 = integrator.f
   integrator.kshortsize = 2
   integrator.fsalfirst = cache.fsalfirst
-  integrator.fsallast = cache.k
+  integrator.fsallast = cache.nlsolver.k
   resize!(integrator.k, integrator.kshortsize)
   integrator.k[1] = integrator.fsalfirst
   integrator.k[2] = integrator.fsallast
@@ -666,7 +666,8 @@ end
 
 function perform_step!(integrator, cache::IRKCCache, repeat_step=false)
   @unpack t,dt,uprev,u,f,p,alg = integrator
-  @unpack tmp,gprev,gprev2,k,f1ⱼ₋₁,f1ⱼ₋₂,f2ⱼ₋₁,du₁,du₂,z,W,atmp,nlsolver = cache
+  @unpack gprev,gprev2,f1ⱼ₋₁,f1ⱼ₋₂,f2ⱼ₋₁,du₁,du₂,atmp,nlsolver = cache
+  @unpack tmp,k,z = nlsolver
   @unpack minm = cache.constantcache
   @unpack f1, f2 = integrator.f
 
@@ -767,7 +768,7 @@ function perform_step!(integrator, cache::IRKCCache, repeat_step=false)
   if isnewton(nlsolver) && integrator.opts.adaptive
     update_W!(integrator, cache, dt, false)
     @.. gprev = dt*0.5*(du₂ - f2ⱼ₋₁) + dt*(0.5 - μs₁)*(du₁ - f1ⱼ₋₁)
-    cache.linsolve(vec(tmp),get_W(nlsolver),vec(gprev),false)
+    nlsolver.linsolve(vec(tmp),get_W(nlsolver),vec(gprev),false)
     calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
     integrator.EEst = integrator.opts.internalnorm(atmp,t)
   end

diff --git a/src/rkc_utils.jl b/src/rkc_utils.jl
@@ -99,7 +99,11 @@ end
 function maxeig!(integrator, cache::OrdinaryDiffEqMutableCache)
   isfirst = integrator.iter == 1 || integrator.u_modified
   @unpack t, dt, uprev, u, f, p, fsalfirst = integrator
-  fz, z, atmp = cache.k, cache.tmp, cache.atmp
+  if cache isa IRKCCache
+    fz, z, atmp = cache.nlsolver.k, cache.nlsolver.tmp, cache.atmp
+  else
+    fz, z, atmp = cache.k, cache.tmp, cache.atmp
+  end
   ccache = cache.constantcache
   maxiter = (typeof(integrator.alg) <: Union{ESERK4,ESERK5,SERK2}) ? 100 : 50
   safe = (typeof(integrator.alg) <: RKCAlgs) ? 1.0 : 1.2