I am in the process of implementing the math of the following paper (Position Based Fluid) with java, but have still trouble following the paper and are confused about some formulas. And yes i dont have any math professional backgrunds at all, so calculus, differential equations are out of my league but still i really want to implement this paper - my entire game concept is based on fluid dynamics, but there are no existing physics engine which fits my requirements.

Here some part of formulas from the right side of page 2:

Can you please help me to break this down - so that a non math-professionals understand it?

Especially in equation (9) what is the sum of k? Sum over all neighbor particles where every neightbor index is called k and the sum is incremented for every result of the formula "gradient p k c i"?

The two vertical lines between some symbols mean scalar right?

The symbol "pk" after the gradient and before the W confuse me a lot.

I already have implemented the density estimator (p i) and the density constraint (C i) and have all 2 kernels (poly 6 and spiky gradient).

Of course there is much more, but i want to do it step by step.

Thanks in regards,

Final

Its a shame that i never had the chance to study calculus and all that math stuff.

K i think i implemented it - but it explodes immediatly when i start it...

Is that correct implemented?:

private float Wpoly6 (Vec2f r, float kH) {
	float rlen = r.length();
	if (rlen > kH || rlen == 0) {
		return 0;
	}
	final float tmp = kH * kH - rlen * rlen;
	final float term = 315.0f / (64.0f * (float)Math.PI * (float)Math.pow(kH, 9));
	return term * tmp * tmp * tmp;
}

private Vec2f gradWspiky (Vec2f r, float kH) {
	float rlen = r.length();
	if (rlen > kH || rlen == 0) {
		return new Vec2f (0f, 0f);
	}
	final float tmp = kH - rlen;
	final float term = 45.0f / ((float)Math.PI * (float)Math.pow(kH, 6) * rlen);
	return r.mult(-term * tmp * tmp);
}

private Vec2f calculateCiGradient(Vec2f position_i, int index) {
	Vec2f accum = new Vec2f(0.0f, 0.0f);
	for (int j : neighborSet[index]) {
		final Vec2f position_j = predictedPositionBuffer[j];
		final Vec2f distance = position_i.sub(position_j);
		accum = accum.add(gradWspiky(distance, kH));
	}
	return accum.mult(invRestDensity);
}

private void calculateLamda(int index) {
	final Vec2f position_i = predictedPositionBuffer[index];

	// Calculate standard SPH density estimator (equation 2)
	float p_i = 0.0f;
	for (int j : neighborSet[index]) {
		final Vec2f position_j = predictedPositionBuffer[j];
		final Vec2f distance = position_i.sub(position_j);
		p_i += Wpoly6(distance, kH);
	}
	
	// Calculate density constraint (equation 1)
	float C_i = (p_i / restDensity) - 1f;
	
	// Calculate C_i gradient sum (equation 9)
	float C_i_gradient, sum_C_i_gradient = 0.0f;
	for (int j : neighborSet[index]) {
		final Vec2f position_j = predictedPositionBuffer[j];
		final Vec2f distance = position_i.sub(position_j);
		
		// Calculate gradient when k = j
		C_i_gradient = gradWspiky(distance, kH).mult(-invRestDensity).length();
		sum_C_i_gradient += C_i_gradient * C_i_gradient;
	}
	
	// Add gradient when k = i (equation 8)
	C_i_gradient = calculateCiGradient(position_i, index).length();
	sum_C_i_gradient += C_i_gradient * C_i_gradient;

	// Calculate lamda (equation 11)
	float sum_C_i = sum_C_i_gradient + epsilon;
	synchronized (lamdaLock) {
		lambdaBuffer[index] = -1.0f / (C_i / sum_C_i);
	}
}

private void calculateDeltaPos(int index) {
	final Vec2f position_i = predictedPositionBuffer[index];
	
	final float lambda_i = lambdaBuffer[index];

	Vec2f delta = new Vec2f(0.0f, 0.0f);
	
	// Calculate delta q (some point inside the smoothing kernel - equation 13)
	final Vec2f d_q = new Vec2f(kH * 0.75f).mult(deltaQ).add(position_i);
	
	// Calculate divider term for s_corr (equation 13)
	final Vec2f d_q_distance = position_i.sub(d_q);
	final float d_q_term = Wpoly6(d_q_distance, kH);
	
	// Calculate delta p i (equation 14)
	for (int j : neighborSet[index]) {
		Vec2f position_j = predictedPositionBuffer[j];
		Vec2f distance = position_i.sub(position_j);

		float k_term = Wpoly6(distance, kH) / d_q_term;
		float s_corr = -pressureK * (float)Math.pow(k_term, pressureN);
		float lambda_j = lambdaBuffer[j];
		Vec2f gradient = gradWspiky(distance, kH);
		delta = delta.add(gradient.mult(lambda_i + lambda_j + s_corr));
	}
	synchronized (posDeltaLock) {
		posDeltaBuffer[index].set(delta.mult(invRestDensity));
	}
}

Here are my properties:

particleRadius = 0.1f;
kH = particleRadius * 4f;
restDensity = 1000f;
invRestDensity = 1.0f / restDensity;
epsilon = 0.01f;
deltaQ = 0.1f * kH;
pressureK = 0.1f;
pressureN = 4;

I still havent found the issue, implementation looks right but which confuses me a lot are the gradient symbol on the right side in equations (7, 8).

To debug it better i added a single step mode to see what causes the instability. The simulation starts with just 4 particles, after the first step it gets a very small delta going towards the center of the origin where the particles are aligned, then on the next step the delta increases much more and after the third step the particles starts to going in the opposite direction.

Also i double checked the neighbor search which looked fine, but for the sake of simplicity i changed it to a O^N method.

Please help me to get this right.

This topic is not closed yet, but i still hope that someone can help with that to get this working correctly.

But in the meantime - i just implemented the "Particle-based Viscoelastic Fluid Simulation" paper by Clavet in just 4 hours and i have a fully working simulation right now.

There are still stuff missing - like the spring calculations and the stickiness thingy which i completely left it out but these will come later.

I found the issue on PBF - it was just some parameter tweaking. So my math was correct :-)

Oh and another things i fixed, the poly6 and spiky kernel formulas was based on 3D vectors - therefore i changed it to 2D instead.

Thread can be closed.