Fix multi-material support when combining models with different internal variables

Author: termi-official

src/discretization/operator.jl

@@ -161,6 +161,12 @@ function AssembledNonlinearOperator(integrator::NonlinearIntegrator, dh::Abstrac
     )
 end
 
+function Base.show(io::IO, cache::AssembledNonlinearOperator)
+    println(io, "AssembledNonlinearOperator:")
+    Base.show(io, typeof(cache.integrator))
+    Base.show(io, MIME"text/plain"(), cache.dh)
+end
+
 getJ(op::AssembledNonlinearOperator) = op.J
 
 function update_linearization!(op::AssembledNonlinearOperator, u::AbstractVector, time)

src/ferrite-addons/quadrature_iterator.jl

@@ -19,9 +19,9 @@ A helper to loop over the quadrature points in some rule or cache with type [`Qu
 struct QuadratureIterator{QR<:QuadratureRule}
     qr::QR
 end
-QuadratureIterator(fqr::FacetQuadratureRule, local_face_idx::Int) = QuadratureIterator(fqr.facet_rules[local_face_idx])
+QuadratureIterator(fqr::FacetQuadratureRule, local_face_idx::Int) = QuadratureIterator(fqr.face_rules[local_face_idx])
 QuadratureIterator(cv::CellValues) = QuadratureIterator(cv.qr)
-QuadratureIterator(fv::FacetValues) = QuadratureIterator(fv.fqr.facet_rules[fv.current_facet[]])
+QuadratureIterator(fv::FacetValues) = QuadratureIterator(fv.fqr.face_rules[fv.current_facet[]])
 
 function Base.iterate(iterator::QuadratureIterator, i = 1)
     i > getnquadpoints(iterator.qr) && return nothing

src/modeling/functions.jl

@@ -93,8 +93,7 @@ end
 
 solution_size(f::QuasiStaticFunction) = ndofs(f.dh)+ndofs(f.lvh)
 function local_function_size(f::QuasiStaticFunction)
-    length(f.lvh.dh.subdofhandlers) == 0 && return 0
-    return sum(Ferrite.n_components.(f.lvh.dh.subdofhandlers[1].field_interpolations); init=0)
+    error("Local function size of QuasiStaticFunction can vary!")
 end
 function default_initial_condition!(u::AbstractVector, f::QuasiStaticFunction)
     fill!(u, 0.0)
@@ -119,3 +118,17 @@ end
 
 gather_internal_variable_infos(model::QuasiStaticModel) = gather_internal_variable_infos(model.material_model)
 gather_internal_variable_infos(model::AbstractMaterialModel) = InternalVariableInfo[]
+
+@unroll function __get_material_model_multi(materials, domains, sdh)
+    idx = 1
+    @unroll for material ∈ materials
+        if first(domains[idx]) ∈ sdh.cellset
+            return material
+        end
+        idx += 1
+    end
+    error("MultiDomainIntegrator is broken: Requested to construct an internal cache for a SubDofHandler which is not associated with the integrator.")
+end
+__get_material_model(model::MultiMaterialModel, sdh) = __get_material_model_multi(model.materials, model.domains, sdh)
+__get_material_model(model::AbstractMaterialModel, sdh) = model
+get_material_model(f::QuasiStaticFunction, sdh) = __get_material_model(f.integrator.volume_model.material_model, sdh)

src/modeling/solid/elements.jl

@@ -117,7 +117,7 @@ function assemble_element!(residualₑ::AbstractVector, uₑ::AbstractVector, ge
     end
 end
 
-struct MultiMaterialModel{MaterialTuple}
+struct MultiMaterialModel{MaterialTuple} <: AbstractMaterialModel
     materials::MaterialTuple
     domains::Vector{OrderedSet{Int}} # These must match the subdofhandler sets and hence be disjoint
     domain_names::Vector{String}     # These must match the subdofhandler sets and hence be disjoint
@@ -174,17 +174,3 @@ end
 function setup_internal_cache(model::MultiMaterialModel, qr::QuadratureRule, sdh::SubDofHandler)
     return setup_internal_cache_multi(model.materials, model.domains, qr, sdh)
 end
-
-@unroll function __get_material_model_multi(materials, domains, sdh)
-    idx = 1
-    @unroll for material ∈ materials
-        if first(domains[idx]) ∈ sdh.cellset
-            return material
-        end
-        idx += 1
-    end
-    error("MultiDomainIntegrator is broken: Requested to construct an internal cache for a SubDofHandler which is not associated with the integrator.")
-end
-__get_material_model(model::MultiMaterialModel, sdh) = __get_material_model_multi(model.materials, model.domains, sdh)
-__get_material_model(model::AbstractMaterialModel, sdh) = model
-get_material_model(f::QuasiStaticFunction, sdh) = __get_material_model(f.integrator.volume_model.material_model, sdh)

src/modeling/solid/materials.jl

@@ -251,9 +251,9 @@ function _solve_local_sarcomere_dQdF(dQdλ, dλdF, λ, F, coefficients, active_t
 end
 
 function solve_local_constraint(F::Tensor{2,dim}, coefficients, material_model::ActiveStressModel, state_cache::GenericFirstOrderRateIndependentCondensationMaterialStateCache, geometry_cache, qp, time) where dim
+    Qflat, Qprevflat = _query_local_state(state_cache, geometry_cache, qp)
     # Early out if any of the previous local solves failed
     if state_cache.local_solver_cache.retcode ∉ (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
-        Qflat, _ = _query_local_state(state_cache, geometry_cache, qp)
         return Qflat, zero(Tensor{4,dim,Float64,4^dim})
     end
 
@@ -310,11 +310,9 @@ function solve_local_constraint(F::Tensor{2,dim}, coefficients, material_model::
         return Q, J
     end
 
-    Qflat, Qprevflat = _query_local_state(state_cache, geometry_cache, qp)
     Q, J = solve_internal_timestep(material_model, state_cache, λ, dλdt, Qflat, Qprevflat)
     # Abort if local solve failed
     if state_cache.local_solver_cache.retcode ∉ (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
-        Qflat, _ = _query_local_state(state_cache, geometry_cache, qp)
         return Qflat, zero(Tensor{4,dim,Float64,4^dim})
     end
     Qflat .= Q
@@ -336,13 +334,14 @@ function solve_local_constraint(F::Tensor{2,dim}, coefficients, material_model::
 end
 
 # Some debug materials
-Base.@kwdef struct LinearMaxwellMaterial{T} <: AbstractMaterialModel
+Base.@kwdef struct LinearMaxwellMaterial{T, sdim} <: AbstractMaterialModel
     E₀::T
     E₁::T
     μ::T
     η₁::T
     ν::T
 end
+LinearMaxwellMaterial(E₀::T, Eₗ::T, μ::T, η₁::T, ν::T) where T = LinearMaxwellMaterial{T,3}(E₀, Eₗ, μ, η₁, ν)
 
 local_function_size(model::QuasiStaticModel) = local_function_size(model.material_model)
 function local_function_size(model::AbstractMaterialModel)
@@ -452,7 +451,7 @@ function stress_function(material::LinearMaxwellMaterial, ε, coefficients, ε
     return E₀ * ℂ ⊡ ε + E₁ * ℂ ⊡ (ε - εᵛ)
 end
 
-function stress_and_tangent(material_model::LinearMaxwellMaterial, F::Tensor{2}, coefficients, εᵛ)
+function stress_and_tangent(material_model::LinearMaxwellMaterial, F::Tensor{2}, coefficients, εᵛ::SymmetricTensor{2})
     ε = symmetric(F - one(F))
     ∂σ∂ε, σ = Tensors.gradient(ε->stress_function(material_model, ε, coefficients, εᵛ), ε, :all)
     return σ, ∂σ∂ε
@@ -463,9 +462,13 @@ function setup_coefficient_cache(m::LinearMaxwellMaterial, qr::QuadratureRule, s
 end
 
 function setup_internal_cache(material_model::LinearMaxwellMaterial, qr::QuadratureRule, sdh::SubDofHandler)
-    return nothing # FIXME what should we do here? :)
+    return EmptyRateIndependentCondensationMaterialStateCache()
 end
 
-function gather_internal_variable_infos(model::LinearMaxwellMaterial)
-    return InternalVariableInfo(:εᵛ, 6) # TODO dimension info
+function gather_internal_variable_infos(model::LinearMaxwellMaterial{T, sdim}) where {T, sdim}
+    if sdim == 3
+        return InternalVariableInfo(:εᵛ, 6)
+    else
+        return InternalVariableInfo(:εᵛ, 4)
+    end
 end

src/solver/nonlinear/multilevel_newton_raphson.jl

@@ -12,6 +12,30 @@ Base.@kwdef mutable struct GenericLocalNonlinearSolverCache{JacobianType, Residu
     retcode::SciMLBase.ReturnCode.T = SciMLBase.ReturnCode.Default
 end
 
+function Base.show(io::IO, cache::GenericLocalNonlinearSolverCache)
+    println(io, "NewtonRaphsonSolverCache:")
+    Base.show(io, cache.params)
+    println(io, "J=$(typeof(cache.J)) with size $(size(cache.J))")
+    println(io, "R=$(typeof(cache.residual)) with size $(size(cache.residual))")
+    println(io, "C=$(typeof(cache.rhs_corrector)) with size $(size(cache.rhs_corrector))")
+    println(io, "outer_tol=$(cache.outer_tol)")
+    println(io, "status=$(cache.retcode)")
+end
+
+function check_local_solve_covergence(local_solver_cache::GenericLocalNonlinearSolverCache)
+    return local_solver_cache.retcode ∉ (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
+end
+function check_local_solve_covergence(local_solver_cache::Tuple)
+    return any(check_local_solve_covergence.(local_solver_cache))
+end
+
+function set_local_solver_tol(local_solver_cache::GenericLocalNonlinearSolverCache, tol)
+    local_solver_cache.outer_tol = tol
+end
+function set_local_solver_tol(local_solver_cache::Tuple, tol)
+    set_local_solver_tol.(local_solver_cache, tol)
+end
+
 """
     MultilevelNewtonRaphsonSolver{T}
 
@@ -29,6 +53,18 @@ struct MultiLevelNewtonRaphsonSolverCache{gCacheType, lCacheType} <: AbstractNon
     local_solver_cache::lCacheType
 end
 
+function Base.show(io::IO, cache::MultiLevelNewtonRaphsonSolverCache)
+    println(io, "MultiLevelNewtonRaphsonSolverCache:")
+    Base.show(io, cache.global_solver_cache)
+    if cache.local_solver_cache isa Tuple
+        for local_solver_cache in cache.local_solver_cache
+            Base.show(io, local_solver_cache)
+        end
+    else
+        Base.show(io, cache.local_solver_cache)
+    end
+end
+
 function nlsolve!(u::AbstractVector, f::AbstractSemidiscreteFunction, mlcache::MultiLevelNewtonRaphsonSolverCache, t)
     cache = mlcache.global_solver_cache
 
@@ -39,21 +75,21 @@ function nlsolve!(u::AbstractVector, f::AbstractSemidiscreteFunction, mlcache::M
     residualnormprev = 0.0
     Θ1prev = length(Θks) > 0 ? first(Θks) : 0.0
     resize!(Θks, 0)
-    mlcache.local_solver_cache.outer_tol = Inf
+    set_local_solver_tol(mlcache.local_solver_cache, Inf)
     while true
         cache.iter += 1
         residual .= 0.0
         @timeit_debug "update operator" update_linearization!(op, residual, u, t)
         # Check if local solve failed
-        if mlcache.local_solver_cache.retcode ∉ (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
+        if check_local_solve_covergence(mlcache.local_solver_cache)
             @debug "Some local newton did not converge. Aborting. ||r|| = $residualnorm" _group=:nlsolve
             return false
         end
         @timeit_debug "elimination" eliminate_constraints_from_linearization!(cache, f)
         linear_solver_cache.isfresh = true # Notify linear solver that we touched the system matrix
 
         residualnorm = residual_norm(cache, f)
-        mlcache.local_solver_cache.outer_tol = residualnorm
+        set_local_solver_tol(mlcache.local_solver_cache, residualnorm^2)
         if residualnorm < cache.parameters.tol && cache.iter > 1 # Do at least two iterations to get a sane convergence estimate
             break
         elseif cache.iter > cache.parameters.max_iter

src/solver/nonlinear/newton_raphson.jl

@@ -27,6 +27,12 @@ mutable struct NewtonRaphsonSolverCache{OpType, ResidualType, T, NewtonType <: N
     iter::Int
 end
 
+function Base.show(io::IO, cache::NewtonRaphsonSolverCache)
+    println(io, "NewtonRaphsonSolverCache:")
+    Base.show(io, cache.parameters)
+    Base.show(io, cache.op)
+end
+
 function setup_solver_cache(f::AbstractSemidiscreteFunction, solver::NewtonRaphsonSolver{T}) where {T}
     @unpack inner_solver = solver
     op = setup_operator(f, solver)

src/solver/time/euler.jl

@@ -180,16 +180,12 @@ end
 function setup_solver_cache(wrapper::BackwardEulerStageAnnotation, solver::AbstractNonlinearSolver)
     _setup_solver_cache(wrapper, wrapper.f, solver)
 end
-@inline function _setup_solver_cache(wrapper::BackwardEulerStageAnnotation, f::QuasiStaticFunction, solver::MultiLevelNewtonRaphsonSolver)
-    @unpack integrator, dh = f
-    @unpack volume_model, face_model = integrator
-    @unpack local_solver, newton = solver
-
-    # TODO add an abstraction layer to autoamte the steps below
-    singleQsize = local_function_size(f)
-    local_solver_cache = GenericLocalNonlinearSolverCache(
+function _setup_local_solver_cache(local_solver::GenericLocalNonlinearSolver, material_model::AbstractMaterialModel)
+    singleQsize = local_function_size(material_model)
+    @debug "Setting up local nonlinear solver with size(Q)=$(singleQsize) for material $(material_model)" _group=:nlsolve
+    return GenericLocalNonlinearSolverCache(
         # Solver parameters
-        solver.local_solver,
+        local_solver,
         # Buffers
         zeros(singleQsize, singleQsize),
         zeros(singleQsize),
@@ -199,6 +195,17 @@ end
         # Local convergence 
         SciMLBase.ReturnCode.Default,
     )
+end
+function _setup_local_solver_cache(local_solver::GenericLocalNonlinearSolver, material_models::MultiMaterialModel)
+    return map(material_model -> _setup_local_solver_cache(local_solver, material_model), material_models.materials)
+end
+@inline function _setup_solver_cache(wrapper::BackwardEulerStageAnnotation, f::QuasiStaticFunction, solver::MultiLevelNewtonRaphsonSolver)
+    @unpack integrator, dh = f
+    @unpack volume_model, face_model = integrator
+    @unpack local_solver, newton = solver
+
+    # TODO add an abstraction layer to autoamte the steps below
+    local_solver_cache = _setup_local_solver_cache(solver.local_solver, f.integrator.volume_model.material_model)
 
     # Extract condensable parts
     Q     = @view wrapper.u[(ndofs(dh)+1):end]
@@ -237,10 +244,13 @@ end
 
     newton_cache = NewtonRaphsonSolverCache(op, residual, newton, inner_cache, T[], 0)
 
-    return MultiLevelNewtonRaphsonSolverCache(
+    cache = MultiLevelNewtonRaphsonSolverCache(
         newton_cache, # setup_solver_cache(G, solver.newton),
         local_solver_cache, #setup_solver_cache(L, solver.local_newton), # FIXME pass
     )
+    @debug "Setting up Multi-Level Newton-Raphson solver." _group=:nlsolve
+    @debug cache _group=:nlsolve
+    return cache
 end
 
 # TODO Refactor the setup into generic parts and use multiple dispatch for the specifics.
@@ -347,7 +357,7 @@ end
             n_ivs_per_qp = local_function_size(material_model)
             return GenericFirstOrderRateIndependentCondensationMaterialStateCache(
                 # Pass the model
-                wrapper.f,
+                material_model,
                 # And some cache to speed up evaluation of f and associated coefficients
                 internal_cache,
                 # Pass global solution info
@@ -356,7 +366,7 @@ end
                 # Current time step length
                 wrapper.Δt,
                 # Local nonlinear solver cache
-                wrapper.local_solver_cache,
+                wrapper.local_solver_cache[idx],
                 # This one holds information about the local dofs inside u and uprev
                 wrapper.lvh,
                 # Buffer for Q and Qprev
@@ -367,13 +377,12 @@ end
         idx += 1
     end
     error("MultiDomainIntegrator is broken: Requested to construct an element cache for a SubDofHandler which is not associated with the integrator.")
-
 end
-function setup_internal_cache_backward_euler_unwrap(wrapper::BackwardEulerStageFunctionWrapper{<:QuasiStaticModel}, material_model::AbstractMaterialModel internal_cache::Union{RateIndependentCondensationMaterialStateCache, RateDependentCondensationMaterialStateCache}, qr::QuadratureRule, sdh::SubDofHandler)
-    n_ivs_per_qp = local_function_size(wrapper.f.material_model)
+function setup_internal_cache_backward_euler_unwrap(wrapper::BackwardEulerStageFunctionWrapper{<:QuasiStaticModel}, material_model::AbstractMaterialModel, internal_cache::Union{RateIndependentCondensationMaterialStateCache, RateDependentCondensationMaterialStateCache}, qr::QuadratureRule, sdh::SubDofHandler)
+    n_ivs_per_qp = local_function_size(material_model)
     return GenericFirstOrderRateIndependentCondensationMaterialStateCache(
         # Pass the model
-        wrapper.f,
+        material_model,
         # And some cache to speed up evaluation of f and associated coefficients
         internal_cache,
         # Pass global solution info
@@ -391,7 +400,7 @@ function setup_internal_cache_backward_euler_unwrap(wrapper::BackwardEulerStageF
     )
 end
 function setup_internal_cache(wrapper::BackwardEulerStageFunctionWrapper{<:QuasiStaticModel}, qr::QuadratureRule, sdh::SubDofHandler)
-    return setup_internal_cache_backward_euler_unwrap(wrapper, setup_internal_cache(wrapper.f.material_model, qr, sdh), qr, sdh)
+    return setup_internal_cache_backward_euler_unwrap(wrapper, wrapper.f.material_model, setup_internal_cache(wrapper.f.material_model, qr, sdh), qr, sdh)
 end
 
 function setup_boundary_cache(wrapper::BackwardEulerStageFunctionWrapper, fqr, sdh)