Equation for velocity verlet

[boxiXia]

Hi Jacob,

I was going through your code on the velocity verlet: https://github.com/jacobaustin123/Titan/blob/951733fc6caa1c40989664fbe7802a8c10bf9215/src/sim.cu#L1152-L1155 It was implemented differently than the one I found in Wikipedia and also in this page, specifically the next position is calcualted based on the previous velocity and acceleration; and in your code the position is updated base on the current velocity and acceleration.

I checked both implementation, your implementation seems to be more stable in larger time step compared to the Wikipedia implementation. Mathematically they are different, can you confirm that your implantation is also correct? Thanks!

[jacobaustin123]

Hi Boxi,

Yes, it's correct, and it agrees with the Wikipedia page (the second algorithm). Note that a{n+1} has to be calculated with respect to the new positions x{n+1}. All we're doing here is re-numbering the lines 1, 2, 3 -> 2, 3, 1. It doesn't work if the acceleration depends on the velocity, but it doesn't here.

[boxiXia]

Hi Jacob, Thanks for your reply! I agree that "a{n+1} has to be calculated with respect to the new positions x{n+1}". The part that I have question is the position update. Please correct me if I'm wrong, for position update, you use: X(t+dt) = V(t+dt)dt + 0.5 a(t+dt) dt^2 This is different form the Wikipedia: X(t+dt) = V(t) dt + 0.5 a(t) dt^2 Thus your version for position update use velocity and acceleration at step {n+1} instead of step {n}.

[jacobaustin123]

I think if you imagine renumber all the steps in the algorithm it makes sense. Step 2 becomes step 1, step 3 becomes step 2, step 1 becomes step 3. First we calculate a_{n+1}. Then velocity is updated using the new and old velocities, and finally the position is updated using the new velocities which are the old velocities under renumbering. Does that make sense? Because before step 1 was the first step, so everything is numbered a{n} even though it's actually a{n+1} because an iteration has passed.

[jacobaustin123]

On this note, I do want to ask your opinion about constraint handling. Right now contacts are handled by inserting stiff springs at contact points that apply repulsive forces. It's also possible to very carefully handle constraints by manually setting positions and velocities such that energy is conserved. For example, we can do a form of Verlet integration

mass.acc = mass.force / mass.m;
Vec temp = mass.pos;
mass.pos += mass.pos - mass.vel + mass.acc * pow(dt, 2);
mass.vel = temp;

if (mass.pos[2] < 0) { // handle a single contact plane at z = 0
    double t = - mass.vel[2] / (mass.pos[2] - mass.vel[2]);
    Vec new_pos = mass.vel + t * (mass.pos - mass.vel);

    Vec vel = mass.pos - mass.vel;
    vel[2] = - vel[2];
    mass.vel = new_pos - t * vel;
    mass.pos = mass.vel + vel;
}

that exactly handles elastic collisions with the ground and has better energy conservation and accuracy than the stiff spring approach. On the other hand, it's more complicated and may end up being slower and easier to make mistakes with.

[boxiXia]

I just check your velocity Verlet with Euler and with Störmer–Verlet. I simulated a model with 3k mass and 40k spring at dt=2e-5 s, spring dampling and contact damping set to 0, the energy deviation is normalized by the total energy at T=0,(implementation), here is the result : Your variation of Verlet seems to have a larger error compared to the Euler and Störmer–Verlet. Could you check energy conservation with a large lattice and see if this is the case?

About constraint handling: setting a very large normal force coefficient for contact constraints seems to help reducing the total energy deviation, the downside is a small time-step required. I do agree applying some form of Verlet integration may help gaining better accuracy, my question is whether this will interfere with the "normal" mass update?

[jacobaustin123]

So I noticed that deviation as well. I'm not 100% sure, but Störmer–Verlet is supposed to be more accurate than velocity Verlet because velocity is defined implicitly, so there is less numerical error. The trend does extend to larger lattices. The reason I didn't implement Störmer–Verlet is that atm we allow the user to set the velocity manually, which doesn't work for Störmer–Verlet.

As for constraints, we would have to restructure the way constraints are applied. Basically, this is how FleX works. Constraints are applied after/while position updates have been applied, based on position-based constraint solvers. That might be a better approach.

[jacobaustin123]

On the point of energy conservation for constraints, here's a comparison for a 20x20x20 bouncing cube using Verlet. One uses spring constraints and the other manually sets the position and velocity to avoid collisions.