Zero pivots in factorization error when computing compliance matrix

[tylercritchfield]

Here's a basic example that leads to a factorization error with zero pivots.

using GXBeam

xaf = [1.0, 0.99669031, 0.98707731, 0.97177624, 0.95120013, 0.92554805, 0.89513234, 0.86042844, 0.82194932, 0.78024605, 0.73589923, 0.68948289, 0.64155561, 0.59265387, 0.5432852, 0.49390561, 0.44491485, 0.39667271, 0.34949257, 0.30376269, 0.25995386, 0.21850111, 0.17980252, 0.14421934, 0.11207729, 0.08364666, 0.05914156, 0.03873701, 0.02260069, 0.01080931, 0.00329342, 0.0, 0.00480487, 0.0111091, 0.02423919, 0.04202841, 0.0645839, 0.09186189, 0.1236919, 0.15984299, 0.20001108, 0.24382591, 0.29085133, 0.34057636, 0.39242047, 0.4457723, 0.49999633, 0.55444002, 0.60844101, 0.66133467, 0.71246137, 0.7611783, 0.80686551, 0.84892722, 0.8867981, 0.91995456, 0.9479178, 0.97026156, 0.98661057, 0.99662097, 1.0]
yaf = [0.0, 0.00049003, 0.00207033, 0.00475426, 0.00825965, 0.01241627, 0.01723294, 0.0226708, 0.02862516, 0.034964, 0.0415088, 0.0480415, 0.05432933, 0.06013032, 0.06520439, 0.06931828, 0.07229087, 0.07399667, 0.0744182, 0.07365704, 0.07178323, 0.06883168, 0.064851, 0.05990381, 0.05406639, 0.04742811, 0.04011465, 0.03230144, 0.02418269, 0.01593784, 0.00791465, 0.0, -0.00634666, -0.0102234, -0.01556325, -0.02006888, -0.02353976, -0.02598343, -0.02742948, -0.02793392, -0.02758763, -0.026503, -0.02481832, -0.02268597, -0.02024403, -0.01760693, -0.01488444, -0.01217351, -0.00956485, -0.00713548, -0.0049519, -0.0030634, -0.00151205, -0.00032506, 0.00048717, 0.00093779, 0.00105862, 0.00091038, 0.00057426, 0.00018251, 0.0]

chord = 0.0863
twist = 17.7*pi/180
pitch_axis = 0.25 / chord

uni = Material(1.75e11, 8.0e9, 8.0e9, 5.0e9, 5.0e9, 5.0e9, 0.3, 0.3, 0.3, 1600.0)
weave = Material(8.5e10, 8.5e10, 8.5e10, 5.0e9, 5.0e9, 5.0e9, 0.1, 0.1, 0.1, 1600.0)
materials = [uni, weave]

thickness = [0.001, 0.01853, 0.001]
orientation = [45.0,0.0,45.0]*pi/180
material_ids = [2;1;2]
layup1 = Layer.(materials[material_ids], thickness, orientation)

segments = [layup1]

xbreak = [0.0, 1.0]
web_location = [] #no webs
webs = []

nodes, elements = afmesh(xaf, yaf, chord, twist, pitch_axis, xbreak, web_location, segments, webs)
cache = initialize_cache(nodes, elements, eltype(chord))

compliance, _, _ = compliance_matrix(nodes, elements; cache)
mass, _ = mass_matrix(nodes, elements)

Here is the stacktrace:

ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.
Stacktrace:
 [1] ldlt!(F::SuiteSparse.CHOLMOD.Factor{Float64}, A::SuiteSparse.CHOLMOD.Sparse{Float64}; shift::Float64, check::Bool)
   @ SuiteSparse.CHOLMOD /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1472
 [2] #ldlt!#11
   @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1497 [inlined]
 [3] ldlt!
   @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SuiteSparse/src/cholmod.jl:1497 [inlined]
 [4] factorize(A::LinearAlgebra.Symmetric{Float64, SparseArrays.SparseMatrixCSC{Float64, Int64}})
   @ SparseArrays /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/linalg.jl:1634
 [5] linearsolve2(AM::SparseArrays.SparseMatrixCSC{Float64, Int64}, B2::SparseArrays.SparseMatrixCSC{Float64, Int64})
   @ GXBeam ~/.julia/dev/GXBeam/src/section.jl:559
 [6] compliance_matrix(nodes::Vector{Node{Float64}}, elements::Vector{MeshElement{Vector{Int64}, Float64}}; cache::GXBeam.SectionCache{Matrix{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Int64}}, gxbeam_order::Bool)
   @ GXBeam ~/.julia/dev/GXBeam/src/section.jl:736
 [7] gxbeam_test2(i::Int64)
   @ Main ~/Dropbox/Research-Ning/projects/Prop_out/src/debugging/gxbeam_test.jl:106
 [8] top-level scope
   @ REPL[30]:1

I'm using GXBeam 0.4.1 and Julia 1.6.6.

I can try to figure out what's going on but from a brief look into it it seems others have gotten this error before (here and here) on older versions of Julia, where it was determined it was a bug in Julia itself that was fixed in subsequent releases. After seeing this I tried it again with Julia 1.8.0 with the same result.

@andrewning did you see this at all during the code development?

[andrewning]

If you see an error in a linear solve, that almost certainly means the problem is ill-posed (not a bug in Julia) - in the case of section analysis probably a messed up mesh. I would recommend always visualizing your mesh while creating: plotmesh(nodes, elements, PyPlot). If you do that in this case you can see it goes crazy. Looking a little closer the problem is your chord (0.0863, and so thickness is an order of magnitude less than that) is super tiny relative to the internal thicknesses (0.01853) so the structure cannot physically fit.

[tylercritchfield]

That makes sense, I hadn't yet given much thought to the thicknesses. Thanks for the feedback.

[tylercritchfield]

Reopening as I've found this error again, on a seemingly feasible geometry. The following function can run two different geometries, where the works keyword distinguishes between the two.

function gxbeam_test(;works=false)

    chord = 0.1
    twist = 0.4456941218448864
    pitch_axis = 0.25 / chord

    xaf = [1.0, 0.99669031, 0.98707731, 0.97177624, 0.95120013, 0.92554805, 0.89513234, 0.86042844, 0.82194932, 0.78024605, 0.73589923, 0.68948289, 0.64155561, 0.59265387, 0.5432852, 0.49390561, 0.44491485, 0.39667271, 0.34949257, 0.30376269, 0.25995386, 0.21850111, 0.17980252, 0.14421934, 0.11207729, 0.08364666, 0.00329342, 0.0, 0.0645839, 0.09186189, 0.1236919, 0.15984299, 0.20001108, 0.24382591, 0.29085133, 0.34057636, 0.39242047, 0.4457723, 0.49999633, 0.55444002, 0.60844101, 0.66133467, 0.71246137, 0.7611783, 0.80686551, 0.84892722, 0.8867981, 0.91995456, 0.9479178, 0.97026156, 0.98661057, 0.99662097, 1.0]
    yaf = [0.0, 0.00049003, 0.00207033, 0.00475426, 0.00825965, 0.01241627, 0.01723294, 0.0226708, 0.02862516, 0.034964, 0.0415088, 0.0480415, 0.05432933, 0.06013032, 0.06520439, 0.06931828, 0.07229087, 0.07399667, 0.0744182, 0.07365704, 0.07178323, 0.06883168, 0.064851, 0.05990381, 0.05406639, 0.04742811,  0.00791465,  0.0, -0.02353976, -0.02598343, -0.02742948, -0.02793392, -0.02758763, -0.026503, -0.02481832, -0.02268597, -0.02024403, -0.01760693, -0.01488444, -0.01217351, -0.00956485, -0.00713548, -0.0049519, -0.0030634, -0.00151205, -0.00032506,  0.00048717,  0.00093779,  0.00105862,  0.00091038,  0.00057426,  0.00018251,  0.0]

    uni = GXBeam.Material(1.75e11, 8.0e9, 8.0e9, 5.0e9, 5.0e9, 5.0e9, 0.3, 0.3, 0.3, 1600.0)
    weave = GXBeam.Material(8.5e10, 8.5e10, 8.5e10, 5.0e9, 5.0e9, 5.0e9, 0.1, 0.1, 0.1, 1600.0)
    materials = [uni, weave]

    if works # no exception
        thickness = [1e-4, 0.00275, 1e-4] # total thickness is 2.95% of the chord
    else # zero-pivot exception
        thickness = [9.99999999999994e-5, 0.002749999999999999, 9.99999999999994e-5]
    end

    orientation = [45.0, 0.0, 45.0]*pi/180
    material_ids = [2,1,2]

    layup = GXBeam.Layer.(materials[material_ids], thickness, orientation)
    segments = [layup]

    xbreak = [0.0, 1.0]
    web_location = []
    webs = []

    nodes, elements = GXBeam.afmesh(xaf, yaf, chord, twist, pitch_axis, xbreak, web_location, segments, webs)
    cache = GXBeam.initialize_cache(nodes, elements, eltype(chord))

    figure(); GXBeam.plotmesh(nodes, elements, PyPlot)

    compliance, _, _ = GXBeam.compliance_matrix(nodes, elements; cache)
    mass, _ = GXBeam.mass_matrix(nodes, elements)

    return compliance, mass
end

Here's a pasted image of the mesh (I've really coursened up my airfoil at the leading edge to ensure the mesh elements weren't overlapping even with some variation in the ply thickness. The example above has a laminate thickness of about 3% of the chord):

The only difference between the two examples is the ply thickness, which for all intents and purposes is the same (the plotted meshes will look identical even if you zoom in as far as you can with PyPlot), but are ever so slightly off enough to cause some havoc. At a closer look, 8 of the 116 nodes are not identical. Here's node 11 compared between the two:

julia> nodes_workstrue[11]
GXBeam.Node{Float64}(0.03618109685895313, 0.10505151034301377)

julia> nodes_worksfalse[11]
GXBeam.Node{Float64}(0.03618109685895313, 0.10505151034301378)

The very last decimal place of the y coordinate is off by one value. The other differing nodes are similar.

I can put a breakpoint in section.jl at line 559 and compare the symmetric-sparse AM matrix (before it's factorized, which is where the exception occurs), and there are some elements that are a bit different. See for instance:

julia> AM_1[5,5]
5.99783700658042e12
julia> AM_2[5,5]
5.997877516217298e12

But all of the elements remain at the same order of magnitude, same sign, etc. and honestly, nothing jumps out at me as being wrong so I'm not sure what to do from there...

[andrewning]

I can't reproduce. Both run fine for me. Again, an error like this almost certainly means you have an ill-conditioned matrix. You have some extremely skewed mesh elements, especially towards the trailing edge. I wouldn't be surprised to see numerical issues with that. Also, the outer plys are so thin relative to your center laminate (with properties that aren't so different) that might not be worth including. But at a minimum I would definitely increase the resolution at the leading and trailing edge (especially trailing edge). Way too coarse.

By the way laminate thickness that is 3% of chord isn't really the most relevant metric. the airfoil thickness is an order of magnitude smaller, and the thickness comes from both sides, so the total laminate thickness is like 60% of the airfoil thickness.

[tylercritchfield]

Thanks. I'll work on finding the right resolution balance and give the thickness problem some more thought.