In this work, we present a hierarchical domain decomposition method for the multi-scale Boltzmann equation based on moment realizability matrices, a concept introduced by Levermore, Morokoff, and Nadiga in \cite{lev-mor-nad-1998}. This criterion is used to dynamically partition the two-dimensional spatial domain into three regimes: the Euler regime, an intermediate kinetic regime governed by the ES-BGK model, and the full Boltzmann regime. The key advantage of this approach lies in the use of Euler equations in regions where the flow is near hydrodynamic equilibrium, the ES-BGK model in moderately non-equilibrium regions where a fluid description is insufficient but full kinetic resolution is not yet necessary, and the full Boltzmann solver where strong non-equilibrium effects dominate, such as near shocks and boundary layers. This allows for both high accuracy and significant computational savings, as the Euler solver and the ES-BGK models are considerably cheaper than the full kinetic Boltzmann model. To ensure accurate and efficient coupling between regimes, we employ asymptotic-preserving (AP) numerical schemes and fast spectral solvers for evaluating the Boltzmann collision operator. Among the main novelties of this work are the use of a full 2D spatial and 3D velocity decomposition, the integration of three distinct physical regimes within a unified solver framework, and a parallelized implementation exploiting CPU multithreading. This combination enables robust and scalable simulation of multiscale kinetic flows with complex geometries.

39 pages, 9 figures, 1 table