We present a higher-order convergent numerical solver for active polar hydrodynamics in three-dimensional domains of arbitrary shape, along with a scalable open-source software implementation for shared- and distributed-memory parallel computers. This enables the computational study of the nonlinear dynamics of out-of-equilibrium materials from first principles. We numerically solve the nonlinear active Ericksen-Leslie hydrodynamic equations of three-dimensional (3D) active nematics using both a meshfree and a hybrid particle-mesh method in either the Eulerian or Lagrangian frame of reference. The solver is validated against a newly derived analytical solution in 3D and implemented using the OpenFPM software library for scalable scientific computing. We then apply the presented method to studying the transition of 3D active polar fluids to spatiotemporal chaos, the emergence of coherent angular motion in a 3D annulus, and chiral vortices in symmetric and asymmetric 3D shapes resembling dividing cells. Overall, this provides a robust and efficient open-source simulation framework for 3D active matter with verified numerical convergence and scalability on parallel computers.