Differentiable physical simulators are proving to be valuable tools for developing data-driven models for computational fluid dynamics (CFD). By enabling end-to-end training of machine learning (ML) models embedded directly within CFD solvers, they allow novel learning algorithms that combine the generalization and efficiency of physics-based simulations with the adaptability of deep learning. In this study, we introduce a differentiable physics framework for embedding deep learning models within a finite element solver for incompressible Navier-Stokes equations, specifically applying this approach to learn a subgrid-scale (SGS) closure with graph neural networks (GNNs) for large eddy simulations (LES). Solver differentiability is implemented through the discrete adjoint method and integrated with the PyTorch ecosystem to enable fully coupled rollout training over many timesteps. Using this framework, we train a SGS closure for a turbulent flow over a three-dimensional backward-facing step (BFS), a canonical separated-flow configuration. The learned GNN-based closure achieves low prediction error and long-term stability, accurately recovering key turbulence statistics and preserving multiscale structures. Generalization capability is evaluated on unseen geometries, including a BFS with increased step height and a channel with a cavity. Even for the cavity flow, which features reattachment on the opposing cavity wall, the learned closure outperforms traditional SGS models in both accuracy and turbulence statistics without retraining. Finally, we demonstrate the practicality of the framework by considering data-limited learning scenarios in which only flow measurements downstream of the step are available during training. Overall, this work establishes a viable pathway toward physically consistent and generalizable data-driven SGS modeling on complex, unstructured meshes through differentiable physics.