Is solving Verlet constraints parallelizable?

Question

I am a computer science student and I came across the Verlet Cloth simulation. An implementation is here. I was interested in how such computation would lend itself to parallelization.

The problem looks like a variant of N-body computation where we could compute constraints on each point and them move the points based on the effect of the constraint. I tried along those lines and ended up with "logical" errors. I am sure my understanding of the physics is incomplete.

Here is what I tried. Assume a constraint C between A and B. Instead of trying to resolve C once and update both A and B, I tried resolving C twice once along with A's constraints and update just A and once along with B's constraints and just update B.

From this link, I am wondering whether the constraint resolution for each point can be done independently of the other, when each point is constrained by multiple other points.

dist = Math.abs(A - B)
diff = (dist - r) / dist
moveBy = 0.5 * diff * dist
A = A + moveBy
B = B - moveBy

Doing this for each constraint works. But the moment, I compute all these "moveBy" values for all constraints together and move the points by the sum of these values, things don't work. This makes me think that the order in which the points are resolved for constraints and the order in which constraints are applied matters.Is that right?

Answer 1

I'm not quite sure if this answers your question, or whether its more of a length comment. If I am completely off-base here, do let me know.

So suppose that A and B are of some type Point, and you call a method of Point, resolve(Point that), which is defined as

dist   = Math.abs(this.x - that.x)
diff   = (dist - r) / dist
moveBy = 0.5 * diff * dist
this.x = this.x + moveBy
that.x = that.x - moveBy

I'm not quite sure if I understood you correctly, but are you saying that you're calling this method for both A and B, i.e. A.resolve(B) and B.resolve(A)? Don't. moveBy is the same in both calls, but this and that apply it with different signs, so it would cancel out. This is to say that the order matters. Note that just calling A.resolve(B) would work.

Right, I wanted to cover all bases, and I don't think you actually did what I proposed above. Perhaps, instead, what you have is the following for resolve(Point that),

dist = Math.abs(this.x - that.x)
diff = (dist - r) / dist
this.moveBy = 0.5 * diff * dist

and it is only later in another function that you actually apply the this.x += this.moveBy. This is almost ok. Not quite, though. The reason is quite simple: Math.abs. Let's look at the equations of harmonic motion. Here's the energy of the bond between particles A and B: $$E = \frac{k}{2}(\|\mathbf{r}_\text{A} - \mathbf{r}_\text{B}\| - r_0)^2$$ Using Newton's equations of motion, $$m_\text{A}\ddot{\mathbf{r}}_\text{A} = - \frac{\partial E}{\partial \mathbf{r}_\text{A}} = -k(\|\mathbf{r}_\text{A} - \mathbf{r}_\text{B}\|-r_0)\frac{\mathbf{r}_\text{A} - \mathbf{r}_\text{B}}{\|\mathbf{r}_\text{A} - \mathbf{r}_\text{B}\|}$$

How does this compare with the code above? The code does not preserve directionality (because of the absolute value in the code), and thus all the forces are to one specific direction (positive $x$, $y$ axes).

If there are more particles, the equations don't really get any more complicated: on the left hand side we have the same thing as before, and on the right more force terms. This is to say that in the Eulerian integration scheme that you are most probably using, you can first compute all the forces, and from those, all the velocities, and from those, all the displacements. Computing the force for one particle ought to be completely independent of the computations done at some other particle (well, independent is a poor choice of words, as half the forces are the same, but with the opposite direction as per Newton's third law, but in the sense of parallel computation they do not affect one another).