Summary: Coupling recent imaging capabilities with microstructural finite element ($μ$FE) analysis offers a powerful tool to determine bone stiffness and strength. It shows high potential to improve the individual fracture risk prediction, a tool much needed in the diagnosis and treatment of osteoporosis, that is, according to the World Health Organization (WHO), second only to cardiovascular disease as a leading health-care problem. We adapted a multilevel preconditioned conjugate gradient method to solve the very large voxel models that arise in the $μ$FE bone structure analysis. The intricate microstructure properties of bone lead to sparse matrices with billions of rows, thus rendering this application to be an ideal candidate for massively parallel architectures such as the BG/L Supercomputer. In this work we present our progress as well as the challenges we were able to identify in our quest to achieve scalability to thousands of BG/L cores.