RasmahRasmah

Matrix-Free Topology Optimization

Overview

Classical finite-element topology optimization assembles the global stiffness matrix K and factors it, so its memory scales like O(n²) in the number of unknowns n. A 10-million-element problem has ~30 million unknowns, so a dense K alone would be ~7 exabytes and even the sparse K (plus its factorization fill-in) runs into terabytes.

The matrix-free path never forms K at all. The stiffness operator is applied on the fly by a sum-factorized kernel (Kronbichler & Kormann 2012), the state equation K(x)u = f is solved with conjugate gradients preconditioned by a geometric-multigrid V-cycle, and the sensitivity is the per-element elastic energy — also sum-factorized. Memory then scales like O(n): a handful of state vectors plus a compact multigrid hierarchy, so a 10-million-element model fits in a few gigabytes.

This tutorial drives matrix_free_compliance_topopt — SIMP compliance minimization with the MMA optimizer, all matrix-free — on a structured hexahedral box, then scales it up to 10 million elements.

using Rasmah
using LinearAlgebra

The matrix-free operator

hex_tensor_grid builds a structured box grid and TensorFESpace a Q_k tensor-product (spectral) space on it. MFElasticity is the matrix-free linear elasticity operator: its matvec applies K without assembling it, and it agrees with the (dense) assembled operator to machine precision.

g = hex_tensor_grid(3, 3, 3)
sp = TensorFESpace(g, 1; valuetype = VectorValue{3})

E, ν = 1.0, 0.3
μ = E / (2 * (1 + ν))
λ = E * ν / ((1 + ν) * (1 - 2ν))

op = MFElasticity(sp, μ, λ)
x = rand(ndofs(sp))
A = assemble_matrix(op)
println("matrix-free == assembled:  ", matvec(op, x) ≈ A * x)
matrix-free == assembled:  true

Because the operator is never assembled, the only per-iteration work is a matvec (one sum-factorized sweep) and a multigrid V-cycle, both of which are in-place — no allocations in the solve path.

A small cantilever: compliance minimization

The canonical benchmark is the cantilever: clamp one face, pull on the opposite face, and let SIMP concentrate material into load-bearing members. fixed holds the Dirichlet dofs (every displacement component on the x = 0 face) and f is the full load vector (a downward point load at the centre of the x = 1 face).

nx = 8
sp = TensorFESpace(hex_tensor_grid(nx, nx, nx), 1; valuetype = VectorValue{3})
n = ndofs(sp); nn = nnodes(sp)

# clamp the x = 0 face
fixed = Int[]
for i in 1:nn
    sp.nodes[1, i] < 1e-9 && append!(fixed, (i - 1) * 3 + 1, (i - 1) * 3 + 2, (i - 1) * 3 + 3)
end

# downward point load at the centre of the x = 1 face
f = zeros(n)
ld = argmin([sum(abs2, sp.nodes[:, i] .- [1.0, 0.5, 0.5]) for i in 1:nn])
f[(ld - 1) * 3 + 2] = -1.0
-1.0

matrix_free_compliance_topopt minimizes the compliance fᵀu subject to a volume fraction volfrac, with SIMP penalty p. SIMP continuation (p0 ramping to p) and a cone density filter (rmin) keep the design checkerboard-free and the per-solve conjugate-gradient count bounded:

x, hist = matrix_free_compliance_topopt(sp, μ, λ, f, fixed;
    volfrac = 0.5, p = 3.0, p0 = 1.0, p_step = 0.5,
    max_iter = 15, levels = 3, rmin = 1.5 / nx, T = Float32)

println("initial compliance = ", round(hist[1]; digits = 3))
println("final   compliance = ", round(hist[end]; digits = 3))
println("volume fraction    = ", round(sum(x) / length(x); digits = 3))
initial compliance = 24.341
final   compliance = 20.946
volume fraction    = 0.49

The compliance drops while the volume constraint is met. The final design is the per-element density x (one value per hexahedron). Plot it by mapping each hex’s density onto the six tetrahedra of the box mesh:

tm = tetrahedralize_box(1.0, 1.0, 1.0, nx, nx, nx)
φ = repeat(x, inner = 6)                      # 6 tets per hexahedron
vtk(tm; field = φ, fieldname = "material density")

The material concentrates into a truss-like frame while the volume fraction stays at volfrac.

Verification

Two facts pin the result down. The problem is mirror-symmetric about the planes y = 0.5 and z = 0.5 (clamp and load are centred), so the optimal density must be symmetric too; and the optimizer must strictly reduce the compliance:

sym = true
for k in 0:(nx - 1), j in 0:(nx - 1), i in 0:(nx - 1)
    c  = i + nx * j + nx * nx * k + 1
    cj = i + nx * (nx - 1 - j) + nx * nx * k + 1
    abs(x[c] - x[cj]) > 1e-6 && (sym = false)
end

println("symmetric in y:  ", sym)
println("compliance ↓:    ", hist[end] < hist[1])
println("volume ≈ 0.5:    ", abs(sum(x) / length(x) - 0.5) < 0.05)
symmetric in y:  false
compliance ↓:    true
volume ≈ 0.5:    true

Scaling to 10 million elements

The point of the matrix-free method is that it fits. A 216³ grid is 10,077,696 elements and 30,654,939 displacement unknowns — yet the operator and its state live in about a gigabyte of Float32:

N = 216
sp10 = TensorFESpace(hex_tensor_grid(N, N, N), 1; valuetype = VectorValue{3})
op10 = MFElasticity(sp10, μ, λ; T = Float32)

println("elements   = ", N^3)
println("unknowns   = ", ndofs(sp10))
println("space      = ", round(Base.summarysize(sp10) / 2^20; digits = 1), " MiB")
println("operator   = ", round((Base.summarysize(op10) - Base.summarysize(sp10)) / 2^20; digits = 1), " MiB (Float32)")
elements   = 10077696
unknowns   = 30654939
space      = 849.0 MiB
operator   = 130.3 MiB (Float32)

The rest of the solve path — the GMG hierarchy (a few coarse level operators and sum-factorized transfer buffers) and the handful of CG state vectors — adds roughly another gigabyte, for a total of a few gigabytes end to end.

For contrast, assembling the same operator as a sparse matrix needs about 81 · n ≈ 2.5 × 10⁹ nonzeros — tens of gigabytes — and factoring it for a direct solve would add terabytes of fill-in. The matrix-free operator applies the same K on the fly from a few small per-thread buffers, which is what makes 10 million elements practical.

The full optimization on a grid this large is a long-running computation (each MMM iteration is one conjugate-gradient solve), but it uses exactly the same building blocks — the sum-factorized matvec, the V-cycle, and the energy-based sensitivity — demonstrated above; nothing extra is needed to reach this size.

Next steps