Differentiable multiphysics and topology optimization

Overview

The weak-form layer is differentiable: assemble and solve with automatic differentiation, so CAD → mesh → weak form → solve → objective gradients flow end to end. This tutorial shows compliance gradients and density-based topology optimization.

using Rasmah, LinearAlgebra

E, ν = 1e3, 0.3
C = elasticity_matrix(E, ν)

6×6 FixedSizeArrays.FixedSizeArray{Float64, 2, Memory{Float64}}:
 1346.15    576.923   576.923    0.0      0.0      0.0
  576.923  1346.15    576.923    0.0      0.0      0.0
  576.923   576.923  1346.15     0.0      0.0      0.0
    0.0       0.0       0.0    384.615    0.0      0.0
    0.0       0.0       0.0      0.0    384.615    0.0
    0.0       0.0       0.0      0.0      0.0    384.615

Differentiable compliance

Differentiate compliance w.r.t. the box width through mesh generation, assembly, and the solve:

function wf_compliance(w)
    m = tetrahedralize_box(w, 1.0, 1.0, 2, 2, 2)
    n = size(m.nodes, 2)
    fixed = Int[]
    for i in 1:n
        m.nodes[1, i] < 1e-9 && append!(fixed, 3i - 2, 3i - 1, 3i)
    end
    sp = FESpace(m, ReferenceFE(3, 1); valuetype = VectorValue{3})
    a = BilinearForm(sp, sp, elasticity_integrand(C))
    loads = zeros(ndofs(sp))
    for i in 1:n
        m.nodes[1, i] > 0.9 * w && (loads[3i - 2] = 1.0 / 9)
    end
    K = assemble_matrix(a)
    free = setdiff(1:ndofs(sp), fixed)
    u = cholesky(K[free, free]) \ loads[free]
    return dot(loads[free], u)
end

g_ad = derivative(wf_compliance, 1.0)
g_fd = (wf_compliance(1.0 + 1e-6) - wf_compliance(1.0 - 1e-6)) / 2e-6
println("∇compliance: AD = ", round(g_ad; digits = 6), ", FD = ", round(g_fd; digits = 6))

∇compliance: AD = 0.000911, FD = 0.000911

The assembled operator solves through the native sparse factorization, so the whole chain — mesh generation, assembly, and solve — is differentiable; compliance_gradient is the vector form over several parameters.

Density-based topology optimization (SIMP)

ComplianceTopOpt minimizes compliance on a fixed mesh by scaling each element’s stiffness E → φᵖE; the density field itself comes from an SDF, so the design is differentiable:

m = tetrahedralize_box(1.0, 1.0, 1.0, 10, 10, 10)
n = size(m.nodes, 2)
fixed = Int[]; pull = Int[]
for i in 1:n
    x = m.nodes[:, i]
    if x[1] ≈ 0.0
        append!(fixed, 3i - 2, 3i - 1, 3i)
    elseif x[1] ≈ 1.0
        push!(pull, i)
    end
end
loads = zeros(3n)
for i in pull
    loads[3i - 2] = 1.0 / length(pull)
end

φ, hist = solve(ComplianceTopOpt(), m, Material(youngs_modulus=E, poisson_ratio=ν);
                fixed = fixed, loads = loads, volfrac = 0.5, max_iter = 40)
println("compliance trace: ", round(ustrip(hist[end]); digits = 4),
        ", volume fraction = ", round(sum(φ) / size(m.elements, 2); digits = 3))

compliance trace: 0.008, volume fraction = 0.5

Compliant mechanisms (topopt_mechanism, see multiphysics-mechanisms) use the same machinery with a mutual-energy objective, and heat conduction has the analogous ThermalComplianceTopOpt.