I thought I'd got this finished, but in its current state I get an error in the tests:

HDF5 vs Memory: Test Failed at C:\Users\John\.julia\dev\Luna\test\test_output.jl:152
  Expression: read((file["stats"])["ω0"]) == (mem.data["stats"])["ω0"]

This works fine on current master. I'm struggling to understand how this occurs. All other tests pass, including the few new ones I added, apart from the frequency centroid.

The problem isn't just with ω0, it's all the stats arrays. @test read(file["stats"]["z"]) == mem.data["stats"]["z"] also fails. The reason is that you inadvertently fixed a bug in HDF5Output but not in MemoryOutput: before we were (incorrectly) saving the final propagation point even if its z-value was actually beyond the end of the propagation. e.g. in test_output the final step of the propagation propagates to 0.05010315819460643 so very slightly beyond the fibre length. Your change to HDF5Output fixes this by appending to the temporary statistics cash after the saving to disk, so for the last propagation step the statistics simply never get saved. For MemoryOutput this doesn't work since the statistics results always get added immediately. Therefore read(file["stats"]["z"]) has one fewer elements than mem.data["stats"]["z"]. @test read(file["stats"]["ω0"]) == mem.data["stats"]["ω0"][1:end-1] passes.

Just looking into how best to fix this.

As an additional issue, your added check:

    if ts != t
        append_stats!(o, o.statsfun(y, t, dt))
    end

doesn't actually do anything--append_stats! is still always called, because at the end of the loop, ts is always 0 (the loop runs until o.save_cond returns (False, 0)

(and yes I realised I suggested this change)

. @test read(file["stats"]["z"]) == mem.data["stats"]["z"] also fails.

That's weird. If I comment out just the ω0 tests everything passes for me. (At least I thought it did!)

That's because the test for the z stat isn't in the test file. I added it just now while debugging.

This can be fixed (and probably should be) by forcing the RK45 solver to finish precisely tmax rather than wherever it wants to, i.e. adding

# make sure we don't propagate beyond the end
if (s.tn + s.dtn) > tmax
    s.dtn = tmax-s.tn
end

in the main loop in RK45.solve() and changing the while loop condition to

while s.tn < tmax
    [...]

instead of <=.
I think this makes sense anyway. For constant waveguides this makes very little difference, but for gradient/taper simulations, propagating past the end of the waveguide very subtly changes the dynamics because suddenly the pressure/core size is fixed again. Since the last slice is usually the most important, we also gain a tiny bit by using the 5th order accurate solution instead of the 4th order accurate interpolant.

I have your branch checked out so could push it to this PR. Just running the whole test suite.

OK, just remember to check for too small step size, which can also cause numerical issues.

I think this might need a bit more thought actually... interleaving the interpolated stats with the directly calculated ones creates weird results:

Well, that was the point. The interpolated field is very slightly different from the non-interpolated one.

The best fix (IMO) is to make the propghataor step to each save position, the old fashioned way, for a very slight loss of performance (probably negligible) but increase in precision.

Sure we can think about that. It's a much bigger overhaul of the Luna internals though, so definitely not going to happen quickly.

Why was this closed? I suppose we could open an issue to cover the same point. But this was a useful reminder to make the propagator stop at exactly the saved steps and not use interpolation. I remember in fnfep the interpolation (done in the same way as in Luna) didn't help much for any decently complex simulation anyway.

See #159, no idea why