Advancements in computing technology have revolutionized the efficiency and cost-effectiveness of realistic reactor core simulations. The discrete ordinates ( \(S_N\) ) method is a widely adopted approach for numerically solving the Boltzmann transport equation (BTE), which describes neutron distribution in nuclear reactors. Recently, the MT-3000, a novel multizone heterogeneous architecture designed for high-performance computing, has been developed, offering a peak double-precision performance of 11.6 TFLOPS at 1.2 GHz. In this study, we propose an efficient four-level heterogeneous \(S_N\) parallel algorithm for structured hexahedral grids on the MT-3000 system. The algorithm incorporates a two-level KBA strategy with an optimized communication scheme among the MT-3000’s acceleration cores and employs a software caching technique to reduce memory access latency. Numerical experiments reveal that our algorithm achieves 1.68 TFLOPS on a single MT-3000 chip, representing 14.5% of its peak performance. The heterogeneous parallel algorithm demonstrates strong scaling efficiency, consistently exceeding 50% across 4 to 256 MT-3000 systems. Furthermore, we develop a performance model that closely aligns with the experimental results.