Classical finite-element solvers as differentiable graph programs on GPUs

Abstract Conjugate heat transfer in turbomachinery, heat exchangers, and energy systems couples buoyant flow to the surrounding solid on complex geometries. The cost of simulating it is dominated by large sparse linear solves on CPUs. We recast a mixed finite element discretization on unstructured meshes, and its multigrid pressure solve, as sparse graph layers within an automatic differentiation framework. Nothing is trained: weights are entries of the assembled operator, and manufactured solutions recover theoretical convergence orders. Reverse-mode sensitivities match central finite differences; one gradient costs 5.6 forward solves with unrolling and 4.9 with implicit differentiation, nearly constant from one to 32 parameters. Implicit differentiation through every linear solve substantially reduces truncation-induced sensitivity error and makes memory independent of forward iteration count. Against three established solvers, no single code is fastest; the ranking depends on the regime, the size of the problem, and the accuracy demanded. With the discretization, mesh, time step, and tolerances held fixed, the same problem runs up to 181 times faster per step on one GPU than on 32 CPU cores, and the pressure solve alone runs 18.7 times faster. The graph representation depends on operator sparsity rather than the governing equations, providing a route to other assembled discretizations.