# Is solving Verlet constraints parallelizable?

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?

• Would Computational Science be a better home for this question? – Qmechanic Dec 9 '14 at 11:40
• In general, you have to deal with all the constraints simultaneously. See en.wikipedia.org/wiki/Constraint_algorithm for details. Two methods that are particularly amenable to parallelization are M-SHAKE and P-LINCS. – Juan M. Bello-Rivas Dec 9 '14 at 13:40
• @JuanM.BelloRivas Many of the algorithms have limitations, particularly when connectivity is high or there are loops. As I understand it, the cloth simulation is little else but loops: Do M-SHAKE and P-LINCS handle these gracefully? I should also say that having run the code (didn't look at the source, though), it seems that the constraints are just harmonic bonds i.e. the meaning is different from constraints in molecular simulations (which are stiff, usually through Lagrange multipliers). – alarge Dec 9 '14 at 17:58
• @Juan M Bello Rivas, what does the simultaneously mean. Because, at least from the code, each constraint alters the position of points and this constrained position is what the next constraint applies on. Or is there something I am missing. – Sruti Dec 9 '14 at 20:25
• @alarge SHAKE (and M-SHAKE) do handle loops and I don't think it would be a problem for LINCS. I see your point, though. What the OP is referring to is commonly called a restraint not a constraint in the MD literature. – Juan M. Bello-Rivas Dec 10 '14 at 0:35

## 1 Answer

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