[ANN] LowLevelFEM.jl v1.13.28: help test multithreaded FEM assembly

I have recently added multithreaded assembly of bilinear and linear forms to the new stable version of LowLevelFEM.jl.

The implementation appears to work correctly, but at the moment I only have access to an older laptop with a 2-core / 4-thread Intel i5-7300U. On this machine, matrix assembly scales well from one to two threads, but—as expected—Hyper-Threading provides only a small additional improvement.

I would therefore be grateful if anyone with a more recent or higher-core-count CPU could run a small benchmark and share the output here.

The benchmark:

  • generates a fixed 3D finite element mesh;
  • measures bilinear-form (matrix) assembly;
  • measures linear-form (vector) assembly;
  • runs separate Julia processes with 1, 2, 4, … threads;
  • reports the median runtime and allocations;
  • prints checksums to verify that the results remain consistent.

How to run the benchmark

  1. Install or update LowLevelFEM.jl from the Julia REPL:
using Pkg
Pkg.add("LowLevelFEM")
Pkg.update("LowLevelFEM")
  1. Download the benchmark script:

Download benchmark_multithread.jl

  1. Open a terminal in the folder containing the downloaded file and run:
julia benchmark_multithread.jl
  1. Copy the complete terminal output and paste it into your reply.

No Git installation or repository checkout is required.

The first run may take some time because Julia needs to install and precompile the dependencies. Compilation is excluded from the reported assembly times.

Please copy the complete terminal output into your reply. The script automatically reports the relevant system information, including the CPU model, number of logical threads, Julia version, memory, mesh size, and LowLevelFEM version.

On my machine, the matrix-assembly results for approximately 290,000 tetrahedral elements were:

Threads: 1    Matrix assembly: 6.989 s
Threads: 2    Matrix assembly: 4.606 s
Threads: 4    Matrix assembly: 4.541 s

This corresponds to speedups of approximately (1.52\times) with two threads and (1.54\times) with four threads on a two-core CPU.

I am particularly interested in results from CPUs with 4, 8, 16, or more physical cores. Reports from Windows and macOS would also be very useful.

Thank you for helping test it!

ERROR: The following 1 direct dependency failed to precompile:

LowLevelFEM

Failed to precompile LowLevelFEM [6171b9fb-adbf-4751-adb9-5faded75de07] to "/home/ofekb/.julia/compiled/v1.12/LowLevelFEM/jl_LmRsSR".
ERROR: LoadError: InitError: could not load library "/home/ofekb/.julia/artifacts/fe56f62d30a5f4c8effc9aac1f65c55e839242e1/lib/libarpack.so"
libquadmath.so.0: cannot open shared object file: No such file or directory
Stacktrace:
  [1] #dlopen#3
    @ ./libdl.jl:120 [inlined]
  [2] dlopen(s::String, flags::UInt32)
    @ Base.Libc.Libdl ./libdl.jl:119
  [3] macro expansion
    @ ~/.julia/packages/JLLWrappers/Kp7TC/src/products/library_generators.jl:63 [inlined]
  [4] __init__()
    @ Arpack_jll ~/.julia/packages/Arpack_jll/g8D7r/src/wrappers/x86_64-linux-gnu-libgfortran5.jl:10
  [5] run_module_init(mod::Module, i::Int64)
    @ Base ./loading.jl:1445
  [6] register_restored_modules(sv::Core.SimpleVector, pkg::Base.PkgId, path::String)
    @ Base ./loading.jl:1433
  [7] _include_from_serialized(pkg::Base.PkgId, path::String, ocachepath::String, depmods::Vector{Any}; register::Bool)
    @ Base ./loading.jl:1321
  [8] _include_from_serialized
    @ ./loading.jl:1276 [inlined]
  [9] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String, build_id::UInt128, stalecheck::Bool; reasons::Dict{String, Int64}, DEPOT_PATH::Vector{String})
    @ Base ./loading.jl:2117
 [10] _require_search_from_serialized
    @ ./loading.jl:2011 [inlined]
 [11] __require_prelocked(pkg::Base.PkgId, env::String)
    @ Base ./loading.jl:2629
 [12] _require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:2495
 [13] macro expansion
    @ ./loading.jl:2423 [inlined]
 [14] macro expansion
    @ ./lock.jl:376 [inlined]
 [15] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2388
 [16] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2364
 [17] top-level scope
    @ ~/.julia/packages/Arpack/FCvNd/src/Arpack.jl:8
 [18] include(mod::Module, _path::String)
    @ Base ./Base.jl:306
 [19] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::String)
    @ Base ./loading.jl:3028
 [20] top-level scope
    @ stdin:5
 [21] eval(m::Module, e::Any)
    @ Core ./boot.jl:489
 [22] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:2874
 [23] include_string
    @ ./loading.jl:2884 [inlined]
 [24] exec_options(opts::Base.JLOptions)
    @ Base ./client.jl:315
 [25] _start()
    @ Base ./client.jl:550
during initialization of module Arpack_jll
in expression starting at /home/ofekb/.julia/packages/Arpack/FCvNd/src/Arpack.jl:2
in expression starting at stdin:5
ERROR: LoadError: Failed to precompile Arpack [7d9fca2a-8960-54d3-9f78-7d1dccf2cb97] to "/home/ofekb/.julia/compiled/v1.12/Arpack/jl_ogoD6X".
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:44
  [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool; flags::Cmd, cacheflags::Base.CacheFlags, reasons::Dict{String, Int64}, loadable_exts::Nothing)
    @ Base ./loading.jl:3315
  [3] (::Base.var"#__require_prelocked##0#__require_prelocked##1"{Base.PkgId, String, Dict{String, Int64}})()
    @ Base ./loading.jl:2683
  [4] mkpidlock(f::Base.var"#__require_prelocked##0#__require_prelocked##1"{Base.PkgId, String, Dict{String, Int64}}, at::String, pid::Int32; kwopts::@Kwargs{stale_age::Int64, wait::Bool})
    @ FileWatching.Pidfile /nix/store/i1srz1q8m3cw26qc2z9hlizc53p7nba6-julia-1.12.6/share/julia/stdlib/v1.12/FileWatching/src/pidfile.jl:96
  [5] #mkpidlock#7
    @ /nix/store/i1srz1q8m3cw26qc2z9hlizc53p7nba6-julia-1.12.6/share/julia/stdlib/v1.12/FileWatching/src/pidfile.jl:91 [inlined]
  [6] trymkpidlock(::Function, ::Vararg{Any}; kwargs::@Kwargs{stale_age::Int64})
    @ FileWatching.Pidfile /nix/store/i1srz1q8m3cw26qc2z9hlizc53p7nba6-julia-1.12.6/share/julia/stdlib/v1.12/FileWatching/src/pidfile.jl:117
  [7] #invokelatest_gr#232
    @ ./reflection.jl:1297 [inlined]
  [8] invokelatest_gr
    @ ./reflection.jl:1289 [inlined]
  [9] maybe_cachefile_lock(f::Base.var"#__require_prelocked##0#__require_prelocked##1"{Base.PkgId, String, Dict{String, Int64}}, pkg::Base.PkgId, srcpath::String; stale_age::Int64)
    @ Base ./loading.jl:3886
 [10] maybe_cachefile_lock
    @ ./loading.jl:3883 [inlined]
 [11] __require_prelocked(pkg::Base.PkgId, env::String)
    @ Base ./loading.jl:2669
 [12] _require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:2495
 [13] macro expansion
    @ ./loading.jl:2423 [inlined]
 [14] macro expansion
    @ ./lock.jl:376 [inlined]
 [15] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2388
 [16] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2364
 [17] top-level scope
    @ ~/.julia/packages/LowLevelFEM/1cgQE/src/LowLevelFEM.jl:6
 [18] include(mod::Module, _path::String)
    @ Base ./Base.jl:306
 [19] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::Nothing)
    @ Base ./loading.jl:3028
 [20] top-level scope
    @ stdin:5
 [21] eval(m::Module, e::Any)
    @ Core ./boot.jl:489
 [22] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:2874
 [23] include_string
    @ ./loading.jl:2884 [inlined]
 [24] exec_options(opts::Base.JLOptions)
    @ Base ./client.jl:315
 [25] _start()
    @ Base ./client.jl:550

i have tried to run the benchmark but it failed at precompile. is it because of NixOS or is it an actual bug?

Since I dual boot, here is for both Linux and Windows:

Linux:

========================================================================
Threads:              1
Julia:                1.12.6
OS / architecture:    Linux / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         14.9 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    /home/bjarke/.julia/packages/LowLevelFEM/1cgQE/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51919
Elements:             289960
Matrix assembly:      3.706862 s, 10279.1 MiB
Vector assembly:      0.963317 s, 672.3 MiB
Matrix size:          (155757, 155757)
Matrix nonzeros:      6820659
Matrix checksum:      5.9626472648233175e-9
Vector checksum:      6.000000000000011

========================================================================
Threads:              2
Julia:                1.12.6
OS / architecture:    Linux / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         14.9 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    /home/bjarke/.julia/packages/LowLevelFEM/1cgQE/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51919
Elements:             289960
Matrix assembly:      2.279432 s, 10270.3 MiB
Vector assembly:      0.636823 s, 672.3 MiB
Matrix size:          (155757, 155757)
Matrix nonzeros:      6820659
Matrix checksum:      5.9626472648233175e-9
Vector checksum:      6.000000000000011

========================================================================
Threads:              4
Julia:                1.12.6
OS / architecture:    Linux / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         14.9 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    /home/bjarke/.julia/packages/LowLevelFEM/1cgQE/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51919
Elements:             289960
Matrix assembly:      1.692643 s, 10270.3 MiB
Vector assembly:      0.540016 s, 672.4 MiB
Matrix size:          (155757, 155757)
Matrix nonzeros:      6820659
Matrix checksum:      5.9626472648233175e-9
Vector checksum:      6.000000000000011

========================================================================
Threads:              8
Julia:                1.12.6
OS / architecture:    Linux / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         14.9 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    /home/bjarke/.julia/packages/LowLevelFEM/1cgQE/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51919
Elements:             289960
Matrix assembly:      1.414195 s, 10270.4 MiB
Vector assembly:      0.483964 s, 672.4 MiB
Matrix size:          (155757, 155757)
Matrix nonzeros:      6820659
Matrix checksum:      5.9626472648233175e-9
Vector checksum:      6.000000000000011

========================================================================
Threads:              16
Julia:                1.12.6
OS / architecture:    Linux / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         14.9 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    /home/bjarke/.julia/packages/LowLevelFEM/1cgQE/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51919
Elements:             289960
Matrix assembly:      1.352271 s, 10270.5 MiB
Vector assembly:      0.49996 s, 672.4 MiB
Matrix size:          (155757, 155757)
Matrix nonzeros:      6820659
Matrix checksum:      5.9626472648233175e-9
Vector checksum:      6.000000000000011

Windows:

========================================================================
Threads:              1
Julia:                1.12.6
OS / architecture:    NT / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         15.3 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    C:\Users\bjarke\.julia\packages\LowLevelFEM\1cgQE\src\LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51907
Elements:             289940
Matrix assembly:      3.405405 s, 10278.4 MiB
Vector assembly:      0.941553 s, 672.3 MiB
Matrix size:          (155721, 155721)
Matrix nonzeros:      6820007
Matrix checksum:      -8.731149137020111e-10
Vector checksum:      6.000000000000011

========================================================================
Threads:              2
Julia:                1.12.6
OS / architecture:    NT / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         15.3 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    C:\Users\bjarke\.julia\packages\LowLevelFEM\1cgQE\src\LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51907
Elements:             289940
Matrix assembly:      2.333201 s, 10269.6 MiB
Vector assembly:      0.560898 s, 672.3 MiB
Matrix size:          (155721, 155721)
Matrix nonzeros:      6820007
Matrix checksum:      -8.731149137020111e-10
Vector checksum:      6.000000000000011

========================================================================
Threads:              4
Julia:                1.12.6
OS / architecture:    NT / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         15.3 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    C:\Users\bjarke\.julia\packages\LowLevelFEM\1cgQE\src\LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51907
Elements:             289940
Matrix assembly:      1.674898 s, 10269.6 MiB
Vector assembly:      0.494225 s, 672.3 MiB
Matrix size:          (155721, 155721)
Matrix nonzeros:      6820007
Matrix checksum:      -8.731149137020111e-10
Vector checksum:      6.000000000000011

========================================================================
Threads:              8
Julia:                1.12.6
OS / architecture:    NT / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         15.3 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    C:\Users\bjarke\.julia\packages\LowLevelFEM\1cgQE\src\LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51907
Elements:             289940
Matrix assembly:      1.46257 s, 10269.7 MiB
Vector assembly:      0.426041 s, 672.3 MiB
Matrix size:          (155721, 155721)
Matrix nonzeros:      6820007
Matrix checksum:      -8.731149137020111e-10
Vector checksum:      6.000000000000011

========================================================================
Threads:              16
Julia:                1.12.6
OS / architecture:    NT / x86_64
CPU:                  AMD Ryzen 7 8845HS w/ Radeon 780M Graphics
Logical CPU threads:  16
Total memory:         15.3 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    C:\Users\bjarke\.julia\packages\LowLevelFEM\1cgQE\src\LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51907
Elements:             289940
Matrix assembly:      1.317764 s, 10269.8 MiB
Vector assembly:      0.41528 s, 672.3 MiB
Matrix size:          (155721, 155721)
Matrix nonzeros:      6820007
Matrix checksum:      -8.731149137020111e-10
Vector checksum:      6.000000000000011

Thanks for trying the benchmark and reporting this.

This is not a bug in the benchmark or in the new multithreaded assembly. The failure occurs earlier, while LowLevelFEM is being loaded: Arpack_jll cannot find libquadmath.so.0.

Since your Julia installation comes from the Nix store, this appears to be a NixOS library-visibility/packaging issue. LowLevelFEM currently imports Arpack during startup, so an Arpack loading failure prevents the whole package from loading even though this benchmark does not use the eigensolver.

The simplest workaround would probably be to run the benchmark using Julia’s official generic Linux binary. Alternatively, libquadmath needs to be made available to Julia in the Nix environment.

I will also consider whether Arpack can be loaded only when the eigensolver is actually used, so that an optional feature cannot prevent the rest of LowLevelFEM from starting.

Thank you for testing it on both Linux and Windows! The results are very useful, and it is good to see that the multithreaded assembly behaves consistently on both systems.

On macOS

========================================================================
Threads:              1
Julia:                1.12.6
OS / architecture:    Darwin / aarch64
CPU:                  Apple M4
Logical CPU threads:  4
Total memory:         16.0 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    /Users/mat/.julia/packages/LowLevelFEM/1cgQE/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51809
Elements:             289338
Matrix assembly:      1.965245 s, 10304.6 MiB
Vector assembly:      0.49342 s, 677.7 MiB
Matrix size:          (155427, 155427)
Matrix nonzeros:      6806547
Matrix checksum:      4.849425749853253e-9
Vector checksum:      6.000000000000011

========================================================================
Threads:              2
Julia:                1.12.6
OS / architecture:    Darwin / aarch64
CPU:                  Apple M4
Logical CPU threads:  4
Total memory:         16.0 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    /Users/mat/.julia/packages/LowLevelFEM/1cgQE/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51809
Elements:             289338
Matrix assembly:      1.246643 s, 10297.6 MiB
Vector assembly:      0.336808 s, 718.7 MiB
Matrix size:          (155427, 155427)
Matrix nonzeros:      6806547
Matrix checksum:      4.849425749853253e-9
Vector checksum:      6.000000000000011

========================================================================
Threads:              4
Julia:                1.12.6
OS / architecture:    Darwin / aarch64
CPU:                  Apple M4
Logical CPU threads:  4
Total memory:         16.0 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.28
Package tree hash:    /Users/mat/.julia/packages/LowLevelFEM/1cgQE/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51809
Elements:             289338
Matrix assembly:      0.910792 s, 10297.7 MiB
Vector assembly:      0.274606 s, 718.7 MiB
Matrix size:          (155427, 155427)
Matrix nonzeros:      6806547
Matrix checksum:      4.849425749853253e-9
Vector checksum:      6.000000000000011

Thank you for testing it on Apple Silicon! A 2.16× speedup for matrix assembly on four threads is a very useful result, and it is good to see that the implementation works consistently on macOS as well.

Ran the benchmark on a bigger machine, results below.

CPU:                  2 x AMD EPYC 9354 32-Core Processor (64 physical cores)
Total memory:         755.6 GiB
OS / architecture:    Linux / x86_64
Julia:                1.12.6
LowLevelFEM:          v1.13.28
BLAS threads:         1

The measurement loop and problem setup are identical to benchmark_multithread.jl, run with 1–32 threads. I also ran a finer mesh (h = 0.0125) in addition to the default.

h = 0.025 (51,919 nodes / 289,960 elements / 155,757 dofs):

Threads Matrix assembly Speedup Vector assembly Speedup
1 4.53 s 1.0× 1.11 s 1.0×
2 2.69 s 1.7× 0.72 s 1.5×
4 2.40 s 1.9× 0.68 s 1.6×
8 1.75 s 2.6× 0.58 s 1.9×
16 1.56 s 2.9× 0.54 s 2.1×
32 1.35 s 3.4× 0.53 s 2.1×

h = 0.0125 (385,009 nodes / 2,276,675 elements / 1,155,027 dofs):

Threads Matrix assembly Speedup Vector assembly Speedup
1 37.95 s 1.0× 9.62 s 1.0×
2 25.08 s 1.5× 7.35 s 1.3×
4 20.00 s 1.9× 5.84 s 1.6×
8 15.47 s 2.5× 5.20 s 1.8×
16 15.29 s 2.5× 5.33 s 1.8×
32 13.34 s 2.8× 5.41 s 1.8×
For my own curiousity I compared with Ferrite.jl

I also ran the equivalent assembly with Ferrite.jl on the same machine and the exact same mesh (same gmsh-generated .msh file), material parameters, and quadrature.

LowLevelFEM builds a new sparse matrix on every call, while Ferrite assembles into a pre-allocated sparsity pattern. The numbers in the table is just the re-assembly and the one-time setup costs for Ferrite are included below the tables.

Matrix assembly, same mesh, h = 0.025 (155,757 dofs):

Threads LowLevelFEM Ferrite Ferrite speedup (vs 1 thread)
1 4.53 s 2.03 s 1.0×
2 2.69 s 1.17 s 1.7×
4 2.40 s 0.59 s 3.4×
8 1.75 s 0.40 s 5.1×
16 1.56 s 0.21 s 9.5×
32 1.35 s 0.13 s 16.1×

Matrix assembly, same mesh, h = 0.0125 (1,155,027 dofs):

Threads LowLevelFEM Ferrite Ferrite speedup (vs 1 thread)
1 37.95 s 17.12 s 1.0×
2 25.08 s 8.33 s 2.1×
4 20.00 s 4.18 s 4.1×
8 15.47 s 2.40 s 7.1×
16 15.29 s 1.47 s 11.6×
32 13.34 s 0.85 s 20.2×

Assembling the matrix and the vector in the same loop (the idiomatic way in Ferrite) costs essentially the same as the matrix assemble alone so the vector comes along for free when both are needed.

One-time setup costs for Ferrite at h = 0.0125:

  • Sparse matrix allocation: ~4–5s (serial process so doesn’t change with number of threads)
  • Grid coloring: ~16s serially and 5.2s at 32 threads

For a single assembly, serial Ferrite (17 s assembly + 16 s coloring + 5 s pattern allocation = 38 s) compares with LowLevelFEM. For something like 10 timesteps or Newton iterations Ferrite’s speedup from 1 to 32 threads would roughly be 11× (192 s → 18 s, including the one-time setup) and LowLevelFEM speedup roughly 2.8× (380 s → 133 s). For 100 steps/iterations the setup costs amortize away and Ferrite reaches roughly 18× (1733 s → 94 s) while LowLevelFEM stays at 2.8× (3795 s → 1334 s). Maybe there is a way to pre-allocate in LowLevelFEM too though, I don’t know.

(Disclaimer: I used Claude to put together this summary from the benchmark data. Happy to share the benchmark scripts if you are curious.)

Thank you — this is an extremely useful comparison. Preallocating the CSC sparsity pattern and updating only nzval is indeed possible in LowLevelFEM, and it is actually the next assembly optimization I am working on; your results clearly show how important this will be for nonlinear and transient problems.

Ferrite’s multithreaded scaling is very impressive. I would also be very interested in the benchmark scripts, especially to ensure that I reproduce the Ferrite setup and comparison correctly.

See LowLevelFEM vs Ferrite benchmarking · GitHub

The latest LowLevelFEM version no longer relies on Julia’s built-in sparse matrix constructor from triplets. It uses a CSC-based assembly approach instead, which has nearly doubled the assembly speed, depending on the problem size.

The CSC pattern can also be reused for subsequent matrix assemblies, so it does not have to be recomputed as long as the mesh remains unchanged.

This functionality is included in LowLevelFEM v1.13.29. You can try it by downloading the updated benchmark script again from the link:

(second script reuses the CSC pattern)

Here is the output for the second script on macOS

========================================================================
Threads:              1
Julia:                1.12.6
OS / architecture:    Darwin / aarch64
CPU:                  Apple M4
Logical CPU threads:  4
Total memory:         16.0 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.29
Package tree hash:    /Users/mat/.julia/packages/LowLevelFEM/Ua7Gl/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51809
Elements:             289338
Matrix assembly:      1.019934 s, 12.7 MiB
Vector assembly:      0.501133 s, 682.6 MiB
Matrix size:          (155427, 155427)
Matrix nonzeros:      6807141
Matrix checksum:      1.209351049965699e10
Vector checksum:      6.000000000000011

========================================================================
Threads:              2
Julia:                1.12.6
OS / architecture:    Darwin / aarch64
CPU:                  Apple M4
Logical CPU threads:  4
Total memory:         16.0 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.29
Package tree hash:    /Users/mat/.julia/packages/LowLevelFEM/Ua7Gl/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51809
Elements:             289338
Matrix assembly:      0.622872 s, 64.6 MiB
Vector assembly:      0.338922 s, 708.1 MiB
Matrix size:          (155427, 155427)
Matrix nonzeros:      6807141
Matrix checksum:      1.209351049965699e10
Vector checksum:      6.000000000000011

========================================================================
Threads:              4
Julia:                1.12.6
OS / architecture:    Darwin / aarch64
CPU:                  Apple M4
Logical CPU threads:  4
Total memory:         16.0 GiB
BLAS threads:         1
LowLevelFEM version:  1.13.29
Package tree hash:    /Users/mat/.julia/packages/LowLevelFEM/Ua7Gl/src/LowLevelFEM.jl
Mesh size parameter:  0.025
Repetitions:          5
Nodes:                51809
Elements:             289338
Matrix assembly:      0.458384 s, 263.3 MiB
Vector assembly:      0.266647 s, 762.9 MiB
Matrix size:          (155427, 155427)
Matrix nonzeros:      6807141
Matrix checksum:      1.209351049965699e10
Vector checksum:      6.000000000000011