Refactor cache

src/MaterialModels.jl

@@ -33,6 +33,12 @@ Defines the type of strain measure the a material uses, i.e Deformation gradient
 """
 abstract type StrainMeasure end
 
+"""
+    AbstractCache
+Stores matrices, vectors etc. to avoid re-allcating memory each time the material routine is called.
+"""
+abstract type AbstractCache end
+
 """
     material_response(m::AbstractMaterial, Δε::SymmetricTensor{2,3}, state::AbstractMaterialState, Δt; cache, options)
 
@@ -64,7 +70,7 @@ When multithreading is used, each threads needs its own cache.
 Returns `nothing` for materials that don't need a cache.
 """
 function get_cache(::AbstractMaterial)
-    return nothing
+    nothing 
 end
 """
     update_cache!(cache::OnceDifferentiable, f)

src/Plastic.jl

@@ -48,14 +48,18 @@ 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})
 
+struct PlasticCache{T<:OnceDifferentiable} <: AbstractCache
+    nlsolve_cache::T
+end
+
 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.
     f(r_vector, x_vector) = vector_residual!(((x)->MaterialModels.residuals(x, m, state, zero(SymmetricTensor{2,3}))), r_vector, x_vector, m)
     v_cache = Vector{Float64}(undef, get_n_scalar_equations(m))
     cache = NLsolve.OnceDifferentiable(f, v_cache, v_cache; autodiff = :forward)
-    return cache
+    return PlasticCache(cache)
 end
 
 # constitutive driver operates in 3D, so these can always be 3D
@@ -122,28 +126,30 @@ 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
 
+    nlsolve_cache = cache.nlsolve_cache
+
     σ_trial = m.Eᵉ ⊡ (ε - state.εᵖ)
     Φ = sqrt(3/2)*norm(dev(σ_trial-state.α)) - m.σ_y - state.κ
 
     if Φ <= 0
         return σ_trial, m.Eᵉ, state
     else
         # 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)
+        # nlsolve_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!(nlsolve_cache, f)
         # initial guess
         x0 = ResidualsPlastic(σ_trial, state.κ, state.α, state.μ)
         # convert initial guess to vector
-        tomandel!(cache.x_f, x0)
+        tomandel!(nlsolve_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(nlsolve_cache, nlsolve_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 = nlsolve_cache.DF
             inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx))
             ∂σ∂ε = inv_J_σσ ⊡ m.Eᵉ
             return x.σ, ∂σ∂ε, PlasticState(εᵖ, x.κ, x.α, x.μ)

src/wrappers.jl

@@ -28,7 +28,7 @@ function material_response(
     Δε::AbstractTensor{2,d,T},
     state::AbstractMaterialState,
     Δt = nothing;
-    cache = nothing,
+    cache = get_cache(m),
     options = Dict{Symbol, Any}(),
     ) where {d,T}
 
@@ -45,7 +45,7 @@ function material_response(
     Δε::AbstractTensor{2,d,T},
     state::AbstractMaterialState,
     Δt = nothing;
-    cache::Union{Any, Nothing} = nothing, #get_cache(m), #TODO: create AbstractCache type
+    cache =  get_cache(m),
     options = Dict{Symbol, Any}(),
     ) where {d, T}
 
@@ -113,4 +113,4 @@ get_nonzero_indices(::PlaneStress{2}, ::Tensor{2,3}) = [1, 2, 6, 9] # for voigt/
 get_zero_indices(::UniaxialStress{1}, ::SymmetricTensor{2,3}) = collect(2:6) # for voigt/mandel format, do not use on tensor data!
 get_nonzero_indices(::UniaxialStress{1}, ::SymmetricTensor{2,3}) = [1] # for voigt/mandel format, do not use on tensor data!
 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!
+get_nonzero_indices(::UniaxialStress{1}, ::Tensor{2,3}) = [1] # for voigt/mandel format, do not use on tensor data!