Abstract: The random ray method is a promising approach for solving three-dimensional heterogeneous neutron transport problems. It samples rays randomly in the spatial and angular phase space and generates ray trajectories during the transport sweep. Therefore, it avoids the storage of a large number of predefined ray tracks and angular fluxes. This feature provides favorable memory efficiency and makes the method suitable for large-scale three-dimensional reactor calculations. However, the fission source generally converges slowly in eigenvalue calculations. A large number of inactive generations is often required before active tallies can be collected. As a result, the source-convergence stage may account for a considerable part of the total computational cost. This problem limits the practical efficiency of the random ray method in high-fidelity reactor analysis. The objective of this study is to accelerate fission-source convergence in three-dimensional random ray transport calculations while maintaining the conservation consistency between the random ray transport solution and the coarse-mesh equation. A consistent coarse-mesh finite difference, or consistent CMFD, acceleration method was developed in this work. A coarse-mesh neutron balance equation was constructed by using quantities tallied during the random ray transport calculation. In addition to the neutron net currents across coarse-mesh surfaces, neutron current contributions associated with the starting and ending points of active ray segments were introduced. In the random ray method, each active ray segment has a finite length. Its starting and ending points do not generally coincide with coarse-mesh boundaries. If these endpoint contributions are neglected, the neutron balance represented by the coarse-mesh equation may differ from that represented by the random ray transport calculation. This inconsistency may weaken the effectiveness of the CMFD acceleration and may also affect the stability of the feedback process. To address this problem, the balance residual caused by the endpoints of active ray segments was evaluated from the random ray tally information. The endpoint balance residual was then converted into a pseudo-absorption correction term in the coarse-mesh equation. With this treatment, the coarse-mesh equation reproduced the neutron conservation relation of the random ray transport calculation more consistently. The equation retained the conventional coarse-mesh balance form while accounting for the additional balance contribution introduced by the finite active ray segments. A statistical-window strategy was also adopted to reduce the influence of statistical fluctuations on the coarse-mesh equation. Each statistical window consisted of several consecutive inactive generations. During each window, the coarse-mesh fluxes, surface net currents, and endpoint balance quantities were accumulated. The CMFD equation was not reconstructed and solved after every individual generation. Instead, the tally quantities averaged over the statistical window were used to construct the consistent CMFD equation after the window was completed. The solution of the coarse-mesh equation was then used to update the fission source for the subsequent random ray transport calculation. This treatment reduced the influence of generation-wise statistical fluctuations on the CMFD coefficients and feedback correction. The proposed method was implemented in the NECP-MCX code. Numerical verification was performed using the C5G7-3D benchmark and the C5G7-3D Extension benchmark. These problems contain representative three-dimensional heterogeneous reactor configurations and provide suitable test cases for evaluating both calculation accuracy and acceleration performance. Calculations with and without consistent CMFD acceleration were compared under equivalent computational conditions. The eigenvalue and spatial power results were examined to confirm that the acceleration procedure did not introduce a noticeable loss of accuracy. The computational time required for the inactive generations was used to evaluate the acceleration performance. The numerical results show that the proposed consistent CMFD method effectively accelerates fission-source convergence while maintaining calculation accuracy. For the C5G7-3D benchmark, the maximum speedup in the computational time of the inactive generations reaches 5.67. For the C5G7-3D Extension benchmark, the speedup ranges from 2.31 to 2.49. The accelerated solutions remain consistent with the corresponding unaccelerated solutions in the quantities used to assess calculation accuracy. These results demonstrate that the method can substantially reduce the computational cost of the source-convergence stage for different three-dimensional benchmark configurations. In conclusion, the proposed consistent CMFD method provides an effective acceleration strategy for three-dimensional random ray transport calculations. The introduction of endpoint current contributions accounts for the neutron balance effect of finite active ray segments. The conversion of the endpoint balance residual into a pseudo-absorption correction term maintains conservation consistency between the coarse-mesh equation and the random ray transport calculation. The statistical-window strategy reduces the influence of statistical fluctuations on the CMFD feedback process. The benchmark results confirm that the method maintains calculation accuracy, accelerates fission-source convergence, and reduces the computational cost of inactive generations. Therefore, the method has good potential for efficient high-fidelity three-dimensional heterogeneous reactor transport analysis based on the random ray method.