Gradient-based algorithms used in aerodynamic optimization usually rely on the adjoint method to compute the required gradients. The continuous adjoint method offers a clear insight into the adjoint equations and their boundary conditions, but discretization schemes may significantly affect the accuracy of gradients. On the other hand, the discrete adjoint method computes sensitivities consistent with the discretized flow equations, with a higher memory footprint though. This work bridges the gap between the two adjoint variants by proposing consistent discretization schemes for the (continuous) adjoint PDEs and the adjoint boundary conditions; these are developed and justified, based on the discrete adjoint equations that have also been derived and used as reference. This article is dealing with inviscid flow models for compressible flows; the Roe’s upwind scheme is used to discretize the inviscid fluxes. Development has been made on the in-house GPU-accelerated PUMA code, which employs a multi-stage Runge-Kutta scheme for the integration of the governing equations in pseudo-time and Mixed Precision Arithmetics to reduce the memory footprint without affecting code’s accuracy. The developed continuous adjoint with the proposed discretization schemes replicates the outcome of discrete adjoint, avoiding its memory-related flaws. The capabilities of the new (consistent) continuous adjoint are demonstrated in 2D/3D cases. Initially, the accuracy of the sensitivity derivatives computed based on the “classical” and consistent adjoint schemes running on different meshes (from coarse to fine), as well as its impact on the optimization convergence are investigated. An isolated airfoil is selected for the parametric study. The sensitivities of the consistent adjoint are in excellent agreement with finite differences. Once the gradient accuracy has been verified, in the last part of this paper, the constrained shape optimization of an isolated wing is carried out.