reuse_chol -> reuse_fact and avoid try catch in direct solver (#171)

* reuse_chol -> reuse_fact and avoid try catch

* bump version

* fix typo

* formatting

* remove ldiv!

* lower bound NonconvexPercival

* split the tests some more

.github/workflows/CI.yml

@@ -25,6 +25,7 @@ jobs:
           - Extended_Tests
           - Examples_1
           - Examples_2
+          - Examples_3
           - WCSMO14_1
           - WCSMO14_2
     steps:

Project.toml

@@ -1,7 +1,7 @@
 name = "TopOpt"
 uuid = "53a1e1a5-51bb-58a9-8a02-02056cc81109"
 authors = ["mohamed82008 <[email protected]>", "yijiangh <[email protected]>"]
-version = "0.8.3"
+version = "0.9.0"
 
 [deps]
 AbstractDifferentiation = "c29ec348-61ec-40c8-8164-b8c60e9d9f3d"
@@ -65,7 +65,7 @@ MappedArrays = "0.4"
 NearestNeighbors = "0.4"
 Nonconvex = "2"
 NonconvexMMA = "1"
-NonconvexPercival = "0.1"
+NonconvexPercival = "0.1.4"
 NonconvexSemidefinite = "0.1.7"
 Parameters = "0.12"
 Preconditioners = "0.3, 0.4, 0.5, 0.6"

src/FEA/direct_displacement_solver.jl

@@ -49,13 +49,12 @@ function (s::DirectDisplacementSolver{T})(
     ::Type{Val{safe}}=Val{false},
     ::Type{newT}=T;
     assemble_f=true,
-    reuse_chol=false,
+    reuse_fact=false,
     rhs=assemble_f ? s.globalinfo.f : s.rhs,
     lhs=assemble_f ? s.u : s.lhs,
     kwargs...,
 ) where {T,safe,newT}
     globalinfo = s.globalinfo
-    N = size(globalinfo.K, 1)
     assemble!(
         globalinfo,
         s.problem,
@@ -75,40 +74,24 @@ function (s::DirectDisplacementSolver{T})(
         end
     end
     nans = false
-    if !reuse_chol
-        try
-            if T === newT
-                if s.qr
-                    globalinfo.qrK = qr(K.data)
-                else
-                    globalinfo.cholK = cholesky(Symmetric(K))
-                end
-            else
-                if s.qr
-                    globalinfo.qrK = qr((newT.(K)).data)
-                else
-                    globalinfo.cholK = cholesky(Symmetric(newT.(K)))
-                end
-            end
-        catch err
-            lhs .= T(NaN)
-            nans = true
-        end
-    end
-    if !nans
-        if T === newT
-            if s.qr
-                lhs .= globalinfo.qrK \ rhs
-            else
-                lhs .= globalinfo.cholK \ rhs
-            end
+    if !reuse_fact
+        newK = T === newT ? K : newT.(K)
+        if s.qr
+            globalinfo.qrK = qr(newK.data)
         else
-            if s.qr
-                lhs .= globalinfo.qrK \ newT.(rhs)
+            cholK = cholesky(Symmetric(K); check=false)
+            if issuccess(cholK)
+                globalinfo.cholK = cholK
             else
-                lhs .= globalinfo.cholK \ newT.(rhs)
+                @warn "The global stiffness matrix is not positive definite. Please check your boundary conditions."
+                lhs .= T(NaN)
+                nans = true
             end
         end
     end
+    nans && return nothing
+    new_rhs = T === newT ? rhs : newT.(rhs)
+    fact = s.qr ? globalinfo.qrK : globalinfo.cholK
+    lhs .= fact \ new_rhs
     return nothing
 end

src/Functions/block_compliance.jl

@@ -143,7 +143,7 @@ function compute_approx_bc(bc, F, V, Y)
     bc.method.sample_once || bc.method.sample_method(V)
     for i in 1:nv
         @views mul!(solver.rhs, F, V[:, i])
-        solver(; assemble_f=false, reuse_chol=(i > 1))
+        solver(; assemble_f=false, reuse_fact=(i > 1))
         invKFv = solver.lhs
         Y[:, i] .= invKFv
         temp = F' * invKFv
@@ -161,7 +161,7 @@ function compute_jtvp!_bc(out, bc, method::DiagonalEstimation, w)
         temp .= 0
         #q_i = K^-1 F (w .* v_i)
         @views mul!(solver.rhs, F, w .* V[:, i])
-        solver(; assemble_f=false, reuse_chol=(i > 1))
+        solver(; assemble_f=false, reuse_fact=(i > 1))
         Q[:, i] = solver.lhs
         #<q_i, dK/dx_e, y_i>
         @views compute_inner(temp, Q[:, i], Y[:, i], solver)

src/Functions/displacement.jl

@@ -79,7 +79,7 @@ function ChainRulesCore.rrule(dp::Displacement, x::PseudoDensities)
         else
             solver.rhs .= Δ
         end
-        solver(; reuse_chol=true, assemble_f=false)
+        solver(; reuse_fact=true, assemble_f=false)
         dudx_tmp .= 0
         for e in 1:length(x.x)
             _, dρe = get_ρ_dρ(x.x[e], penalty, xmin)

src/Functions/mean_compliance.jl

@@ -75,7 +75,7 @@ function compute_exact_ec(ec, x, grad, F, n)
     grad .= 0
     for i in 1:size(F, 2)
         @views solver.rhs .= F[:, i]
-        solver(; assemble_f=false, reuse_chol=(i > 1))
+        solver(; assemble_f=false, reuse_fact=(i > 1))
         u = solver.lhs
         obj += compute_compliance(
             cell_comp, grad_temp, cell_dofs, Kes, u, black, white, varind, x, penalty, xmin
@@ -104,7 +104,7 @@ function compute_approx_ec(ec, x, grad, F, V, n)
     ec.method.sample_once || ec.method.sample_method(V)
     for i in 1:nv
         @views mul!(solver.rhs, F, V[:, i])
-        solver(; assemble_f=false, reuse_chol=(i > 1))
+        solver(; assemble_f=false, reuse_fact=(i > 1))
         invKFv = solver.lhs
         obj += compute_compliance(
             cell_comp,

src/Functions/truss_stress.jl

@@ -81,7 +81,7 @@ end
 #     u = dp(x)
 #     return u, Δ -> begin # v
 #         solver.rhs .= Δ
-#         solver(reuse_chol = true, assemble_f = false)
+#         solver(reuse_fact = true, assemble_f = false)
 #         dudx_tmp .= 0
 #         for e in 1:length(x)
 #             _, dρe = get_ρ_dρ(x[e], penalty, xmin)

test/runtests.jl

@@ -33,18 +33,21 @@ if GROUP == "All" || GROUP == "Examples_1"
     @safetestset "CSIMP example" begin
         include("examples/csimp.jl")
     end
+end
+
+if GROUP == "All" || GROUP == "Examples_2"
     @safetestset "Global stress example" begin
         include("examples/global_stress.jl")
     end
     @safetestset "Local stress example" begin
         include("examples/local_stress.jl")
     end
+end
+
+if GROUP == "All" || GROUP == "Examples_3"
     @safetestset "More examples" begin
         include("examples/test_examples.jl")
     end
-end
-
-if GROUP == "All" || GROUP == "Examples_2"
     @safetestset "Neural network example" begin
         include("examples/neural.jl")
     end