A continuous adjoint method is developed and implemented within the OpenFOAM framework for the gradient-based optimization of cavitating flows modeled using a barotropic homogeneous equilibrium model. The primal solver is formulated using the preconditioned Euler equations with a pseudo-compressibility treatment, allowing a compressible-type numerical framework to be applied to cavitating flows characterized by strong numerical stiffness due to the coexistence of highly compressible mixture regions and nearly incompressible liquid phases. The cell-centered Jameson–Schmidt–Turkel scheme \cite{Jameson}, for unstructured meshes \cite{Stolcis}, solving for pressure and velocity (i.e. primitive flow variables) is used. The adjoint equations are derived with respect to the same set of primitive variables. Consistency of the discretization is achieved using the Think-Discrete Do-Continuous (TDDC) approach, originally proposed in \cite{Kontou} for a single-phase solver, and a vertex-centered finite-volume formulation with the Roe scheme is used. In the TDDC adjoint, the discretized primal residuals are first hand-differentiated and, based on them, discretization of the flux schemes of the continuous adjoint PDEs are derived. This preserves the physical interpretability of continuous adjoints while ensuring accurate/consistent gradient computation. The proposed adjoint solver computes sensitivity derivatives which are in excellent agreement with finite-difference results while maintaining a minimal memory footprint. The method is demonstrated in two hydrofoil cases targeting minimum cavitation. In both cases, a volumetric B-spline shape parameterization is employed, and the computed gradients effectively support the minimization. The first author has received funding from the EU Horizon 2021 research and innovation program under the Marie Sklodowska Curie Grant Agreement No. 101072851 (MFLOPS).