Update to new MMB diffinterface

Project.toml

@@ -21,10 +21,12 @@ Newton = "0.2"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+MaterialModelsTesting = "882b014b-b96c-4115-8629-e17fb35110d2"
 
 [targets]
-test = ["Test", "ForwardDiff", "FiniteDiff"]
+test = ["Test", "ForwardDiff", "FiniteDiff", "MaterialModelsTesting"]
 
 [sources]
-MaterialModelsBase = {rev = "main", url = "https://github.com/knutam/MaterialModelsBase.jl"}
-Newton = {rev = "main", url = "https://github.com/knutam/Newton.jl"}
+MaterialModelsBase = {url = "https://github.com/knutam/MaterialModelsBase.jl"}
+Newton = {url = "https://github.com/knutam/Newton.jl"}
+MaterialModelsTesting = {url = "https://github.com/KnutAM/MaterialModelsTesting.jl"}

docs/Project.toml

@@ -8,5 +8,5 @@ Newton = "83aa5b51-0588-403c-85e4-434ec185aae7"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 
 [sources]
-MaterialModelsBase = {rev = "main", url = "https://github.com/knutam/MaterialModelsBase.jl"}
-Newton = {rev = "main", url = "https://github.com/knutam/Newton.jl"}
+MaterialModelsBase = {url = "https://github.com/knutam/MaterialModelsBase.jl"}
+Newton = {url = "https://github.com/knutam/Newton.jl"}

docs/src/devdocs.md

@@ -15,6 +15,4 @@ compute_stress
 compute_tangent
 compute_stress_and_tangent
 static_vector
-get_resid_eltype
-get_num_unknowns
 ```

src/CrystalPlasticity.jl

@@ -20,7 +20,7 @@ The parameters for this model are,
 See [Meyer (2020)](https://doi.org/10.1016/j.ijsolstr.2020.04.037) for a description of this model
 in a finite strain framework. 
 """
-@kwdef struct CrystalPlasticity{C, E, T, OS} <: AbstractMaterial
+@kwdef struct CrystalPlasticity{C, E, T, OS, T_tol} <: AbstractMaterial
     crystal::C          # E.g. BCC, FCC, etc.
     elastic::E          # Elastic definition
     yield::T            # Initial yield limit
@@ -33,6 +33,8 @@ in a finite strain framework.
     Hkin::T # Kinematic hardening modulus 
     β∞::T   # Kinematic saturation stress
     overstress::OS      # Overstress function
+    maxiter::Int = 100
+    tolerance::T_tol = sqrt(eps(typeof(q)))
 end
 
 struct CrystalPlasticityState{T, nslip} <: MMB.AbstractMaterialState
@@ -96,7 +98,7 @@ function MMB.material_response(mat::CrystalPlasticity, ϵ::SymmetricTensor{2,3},
         Δγ0 = zero(typeof(τ_trial_red))
         ν_trial = sign.(τ_trial_red)
         rf(x) = residual(mat, x, old, ϵ, ν_trial, Δt)
-        Δγ, ∂r∂x, converged = newtonsolve(rf, Δγ0)
+        Δγ, ∂r∂x, converged = newtonsolve(rf, Δγ0; maxiter = mat.maxiter, tol = mat.tolerance)
         if converged
             σ = calculate_stress(mat, Δγ, old, ϵ, ν_trial)
             # R(Δγ(ϵ), ϵ) = 0 => ∂R/∂x ∂γ/∂ϵ = - ∂R/∂ϵ

src/Elastic.jl

@@ -55,7 +55,7 @@ struct LinearElastic{T, case, N} <: AbstractMaterial
     p::SVector{N,T}
 end
 
-MMB.get_params_eltype(::LinearElastic{T}) where T = T
+MMB.get_vector_eltype(::LinearElastic{T}) where T = T
 
 # General symmetry
 LinearElastic{:general}(C::SymmetricTensor) = LinearElastic(C)
@@ -94,7 +94,7 @@ end
 calculate_stress(m::LinearElastic, ϵ::SymmetricTensor) = m.C⊡ϵ
 
 # Functions for conversion between material and parameter vectors
-MMB.get_num_params(::LinearElastic{<:Any,<:Any,N}) where{N} = N
+MMB.get_vector_length(::LinearElastic{<:Any,<:Any,N}) where{N} = N
 
 function MMB.fromvector(v::AbstractVector, ::LinearElastic{<:Any, :isotropic}; offset=0)
     return LinearElastic{:isotropic}(;E=v[offset+1], ν=v[offset+2])
@@ -115,12 +115,10 @@ function MMB.differentiate_material!(deriv::MaterialDerivatives{T}, m::LinearEla
     tomandel!(deriv.dσdϵ, dσdϵ)
     p = tovector(m)
     σ_from_param(p_vector) = tomandel(calculate_stress(fromvector(p_vector, m), ϵ))
-    ForwardDiff.jacobian!(deriv.dσdp,σ_from_param,p)
+    ForwardDiff.jacobian!(deriv.dσdp, σ_from_param, p)
 
-    fill!(deriv.dσdⁿs,  zero(T))
-    fill!(deriv.dsdϵ,  zero(T))
-    fill!(deriv.dsdp,  zero(T))
-    fill!(deriv.dsdⁿs,  zero(T))
+    fill!(deriv.dsdϵ, zero(T))
+    fill!(deriv.dsdp, zero(T))
 end
 
 # Show methods

src/ExtraOutputs.jl

@@ -10,25 +10,27 @@ mutable struct DiffOutputHelper{TX, T} <: AbstractExtraOutput
     const ∂R∂ⁿsᴹ::Matrix{T}
     const ∂R∂pᴹ::Matrix{T}
     const ∂s∂Xᴹ::Matrix{T}
-    const ∂s∂ⁿsᴹ::Matrix{T}
     const ∂s∂pᴹ::Matrix{T}
+    const dσdⁿsᴹ::Matrix{T}
+    const dsdⁿsᴹ::Matrix{T}
 end
 
 function DiffOutputHelper(m::AbstractMaterial)
-    T = MMB.get_params_eltype(m)
+    T = MMB.get_vector_eltype(m)
     X = initial_guess(m, initial_material_state(m), zero(SymmetricTensor{2,3}))
-    NR = get_num_unknowns(X)
+    NR = get_vector_length(X)
     Nσ = get_num_tensorcomponents(m)
     Ns = get_num_statevars(m)
-    Np = get_num_params(m)
+    Np = get_vector_length(m)
     dRdX_invᴹ  = zeros(T, NR, NR)
     ∂R∂ϵᴹ  = zeros(T, NR, Nσ)
     ∂R∂ⁿsᴹ = zeros(T, NR, Ns)
     ∂R∂pᴹ  = zeros(T, NR, Np)
     ∂s∂Xᴹ  = zeros(T, Ns, NR)
-    ∂s∂ⁿsᴹ = zeros(T, Ns, Ns)
     ∂s∂pᴹ  = zeros(T, Ns, Np)
-    return DiffOutputHelper(X, dRdX_invᴹ, false, ∂R∂ϵᴹ, ∂R∂ⁿsᴹ, ∂R∂pᴹ, ∂s∂Xᴹ, ∂s∂ⁿsᴹ, ∂s∂pᴹ)
+    dσdⁿsᴹ = zeros(T, Nσ, Ns)
+    dsdⁿsᴹ = zeros(T, Ns, Ns)
+    return DiffOutputHelper(X, dRdX_invᴹ, false, ∂R∂ϵᴹ, ∂R∂ⁿsᴹ, ∂R∂pᴹ, ∂s∂Xᴹ, ∂s∂pᴹ, dσdⁿsᴹ, dsdⁿsᴹ)
 end
 
 update_extras!(extras::DiffOutputHelper) = (extras.updated = false)
@@ -37,4 +39,4 @@ function update_extras!(extras::DiffOutputHelper, X, dRdX_invᴹ)
     extras.X = X
     extras.dRdX_invᴹ .= dRdX_invᴹ
     extras.updated = true
-end
+end

src/FiniteStrainPlastic.jl

@@ -28,24 +28,24 @@ m = FiniteStrainPlastic(elastic = CompressibleNeoHooke(G=80.e3, K=160.e3),
 ```
 
 """
-struct FiniteStrainPlastic{E,YLD,IsoH,KinH,OS} <: AbstractMaterial
+struct FiniteStrainPlastic{E,YLD,IsoH,KinH,OS,T} <: AbstractMaterial
     elastic::E          # Elastic definition
     yield::YLD          # Initial yield limit
     isotropic::IsoH     # Tuple of isotropic hardening definitions
     kinematic::KinH     # Tuple of kinematic hardening definitions
     overstress::OS      # Overstress function
+    maxiter::Int        # Maximum number of iterations for local problem (default 10)
+    tolerance::T        # Absolute tolerance for local problem (default √eps(T))
 end
-function FiniteStrainPlastic(;elastic::AbstractHyperElastic, yield, isotropic=nothing, kinematic=nothing, overstress=RateIndependent())
+function FiniteStrainPlastic(;elastic::AbstractHyperElastic, yield, isotropic=nothing, kinematic=nothing, overstress=RateIndependent(), maxiter = 100, tolerance = sqrt(eps(MMB.get_vector_eltype(elastic))))
     iso = maketuple_or_nothing(isotropic)
     kin = maketuple_or_nothing(kinematic)
     yld = default_yield_criteria(yield)
-    return FiniteStrainPlastic(elastic, yld, iso, kin, overstress)
+    return FiniteStrainPlastic(elastic, yld, iso, kin, overstress, maxiter, tolerance)
 end
 
-
-MMB.get_params_eltype(m::FiniteStrainPlastic) = typeof(initial_yield_limit(m.yield))
-
-MMB.get_tensorbase(::FiniteStrainPlastic) = Tensor{2,3}
+MMB.get_vector_eltype(m::FiniteStrainPlastic) = typeof(initial_yield_limit(m.yield))
+MMB.get_tensorbase(::FiniteStrainPlastic) = Tensor{2, 3}
 
 # Definition of material state
 struct FiniteStrainPlasticState{NKin,NIso,TFp,Tκ<:NTuple{NIso},TFk<:NTuple{NKin}} <: AbstractMaterialState
@@ -55,7 +55,7 @@ struct FiniteStrainPlasticState{NKin,NIso,TFp,Tκ<:NTuple{NIso},TFk<:NTuple{NKin
 end
 
 function MMB.initial_material_state(m::FiniteStrainPlastic)
-    T = MMB.get_params_eltype(m)
+    T = MMB.get_vector_eltype(m)
     I2 = one(Tensor{2,3,T})
     FiniteStrainPlasticState(
         I2,                                  # Fp
@@ -78,7 +78,7 @@ function FiniteStrainPlasticResidual(Mred, Δλ, κ::Tuple, Mk::Tuple)
     return FiniteStrainPlasticResidual(Mred, Δλ, map(x->convert(Tκ, x), κ), map(x->convert(TMk, x), Mk))
 end
 
-function get_num_unknowns(::FiniteStrainPlasticResidual{NKin_m1,NIso}) where{NKin_m1,NIso}
+function MMB.get_vector_length(::FiniteStrainPlasticResidual{NKin_m1,NIso}) where{NKin_m1,NIso}
     return 10 + NIso + 9*NKin_m1
 end
 
@@ -98,22 +98,12 @@ function MMB.fromvector(v::AbstractVector{T}, ::FiniteStrainPlasticResidual{NKin
     return FiniteStrainPlasticResidual(Mred, Δλ, κ, Mk)
 end
 
-function MMB.allocate_material_cache(m::FiniteStrainPlastic)
-    T = MMB.get_params_eltype(m)
-    s = initial_material_state(m)
-    F = one(Tensor{2,3})
-    x = initial_guess(m, s, F)
-    xv = Vector{T}(undef, get_num_unknowns(x))
-    rf!(r_vector, x_vector) = vector_residual!((x)->residual(x, m, old, F, zero(T)), r_vector, x_vector, x)
-    return NewtonCache(xv)
-end
-
 function mandel_stress(m_el::AbstractHyperElastic, Fe::Tensor{2,3})
     Ce = tdot(Fe)
     return Ce ⋅ compute_stress(m_el, Ce)
 end
 
-function MMB.material_response(m::FiniteStrainPlastic, F::Tensor{2,3}, old::FiniteStrainPlasticState, Δt, cache, extras)
+function MMB.material_response(m::FiniteStrainPlastic, F::Tensor{2,3}, old::FiniteStrainPlasticState, Δt, _, extras)
     Fpinv = inv(old.Fp)
     Fe = F ⋅ Fpinv
     M = mandel_stress(m.elastic, Fe)
@@ -127,7 +117,7 @@ function MMB.material_response(m::FiniteStrainPlastic, F::Tensor{2,3}, old::Fini
     else
         x0 = initial_guess(m, old, M)
         rf(x) = tovector(SVector, residual(fromvector(x, x0), m, old, F, Δt))
-        x_vector, ∂R∂X, converged = newtonsolve(rf, tovector(SVector, x0))
+        x_vector, ∂R∂X, converged = newtonsolve(rf, tovector(SVector, x0); maxiter = m.maxiter, tol = m.tolerance)
 
         if converged
             x_sol = fromvector(x_vector, x0)
@@ -152,7 +142,7 @@ function MMB.material_response(m::FiniteStrainPlastic, F::Tensor{2,3}, old::Fini
 
             return P, dPdF, new
         else
-            throw(MMB.NoLocalConvergence("$(typeof(m)): newtonsolve! didn't converge, F = ", F))
+            throw(MMB.NoLocalConvergence("$(typeof(m)): newtonsolve didn't converge, F = ", F))
         end
     end
 end
@@ -243,49 +233,49 @@ end
 
 #=
 # Functions for conversion between material and parameter vectors
-function MMB.get_num_params(m::Plastic)
+function MMB.get_vector_length(m::Plastic)
     return sum((
-        get_num_params(m.elastic),
-        get_num_params(m.yield),
+        get_vector_length(m.elastic),
+        get_vector_length(m.yield),
         sum(get_num_params, m.isotropic),
         sum(get_num_params, m.kinematic),
-        get_num_params(m.overstress)))
+        get_vector_length(m.overstress)))
 end
 
 function fromvectortuple(v::AbstractVector, materialtuple; offset=0)
     n = 0
     mout = map(materialtuple) do mt
         m = fromvector(v, mt; offset=(offset+n))
-        n += get_num_params(m)
+        n += get_vector_length(m)
         m
     end
     return mout, n
 end
 
 function MMB.fromvector(v::AbstractVector, m::Plastic; offset=0)
     i = offset
-    elastic = fromvector(v, m.elastic, offset=i); i += get_num_params(m.elastic)
-    yield = fromvector(v, m.yield, offset=i); i += get_num_params(m.yield)
+    elastic = fromvector(v, m.elastic, offset=i); i += get_vector_length(m.elastic)
+    yield = fromvector(v, m.yield, offset=i); i += get_vector_length(m.yield)
 
     isotropic, nisoparam = fromvectortuple(v, m.isotropic, offset=i); i+=nisoparam
     kinematic, nkinparam = fromvectortuple(v, m.kinematic, offset=i); i+=nkinparam
 
-    overstress = fromvector(v, m.overstress, offset=i); # i += get_num_params(m.overstress)
+    overstress = fromvector(v, m.overstress, offset=i); # i += get_vector_length(m.overstress)
 
     return Plastic(;elastic, yield, isotropic, kinematic, overstress)
 end
 
 function MMB.tovector!(v::AbstractVector, m::Plastic; offset=0)
     i = offset
-    tovector!(v, m.elastic; offset=i); i += get_num_params(m.elastic)
-    tovector!(v, m.yield; offset=i); i += get_num_params(m.yield)
+    tovector!(v, m.elastic; offset=i); i += get_vector_length(m.elastic)
+    tovector!(v, m.yield; offset=i); i += get_vector_length(m.yield)
     for iso in m.isotropic
-        tovector!(v, iso;offset=i); i+= get_num_params(iso)
+        tovector!(v, iso;offset=i); i+= get_vector_length(iso)
     end
     for kin in m.kinematic
-        tovector!(v, kin; offset=i); i+= get_num_params(kin)
+        tovector!(v, kin; offset=i); i+= get_vector_length(kin)
     end
-    tovector!(v, m.overstress; offset=i); # i += get_num_params(m.overstress)
+    tovector!(v, m.overstress; offset=i); # i += get_vector_length(m.overstress)
     return v
 end