The Non-Orthogonal Multiple Access (NOMA) technique has enormous potential for wireless communications in the fifth generation (5G) and beyond. Researchers have recently become interested in the combination of NOMA and cooperative relay. Even though geometric-based stochastic channel models (GBSM) have been found to provide better, practical, and realistic channel properties of massive multiple-input multipleoutput (mMIMO) systems, the assessment of Cooperative Relay NOMA (CR-NOMA) with mMIMO system is largely based on correlated-based stochastic channel model (CBSM). We believe that this is a result of computational difficulties. Again, not many discussions have been done in academia about how well CR-NOMA systems perform when large antenna transmitters with the GBSM channel model are used. As a result, it is critical to investigate the mMIMO CR-NOMA system with the GBSM channel model that takes into account channel parameters such as path loss, delay profile, and tilt angle. Moreover, the coexistence of large antenna transmitters and coding methods requires additional research. In this research, we propose a twostage, three-dimension (3D) GBSM mMIMO channel model from the 3GPP, in which the transmitter is modelled as a cylindrical array (CA) to investigate the efficiency of CR-NOMA. By defining antenna elements placement vectors using the actual dimensions of the antenna array and incorporating them into the threedimension (3D) channel model, we were able to increase the analytical tractability of the 3D GBSM. Bit-error rates, achievable rates, and outage probabilities (OP) are investigated utilizing the decode-and-forward (DF) coding method: the results are compared with that of a system using the CBSM channel model. Despite the computational difficulties of the proposed GBSM system, there is no difference in performance between CBSM and GBSM.