Erroneous results from the examples

[EdmundBenedict]

I'm sorry to continuously bother you with problems I've had with running the code, but the results I get from these examples don't seem correct, and I cannot figure out what could be wrong.
I have been comparing the results to the graphs in 170215-merging-flux.pdf. In the case of the single rope in cylindrical coordinates, the value of the central flux is too low and it changes in the wrong direction:

When I plot the parallel current in X and Z, the initial X distribution seems to have the correct shape and size, but it does not decrease by enough at t=5s (it barely changes at all when it should have been spreading out along X) and then the distribution at t=500s actually has a lower peak, and doesn't spread by enough.

The distribution of the parallel current along Z is worse, however: it has large negative steps at the edge of the distribution. Similar to that along X, the distribution does not change greatly at t=5s, and the peak is too low at t=500s. In this case the distribution should have become more narrow as the simulation progressed.

I had thought that the example outlined in 170215-merging-flux.pdf had different parameters (since Bxy is set to 1T rather than 0.5T as it was in the document) and so I have tried to vary the parameters to recover those graphs.
I found that if I changed Bxy to 0.5T and multiplied the Jpar distribution by a factor of -2, I could recover something more correct for the change in flux over time (although the actual change is ~2x what it should be):

And in this case the Jpar distributions are inverted and are around double what they should be:

If I reduce the factor I multiply Jpar to -1, the flux has a more correct change in value over the simulation, but its starting value is around half of what it should be. Similarly, the Jpar distributions' starting peak values are of the correct magnitude but the distributions are still inverted.
If I remove the Bxy scaling from the dx variable, so that it is equal to

dx = Lx / nx

instead of

dx = Bxy * Lx / nx

I find that at least one of the current distributions gets narrower as the calculation proceeds, but in this case it is the wrong one, the current distribution gets narrower along X and wider along Z rather than vice versa:

I have also had difficulty in carrying out the merging example. I found that the merging process itself occurred much too quickly and the change in central flux is once again much too small and in the wrong direction:

This can also be seen in the plots of the parallel current, in which the merging can be seen to occur too early, result in current distributions that are too small, and have a post merging distribution that is quite different (i.e. the the ring of zero current around the centre):

I found that if I increased the current, changed the resistivity and viscosity to the values mentioned in Stanier (2013) and decreased the grid size to that in the example I could get a change in centroid flux that is closer to what is shown there, but the merging still occurs much too quickly:

Do you know what I am doing wrong here? Or what the issue with my code could be? I have tried adjusting everything but I cannot seem to recover the results from the examples.
I am using the BOUT-dev version, with the configuration:

configure: -------------------------
configure: Configuration summary
configure: -------------------------
configure: PETSc support : no (has SUNDIALS: no)
configure: SLEPc support : no
configure: IDA support : yes
configure: CVODE support : yes
configure: ARKODE support : yes
configure: FFTW support : yes
configure: NetCDF support : yes
configure: Parallel-NetCDF support : no
configure: HDF5 support : yes (parallel: no)
configure: MUMPS support : no
configure: Lapack support : yes
configure: Scorep support : no
configure: OpenMP support : no (schedule: )
configure: Natural language support: yes (path: ${datarootdir}/locale)
configure: System UUID generator : no
configure:

Although I have also compiled it with the PetSc library, it did not seem to make any difference.
Thank you, and I am sorry again for the trouble.

[EdmundBenedict]

I went back to the resistive single current filament to see if I could get it to match the given results by tweaking the parameters. I found that I could get the parallel current density profiles in their initial and final stages to match much more closely if I changed the Z-spacing from
dz = Lz / nz
to
dz = Lz * 2*pi / nz
but this did result in the density profile having odd shapes for the t=5 stage:

I also found in this case, that if I ran the calculation for long enough, it could match the magnitude of the change in flux, although the starting value and the sign of the change were still incorrect:

Finally, I decided to try scaling the Jpar distribution (as this seemed to be the only parameter that made a difference, changing Bxy or J had no effect) to see if I could recover the correct change in flux. I found that I could get a very similar plot, but the change in flux was slightly too small:

and of course the parallel current density distributions are now incorrect:

I wanted to check, are the flux evolution and parallel current density graphs definitely from the same calculation? I just cannot find how to recover that starting flux and its change without having negative current densities.
Thank you.