Geode

Geode is a single C++20 engine that runs the entire generative-design loop on one machine: SDF geometry voxelized onto one grid, co-simulated across thermal, fluid, structural, electromagnetic and acoustic domains, then pushed through gradient-based topology optimization under real LPBF and CNC manufacturability constraints, with no geometry shuttling between separate simulation and optimization tools.

Topology optimization inverts the usual workflow: instead of drawing a part and validating it with simulation, you state the physics (where it's bolted, where the load enters, how loud it may be, how much material you get) and the shape becomes the solver's output. Commercial tools do this, but as cloud black boxes. Geode is my own version of the loop, in code readable to the last intrinsic, and the hardest training ground I could design for myself: adjoint math, multiphysics solvers, task graphs, SIMD kernels, cache and memory discipline, numerical conditioning, the whole low-level stack in one codebase. It also served as a testbed for how far LLM-assisted coding holds up inside a genuinely complex system, and it isn't only an exercise: its first 3D-print test was a drone wing optimized for stability and silence at once.

Geode is a C++20 computational engineering engine (~11k LoC of first-party engine code, ~300 GoogleTest cases) that runs the full generative loop: implicit geometry (signed-distance functions) -> voxelization -> multiphysics simulation -> constraint evaluation -> adjoint sensitivities -> SIMP density update -> repeat until the shape converges. Six physics domains are wired into the pipeline: thermal diffusion (FD), fluid (Stokes / Navier-Stokes / RANS), linear elasticity (voxel FE), electrostatics, magnetostatics, and Helmholtz acoustics, each with a matching sensitivity path so any of them can drive the optimizer. Geometry in and out goes through PicoGK, LEAP71's open voxel kernel: SDF voxelization, marching-cubes meshing, an OpenGL viewer, plus OpenVDB export for everything downstream.

The user-facing API is a fluent builder. This is the entire definition of a drone wing that is simultaneously stiff under lift load and quiet at the rotor's blade-pass frequency:

auto model =
  Geode::Engine::GeodeModel::Create("drone_wing")
    .Material("Ti64")
    .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();

Geode::Engine::CEMRunner runner;
runner.SetModel(std::move(model));
runner.SetOptimizer(std::make_unique<Geode::Optimization::SIMP>());
runner.SetResolution(80, 24, 6, 0.0025f);
runner.Build();
runner.Run();   // 150 iterations later: a wing nobody drew

Everything below is about what happens inside Run(), and about why the hard part turned out not to be speed, but keeping fifty chained numerical systems alive for 150 iterations.

01How it works

The iteration as a DAG, built once

Each optimization iteration is a task graph: Voxelize -> Material -> {Thermal, Fluid, Structural, [Electro, Magneto, Acoustic]} -> Objective -> Sensitivity -> Optimize. The Pipeline class owns tasks as name-keyed std::functions with string dependencies, validates the graph itself (Kahn's algorithm, so a cycle or unknown dependency is a build error with a name in the log, not a Taskflow assertion), then lowers it onto tf::Taskflow. The DAG is constructed once in CEMRunner::Build(); RunUntilConverged() just re-executes it. Physics domains that don't share data run concurrently on the executor; domains the model never mentions (no acoustic BC, no acoustic objective) are never added to the graph at all, so you don't pay for physics you didn't ask for.

One memory discipline for everything

All field data lives in FieldBuffer<T>: a move-only, 64-byte-aligned (std::aligned_alloc), flat SoA allocation. No Eigen in the core, no std::vector<Vec3> anywhere hot: displacement is three separate scalar arrays, interleaved only as 3*i + dof inside the FE solver. Scalar is float throughout: half the memory traffic and twice the SIMD lanes of double, a choice that later extracted its price in conditioning, as covered later in "Where it went sideways." On top of that sits SIMD::VectorOps: compile-time dispatched AVX-512 / AVX2 / scalar kernels for the CG primitives (Dot, Axpy with FMA, Fill, Scale) using aligned loads, and 7-point stencil kernels (Laplacian3D, Jacobi3D) using unaligned loads because a stencil's neighbors are never all aligned at once.

The FE solver: never assemble K

The structural solver is where ~90% of wall time goes, and it never builds a stiffness matrix. Elements are 8-node hexes on the voxel grid, so every element has the same geometry. One 24x24 unit stiffness Ke0 is integrated once (full 2x2x2 Gauss; why not 1-point is covered later, in "Where it went sideways") and each element's contribution is just (E[e]/Eref) * Ke0. That's 576 floats, 2.3 KB, permanently resident in L1. The solve is preconditioned conjugate gradient (Jacobi/diagonal preconditioner, assembled from the same Ke0 diagonal), and the operator application is gather -> 24x24 mat-vec -> scatter, element by element.

The scatter is the parallelization problem: adjacent elements share nodes, so naive parallel accumulation races. The regular grid gives a closed-form answer, 8-color graph coloring:

// 8-color graph coloring for structured hex meshes.
// color = (i%2) + 2*(j%2) + 4*(k%2)
// Elements of the same color are >=2 apart in every axis -> share no nodes
// -> the scatter step is race-free and can run in parallel.
//
// One persistent parallel region wraps all 8 color passes; explicit barriers
// separate passes instead of 8x fork/join (avoids ~40k fork/joins per FE solve).
#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 fully scattered before N+1 reads any node
  }
}

Two details matter more than the coloring itself. First, the parallel region is opened once around all eight passes with explicit barriers between them; opening it per color would mean 8 fork/joins per mat-vec x ~5k CG iterations ~ 40,000 thread-team spinups per solve. Second, the inner 24x24 mat-vec is where 24 turns out to be a gift: 24 = 3 x 8 maps exactly onto three AVX2 FMA operations per row, with no tail loop, no masking, and no padding:

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);
}

The optimizer: SIMP with the failure modes designed out

The density update is classic SIMP with an Optimality-Criteria step: penalized stiffness E(ρ) = E_min + ρ³(E₀ - E_min), volume constraint enforced by bisecting the Lagrange multiplier, and a cone-weighted sensitivity filter to kill checkerboard patterns (the filter is embarrassingly parallel: each output reads a neighborhood and writes one cell, so it's a single omp parallel for collapse(2)). Around that textbook core sit the decisions that made it actually converge, each recorded next to the code it protects:

• Start at ρ = 0.5, not 1.0. A solid start forces the first OC step to shed 60% of material before the sensitivity field means anything, locking in boundary-peeled topology. A neutral start lets structure form from the inside out.
• Move limit 0.1 so the topology changes slower than the FE conditioning degrades.
• *Sensitivity = per-element strain energy `u_eᵀKe0u_e, not ‖u‖²`.* The wrong (and tempting) formula scores the free tip (large displacement, zero stress) as the most important material in the model, and the clamped root as worthless. It inverts the physics.
• Cold-start CG every iteration. Warm-starting from last iteration's non-converged displacement compounds error across the outer loop.
• Uniform speed of sound in the acoustic domain. SIMP-penalizing c sends k²h² >> 6 in void cells and the Helmholtz Jacobi sweep (denominator 6 - k²h²) diverges; void fraction controls absorption instead.

Multiphysics objectives and their adjoints

Every objective the builder accepts has a sensitivity path: structural compliance uses the standard SIMP adjoint dC/dρ = -p*ρ^{p-1}(E₀-E_min)*u_eᵀKe0*u_e (self-adjoint, so it's free once you have the strain energies); thermal, electrostatic, and magnetostatic sensitivities are derived analytically from their respective energy functionals; fluid and acoustic use dedicated adjoint modules (the acoustic one weights by distance to the observer microphone). Manufacturing constraints (minimum feature size and maximum overhang angle for laser powder-bed fusion) enter the same gradient accumulation as the physics, so printability shapes the topology rather than rejecting it afterward.

02Trade-offs

03Benchmarks

The headline optimization is ApplyStiffness, the PCG operator that dominates the loop. Measured on the drone-wing case (80x24x6 grid, 11,520 elements, 34,560 structural DOFs), Ryzen 5 7600X (6C/12T), same convergence tolerances before and after:

That's an 8x wall-time speedup at 704% CPU utilization, i.e. ~59% parallel efficiency across 12 hardware threads, which is about what a gather/scatter-bound kernel deserves: the 24x24 FMA block is compute, but the node gather and colored scatter are memory traffic, and six cores share the memory controller. The three contributions stack roughly as: SIMD on the mat-vec (per-element arithmetic), coloring (parallelism without atomics), and the persistent parallel region (recovering the time the first parallel version spent forking). A conditioning fix rode along in the same commit and was worth as much as the hardware work: raising the SIMP stiffness floor from 1e-3 absolute to 0.1*E₀ cut the global condition number from ~6.4x10⁶ to ~6.4x10⁵, which is the difference between PCG converging in the iteration budget and the solver silently returning MaxIterationsReached garbage for the optimizer to chase.

04Where it went sideways

The NaN that impersonated a threading bug. The optimizer appeared to stall, and the pipeline appeared to run single-threaded; every symptom said scheduler. The actual culprit: the OC update computes xNew = ρ*sqrt(B) where B = -sensitivity/λ. A positive sensitivity (which multi-objective weighting can legitimately produce) made B negative, sqrt returned NaN, every NaN comparison in the bisection went the same way, and the loop ground through 150+ float-underflow iterations instead of ~10, per voxel, per step. Nothing crashed. NaN doesn't crash; it makes healthy code slow and weird somewhere else. The fix is a one-line clamp; the lesson is permanent: in iterative numerics, treat "it's mysteriously slow" as "something upstream is NaN" until proven otherwise.

The sensitivity that optimized the part backwards. The first compliance gradient used E*‖u‖² per cell. Plausible-looking (stiffness times deformation) and exactly inverted: displacement is largest at the free tip where the material does nothing, and zero at the clamped root where the bending moment peaks. The optimizer dutifully reinforced the tip and hollowed the root. The correct quantity is per-element strain energy u_eᵀ*Ke0*u_e. No amount of profiling finds this class of bug; only knowing what the field should look like does, which is why the engine grew LogHealth summaries (min/mean/max/NaN-count per stage per iteration) as a first-class feature.

The stiffness matrix that wasn't. Two structural failures hid in the element formulation itself. The original "FE solver" applied a 6-point Laplacian stencil: a scalar heat operator wearing an elasticity costume. It produced smooth, physical-looking, wrong displacement fields, and had to be replaced with real hexahedral assembly. Then the real element, integrated at 1 Gauss point for speed, brought 12 hourglass modes, zero-energy deformations that make K indefinite, which CG (whose convergence theory starts at "K is SPD") punished with non-convergence rather than an error message. Full 2x2x2 integration is exact for trilinear hexes and closed the case. Iterative solvers don't diagnose your discretization; they just don't converge.

The two-language architecture. Geode originally shipped as a C++ engine behind a flat extern "C" API (geode_add_objective(engine, "Acoustic", "MinSPL", 3.0f)...) with a C# frontend driving PicoGK. It worked, and it was miserable to debug: every investigation crossed a P/Invoke boundary where types flatten to strings and floats, and the marshaling layer was pure maintenance mass. The C# layer was deleted and replaced by the GeodeModel fluent builder in native C++; the string-keyed rule API it had forced into existence stayed, and turned out to be a good extensibility seam anyway. Bonus lesson from the same integration: linking against prebuilt PicoGK meant reconciling three different Boost sonames on one machine; the repo's justfile still carries the init-boost-shims target as a scar.

05What I'd do better next time

The mental model I started with was "make the inner loop fast." The model I ended with is: a topology optimizer is a machine for driving its own linear systems toward ill-conditioning, and everything else is secondary. As SIMP sharpens the design, stiffness contrast between solid and void explodes, κ(K) climbs, and PCG iteration counts, not FLOPs, set the wall time. Half the entries in "Where it went sideways" are this one phenomenon in different clothes. The performance work and the conditioning work turned out to be the same work.

That reframing orders the roadmap. First, geometric multigrid as the CG preconditioner: on a structured voxel grid, coarsening is trivial by construction, and it's the principled cure for both the grid-size and contrast scaling; it would also let the 0.1*E₀ stiffness floor (a solvability hack that makes voids ten times too stiff) retreat back toward physical values. Second, the GPU port: gather / 24x24-FMA / colored-scatter is embarrassingly close to an ideal CUDA kernel, and the CUDABackend stub already reserves the seam. Third, honesty about the weaker adjoints: structural, thermal, electrostatic and magnetostatic sensitivities are derived from their energy functionals, but the fluid and acoustic paths are engineered approximations, and a proper discrete-adjoint framework is what would let me trust a Pareto front between them. And the piece of this codebase I lean on most isn't a solver: it's the per-stage field-health logging, which is the only reason bugs of the "physically plausible but wrong" class were findable at all.

06Further reading

• Geode. Engine source, the drone-wing / heatsink / bracket examples, and the ~300-case test suite (repository not public yet).
• PicoGK (LEAP71). The open-source voxel geometry kernel Geode uses for SDF voxelization, marching-cubes meshing, and visualization.
• Taskflow. The C++ task-graph runtime scheduling Geode's per-iteration DAG.
• OpenVDB. The sparse volume format Geode exports optimized density fields to.
• Bendsøe, M. P. (1989). "Optimal shape design as a material distribution problem." Structural Optimization, 1(4), 193-202. The paper that introduced density-based (SIMP-style) topology optimization.
• Sigmund, O. (2001). "A 99 line topology optimization code written in Matlab." Structural and Multidisciplinary Optimization, 21(2), 120-127. The canonical SIMP + Optimality-Criteria formulation Geode's optimizer generalizes to 3D multiphysics.
• Andreassen, E., Clausen, A., Schevenels, M., Lazarov, B. S., & Sigmund, O. (2011). "Efficient topology optimization in MATLAB using 88 lines of code." Structural and Multidisciplinary Optimization, 43(1), 1-16. A more modern SIMP reference implementation, including the sensitivity filtering Geode's omp parallel for collapse(2) filter is based on.
• Bendsøe, M. P., & Sigmund, O. (2003). Topology Optimization: Theory, Methods, and Applications. Springer. The theory behind penalization, filtering, and checkerboard control.
• Flanagan, D. P., & Belytschko, T. (1981). "A uniform strain hexahedron and quadrilateral with orthogonal hourglass control." International Journal for Numerical Methods in Engineering, 17(5), 679-706. The classic treatment of the hourglass-mode failure Geode's under-integrated element hit (see "Where it went sideways").
• Langelaar, M. (2017). "An additive manufacturing filter for topology optimization of print-ready designs." Structural and Multidisciplinary Optimization, 55(3), 871-883. Manufacturability filtering for LPBF overhang and minimum feature size, the same class of constraint Geode's MinFeatureSize enforces.
• Briggs, W. L., Henson, V. E., & McCormick, S. F. (2000). A Multigrid Tutorial (2nd ed.). SIAM. Background for the geometric multigrid preconditioner planned in "What I'd do better next time."
• Saad, Y. (2003). Iterative Methods for Sparse Linear Systems (2nd ed.). SIAM. Preconditioned conjugate gradient and the conditioning theory behind Geode's PCG solver.