LowLevelFEM vs Ferrite benchmarking

ferrite_benchmark.jl

# Minimal standalone Ferrite assembly benchmark (see ../README.md for the full setup).
# Usage: julia --threads=N ferrite_benchmark.jl
# Requires: Ferrite, FerriteGmsh, OhMyThreads

using Ferrite, FerriteGmsh, OhMyThreads, SparseArrays, LinearAlgebra

const h = 0.025
const repeats = 5

function make_grid(h)
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.model.add("benchmark")
    volume = gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)
    gmsh.model.occ.synchronize()
    group = gmsh.model.addPhysicalGroup(3, [volume])
    gmsh.model.setPhysicalName(3, group, "body")
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)
    gmsh.model.mesh.generate(3)
    grid = FerriteGmsh.togrid()
    gmsh.finalize()
    return grid
end

function material_stiffness()
    E = 2.0e5; ν = 0.3
    λ = E * ν / ((1 + ν) * (1 - 2ν)); μ = E / (2(1 + ν))
    δ(i, j) = i == j ? 1.0 : 0.0
    return SymmetricTensor{4, 3}() do i, j, k, l
        λ * δ(i, j) * δ(k, l) + μ * (δ(i, k) * δ(j, l) + δ(i, l) * δ(j, k))
    end
end

function assemble_cell!(Ke, fe, cv, C, b)
    for q in 1:getnquadpoints(cv)
        dΩ = getdetJdV(cv, q)
        for i in 1:getnbasefunctions(cv)
            fe === nothing || (fe[i] += (shape_value(cv, q, i) ⋅ b) * dΩ)
            if Ke !== nothing
                ∇δui = shape_symmetric_gradient(cv, q, i)
                for j in 1:getnbasefunctions(cv)
                    Ke[i, j] += (∇δui ⊡ C ⊡ shape_symmetric_gradient(cv, q, j)) * dΩ
                end
            end
        end
    end
end

# Grid coloring approach from the Ferrite howto on multithreaded assembly:
# parallelize the cell loop within each color; task local scratch via @local.
function assemble_global!(K, f, dh, colors, cvt, C, b; ntasks = Threads.nthreads())
    K === nothing || (nonzeros(K) .= 0)
    f === nothing || fill!(f, 0)
    n = ndofs_per_cell(dh)
    for color in colors
        scheduler = OhMyThreads.DynamicScheduler(; ntasks)
        OhMyThreads.@tasks for cellidx in color
            @set scheduler = scheduler
            @local scratch = (;
                cc = CellCache(dh), cv = copy(cvt),
                Ke = K === nothing ? nothing : zeros(n, n),
                fe = f === nothing ? nothing : zeros(n),
                asm = K === nothing ? nothing :
                    start_assemble(K, f === nothing ? Float64[] : f; fillzero = false),
            )
            (; cc, cv, Ke, fe, asm) = scratch
            reinit!(cc, cellidx); reinit!(cv, cc)
            Ke === nothing || fill!(Ke, 0)
            fe === nothing || fill!(fe, 0)
            assemble_cell!(Ke, fe, cv, C, b)
            if asm !== nothing
                fe === nothing ? assemble!(asm, celldofs(cc), Ke) :
                    assemble!(asm, celldofs(cc), Ke, fe)
            else
                assemble!(f, celldofs(cc), fe)  # safe: colors don't share dofs
            end
        end
    end
end

function median_timed(fn, repeats)
    stats = map(1:repeats) do _
        GC.gc()
        @timed fn()
    end
    return sort(stats; by = s -> s.time)[cld(repeats, 2)]
end

function main()
    BLAS.set_num_threads(1)
    grid = make_grid(h)
    ip = Lagrange{RefTetrahedron, 1}()^3
    cv = CellValues(QuadratureRule{RefTetrahedron}(3), ip)
    dh = close!(add!(DofHandler(grid), :u, ip))
    coloring = @timed create_coloring(grid)
    pattern = @timed allocate_matrix(dh)
    colors, K = coloring.value, pattern.value
    f = zeros(ndofs(dh))
    C, b = material_stiffness(), Vec{3}((1.0, 2.0, 3.0))

    println("Threads: ", Threads.nthreads(), ", elements: ", getncells(grid),
        ", dofs: ", ndofs(dh))
    println("Coloring:           ", round(coloring.time; digits = 3), " s")
    println("Pattern allocation: ", round(pattern.time; digits = 3), " s")
    for (label, fn) in [
            "Matrix assembly" => () -> assemble_global!(K, nothing, dh, colors, cv, C, b),
            "Vector assembly" => () -> assemble_global!(nothing, f, dh, colors, cv, C, b),
            "Matrix + vector" => () -> assemble_global!(K, f, dh, colors, cv, C, b),
        ]
        fn() # warm-up (exclude compilation)
        s = median_timed(fn, repeats)
        println(label, ":    ", round(s.time; digits = 4), " s, ",
            round(s.bytes / 2^20; digits = 1), " MiB")
    end
    println("Checksums: ", sum(K), ", ", sum(f))
end

main()

llfem_benchmark.jl

# Minimal standalone LowLevelFEM assembly benchmark, essentially the worker part of
# the upstream benchmark_multithread.jl (see ../README.md for the full setup).
# Usage: julia --threads=N llfem_benchmark.jl
# Requires: LowLevelFEM

using LowLevelFEM, SparseArrays, LinearAlgebra

const h = 0.025
const repeats = 5

function make_mesh(h)
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0)
    gmsh.model.add("benchmark")
    volume = gmsh.model.occ.addBox(0, 0, 0, 1, 1, 1)
    gmsh.model.occ.synchronize()
    group = gmsh.model.addPhysicalGroup(3, [volume])
    gmsh.model.setPhysicalName(3, group, "body")
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", h)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", h)
    gmsh.model.mesh.generate(3)
    return length(gmsh.model.mesh.getNodes()[1]),
        sum(length, gmsh.model.mesh.getElements(3)[2])
end

function median_timed(fn, repeats)
    stats = map(1:repeats) do _
        GC.gc()
        @timed fn()
    end
    return sort(stats; by = s -> s.time)[cld(repeats, 2)]
end

function main()
    BLAS.set_num_threads(1)
    try
        nnodes, nelements = make_mesh(h)
        mat = Material("body") # E = 2.0e5, ν = 0.3 defaults
        problem = Problem([mat]; type = :VectorField, dim = 3, field = :u)
        μ, λ = mat.μ, mat.λ
        D = [
            λ + 2μ  λ       λ       0.0  0.0  0.0
            λ       λ + 2μ  λ       0.0  0.0  0.0
            λ       λ       λ + 2μ  0.0  0.0  0.0
            0.0     0.0     0.0     μ    0.0  0.0
            0.0     0.0     0.0     0.0  μ    0.0
            0.0     0.0     0.0     0.0  0.0  μ
        ]
        println("Threads: ", Threads.nthreads(), ", elements: ", nelements,
            ", nodes: ", nnodes)
        assemble_K = () -> ∫(SymGrad(problem) ⋅ D ⋅ SymGrad(problem); multithread = true)
        assemble_f = () -> ∫(problem ⋅ [1.0, 2.0, 3.0]; multithread = true)
        K = assemble_K(); f = assemble_f() # warm-up (exclude compilation)
        for (label, fn) in ["Matrix assembly" => assemble_K, "Vector assembly" => assemble_f]
            s = median_timed(fn, repeats)
            println(label, ": ", round(s.time; digits = 4), " s, ",
                round(s.bytes / 2^20; digits = 1), " MiB")
        end
        println("Checksums: ", sum(K.A), ", ", sum(f.a))
    finally
        gmsh.isInitialized() != 0 && gmsh.finalize()
    end
end

main()

Thank you, this is very helpful. I see that the Ferrite benchmark reuses the allocated sparsity pattern, while the current LowLevelFEM benchmark still goes through the high-level assembly path and rebuilds more of the matrix structure. This is exactly the optimization I am working on next: build the CSC pattern once and update only the numerical values for repeated assemblies.

Cool, I just didn't know how to do that in LowLevelFEM (but it is the default mode in Ferrite so can't really compare the other approach there).

Until recently, I used Julia’s built-in sparse matrix constructor to assemble the matrix from triplets, using sparse(I, J, V, n, m). In LowLevelFEM v1.13.29, I switched to a CSC-based assembly approach.

The new CSC-based implementation approximately doubles the matrix-assembly speed and substantially reduces memory usage. The sparsity pattern can also optionally be constructed once and reused for subsequent assemblies. With four threads, this reusable-pattern version achieves an approximately 2.8× speedup compared with the previous approach.

I am not currently planning to develop coloring or similar low-level assembly techniques. The present performance is already sufficient for my purposes, and my main focus is different from Ferrite’s: I would like to further strengthen the DSL and the multiphysics capabilities of LowLevelFEM.

My main goal is that mathematical expressions from a paper can be translated into executable Julia code with as little syntactic modification as possible, while keeping the code readable for engineers who are not professional programmers.

If you download the benchmark_multithread.jl file from the link provided earlier and update LowLevelFEM to v1.13.29, it will run using the new CSC-based assembly. The benchmark_multithread_2.jl script can be used to test the optional reuse of the CSC sparsity pattern as well.