Numerical modeling of  nonlinear, fully coupled incompressible magnetohydrodynamics problem in three dimensions is usually a challenge work  because of  its complex features.  Accuracy of FEMs for the MHD system is deteriorated by the presence of high Hartmann numbers, interior or boundary layers, shock fronts and complex geometries, which will be frequently encountered in practice.   It  makes the adaptive mesh refinement indispensable for the  approximations by   finite element methods. We consider a mixed finite element method for the numerical discretization of a stationary incompressible magnetohydrodynamics problem in three dimensions with its velocity field is discretized using $H^1$ conforming elements and the magnetic field is approximated by curl-conforming Nédélec elements. Under the assumption that the original model has a unique solution pair, we derive a posteriori error estimates of the incompressible magnetohydrodynamic (MHD) equations with a sharp upper bound. Using these a posteriori error estimates, we construct an adaptive algorithm for computing the solution of 3D magnetohydrodynamics. Numerical experiments are carried out to show the performance of the adaptive finite element method.