From 183eda5a4e4cd9b4f252e602fd88770e788b2a93 Mon Sep 17 00:00:00 2001 From: kimauth Date: Fri, 4 Feb 2022 17:28:40 +0100 Subject: [PATCH] seperate material + stress iteration caches --- .../CrystalViscoPlastic.jl | 12 ++-- .../CrystalViscoPlasticRed.jl | 16 ++--- src/MaterialModels.jl | 1 + src/Plastic.jl | 12 ++-- src/iteration_caches.jl | 63 +++++++++++++++++++ src/wrappers.jl | 21 +------ 6 files changed, 88 insertions(+), 37 deletions(-) create mode 100644 src/iteration_caches.jl diff --git a/src/CrystalViscoPlastic/CrystalViscoPlastic.jl b/src/CrystalViscoPlastic/CrystalViscoPlastic.jl index 8462cf8..9f627be 100644 --- a/src/CrystalViscoPlastic/CrystalViscoPlastic.jl +++ b/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.μ)) diff --git a/src/CrystalViscoPlastic/CrystalViscoPlasticRed.jl b/src/CrystalViscoPlastic/CrystalViscoPlasticRed.jl index eea3862..43012d9 100644 --- a/src/CrystalViscoPlastic/CrystalViscoPlasticRed.jl +++ b/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ᵉ diff --git a/src/MaterialModels.jl b/src/MaterialModels.jl index 3dfff4c..192409e 100644 --- a/src/MaterialModels.jl +++ b/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 diff --git a/src/Plastic.jl b/src/Plastic.jl index 5af01ad..4c76706 100644 --- a/src/Plastic.jl +++ b/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.μ) diff --git a/src/iteration_caches.jl b/src/iteration_caches.jl new file mode 100644 index 0000000..345e67a --- /dev/null +++ b/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 diff --git a/src/wrappers.jl b/src/wrappers.jl index 82bfb6c..45790db 100644 --- a/src/wrappers.jl +++ b/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