Overhaul steepest descent solver for EMS (WIP)

Author: lxmota

examples/ems/awful-cube/awful-cube.yaml

@@ -1,8 +1,7 @@
 type: single
 input mesh file: awful-cube.g
 output mesh file: awful-cube.e
-Exodus output interval: 10
-CSV output interval: 10
+Exodus output interval: 1
 model:
   type: mesh smoothing
   smooth reference: max
@@ -19,12 +18,8 @@ model:
 time integrator:
   type: quasi static
   initial time: 0.0
-  final time: 1.0e-00
-  time step: 1.0e-04
-  maximum time step: 1.0e-4
-  minimum time step: 1.0e-15
-  increase factor: 1.1
-  decrease factor: 0.5
+  final time: 10.0
+  time step: 1.0
 boundary conditions:
   Dirichlet:
     - node set: nsx-
@@ -49,7 +44,12 @@ solver:
   type: steepest descent
   step:  steepest descent
   minimum iterations: 1
-  maximum iterations: 16
+  maximum iterations: 64
   relative tolerance: 1.0e-12
   absolute tolerance: 1.0e-08
-  step length: 1.0e-8
+  step length: 1.0e-3
+  use line search: true
+  line search backtrack factor: 0.5
+  line search decrease factor: 1.0e-04
+  line search maximum iterations: 16
+

examples/single/implicit-dynamic-solid/cube/cube-matrix-free.yaml

@@ -35,8 +35,12 @@ boundary conditions:
 solver:
   type: steepest descent
   step: steepest descent
-  step length: 1.0
+  step length: 1.0e-01
+  use line search: true
+  line search backtrack factor: 0.5
+  line search decrease factor: 1.0e-04
+  line search maximum iterations: 64
   minimum iterations: 1
-  maximum iterations: 1024
+  maximum iterations: 64
   relative tolerance: 1.0e-15
   absolute tolerance: 1.0e-08

examples/single/static-solid/cube/cube-steepest-descent.yaml

@@ -195,12 +195,10 @@ function SolidMechanics(params::Parameters)
     end
     strain_energy = 0.0
     stiffness = spzeros(0, 0)
-    diag_stiffness = Float64[]
     mass = spzeros(0, 0)
     lumped_mass = Float64[]
     body_force = Float64[]
     compute_stiffness = true
-    compute_diag_stiffness = true
     compute_mass = true
     compute_lumped_mass = true
     mesh_smoothing = get(params, "mesh smoothing", false)
@@ -223,14 +221,12 @@ function SolidMechanics(params::Parameters)
         stored_energy,
         strain_energy,
         stiffness,
-        diag_stiffness,
         mass,
         lumped_mass,
         body_force,
         free_dofs,
         time,
         compute_stiffness,
-        compute_diag_stiffness,
         compute_mass,
         compute_lumped_mass,
         failed,
@@ -534,14 +530,6 @@ function add_internal_force!(Fi::MVector{M,T}, grad_op::SMatrix{9,M,T}, stress::
     @einsum Fi[i] += grad_op[j, i] * stress[j] * dV
 end
 
-function add_diag_stiff!(K::MVector{M,T}, grad_op::SMatrix{9,M,T}, C::SMatrix{9,9,T}, dV::T) where {M,T}
-    @inbounds for i in 1:M
-        g = grad_op[:, i]
-        K[i] += dV * dot(g, C * g)
-    end
-    return nothing
-end
-
 function add_lumped_mass!(M::MVector{R,T}, Nξ::SVector{N,T}, density::T, dV::T) where {R,N,T}
     @assert R == 3N
     s = sum(Nξ)
@@ -565,11 +553,9 @@ function compute_flags(model::SolidMechanics, integrator::TimeIntegrator, solver
     is_implicit = is_implicit_dynamic || is_implicit_static
     is_hessian_opt = solver isa HessianMinimizer
     is_matrix_free = solver isa SteepestDescent
-    need_diag_stiffness = is_implicit && is_matrix_free
     need_lumped_mass = is_explicit_dynamic || (is_implicit_dynamic && is_matrix_free)
     need_stiffness = is_implicit && is_hessian_opt
     need_mass = is_dynamic && is_hessian_opt
-    compute_diag_stiffness = need_diag_stiffness && model.compute_diag_stiffness
     compute_lumped_mass = need_lumped_mass && model.compute_lumped_mass
     compute_stiffness = need_stiffness && model.compute_stiffness
     compute_mass = need_mass && model.compute_mass
@@ -580,11 +566,9 @@ function compute_flags(model::SolidMechanics, integrator::TimeIntegrator, solver
         is_implicit,
         is_hessian_opt,
         is_matrix_free,
-        need_diag_stiffness,
         need_lumped_mass,
         need_stiffness,
         need_mass,
-        compute_diag_stiffness,
         compute_lumped_mass,
         compute_stiffness,
         compute_mass,
@@ -598,11 +582,6 @@ function create_threadlocal_arrays(model::SolidMechanics, flags::EvaluationFlags
     energy = zeros(nthreads())
     internal_force = create_threadlocal_coo_vectors(num_dofs)
 
-    diag_stiffness = create_threadlocal_coo_vectors(flags.compute_diag_stiffness ? num_dofs : 0)
-    if flags.compute_diag_stiffness && model.kinematics == Infinitesimal
-        model.compute_diag_stiffness = false
-    end
-
     lumped_mass = create_threadlocal_coo_vectors(flags.compute_lumped_mass ? num_dofs : 0)
     if flags.compute_lumped_mass
         model.compute_lumped_mass = false
@@ -623,19 +602,18 @@ function create_threadlocal_arrays(model::SolidMechanics, flags::EvaluationFlags
     if flags.compute_mass
         model.compute_mass = false
     end
-    return SMThreadLocalArrays(energy, internal_force, diag_stiffness, lumped_mass, stiffness, mass)
+    return SMThreadLocalArrays(energy, internal_force, lumped_mass, stiffness, mass)
 end
 
 function create_element_threadlocal_arrays(num_element_nodes::Int64, flags::EvaluationFlags)
     valN = Val(num_element_nodes)
     energy = zeros(nthreads())
     dofs = create_threadlocal_element_vectors(Int64, valN)
     internal_force = create_threadlocal_element_vectors(Float64, valN)
-    diag_stiffness = create_threadlocal_element_vectors(Float64, flags.compute_diag_stiffness ? valN : Val(0))
     lumped_mass = create_threadlocal_element_vectors(Float64, flags.compute_lumped_mass ? valN : Val(0))
     stiffness = create_threadlocal_element_matrices(Float64, flags.compute_stiffness ? valN : Val(0))
     mass = create_threadlocal_element_matrices(Float64, flags.compute_mass ? valN : Val(0))
-    return SMElementThreadLocalArrays(energy, dofs, internal_force, diag_stiffness, lumped_mass, stiffness, mass)
+    return SMElementThreadLocalArrays(energy, dofs, internal_force, lumped_mass, stiffness, mass)
 end
 
 function reset_element_threadlocal_arrays!(
@@ -651,9 +629,6 @@ function reset_element_threadlocal_arrays!(
     element_arrays_tl.dofs[t] = reshape(3 .* node_indices' .- [2, 1, 0], :)
     element_arrays_tl.energy[t] = 0.0
     fill!(element_arrays_tl.internal_force[t], 0.0)
-    if flags.compute_diag_stiffness == true
-        fill!(element_arrays_tl.diag_stiffness[t], 0.0)
-    end
     if flags.compute_lumped_mass == true
         fill!(element_arrays_tl.lumped_mass[t], 0.0)
     end
@@ -682,10 +657,6 @@ function compute_element_threadlocal_arrays!(
     stress = SVector{9,Float64}(P)
     element_arrays_tl.energy[t] += W * dvol
     add_internal_force!(element_arrays_tl.internal_force[t], grad_op, stress, dvol)
-    if flags.compute_diag_stiffness == true
-        moduli = second_from_fourth(AA)
-        add_diag_stiff!(element_arrays_tl.diag_stiffness[t], grad_op, moduli, dvol)
-    end
     if flags.compute_lumped_mass == true
         add_lumped_mass!(element_arrays_tl.lumped_mass[t], Np, density, dvol)
     end
@@ -709,9 +680,6 @@ function assemble_element_threadlocal_arrays!(
     t = threadid()
     arrays_tl.energy[t] += element_arrays_tl.energy[t]
     assemble!(arrays_tl.internal_force[t], element_arrays_tl.internal_force[t], element_arrays_tl.dofs[t])
-    if flags.compute_diag_stiffness == true
-        assemble!(arrays_tl.diag_stiffness[t], element_arrays_tl.diag_stiffness[t], element_arrays_tl.dofs[t])
-    end
     if flags.compute_lumped_mass == true
         assemble!(arrays_tl.lumped_mass[t], element_arrays_tl.lumped_mass[t], element_arrays_tl.dofs[t])
     end
@@ -729,9 +697,6 @@ function merge_threadlocal_arrays(
 )
     model.strain_energy = sum(arrays_tl.energy)
     model.internal_force = merge_threadlocal_coo_vectors(arrays_tl.internal_force, num_dofs)
-    if flags.compute_diag_stiffness == true
-        model.diag_stiffness = merge_threadlocal_coo_vectors(arrays_tl.diag_stiffness, num_dofs)
-    end
     if flags.compute_lumped_mass == true
         model.lumped_mass = merge_threadlocal_coo_vectors(arrays_tl.lumped_mass, num_dofs)
     end
@@ -787,7 +752,6 @@ function evaluate(model::SolidMechanics, integrator::TimeIntegrator, solver::Sol
                 J = det(F)
                 if J ≤ 0.0 || isfinite(J) == false
                     model.failed = true
-                    model.compute_stiffness = model.compute_diag_stiffness = true
                     model.compute_mass = model.compute_lumped_mass = true
                     norma_log(0, :error, "Non-positive Jacobian detected! This may indicate element distortion.")
                     norma_logf(4, :warning, "det(F) = %.3e", J)

src/model.jl

@@ -34,11 +34,9 @@ struct EvaluationFlags
     is_implicit::Bool
     is_hessian_opt::Bool
     is_matrix_free::Bool
-    need_diag_stiffness::Bool
     need_lumped_mass::Bool
     need_stiffness::Bool
     need_mass::Bool
-    compute_diag_stiffness::Bool
     compute_lumped_mass::Bool
     compute_stiffness::Bool
     compute_mass::Bool
@@ -48,17 +46,15 @@ end
 struct SMThreadLocalArrays{V,M}
     energy::Vector{Float64}
     internal_force::Vector{V}
-    diag_stiffness::Vector{V}
     lumped_mass::Vector{V}
     stiffness::Vector{M}
     mass::Vector{M}
 end
 
-struct SMElementThreadLocalArrays{T,DOFV,IFV,DSV,LMV,SM,MM}
+struct SMElementThreadLocalArrays{T,DOFV,IFV,LMV,SM,MM}
     energy::Vector{T}
     dofs::Vector{DOFV}
     internal_force::Vector{IFV}
-    diag_stiffness::Vector{DSV}
     lumped_mass::Vector{LMV}
     stiffness::Vector{SM}
     mass::Vector{MM}
@@ -78,14 +74,12 @@ mutable struct SolidMechanics <: Model
     stored_energy::Vector{Vector{Float64}}
     strain_energy::Float64
     stiffness::SparseMatrixCSC{Float64,Int64}
-    diag_stiffness::Vector{Float64}
     mass::SparseMatrixCSC{Float64,Int64}
     lumped_mass::Vector{Float64}
     body_force::Vector{Float64}
     free_dofs::BitVector
     time::Float64
     compute_stiffness::Bool
-    compute_diag_stiffness::Bool
     compute_mass::Bool
     compute_lumped_mass::Bool
     failed::Bool

src/model_types.jl

@@ -82,7 +82,6 @@ function SteepestDescent(params::Parameters, model::Model)
     relative_error = 0.0
     value = 0.0
     gradient = zeros(num_dof)
-    lumped_hessian = zeros(num_dof)
     solution = zeros(num_dof)
     initial_norm = 0.0
     converged = false
@@ -105,7 +104,6 @@ function SteepestDescent(params::Parameters, model::Model)
         relative_error,
         value,
         gradient,
-        lumped_hessian,
         solution,
         initial_norm,
         converged,
@@ -496,10 +494,8 @@ function evaluate(integrator::QuasiStatic, solver::MatrixFree, model::SolidMecha
     external_force = model.body_force + model.boundary_force
     if model.inclined_support == true
         solver.gradient = model.global_transform * (model.internal_force - external_force)
-        solver.lumped_hessian = model.global_transform * model.diag_stiffness
     else
         solver.gradient = model.internal_force - external_force
-        solver.lumped_hessian = model.diag_stiffness
     end
     return nothing
 end
@@ -530,32 +526,6 @@ function evaluate(integrator::Newmark, solver::HessianMinimizer, model::SolidMec
     return nothing
 end
 
-function evaluate(integrator::Newmark, solver::MatrixFree, model::SolidMechanics)
-    evaluate(model, integrator, solver)
-    if model.failed == true
-        return nothing
-    end
-    integrator.stored_energy = model.strain_energy
-    β = integrator.β
-    Δt = integrator.time_step
-    inertial_force = model.lumped_mass .* integrator.acceleration
-    kinetic_energy = 0.5 * model.lumped_mass ⋅ (integrator.velocity .* integrator.velocity)
-    integrator.kinetic_energy = kinetic_energy
-    internal_force = model.internal_force
-    external_force = model.body_force + model.boundary_force
-    diag_stiffness = model.diag_stiffness
-    if model.inclined_support == true
-        global_transform = model.global_transform
-        internal_force = global_transform * internal_force
-        external_force = global_transform * external_force
-        diag_stiffness = global_transform * diag_stiffness
-    end
-    solver.lumped_hessian = diag_stiffness + model.lumped_mass / β / Δt / Δt
-    solver.gradient = internal_force - external_force + inertial_force
-    solver.value = model.strain_energy - external_force ⋅ integrator.displacement + kinetic_energy
-    return nothing
-end
-
 function evaluate(integrator::CentralDifference, solver::ExplicitSolver, model::SolidMechanics)
     evaluate(model, integrator, solver)
     if model.failed == true
@@ -607,7 +577,8 @@ function backtrack_line_search(
         copy_solution_source_targets(solver, model, integrator)
         evaluate(integrator, solver, model)
         if model.failed == true
-            return step
+            step_length *= backtrack_factor
+            continue
         end
         resid = solver.gradient[free]
         merit = 0.5 * dot(resid, resid)
@@ -640,7 +611,7 @@ end
 
 function compute_step(integrator::TimeIntegrator, model::SolidMechanics, solver::MatrixFree, _::SteepestDescentStep)
     free = model.free_dofs
-    direction = -solver.gradient[free] ./ solver.lumped_hessian[free]
+    direction = LinearAlgebra.normalize(-solver.gradient[free])
     return backtrack_line_search(integrator, solver, model, direction)
 end
 

src/solver.jl

@@ -60,7 +60,6 @@ mutable struct SteepestDescent <: MatrixFree
     relative_error::Float64
     value::Float64
     gradient::Vector{Float64}
-    lumped_hessian::Vector{Float64}
     solution::Vector{Float64}
     initial_norm::Float64
     converged::Bool