This paper documents a fully turbulent discrete ad joint method for three-dimensional multistage turbomachinery design. The method is based on a duality preserving algorithm and is implemented in the open-source computational fluid dynamics tool SU2. The SU2 Reynolds-averaged Navier–Stokes solver is first extended to treat three dimensional steady turbomachinery flow using a conservative formulation of the mixing-plane coupled to non reflective boundary conditions. The numerical features of the flow solver are automatically inherited by the discrete ad joint solver, ensuring the same convergence rate of the primal solver. The flow solver is then validated against experimental data available for three turbine configurations, namely, a one-and-half axial turbine stage, a transonic radial turbine coupled to a downstream diffuser, and a supersonic mini–organic Rankine cycle radial turbine operating with a fluid made by a heavy molecule. Finally the ad joint-based optimization framework is applied to the concurrent shape optimization of three rows of the axial turbine, demonstrating the advantages deriving from adopting multi row automated design methods in the context of turbomachinery design.