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

Optimize VerletLeapfrog method #2559

Open
wants to merge 4 commits into
base: master
Choose a base branch
from

Conversation

efaulhaber
Copy link

@efaulhaber efaulhaber commented Dec 17, 2024

Checklist

  • Appropriate tests were added (does not require additional tests, method is already tested)
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated (no docs changes necessary, method and API didn't change)
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API (no new docs)

Additional context

The VerletLeapfrog method is currently using the generic Symplectic2Cache implementation:

@muladd function perform_step!(integrator, cache::Symplectic2Cache, repeat_step = false)
@unpack t, dt, f, p = integrator
@unpack a1, a2, b1, b2 = cache.tab
duprev, uprev, _, kuprev = load_symp_state(integrator)
du, u, kdu, ku = alloc_symp_state(integrator)
# update position
@.. broadcast=false u=uprev + dt * b1 * kuprev
# update velocity
f.f1(kdu, duprev, u, p, t)
@.. broadcast=false du=duprev + dt * a1 * kdu
# update position & velocity
tnew = t + a1 * dt
f.f2(ku, du, u, p, tnew)
@.. broadcast=false u=u + dt * b2 * ku
f.f1(kdu, du, u, p, tnew)
@.. broadcast=false du=du + dt * a2 * kdu
f.f1(kdu, du, u, p, tnew)
f.f2(ku, du, u, p, tnew)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 3)
integrator.stats.nf2 += 2
store_symp_state!(integrator, cache, kdu, ku)
end

This contains 3 calls to f1 and 2 calls to f2. The kick-drift-kick scheme of the Leapfrog method only requires a single call to f1 and a single call to f2 in each step, which is the main selling point of the method, making the current implementation basically useless.

I implemented an optimized step for this method, following the VelocityVerlet method, which also has its own step function.

@efaulhaber efaulhaber marked this pull request as ready for review December 17, 2024 11:25
@ChrisRackauckas
Copy link
Member

Yeah that looks right, seems it just got over simplified.

The interpolation should be specialized though so that it goes to linear instead of Hermite then, since the derivative would be off. Or a special centered Hermite would need to be used.

@efaulhaber
Copy link
Author

efaulhaber commented Dec 17, 2024

To be honest, I have no idea what I need to change to make this use a linear interpolation. Could you point me to the code that I need to change?

@ChrisRackauckas
Copy link
Member

I just realized we were doing that a bit bespoke. Set this branch up on top of #2560 with the trait and I think it's good to go.

@efaulhaber
Copy link
Author

Like so?

@@ -171,7 +174,7 @@ end
# v(t+Δt) = v(t) + 1/2*(a(t)+a(t+Δt))*Δt
du = duprev + dt * (half * ku + half * kdu)

OrdinaryDiffEqCore.increment_nf!(integrator.stats, 2)
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
Copy link
Author

Choose a reason for hiding this comment

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

Note that this was a bug with the VelocityVerlet method.

@ChrisRackauckas
Copy link
Member

Rebase onto the latest master and I think this is good to go.

@ChrisRackauckas
Copy link
Member

Test fails.

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.

2 participants