hugo_fiorillo::swe
← projects06 / 07 · infra, system, automation

// project 06

Geode

3D object optimizer against physical laws: heat, fluid flow, structural load, electromagnetics, acoustics.

Physic SimulationHPCC++C#
Length
6 min read

Geode is a single C++20 engine running the whole generative-design loop on one machine: SDF geometry voxelized onto one grid, co-simulated across thermal, fluid, structural, electromagnetic and acoustic domains, then topology-optimized under LPBF and CNC manufacturability constraints. No geometry shuttling between tools.

01Spec

  • C++20, ~11k LoC, ~300 GoogleTest cases.
  • Loop: SDF -> voxelize -> multiphysics -> constraints -> adjoints -> SIMP update -> repeat.
  • Domains: thermal (FD), fluid (Stokes / Navier-Stokes / RANS), elasticity (voxel FE), electrostatics, magnetostatics, Helmholtz acoustics. Each carries a sensitivity path, so any drives the optimizer.
  • Geometry via PicoGK: voxelization, marching cubes, viewer, OpenVDB export.

Whole drone wing, stiff under lift, quiet at blade-pass:

cpp
auto model =
  Geode::Engine::GeodeModel::Create("drone_wing")
    .Fix (lo, { lo.x + spacing, hi.y, hi.z })              // root face clamped
    .Load({ lo.x, lo.y, hi.z - spacing }, hi, { 0, 0, 500.f }) // lift on top skin
    .SoundSource(srcLo, srcHi, 1.f, /*blade-pass*/ 500.f)
    .AcousticObserver(observer)                             // ground-facing mic
    .Minimize("Structural", "MinCompliance", 1.f)
    .Minimize("Acoustic",   "MinSPL",        3.f)
    .MaxVolumeFraction(0.40f)
    .MinFeatureSize(0.005f)                                  // LPBF printability
    .Build();

02DAG and memory

Voxelize -> Material -> {Thermal, Fluid, Structural, [Electro, Magneto, Acoustic]} -> Objective -> Sensitivity -> Optimize.

Pipeline: name-keyed std::functions, string deps, Kahn validation, lowered onto tf::Taskflow. A cycle fails at build with a name in the log, not a Taskflow assertion. Built once, re-executed per iteration. Unmentioned domains never enter the graph.

FieldBuffer<T>: move-only, 64-byte-aligned, flat SoA. No Eigen, no std::vector<Vec3> hot. Displacement is three scalar arrays, interleaved 3*i + dof only in the FE solver. float throughout: half the traffic, twice the lanes, priced later in conditioning. SIMD::VectorOps dispatches AVX-512 / AVX2 / scalar at compile time; CG aligned loads, 7-point stencils unaligned.

03FE solver: never assemble K

~90% of wall time. 8-node hexes on a regular grid share geometry: one 24x24 unit Ke0, full 2x2x2 Gauss, 576 floats, 2.3 KB, L1-resident. Element contribution (E[e]/Eref) * Ke0. PCG, Jacobi preconditioner from the Ke0 diagonal, operator = gather -> mat-vec -> scatter. Adjacent elements share nodes, so naive parallel scatter races.

cpp
// color = (i%2) + 2*(j%2) + 4*(k%2); same-color elements share no nodes.
#pragma omp parallel
{
  for (int color = 0; color < 8; ++color) {
    #pragma omp for schedule(static) nowait
    for (Core::Index idx = 0; idx < colorCount[color]; ++idx) {
      /* gather u_e, Ke0*u_e, scatter: no atomics, no locks */
    }
    #pragma omp barrier   // color N scattered before N+1 reads
  }
}

One parallel region wraps all 8 passes. Per-color fork/join = 8 x ~5k CG iterations = ~40k team spinups per solve.

24 = 3 x 8: three AVX2 FMAs per row, no tail loop, no mask, no padding.

cpp
for (int row = 0; row < 24; ++row) {
  const float* rp = Ke0_ + row * 24;
  __m256 acc = _mm256_mul_ps  (_mm256_loadu_ps(rp),      _mm256_loadu_ps(u_e));
  acc = _mm256_fmadd_ps(_mm256_loadu_ps(rp + 8),  _mm256_loadu_ps(u_e + 8),  acc);
  acc = _mm256_fmadd_ps(_mm256_loadu_ps(rp + 16), _mm256_loadu_ps(u_e + 16), acc);
  Ke_u[row] = scale * hsum256(acc);
}

04Optimizer

SIMP + Optimality Criteria: E(ρ) = E_min + ρ³(E₀ - E_min), volume constraint by bisecting λ, cone-weighted filter against checkerboards. What made it converge:

  • ρ = 0.5 start. A solid start sheds 60% of material before sensitivities mean anything. Locks in boundary-peeled topology.
  • Move limit 0.1. Topology changes slower than conditioning degrades.
  • *Sensitivity = `u_eᵀKe0u_e, not ‖u‖²`.* ‖u‖² scores the free tip highest, the clamped root worthless: inverted physics.
  • Cold-start CG each iteration. Warm starts compound non-converged error outward.
  • Uniform speed of sound. Penalizing c drives k²h² >> 6 in void cells, diverging the Helmholtz Jacobi sweep (denominator 6 - k²h²). Void fraction controls absorption.

Compliance adjoint dC/dρ = -p*ρ^{p-1}(E₀-E_min)*u_eᵀKe0*u_e: self-adjoint, free once strain energies exist. Thermal, electrostatic, magnetostatic from energy functionals. Fluid and acoustic use dedicated modules. Manufacturing constraints enter the same gradient accumulation, so printability shapes topology instead of rejecting it after.

05Trade-offs

ChoseGave upWhy it held
Voxel grid + matrix-free FEBody-fitted meshes, boundary accuracyRegularity is the performance model: shared Ke0, closed-form coloring, no meshing. Staircase washes out in marching cubes
float everywhere~7 significant digits2x lanes, half the traffic, memory-bound solver. Conditioning cost manageable once explicit
One unit Ke0 scaled per elementPer-element ν, non-cubic voxels2.3 KB for the whole mesh, cache-hot. SIMP only modulates E, so scaling is exact
String-keyed objectives/domainsCompile-time physics namesAPI had to survive a flat C ABI. RegisterDomain adds domains without touching an enum. Opposite call to Nott: different layer
OpenMP in solvers, Taskflow across themOne schedulerSolver loops want static data parallelism with barriers. The pipeline wants dependency concurrency
Jacobi-preconditioned CGMultigrid convergenceMatrix-free compatible, 30 lines. Iteration counts grow with grid size and SIMP contrast
Native C++ fluent builderC#/PicoGK frontendEvery bug crossed a P/Invoke boundary flattening types to strings and floats. Its string-keyed rule API stayed as a seam

06Benchmarks

ApplyStiffness dominates. Drone wing, 80x24x6 grid, 11,520 elements, 34,560 DOFs, Ryzen 5 7600X (6C/12T), equal tolerances.

30 SIMP iterationsper FE solve150-iteration full run
Before (scalar, single-threaded)~210 s~7 s~17 min
After (8-color OpenMP + AVX2)26.4 s0.75 s~2.2 min

8x wall time at 704% CPU, ~59% parallel efficiency over 12 threads. The FMA block is compute; gather and colored scatter are memory traffic on a shared controller. Same commit, equal payoff: stiffness floor 1e-3 absolute -> 0.1*E₀ cut κ(K) ~6.4x10⁶ -> ~6.4x10⁵. PCG then converges in budget instead of returning MaxIterationsReached garbage.

07Failure modes hit

  • NaN as a fake threading bug. xNew = ρ*sqrt(B), B = -sensitivity/λ. Multi-objective weighting makes sensitivity positive, so B < 0, sqrt NaN. Every NaN compare in the bisection falls one way; the loop runs 150+ iterations instead of ~10. One-line clamp.
  • Inverted sensitivity. E*‖u‖² reinforced the free tip, hollowed the clamped root. No profiler finds this class. Hence LogHealth: min/mean/max/NaN per stage per iteration.
  • The stiffness matrix that wasn't. First "FE solver" was a 6-point Laplacian: a heat operator in elasticity costume. Real hexes at 1-point Gauss then gave 12 hourglass modes, K indefinite, CG answering non-convergence instead of an error. 2x2x2 is exact for trilinear hexes.
  • Three Boost sonames from linking prebuilt PicoGK. justfile still carries init-boost-shims.

08Next

A topology optimizer drives its own linear systems toward ill-conditioning. SIMP sharpens, solid/void contrast explodes, κ(K) climbs, PCG iteration counts set wall time, not FLOPs.

  1. Geometric multigrid preconditioner. Coarsening is trivial on a structured grid. Fixes grid-size and contrast scaling; the 0.1*E₀ floor then retreats toward physical values.
  2. GPU port. gather / 24x24-FMA / colored-scatter maps to CUDA. CUDABackend stub reserves the seam.
  3. Discrete adjoints. Fluid and acoustic sensitivities are engineered approximations. No trustworthy Pareto front until exact.

09Further reading