Optimization problems governed by partial differential equations are ubiquitous in modern science, engineering, and mathematics. They play a central role in optimal design and control of multiphysics systems, data assimilation, and inverse problems. However, as the complexity of the underlying PDE increases, efficient and robust methods to accurately compute the objective function and its gradient become paramount. To this end, I will present a globally high-order discretization of PDEs and their quantities of interest and the corresponding fully discrete adjoint method for use in a gradient-based PDE-constrained optimization setting. The framework is applied to solve a slew of optimization problems including the design of energetically optimal flapping motions, the design of energy harvesting mechanisms, and data assimilation to dramatically enhance the resolution of magnetic resonance images. In addition, I will demonstrate that the role of optimization in computational physics extends well beyond these traditional design and control problems. I will introduce a new method for the discovery and high-order accurate resolution of shock waves in compressible flows using PDE-constrained optimization techniques. The key feature of this method is an optimization formulation that aims to align discontinuous features of the solution basis with the discontinuities in the solution. The method is demonstrated on a number of one- and two-dimensional transonic and supersonic flow problems. In all cases, the framework tracks the discontinuity closely with curved mesh elements and provides accurate solutions on extremely coarse meshes.