Magnetorheological fluids (MRFs) are smart materials consisting of micro-scale magnetizable particles suspended in a carrier fluid. The rheological properties of a MRF can be changed from a fluid-state to a solid-state upon the application of an external magnetic field. This study reports the development of a particle-level simulation code for magnetic solid spheres moving through an incompressible Newtonian carrier fluid. The numerical algorithm is implemented within an open-source finite-volume solver coupled with an immersed boundary method (FVM-IBM) to perform fully-resolved simulations. The particulate phase of the MRF is modeled using the discrete element method (DEM). The resultant force acting on the particles due to the external magnetic field (i.e., magnetostatic polarization force) is computed based on the Clausius-Mossotti relationship. The fixed and mutual dipole magnetic models are then used to account for the magnetic (MAG) interactions between particles. Several benchmark flows were simulated using the newly-developed FVM-IBM-DEM-MAG algorithm to assess the accuracy and robustness of the calculations. First, the sedimentation of two spheres in a rectangular duct containing a Newtonian fluid is computed without the presence of an external magnetic field, mimicking the so-called drafting-kissing-tumbling (DKT) phenomenon. The numerical results obtained for the DKT case study are verified against published data from the scientific literature. Second, we activate both the magnetostatic polarization and the dipole-dipole forces and resultant torques between the spheres for the DKT case study. Next, we study the robustness of the FVM-IBM-DEM-MAG solver by computing multi-particle chaining (i.e., particle assembly) in a two-dimensional (2D) domain for area volume fractions of 20% (260 particles) and 30% (390 particles) under vertical and horizontal magnetic fields. Finally, the fourth computational experiment describes the multi-particle chaining in a three-dimensional (3D) domain allowing to study fully-resolved MRF simulations of 580 magnetic particles under vertical and horizontal magnetic fields.