With some help I found a major bug. Several indices had a
* instead of
+ when calculating the shift based on the number of passive scalars. After fixing that the constant, no magnetic fields, case will run without producing a non-zero magnetic divergence. The time step does slowly diverge to zero due to some error in (I think) the density/energy interface states. Values of zero are being passed to the HLLD solver which is corrrecting them to 1E-20 as expected but those fields shouldn’t be zero in the first place.