Post-Newtonian (PN) theory provides the analytic foundation for modeling the early inspiral of binary black holes. However, as an asymptotic series, successive PN orders do not necessarily improve agreement with the full nonlinear dynamics. While this has been explored in the extreme-mass-ratio (EMR) limit, comparable-mass systems most relevant to current observations have not been benchmarked as systematically at high PN order. We study the convergence of the PN series for non-spinning and quasi-circular systems by comparing the PN energy flux at future null infinity to a long, high-accuracy numerical relativity (NR) simulation. To enable a gauge-consistent comparison, we place both descriptions in the same BMS frame and calibrate the intrinsic PN parameters by fitting to the NR waveform in the early inspiral. We find that for orbital velocities $v\lesssim 0.45$, higher PN orders continue to reduce the PN–NR flux discrepancy, with EMR terms beyond 4.5PN providing further improvement in the late inspiral. The improvement with PN order is non-monotonic with local extrema around 2.5PN and 4PN. Thus, for the comparable-mass flux comparison considered here, the first local minimum of the finite-order residual does not identify the best PN truncation order among the orders available. As $v$ approaches $\sim 0.5$ near the innermost circular orbit, higher PN orders no longer improve the agreement between NR and PN, indicating a loss of convergence. These results motivate continued high-order PN calculations and clarify the NR accuracy needed to validate them.