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

versions):

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

error=log(abs((E-E2)/E),10)

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 E

_{2}based on the updated positions and velocities and then the error as (E-E

_{2})/E. Both methods should yield similar errors.

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%.

ReplyDeletenotice that your equations for omega have a simple analytical solution (i.e., you don't need to solve them

ReplyDeletenumerically):

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.

Pawel

.. you need to track T (theta) of course, so you can have

ReplyDeleteT = 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!

This comment has been removed by the author.

ReplyDelete[from James]Thanks for the feedback.

ReplyDeleteFor 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.

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

ReplyDelete