The confusion lies in writing the formula as you have. If you instead started with Faradays law,
$$
\frac{\partial \mathbf B}{\partial t}=-\nabla\times\mathbf E\tag{1}
$$
under the assumption of perfect conductivity, (1) becomes
$$
\frac{\partial \mathbf B}{\partial t}=\nabla\times\left(\mathbf u\times\mathbf B\right)\tag{2}
$$
You, of course, have a diffusion term added to (2), which comes from $\mathbf J=\eta\nabla\times\mathbf B$ (cf. this post of mine). Combining these in the steady-state approximation & using some vector calculus, you can get the form you have in your first equation.
In terms of discretizing the PDE, it's easier to work with Eq. (2) above than your first equation, so along the $x$ direction, you'll have
$$
\frac{\partial}{\partial x}\left(\mathbf u\times\mathbf B\right)_x=\frac{\partial}{\partial x}\left(\begin{array}{c}0 \\ v_xB_y-v_yB_x \\ v_xB_z - v_zB_x\end{array}\right)
$$
The $y$ and $z$ directions are pretty similar; applying FVM/FDM to these terms should be straight forward.
Note that, while the $x$ term here is zero (required for satisfying $\nabla\cdot\mathbf B=0$), numerical instabilities can still exist and drive inaccuracies in your model; fixing this requires ''divergence cleaning'' methods. There are some tricks in either the discretization scheme or the algorithm to employ that have been discovered over the years to do this cleaning: Balsara & Kim (2003) and Tóth (2000) cover the many of the schemes available (Tóth's paper is paywalled, but the Balsara & Kim paper is on arXiv).