Significant progress has been made in the last decade to ascertain the core collapse supernova explosion mechanism across progenitor characteristics and to predict supernova observables, such as gravitational wave emission. And progress is accelerating. Yet, all three-dimensional models to date can benefit from increased physics and/or a better treatment of the physics. Three-dimensional general relativistic neutrino-magnetohydrodynamics models with spectral or, better yet, Boltzmann kinetics that include all of the relevant neutrino interactions and a properly constrained nuclear equation of state are needed. These models live in the very challenging, higher dimensions of phase space. And models must be “realizable” – i.e., obey Fermi-Dirac statistics – and conserve lepton number and energy simultaneously. It is a tall order to develop discretizations of the underlying integro-partial differential equations and solution methods for the resultant nonlinear algebraic equations that satisfy all of these requirements. It is an equally tall order to in turn develop optimized code based on them. The variety of outcomes, and the variety of conclusions based on them, in what is arguably a very limited number of multi-physics three-dimensional models produced to date begs for continued model improvements and continued modeling across an increasingly larger population of progenitors. This certainly applies to efforts to model core collapse supernova gravitational wave emission. I will discuss these challenges, as well as our efforts to address them, as we move toward the predictive capability we desire and need.