Computational modeling of physical phenomena has enabled researchers to acquire insight that was previously only observable from costly real-world experiments. The adoption of physics simulations has resulted in numerous advancements in fusion energy, seismic inversion/monitoring, and national defense. However, optimizing simulation studies requires many realizations of the intensive simulations. Manually tweaking control parameters and searching for optimal results can be tedious and inefficient. To tackle such obstacles in simulation studies, we found that differential evolution combined with the covariance matrix adaptation strategy could effectively optimize simulations while simultaneously behaving as an acquisition function to collect samples. The samples collected may be used to train intelligent surrogate models such as a Fourier neural operator (FNO). Once a surrogate is constructed, it could be used to accelerate the sampling of the optimization search space further. This methodology effectively optimized a hydrodynamic simulation modeled by systems of partial differential equations; it may also be extended to simulation optimization in other disciplines.