seperate material + stress iteration caches

src/CrystalViscoPlastic/CrystalViscoPlastic.jl

@@ -59,7 +59,7 @@ end
 Base.zero(::Type{CrystalViscoPlasticState{dim,T,M,S}}) where {dim,T,M,S} = CrystalViscoPlasticState{dim,T,M,S}(zero(SymmetricTensor{2,dim,T,M}), SVector{S,T}(zeros(T, S)), SVector{S,T}(zeros(T, S)), SVector{S,T}(zeros(T, S)))
 initial_material_state(::CrystalViscoPlastic{S}) where S = zero(CrystalViscoPlasticState{3,Float64,6,S})
 
-function get_cache(m::CrystalViscoPlastic)
+function _get_cache(m::CrystalViscoPlastic)
     state = initial_material_state(m)
     # it doesn't actually matter which state and strain step we use here,
     # f is overwritten in the constitutive driver before being used.
@@ -110,23 +110,25 @@ end
 function material_response(m::CrystalViscoPlastic{S}, Δε::SymmetricTensor{2,3,T,6}, state::CrystalViscoPlasticState{3},
     Δt=1.0; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {S, T}
 
+    _cache = solver_cache(cache)
+
     # set the current residual function that depends only on the variables
     f(r_vector, x_vector) = vector_residual!((x->residuals(x,m,state,Δε,Δt)), r_vector, x_vector, m)
-    update_cache!(cache, f)
+    update_cache!(_cache, f)
     # initial guess
     σ_trial = state.σ + m.Eᵉ ⊡ Δε
     x0 = ResidualsCrystalViscoPlastic(σ_trial, state.κ, state.α, state.μ)
     # convert initial guess to vector
-    tomandel!(cache.x_f, x0)
+    tomandel!(_cache.x_f, x0)
     # solve for variables x
     nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton))
     haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method
-    result = NLsolve.nlsolve(cache, cache.x_f; nlsolve_options...)
+    result = NLsolve.nlsolve(_cache, _cache.x_f; nlsolve_options...)
 
     if result.f_converged
         x_mandel = result.zero::Vector{T}
         x = frommandel(ResidualsCrystalViscoPlastic{S}, x_mandel)
-        dRdx = cache.DF
+        dRdx = _cache.DF
         inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx))
         ∂σ∂ε = inv_J_σσ ⊡ m.Eᵉ
         return x.σ, ∂σ∂ε, CrystalViscoPlasticState{3,T,6,S}(x.σ, SVector{S,T}(x.κ), SVector{S,T}(x.α), SVector{S,T}(x.μ))

src/CrystalViscoPlastic/CrystalViscoPlasticRed.jl

@@ -59,7 +59,7 @@ end
 Base.zero(::Type{CrystalViscoPlasticRedState{dim,T,M,S}}) where {dim,T,M,S} = CrystalViscoPlasticRedState{dim,T,M,S}(zero(SymmetricTensor{2,dim,T,M}), SVector{S,T}(zeros(T, S)), SVector{S,T}(zeros(T, S)), SVector{S,T}(zeros(T, S)))
 initial_material_state(::CrystalViscoPlasticRed{S}) where S = zero(CrystalViscoPlasticRedState{3,Float64,6,S})
 
-function get_cache(m::CrystalViscoPlasticRed)
+function _get_cache(m::CrystalViscoPlasticRed)
     state = initial_material_state(m)
     # it doesn't actually matter which state and strain step we use here,
     # f is overwritten in the constitutive driver before being used.
@@ -91,18 +91,20 @@ end
 function material_response(m::CrystalViscoPlasticRed{S}, Δε::SymmetricTensor{2,3,T,6}, state::CrystalViscoPlasticRedState{3},
     Δt=1.0; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where {S, T}
 
+    _cache = solver_cache(cache)
+
     # set the current residual function that depends only on the variables
     f(r_vector, x_vector) = vector_residual!((x->residuals(x,m,state,Δε,Δt)), r_vector, x_vector, m)
-    update_cache!(cache, f)
+    update_cache!(_cache, f)
     # initial guess
     σ_trial = state.σ + m.Eᵉ ⊡ Δε
     x0 = ResidualsCrystalViscoPlasticRed(state.μ)
     # convert initial guess to vector
-    tomandel!(cache.x_f, x0)
+    tomandel!(_cache.x_f, x0)
     # solve for variables x
     nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton))
     haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method
-    result = NLsolve.nlsolve(cache, cache.x_f; nlsolve_options...)
+    result = NLsolve.nlsolve(_cache, _cache.x_f; nlsolve_options...)
 
     if result.f_converged
         x_mandel = result.zero::Vector{T}
@@ -111,11 +113,11 @@ function material_response(m::CrystalViscoPlasticRed{S}, Δε::SymmetricTensor{2
         σ, κ, α = get_state_vars(x.μ, Δε, m, state)
         ############################################
         # tangent computation
-        ∂R∂x = cache.DF
+        ∂R∂x = _cache.DF
         f2(r_vector, x_vector) = vector_residual!((Δε->residuals(x,m,state,Δε,Δt)), r_vector, x_vector, Δε)
         # lucky coincidence that we can reuse the buffers here
-        tomandel!(view(cache.x_f, 1:6), Δε)
-        ∂R∂ε = ForwardDiff.jacobian(f2, cache.F, view(cache.x_f, 1:6)) # ∂R∂ε = ∂R∂Δε
+        tomandel!(view(_cache.x_f, 1:6), Δε)
+        ∂R∂ε = ForwardDiff.jacobian(f2, _cache.F, view(_cache.x_f, 1:6)) # ∂R∂ε = ∂R∂Δε
         ∂X∂ε = -inv(∂R∂x)*∂R∂ε 
         # in this case X is only μ
         ∂σ∂ε = m.Eᵉ

src/MaterialModels.jl

@@ -86,6 +86,7 @@ include("CohesiveMaterials/XuNeedleman.jl")
 
 include("nonlinear_solver.jl")
 include("wrappers.jl")
+include("iteration_caches.jl")
 
 export initial_material_state, get_cache
 export material_response

src/Plastic.jl

@@ -48,7 +48,7 @@ end
 Base.zero(::Type{PlasticState{dim,T,M}}) where {dim,T,M} = PlasticState(zero(SymmetricTensor{2,dim,T,M}), zero(T), zero(SymmetricTensor{2,dim,T,M}), zero(T))
 initial_material_state(::Plastic) = zero(PlasticState{3,Float64,6})
 
-function get_cache(m::Plastic)
+function _get_cache(m::Plastic)
     state = initial_material_state(m)
     # it doesn't actually matter which state and strain step we use here,
     # f is overwritten in the constitutive driver before being used.
@@ -122,6 +122,8 @@ See [NLsolve documentation](https://github.com/JuliaNLSolvers/NLsolve.jl#common-
 function material_response(m::Plastic, ε::SymmetricTensor{2,3,T,6}, state::PlasticState{3},
     Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where T
 
+    _cache = solver_cache(cache)
+
     σ_trial = m.Eᵉ ⊡ (ε - state.εᵖ)
     Φ = sqrt(3/2)*norm(dev(σ_trial-state.α)) - m.σ_y - state.κ
 
@@ -131,19 +133,19 @@ function material_response(m::Plastic, ε::SymmetricTensor{2,3,T,6}, state::Plas
         # set the current residual function that depends only on the variables
         # cache.f = (r_vector, x_vector) -> vector_residual!(((x)->residuals(x,m,state,Δε)), r_vector, x_vector, m)
         f(r_vector, x_vector) = vector_residual!(((x)->residuals(x,m,state,ε)), r_vector, x_vector, m)
-        update_cache!(cache, f)
+        update_cache!(_cache, f)
         # initial guess
         x0 = ResidualsPlastic(σ_trial, state.κ, state.α, state.μ)
         # convert initial guess to vector
-        tomandel!(cache.x_f, x0)
+        tomandel!(_cache.x_f, x0)
         # solve for variables x
         nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton))
         haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method
-        result = NLsolve.nlsolve(cache, cache.x_f; nlsolve_options...)
+        result = NLsolve.nlsolve(_cache, _cache.x_f; nlsolve_options...)
         if result.f_converged
             x = frommandel(ResidualsPlastic, result.zero::Vector{T})
             εᵖ = state.εᵖ + x.μ*(3/(2*sqrt(3/2)*norm(dev(x.σ-x.α)))*dev(x.σ-x.α))
-            dRdx = cache.DF
+            dRdx = _cache.DF
             inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx))
             ∂σ∂ε = inv_J_σσ ⊡ m.Eᵉ
             return x.σ, ∂σ∂ε, PlasticState(εᵖ, x.κ, x.α, x.μ)

src/iteration_caches.jl

@@ -0,0 +1,63 @@
+"""
+AbstractCache
+Stores matrices, vectors etc. to avoid re-allcating memory each time the material routine is called.
+"""
+abstract type AbstractCache end
+
+
+# wraps OnceDiffentiable, can't make this a subtype of AbstractCache without wrapping
+struct SolverCache{C} <: AbstractCache
+    cache::C
+end
+
+struct CacheContainer{C1, C2} <: AbstractCache
+    solver_cache::C1
+    plane_stress_cache::C2
+end
+
+# fallback in case there is no cache defined
+struct PlaneStressCache{TF, TDF, TX}
+    F::TF
+    DF::TDF
+    x_f::TX
+end
+
+# generic fallback, Materials without field σ need to define it
+get_stress_type(state::AbstractMaterialState) = typeof(state.σ)
+
+# fallback for optional cache argument
+function get_plane_stress_cache(m::AbstractMaterial)
+    state = initial_material_state(m)
+    stress_type = get_stress_type(state)
+    T = eltype(stress_type)
+    M = Tensors.n_components(Tensors.get_base(stress_type))
+    cache = PlaneStressCache(Vector{T}(undef,M), Matrix{T}(undef,M,M), Vector{T}(undef,M))
+    return cache
+end
+
+# cache constructors
+_get_cache(::AbstractMaterial) = nothing # fallback, models that require iterative solving have to implement this
+
+function get_cache(m::AbstractMaterial)
+    c = _get_cache(m)
+    isnothing(c) && @warn("$(typeof(m)) doesn't require a cache.")
+    return SolverCache(_get_cache(m))
+end
+
+get_cache(m::AbstractMaterial, dim::AbstractDim) = get_cache(m)
+
+function get_cache(m::AbstractMaterial, dim::Union{PlaneStress, UniaxialStress})
+    solver_cache = _get_cache(m)
+    plane_stress_cache = get_plane_stress_cache(m)
+    if isnothing(solver_cache)
+        return plane_stress_cache
+    else
+        return CacheContainer(solver_cache, plane_stress_cache)
+    end
+end
+
+# accessor functions
+solver_cache(c::SolverCache) = c.cache
+solver_cache(c::CacheContainer) = c.solver_cache
+plane_stress_cache(c::PlaneStressCache) = c
+plane_stress_cache(c::CacheContainer) = c.plane_stress_cache

src/wrappers.jl

@@ -45,7 +45,7 @@ function material_response(
     Δε::AbstractTensor{2,d,T},
     state::AbstractMaterialState,
     Δt = nothing;
-    cache =  get_cache(m),
+    cache =  get_cache(m, dim),
     options = Dict{Symbol, Any}(),
     ) where {d, T}
 
@@ -91,22 +91,3 @@ get_nonzero_indices(::UniaxialStress{1}, ::SymmetricTensor{2,3}) = [1] # for voi
 get_zero_indices(::UniaxialStress{1}, ::Tensor{2,3}) = collect(2:9) # for voigt/mandel format, do not use on tensor data!
 get_nonzero_indices(::UniaxialStress{1}, ::Tensor{2,3}) = [1] # for voigt/mandel format, do not use on tensor data!
 
-# fallback in case there is no cache defined
-struct PlaneStressCache{TF, TDF, TX}
-    F::TF
-    DF::TDF
-    x_f::TX
-end
-
-# generic fallback, Materials without field σ need to define it
-get_stress_type(state::AbstractMaterialState) = typeof(state.σ)
-
-# fallback for optional cache argument
-function get_cache(m::AbstractMaterial)
-    state = initial_material_state(m)
-    stress_type = get_stress_type(state)
-    T = eltype(stress_type)
-    M = Tensors.n_components(Tensors.get_base(stress_type))
-    cache = PlaneStressCache(Vector{T}(undef,M), Matrix{T}(undef,M,M), Vector{T}(undef,M))
-    return cache
-end