Primordial magnetic fields (PMFs) may explain observations of magnetic fields on extragalactic scales. They are most cleanly constrained by observations of details of the cosmic microwave background radiation (CMB) Their effects on cosmic recombination may even be at the heart of the resolution of the Hubble tension. We present an in-detail analysis of the effects of PMFs on cosmic recombination taking into account of all so far known relevant physical processes. To this end we extend the public magneto-hydrodynamic code ENZO with a new cosmic recombination routine, Monte-Carlo simulations of Lyman-\(\alpha\) photon transport, and a Compton drag term in the baryon momentum equation. The resulting code allows us to predict the impact of PMFs on the cosmic ionization history and the clumping of baryons during cosmic recombination. We study the specific case of non-helical PMFs with a Batchelor spectrum. Our results identify the importance of mixing of Lyman-\(\alpha\) photons between overdense- and underdense- regions for small PMF strength. This mixing shows an attractor to the fully mixed case and speeds up recombination beyond the speed-up due to clumping. It leads to enhanced Silk damping which is strongly constrained by CMB observations. We also show that the increase in the ionization fraction at low redshift by hydrodynamic baryon heating due to PMF dissipation is completely compensated by the faster recombination from baryon clumping. We describe and explain the significant differences of these 3D simulations over earlier three-zone models.