Tuesday, June 22, 2010

[from James]

I've discovered that my integrator isn't working properly; I've tried three different formulations of the leapfrog method, with different initial conditions, but with all of them my error keeps growing. Also, the program gradually slows down, even with only 1 particle, and after about 50,000 force calculations it slows to a crawl (and the relative error climbs close to unity).

leapfrog version 1:

R = r + 0.5*h*v # drift for 1/2 step
T = theta + 0.5*h*omega # drift for 1/2 step

dv = R*omega*omega - GM/(R*R) # acceleration in r
dw = -2.0*omega*v/R # acceleration in theta

v = v + h*dv # kick for full step
omega = omega + h*dw # kick for full step

r = R + 0.5*h*v # drift for remaining 1/2 step
theta = T + 0.5*h*omega # drift for remaining 1/2 step

dE = h*(GM*v/(r*r) + v*dv + r*v*omega*omega + r*r*omega*dw)
error = log(abs(dE/E),10)

Version 2 (this one does an initial 1/2 kick before the loop starts):

# incremental acceleration
dv = L*L/(r*r*r) - GM/(r*r)
dw = -2.0*omega*v/r

# drift for full step
r = r + h*v
theta = theta + h*omega

# Kick for full step - note that because of the initial 1/2 kick, this
# will put the velocities 1/2 step ahead of the positions
v = v + h*dv
omega = omega + h*dw

# compute total energy, as a check:
E2 = 0.5*(v*v + r*r*omega*omega) - GM/r
error = log(abs((E-E2)/E),10)

Version 3 (this one modifies the equations to allow full time steps over
the same time interval, but is mathematically equivalent to the other

dv1 = L*L/(r*r*r) - GM/(r*r) # incremental acceleration in r at t
dw1 = -2.0*omega*v/r

r = r + h*v + 0.5*h*h*dv1
theta = theta + h*omega + 0.5*h*h*dw1

dv2 = L*L/(r*r*r) - GM/(r*r) # incremental acceleration in r at t+dt
dw2 = -2.0*omega*v/r

v = v + 0.5*h*(dv1 + dv2)
omega = omega + 0.5*h*(dw1 + dw2)

E2 = 0.5*(v*v + r*r*omega*omega) - GM/r

You probably noticed that I calculated the error in total energy differently in the first method compared to the second two; the first one calculates dE and then the error as dE/E; the second way is to calculate E2 based on the updated positions and velocities and then the error as (E-E2)/E. Both methods should yield similar errors.


  1. Maybe its slowing down due to a memory issue. If you can't identify it try downloading Parallel Studio from intel or nSight from nvidia. You can use both actually. Parallel Studio is good for analyzing memory bottlenecks and to help developers implement parallel code. nSight does similar things but does not have all the bells and whistles. They both work with CUDA and OMP. Anyway I don't know whats wrong with your integrator. I'm having no such problems. If your motion is smooth then its probably not the integrator but rather the equations for motion that is causing error. If you are using the derived r.. that I did or the one from BT then you you should know that the r(vector).. already takes into account centripetal acceleration. The r(scalar).. that you sub in should only be force due to gravity and tau. In my sim it dosen't maintain perfect angular momentum conservation but fluctuates by about +-50%.

  2. notice that your equations for omega have a simple analytical solution (i.e., you don't need to solve them
    rewriting your equation thus (W = omega for short,
    v = dR/dt as always)

    dW/dt = -2 W v /R
    d(ln W)/dt = -2 d(ln R)/dt
    ln W(t) = -2 ln R(t) + const
    W = L/R^2, where L=const (L=specific ang. mom. = R^2 W)

    this is so simple because we don't have any torques in our problem (any non-radial forces). I'd use that to simplify the numerics by a factor of 2 or so.


  3. .. you need to track T (theta) of course, so you can have
    T = T + h*W equation there, but it's just one not two equations.. and it doesn't need to be done in two half-steps either!

  4. This comment has been removed by the author.

  5. [from James]Thanks for the feedback.
    For r'' I'm following B&T, using equations 3.7a and B.24 (I also re-derived it from scratch just to make sure I understood the derivation); the equation for scalar r'' includes force and r*omega^2 terms (or equivalently force and L^2/r^3 terms).

    For omega, I had switched to simply omega=L/r^2 (after re-calculating r, of course), but the calculation of dE included a theta'' term, and the third leapfrog algorithm required theta'' for theta.

  6. [from James] Woo! Got the error under control. Now on to CUDA (or rather PyCUDA).