Initializing improved differentiation interface

.github/workflows/CI.yml

@@ -19,7 +19,7 @@ jobs:
       matrix:
         version:
           - '1'
-          - '1.8'
+          - '1.11'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -4,15 +4,20 @@ authors = ["Knut Andreas Meyer and contributors"]
 version = "0.2.3"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 
 [compat]
-julia = "1.8"
+julia = "1.11"
 
 [extras]
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestMaterials = "aa4f4413-326f-4ddd-a53e-c58801ebfeaf"
+FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
+
+[sources]
+TestMaterials = {path = "test/TestMaterials"}
 
 [targets]
-test = ["Test", "ForwardDiff"]
+test = ["Test", "FiniteDiff", "TestMaterials"]

docs/make.jl

@@ -19,6 +19,7 @@ makedocs(;
         "Conversion" => "conversion.md",
         "Differentiation" => "differentiation.md",
         "Implementation" => "implementing.md",
+        "Devdocs" => "devdocs.md"
     ],
 )
 

docs/src/devdocs.md

@@ -0,0 +1,8 @@
+# Internal documentation
+The following documentation returns to internals of the package and may change
+at any time.
+
+```@docs
+MaterialModelsBase.stress_controlled_indices
+MaterialModelsBase.strain_controlled_indices
+```

src/MaterialModelsBase.jl

@@ -1,5 +1,6 @@
 module MaterialModelsBase
 using Tensors, StaticArrays
+using ForwardDiff: ForwardDiff
 import Base: @kwdef
 
 # General interface for <:AbstractMaterial
@@ -25,7 +26,8 @@ export update_stress_state!                                 # For nonzero stress
 export tovector, tovector!, fromvector                      # Convert to/from `AbstractVector`s
 export get_num_tensorcomponents, get_num_statevars          # Information about the specific material
 export get_num_params, get_params_eltype                    # 
-export MaterialDerivatives, differentiate_material!         # Differentiation routines
+export MaterialDerivatives, StressStateDerivatives          # Derivative collections
+export differentiate_material!                              # Differentiation routines
 export allocate_differentiation_output                      # 
 
 abstract type AbstractMaterial end
@@ -153,8 +155,8 @@ abstract type AbstractExtraOutput end
 struct NoExtraOutput <: AbstractExtraOutput end
 
 include("vector_conversion.jl")
-include("differentiation.jl")
 include("stressiterations.jl")
+include("differentiation.jl")
 include("ErrorExceptions.jl")
 
 end

src/differentiation.jl

@@ -1,19 +1,31 @@
 struct MaterialDerivatives{T}
     dσdϵ::Matrix{T}
-    dσdⁿs::Matrix{T}
+    #dσdⁿs::Matrix{T}
     dσdp::Matrix{T}
     dsdϵ::Matrix{T}
-    dsdⁿs::Matrix{T}
+    #dsdⁿs::Matrix{T}
     dsdp::Matrix{T}
 end
 
+function Base.getproperty(d::MaterialDerivatives, key::Symbol)
+    if key === :dσdⁿs || key === :dsdⁿs
+        error("You are probably assuming MaterialModelsBase v0.2 behavior for differentiation")
+    else
+        @inline getfield(d, key)
+    end
+end
+
 """
     MaterialDerivatives(m::AbstractMaterial)
 
 A struct that saves all derivative information using a `Matrix{T}` for each derivative,
 where `T=get_params_eltype(m)`. The dimensions are obtained from `get_num_tensorcomponents`, 
-`get_num_statevars`, and `get_num_params`. The values should be updated in `differentiate_material!`
-by direct access of the fields, where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current 
+`get_num_statevars`, and `get_num_params`. `m` must support `tovector` and `fromvector`, while 
+the output of `initial_material_state` must support `tovector`, and in addition the element type
+of `tovector(initial_material_state(m))` must respect the element type in `tovector(m)` for any `m`.
+
+The values should be updated in `differentiate_material!` by direct access of the fields, 
+where `σ` is the stress, `ϵ` the strain, `s` and `ⁿs` are the current 
 and old state variables, and `p` the material parameter vector.
 
 * `dσdϵ`
@@ -29,13 +41,12 @@ function MaterialDerivatives(m::AbstractMaterial)
     n_tensor = get_num_tensorcomponents(m)
     n_state = get_num_statevars(m)
     n_params = get_num_params(m)
+    dsdp = ForwardDiff.jacobian(p -> tovector(initial_material_state(fromvector(p, m))), tovector(m))
     return MaterialDerivatives(
         zeros(T, n_tensor, n_tensor),  # dσdϵ
-        zeros(T, n_tensor, n_state),   # dσdⁿs
         zeros(T, n_tensor, n_params),  # dσdp
         zeros(T, n_state, n_tensor),   # dsdϵ
-        zeros(T, n_state, n_state),    # dsdⁿs
-        zeros(T, n_state, n_params)    # dsdp
+        dsdp
         )
 end
 
@@ -66,4 +77,107 @@ allocate_differentiation_output(::AbstractMaterial) = NoExtraOutput()
 Calculate the derivatives and save them in `diff`, see
 [`MaterialDerivatives`](@ref) for a description of the fields in `diff`.
 """
-function differentiate_material! end
+function differentiate_material! end
+
+struct StressStateDerivatives{T}
+    mderiv::MaterialDerivatives{T}
+    dϵdp::Matrix{T}
+    dσdp::Matrix{T}
+    ϵindex::SMatrix{3, 3, Int} # To allow indexing by (i, j) into only
+    σindex::SMatrix{3, 3, Int} # saved values to avoid storing unused rows.
+    # TODO: Reduce the dimensions, for now all entries (even those that are zero) are stored.
+end
+
+function StressStateDerivatives(::AbstractStressState, m::AbstractMaterial)
+    mderiv = MaterialDerivatives(m)
+    np = get_num_params(m)
+    nt = get_num_tensorcomponents(m) # Should be changed to only save non-controlled entries
+    dϵdp = zeros(nt, np)
+    dσdp = zeros(nt, np)
+    # Should be changed to only save non-controlled entries
+    vo = Tensors.DEFAULT_VOIGT_ORDER[3]
+    if nt == 6
+        index = SMatrix{3, 3, Int}(min(vo[i, j], vo[j, i]) for i in 1:3, j in 1:3)
+    else
+        index = SMatrix{3, 3, Int}(vo)
+    end
+    return StressStateDerivatives(mderiv, dϵdp, dσdp, index, index)
+end
+
+"""
+    differentiate_material!(ssd::StressStateDerivatives, stress_state, m, args...)
+
+For material models implementing `material_response(m, args...)` and `differentiate_material!(::MaterialDerivatives, m, args...)`,
+this method will work automatically by
+1) Calling `σ, dσdϵ, state = material_response(stress_state, m, args...)` (except that `dσdϵ::FourthOrderTensor{dim = 3}` is extracted)
+2) Calling `differentiate_material!(ssd.mderiv::MaterialDerivatives, m, args..., dσdϵ::FourthOrderTensor{3})`
+3) Updating `ssd` according to the constraints imposed by the `stress_state`.
+
+For material models that directly implement `material_response(stress_state, m, args...)`, this function should be overloaded directly
+to calculate the derivatives in `ssd`. Here the user has full control and no modifications occur automatically, however, typically the 
+(total) derivatives `ssd.dσdp`, `ssd.dϵdp`, and `ssd.mderiv.dsdp` should be updated. 
+"""
+function differentiate_material!(ssd::StressStateDerivatives, stress_state::AbstractStressState, m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where {N}
+    σ_full, dσdϵ_full, state, ϵ_full = stress_state_material_response(stress_state, m, ϵ, args...)
+    differentiate_material!(ssd.mderiv, m, ϵ_full, args..., dσdϵ_full)
+
+    if isa(stress_state, NoIterationState)
+        copy!(ssd.dσdp, ssd.mderiv.dσdp)
+        fill!(ssd.dϵdp, 0)
+    else
+        sc = stress_controlled_indices(stress_state, ϵ)
+        ec = strain_controlled_indices(stress_state, ϵ)
+        dσᶠdϵᶠ_inv = inv(get_unknowns(stress_state, dσdϵ_full)) # f: unknown strain components solved for during stress iterations
+        ssd.dϵdp[sc, :] .= .-dσᶠdϵᶠ_inv * ssd.mderiv.dσdp[sc, :]
+        ssd.dσdp[ec, :] .= ssd.mderiv.dσdp[ec, :] .+ ssd.mderiv.dσdϵ[ec, sc] * ssd.dϵdp[sc, :]
+        ssd.mderiv.dsdp .+= ssd.mderiv.dsdϵ[:, sc] * ssd.dϵdp[sc, :]
+    end
+    return reduce_tensordim(stress_state, σ_full), reduce_stiffness(stress_state, dσdϵ_full), state, ϵ_full
+end
+
+"""
+    stress_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int}
+
+Get the `N` indices that are stress-controlled in `stress_state`. The tensor input is used to 
+determine if a symmetric or full tensor is used. 
+"""
+function stress_controlled_indices end
+
+"""
+    strain_controlled_indices(stress_state::AbstractStressState, ::AbstractTensor)::SVector{N, Int}
+
+Get the `N` indices that are strain-controlled in `stress_state`. The tensor input is used to 
+determine if a symmetric or full tensor is used. 
+"""
+function strain_controlled_indices end
+
+# NoIterationState
+stress_controlled_indices(::NoIterationState, ::AbstractTensor) = SVector{0,Int}()
+strain_controlled_indices(::NoIterationState, ::SymmetricTensor) = @SVector([1, 2, 3, 4, 5, 6])
+strain_controlled_indices(::NoIterationState, ::Tensor) = @SVector([1, 2, 3, 4, 5, 6, 7, 8, 9])
+
+# UniaxialStress
+stress_controlled_indices(::UniaxialStress, ::SymmetricTensor) = @SVector([2, 3, 4, 5, 6])
+stress_controlled_indices(::UniaxialStress, ::Tensor) = @SVector([2, 3, 4, 5, 6, 7, 8, 9])
+strain_controlled_indices(::UniaxialStress, ::AbstractTensor) = @SVector([1])
+
+# UniaxialNormalStress
+stress_controlled_indices(::UniaxialNormalStress, ::AbstractTensor) = @SVector([2,3])
+strain_controlled_indices(::UniaxialNormalStress, ::SymmetricTensor) = @SVector([1, 4, 5, 6])
+strain_controlled_indices(::UniaxialNormalStress, ::Tensor) = @SVector([1, 4, 5, 6, 7, 8, 9])
+
+# PlaneStress 12 -> 6, 21 -> 9
+stress_controlled_indices(::PlaneStress, ::SymmetricTensor) = @SVector([3, 4, 5])
+stress_controlled_indices(::PlaneStress, ::Tensor) = @SVector([3, 4, 5, 7, 8])
+strain_controlled_indices(::PlaneStress, ::SymmetricTensor) = @SVector([1, 2, 6])
+strain_controlled_indices(::PlaneStress, ::Tensor) = @SVector([1, 2, 6, 9])
+
+# GeneralStressState
+stress_controlled_indices(ss::GeneralStressState{Nσ}, ::AbstractTensor) where Nσ = controlled_indices_from_tensor(ss.σ_ctrl, true, Val(Nσ))
+function strain_controlled_indices(ss::GeneralStressState{Nσ,TT}, ::AbstractTensor) where {Nσ,TT}
+    N = Tensors.n_components(Tensors.get_base(TT)) - Nσ
+    return controlled_indices_from_tensor(ss.σ_ctrl, false, Val(N))
+end
+function controlled_indices_from_tensor(ctrl::AbstractTensor, return_if, ::Val{N}) where N
+    return SVector{N}(i for (i, v) in pairs(tovoigt(SVector, ctrl)) if v == return_if)
+end

src/stressiterations.jl

@@ -25,7 +25,8 @@ See also [`ReducedStressState`](@ref).
 material_response(::AbstractStressState, ::AbstractMaterial, args...)
 
 @inline function material_response(stress_state::AbstractStressState, m::AbstractMaterial, args::Vararg{Any,N}) where N
-    return reduced_material_response(stress_state, m, args...)
+    σ, dσdϵ, state, ϵ_full = stress_state_material_response(stress_state, m, args...)
+    return reduce_tensordim(stress_state, σ), reduce_stiffness(stress_state, dσdϵ), state, ϵ_full
 end
 
 update_stress_state!(::AbstractStressState, σ) = nothing
@@ -168,7 +169,7 @@ get_tolerance(ss::UniaxialNormalStress) = get_tolerance(ss.settings)
 get_max_iter(ss::UniaxialNormalStress) = get_max_iter(ss.settings)
 
 """
-    GeneralStressState(σ_ctrl::AbstractTensor{2,3,Bool}, σ::AbstractTensor{2,3,Bool}; kwargs...)
+    GeneralStressState(σ_ctrl::AbstractTensor{2,3,Bool}, σ::AbstractTensor{2,3}; kwargs...)
 
 Construct a general stress state controlled by `σ_ctrl` whose component is `true` if that 
 component is stress-controlled and `false` if it is strain-controlled. If stress-controlled,
@@ -252,15 +253,15 @@ reduce_tensordim(::Val{dim}, a::SymmetricTensor{2}) where dim = SymmetricTensor{
 reduce_tensordim(::Val{dim}, A::SymmetricTensor{4}) where dim = SymmetricTensor{4,dim}((i,j,k,l)->A[i,j,k,l])
 
 
-function reduced_material_response(stress_state::NoIterationState, 
+function stress_state_material_response(stress_state::NoIterationState, 
         m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where N
 
     ϵ_full = expand_tensordim(stress_state, ϵ)
     σ, dσdϵ, new_state = material_response(m, ϵ_full, args...)
-    return reduce_tensordim(stress_state, σ), reduce_tensordim(stress_state, dσdϵ), new_state, ϵ_full
+    return σ, dσdϵ, new_state, ϵ_full
 end
 
-function reduced_material_response(stress_state::IterationState, 
+function stress_state_material_response(stress_state::IterationState, 
         m::AbstractMaterial, ϵ::AbstractTensor, args::Vararg{Any,N}) where N
 
     # Newton options
@@ -274,8 +275,7 @@ function reduced_material_response(stress_state::IterationState,
         σ_mandel = get_unknowns(stress_state, σ_full)
 
         if norm(σ_mandel) < tol
-            dσdϵ = reduce_stiffness(stress_state, dσdϵ_full)
-            return reduce_tensordim(stress_state, σ_full), dσdϵ, new_state, ϵ_full
+            return σ_full, dσdϵ_full, new_state, ϵ_full
         end
 
         dσdϵ_mandel = get_unknowns(stress_state, dσdϵ_full)
@@ -284,14 +284,22 @@ function reduced_material_response(stress_state::IterationState,
     throw(NoStressConvergence("Stress iterations with the NewtonSolver did not converge"))
 end
 
-reduce_stiffness(::State3D, dσdϵ_full::AbstractTensor{4,3}) = dσdϵ_full
+reduce_stiffness(ss::State3D, dσdϵ_full) = reduce_stiffness(Val(3), ss, dσdϵ_full)
+reduce_stiffness(ss::State2D, dσdϵ_full) = reduce_stiffness(Val(2), ss, dσdϵ_full)
+reduce_stiffness(ss::State1D, dσdϵ_full) = reduce_stiffness(Val(1), ss, dσdϵ_full)
 
-function reduce_stiffness(stress_state, dσdϵ_full::AbstractTensor{4,3})
+reduce_stiffness(::Val{3}, ::State3D, dσdϵ_full::AbstractTensor{4,3}) = dσdϵ_full
+
+function reduce_stiffness(::Union{Val{1}, Val{2}}, stress_state::IterationState, dσdϵ_full::AbstractTensor{4,3})
     ∂σᶠ∂ϵᶠ, ∂σᶠ∂ϵᶜ, ∂σᶜ∂ϵᶠ, ∂σᶜ∂ϵᶜ = extract_substiffnesses(stress_state, dσdϵ_full)
     dσᶜdϵᶜ = ∂σᶜ∂ϵᶜ - ∂σᶜ∂ϵᶠ * (∂σᶠ∂ϵᶠ \ ∂σᶠ∂ϵᶜ)
     return convert_stiffness(dσᶜdϵᶜ, stress_state, dσdϵ_full)
 end
 
+function reduce_stiffness(::Union{Val{1}, Val{2}}, stress_state::NoIterationState, dσdϵ_full::AbstractTensor{4,3})
+    return reduce_tensordim(stress_state, dσdϵ_full)
+end
+
 convert_stiffness(dσᶜdϵᶜ::SMatrix{1,1}, ::State1D, ::SymmetricTensor) = frommandel(SymmetricTensor{4,1}, dσᶜdϵᶜ)
 convert_stiffness(dσᶜdϵᶜ::SMatrix{1,1}, ::State1D, ::Tensor) = frommandel(Tensor{4,1}, dσᶜdϵᶜ)
 convert_stiffness(dσᶜdϵᶜ::SMatrix{3,3}, ::State2D, ::SymmetricTensor) = frommandel(SymmetricTensor{4,2}, dσᶜdϵᶜ)
@@ -411,17 +419,17 @@ function get_full_tensor(state::GeneralStressState{Nσ}, ::TT, v::SVector{Nσ,T}
 end
 
 function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{2,3,T}) where {Nσ, T}
-    shear_factor = T(√2)
+    shear_factor = a isa SymmetricTensor ? T(√2) : one(T)
     s(i,j) = i==j ? one(T) : shear_factor
     f(c) = ((i,j) = c; s(i,j)*(a[i,j]-state.σ[i,j]))
-    return SVector{Nσ}((f(c) for c in state.σm_inds))
+    return SVector{Nσ,T}((f(c) for c in state.σm_inds))
 end
 
 function get_unknowns(state::GeneralStressState{Nσ}, a::AbstractTensor{4,3,T}) where {Nσ,T}
     shear_factor = T(√2)
     s(i,j) = i==j ? one(T) : shear_factor
     f(c1,c2) = ((i,j) = c1; (k,l) = c2; a[i,j,k,l]*s(i,j)*s(k,l))
-    return SMatrix{Nσ,Nσ}((f(c1,c2) for c1 in state.σm_inds, c2 in state.σm_inds))
+    return SMatrix{Nσ,Nσ,T}((f(c1,c2) for c1 in state.σm_inds, c2 in state.σm_inds))
 end
 
 # Extract each part of the stiffness tensor as SMatrix

src/vector_conversion.jl

@@ -83,6 +83,8 @@ Puts the Mandel components of `a` into the vector `v`.
 """
 function tovector! end
 
+tovector!(v::AbstractVector, ::NoMaterialState; kwargs...) = v
+
 # Tensors.jl implementation
 function tovector!(v::AbstractVector, a::SecondOrderTensor; offset = 0)
     return tomandel!(v, a; offset)
@@ -107,6 +109,8 @@ Create a tensor with shape of `TT` with entries from the Mandel components in `v
 """
 function fromvector end
 
+fromvector(::AbstractVector, ::NoMaterialState; kwargs...) = NoMaterialState()
+
 # Tensors.jl implementation
 function fromvector(v::AbstractVector, ::TT; offset = 0) where {TT <: SecondOrderTensor}
     return frommandel(Tensors.get_base(TT), v; offset)