With the algorithm’s suitability for exploiting current petascale and next-generation exascale supercomputers, stable and structure-preserving properties are necessary to develop predictive computational tools. In this dissertation, summation-by-parts (SBP) operators and a new relaxation Runge–Kutta (RRK) scheme are used to construct mimetic and structure-preserving full discretization for non-reactive compressible computational fluid dynamics (CFD) and reaction-diffusion models. In the first chapter, we provide the necessary background and a literature survey that forms the basis of this dissertation. Next, we provide a short overview of entropy stability for general conservation laws. The second chapter covers the analysis of the Eulerian model for compressible and heat-conducting flows. We provide the necessary background of the new system of parabolic partial differential equation (PDE). Then, we present the entropy stability analysis of the model at the continuous level. Subsequently, using the SBP, we construct an entropy-stable discretization of any order for unstructured grids with tensor-product elements. The third chapter discusses the implementation of RRK methods. We start by reviewing the RRK scheme constructed to guarantee conservation or stability with respect to any inner-product norm. Then, we present the extension and generalization of RRK schemes to general convex functionals and their application to compressible fluid flow problems. The final chapter demonstrates the far-reaching capabilities of the SBP operators and RRK schemes presenting the development of a novel fully discrete Lyapunov stable discretization for reaction models with spatial diffusion. Finally, we conclude this dissertation with an overview of our achievements and future research directions.