Modeling three-dimensional magnetohydrostatic equilibria in stellarators is of paramount importance in the design and conduction of fusion experiments. The workhorse methods of the community at this point (VMEC and DESC) simplify the search for equilibria to configurations with nested flux surfaces. It is known, however, that relevant configurations in experiments feature magnetic islands and chaotic regions, violating this assumption. Magnetic relaxation codes provide an opportunity to solve for more general equilibrium solutions. We introduce a structure-preserving magnetic-relaxation solver that does not require a nested flux surface assumption. The novelty of MRX is that it combines several crucial features for the first time. (i) Through structure-preserving mixed finite-elements, we retain central features of the relaxation problem in the discrete setting, such as div B = 0 to machine precision, guaranteed energy-dissipation and helicity preservation. (ii) The use of B-Spline bases and tools from the isogeometric analysis field allows us to use high-order basis functions that exhibit rapid convergence for regular solutions as well as non-uniform meshes. (iii) The code is pure Python using the JAX ecosystem, making it very easy to install, run, and extend as well as performant on GPUs. It is also fully differentiable for future inverse design and optimization applications. The proposed magnetic relaxation solver is tested in several stellarator geometries at low and high values of beta. Future work will address the integration of this code for 3D equilibrium optimization for modeling magnetic islands and chaos in stellarator fusion devices.