Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Interpolated stats #203

Closed
wants to merge 9 commits into from
Closed

Interpolated stats #203

wants to merge 9 commits into from

Conversation

jtravs
Copy link
Contributor

@jtravs jtravs commented Apr 20, 2021

This is an attempt to fix #202
It saves stats at every natural step point (as previously) but also at every saved grid point based on the interpolated field.
Roughly tested so far.

@jtravs jtravs requested a review from chrisbrahms April 20, 2021 12:55
Copy link
Collaborator

@chrisbrahms chrisbrahms left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good, this is more or less how I thought to do it as well. But the initialisation for HDF5Outputs isn't quite right I think (unless I've missed something?)

src/Output.jl Outdated Show resolved Hide resolved
src/Output.jl Outdated Show resolved Hide resolved
@jtravs
Copy link
Contributor Author

jtravs commented Sep 13, 2021

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.

@chrisbrahms
Copy link
Collaborator

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.

@chrisbrahms
Copy link
Collaborator

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)

@chrisbrahms
Copy link
Collaborator

(and yes I realised I suggested this change)

@jtravs
Copy link
Contributor Author

jtravs commented Sep 13, 2021

. @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!)

@chrisbrahms
Copy link
Collaborator

chrisbrahms commented Sep 13, 2021

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

@chrisbrahms
Copy link
Collaborator

chrisbrahms commented Sep 14, 2021

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.

@jtravs
Copy link
Contributor Author

jtravs commented Sep 14, 2021

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

@chrisbrahms
Copy link
Collaborator

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

@jtravs
Copy link
Contributor Author

jtravs commented Sep 14, 2021

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

@jtravs
Copy link
Contributor Author

jtravs commented Sep 14, 2021

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.

@chrisbrahms
Copy link
Collaborator

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

@chrisbrahms chrisbrahms closed this Aug 2, 2023
@jtravs
Copy link
Contributor Author

jtravs commented Aug 3, 2023

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.

@chrisbrahms
Copy link
Collaborator

See #159, no idea why

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Save stats for z=0
2 participants