I am trying to implement this paper for viscoelastic fluid simulation - http://www.ligum.umontreal.ca/Clavet-2005-PVFS/pvfs.pdf - to run on the GPU (using iOS Metal). The paper uses a particle-based method, with spring interactions between local particles to model elasticity. I am having trouble calculating the spring interactions efficiently. Pseudocode from the paper:
foreach neighbor pair i j, (i < j)
q ← ri j /h
if q < 1
if there is no spring i j
add spring i j with rest length h
// tolerable deformation = yield ratio * rest length d←γLij
if rij > L+d // stretch
Lij ←Lij+∆tα(rij−L−d)
else if rij < L−d // compress
Lij ←Lij−∆tα(L−d−rij) foreach spring i j
if Li j > h
remove spring i j
where i and j are particles; rij = |rij|, rij = positionj − positioni and h is the interaction radius; Lij is the rest length;
There are a variable number of springs, and each particle may have many, or no springs, depending on neighbouring particles. I cannot see an efficient way of storing or calculating the spring interactions on the GPU. If I have a buffer with, say, max. 10 springs allocated per particle with spring defined as follows:
struct Spring {
int i; // index of particle i
int j; // index of particle j
float rest_distance = 0.0;
};
Then I could loop through each particle's springs (with some kind of max. value written for no spring / removal of spring) but I think accessing, writing and removing the springs would be very slow and I would have to use atomic operations to sync writing springs across threads. The rest length would need to be preserved + updated on each frame / iteration.
Is there a better way to store and loop through a variable number of springs like this, or is this just not suited to optimisation on GPU?
Thank you for any help!