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 as cloud black boxes. Geode is my own version, in code readable to the last intrinsic: adjoint math, multiphysics solvers, task graphs, SIMD kernels, cache and memory discipline, numerical conditioning, the whole low-level stack in one codebase, and a testbed for how far LLM-assisted coding holds up in a genuinely complex system. 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) running the full generative loop: implicit geometry (SDF) -> voxelization -> multiphysics simulation -> constraint evaluation -> adjoint sensitivities -> SIMP density update -> repeat until the shape converges. Six physics domains are wired in: 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 goes through PicoGK, LEAP71's open voxel kernel: SDF voxelization, marching-cubes meshing, an OpenGL viewer, plus OpenVDB export.

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 why the hard part wasn't 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 built once in CEMRunner::Build(); RunUntilConverged() just re-executes it. Domains that don't share data run concurrently; domains the model never mentions are never added to the graph, 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 (see "Where it went sideways"). On top 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 since a stencil's neighbors are never all aligned at once.

The FE solver: never assemble K

The structural solver takes ~90% of wall time, and 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; see "Where it went sideways" for why not 1-point) and each element's contribution is (E[e]/Eref) * Ke0: 576 floats, 2.3 KB, permanently resident in L1. The solve is preconditioned conjugate gradient (Jacobi/diagonal preconditioner from the same Ke0 diagonal); 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 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 opens once around all eight passes with explicit barriers between them; per-color opening 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: 24 = 3 x 8 maps exactly onto three AVX2 FMA operations per row, no tail loop, no masking, 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

Classic SIMP with an Optimality-Criteria step: penalized stiffness E(ρ) = E_min + ρ³(E₀ - E_min), volume constraint by bisecting the Lagrange multiplier, cone-weighted sensitivity filter against checkerboard patterns (embarrassingly parallel: a single omp parallel for collapse(2)). Around that core sit the decisions that made it converge:

• 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 tempting-but-wrong formula scores the free tip (large displacement, zero stress) as most important 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, free once you have strain energies); thermal, electrostatic, and magnetostatic sensitivities are derived analytically from their energy functionals; fluid and acoustic use dedicated adjoint modules (acoustic weights by distance to the observer microphone). Manufacturing constraints (minimum feature size, maximum overhang angle for LPBF) enter the same gradient accumulation as the physics, so printability shapes the topology instead of rejecting it after.

02Trade-offs

03Benchmarks

The headline optimization is ApplyStiffness, the PCG operator dominating 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/after:

An 8x wall-time speedup at 704% CPU utilization, ~59% parallel efficiency across 12 hardware threads, about what a gather/scatter-bound kernel deserves: the 24x24 FMA block is compute, but 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, coloring (parallelism without atomics), and the persistent parallel region (recovering 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⁵, the difference between PCG converging in budget and silently returning MaxIterationsReached garbage.

04Where it went sideways

The NaN that impersonated a threading bug. The optimizer stalled and the pipeline looked 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+ iterations instead of ~10. Nothing crashed; NaN makes healthy code slow and weird somewhere else. The fix is a one-line clamp; the lesson: treat "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, and exactly inverted. Displacement is largest at the free tip where material does nothing, zero at the clamped root where the bending moment peaks. The optimizer reinforced the tip and hollowed the root. The correct quantity is per-element strain energy u_eᵀ*Ke0*u_e. No profiling finds this bug class; 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).

The stiffness matrix that wasn't. The original "FE solver" applied a 6-point Laplacian stencil, a scalar heat operator wearing an elasticity costume: smooth, physical-looking, wrong displacement fields, 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. 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 was miserable to debug: every investigation crossed a P/Invoke boundary where types flatten to strings and floats. The C# layer was deleted, replaced by the GeodeModel fluent builder in native C++; the string-keyed rule API it forced into existence stayed, and turned out to be a good extensibility seam anyway. Bonus lesson: 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: 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 of "Where it went sideways" is this phenomenon in different clothes.

That reframing orders the roadmap. First, geometric multigrid as the CG preconditioner: coarsening is trivial on a structured voxel grid, and it's the principled cure for both grid-size and contrast scaling; it would let the 0.1*E₀ stiffness floor retreat toward physical values. Second, the GPU port: gather / 24x24-FMA / colored-scatter is 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, 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.