Thermoacoustic instabilities, which arise due to the interaction between flames and acoustics, are sensitive to small changes to the system parameters. In this paper, we apply adjoint-based shape optimization to a 2D finite element Helmholtz solver to find accurately and inexpensively the shape changes that most stabilize a 2D thermoacoustic system in the linear regime. We examine two cases: a Rijke tube and a turbulent swirl combustor. Both systems exhibit an unstable longitudinal mode and we suppress the instability by slightly modifying the geometry. In the case of the turbulent swirl combustor, the sensitivities are higher in the plenum and in the burner than in the combustion chamber, mainly due to the effect of the mean temperature. In the cooler regions, the local wavelength is shorter, which means that geometry changes of a given distance have more influence than they do where the local wavelength is longer. This is the first time adjoint-based shape optimization is applied to 2D Helmholtz solvers in thermoacoustics, after being previously applied to low-order thermoacoustic networks. But Helmholtz solvers have an intrinsic advantage: they can handle complex geometries. The easy scalability of this method to complex 3D geometries makes this tool a strong candidate for the iterative design of thermoacoustically stable combustors.