Water column loses volume without increase in density with WCSPH

[efaulhaber]

I noticed that the initialized water column in dam break examples with WCSPH collapses without an increase in density, as if the particles are initialized too far apart for the corresponding density.
I reproduced that behavior with a simple "fluid at rest" example in 2D.
The fluid is initialized like this. I added the gray box for reference.

When running the simulation, the fluid collapses and comes to rest at about 80% of the initial volume, but without a corresponding increase in density:

I also noticed that I need a very small CFL number of 0.01 for stability. With CFL = 0.1, I get these little explosions in the fluid:

What is going on here? Is there something wrong with my setup?

{
	"Configuration":
	{
		"pause": true,
		"sim2D": true,
		"particleRadius": 0.005,
		"colorMapType": 1,
		"numberOfStepsPerRenderUpdate": 4,
		"density0": 1000,
		"simulationMethod": 0,
		"gravitation": [0,-9.81,0],
		"cflMethod": 1,
		"cflFactor": 0.01,
		"cflMaxTimeStepSize": 0.005,
		"maxIterations": 100,
		"maxError": 0.1,
		"maxIterationsV": 100,
		"maxErrorV": 0.1,
		"stiffness": 100000,
		"exponent": 7,
		"velocityUpdateMethod": 0,
		"enableDivergenceSolver": true,
		"boundaryHandlingMethod": 0,
		"enableVTKExport": true,
		"dataExportFPS": 50,
		"enableRigidBodyVTKExport": true,
		"particleAttributes": "velocity;density"
	},
	"Materials": [
		{
			"id": "Fluid",
			"surfaceTension": 0.0,
			"surfaceTensionMethod": 0,
			"viscosity": 0.02,
			"viscosityMethod": 1,
			"vorticityMethod": 0,
			"vorticity": 0.1,
			"viscosityOmega": 0.0,
			"inertiaInverse": 0.5,
			"colorMapType": 1
		}
	],
	"RigidBodies": [
		{
			"geometryFile": "../models/UnitBox.obj",
			"translation": [1,1.0,0],
			"rotationAxis": [1, 0, 0],
			"rotationAngle": 0,
			"scale": [2, 2, 1],
			"color": [0.1, 0.4, 0.6, 1.0],
			"isDynamic": false,
			"isWall": true,
			"mapInvert": true,
			"mapThickness": 0.0,
			"mapResolution": [30,30,20]
		}
	],
	"FluidBlocks": [
		{
			"denseMode": 0,
			"start": [0, 0.0, -1],
			"end": [1.99, 1.0, 1],
			"translation": [0.005, 0.005, 0.0],
			"scale": [1, 1, 1]
		}
	]
}

[janbender]

Weakly Compressible SPH (WCSPH) does not enforce incompressibility. It just tries to counteract compression by a force which depends on a stiffness parameter. Since the method is explicit, it is only conditionally stable. So if you want to have less compression using WCSPH, you have to increase the stiffness parameter. However, since it is an explicit method, it will tend to get unstable when using a high stiffness. In this case you have to take a small time step size to get stable results. Therefore, WCSPH is finally much slower than any implicit method for almost incompressible fluids (which has been shown in several publications). Btw. just setting the CFL factor to a small value does not help since a fluid which is not in motion will always take the largest possible time step size.

[efaulhaber]

Thanks a lot for the quick reply. That explains the instability part. But the 20% loss in volume can't be explained by the weak compressibility. I used a stiffness of 100,000. To my understanding, that should limit density fluctuations to about 1-2%. But that doesn't explain why the fluid loses 20% of volume without a corresponding change in density (I'd expect the density to be around 1250 then).

I also just checked the incompressible schemes as well. With these, the fluid also collapses to about 80% of the initial volume. Why is that? The denseMode can't be the solution, right? There must be a way to initialize a fluid with the desired volume in a rectangular grid, like with denseMode 0, right?

[janbender]

One explanation for the volume loss is of course that the regular samping at the beginning is not dense at all. This means in the initial configuration the particles have already a density of less then the rest density. Moreover, I just remind that long ago we reduced the initial density by a factor of 0.8 because otherwise the denseMode 2 leads to an initial positive density and therefore to pressure forces. I will put it on my todo list to remove this factor since it is misleading. You can do it yourself in your code:

SPlisHSPlasH/SPlisHSPlasH/FluidModel.cpp

Line 258 in a081124

m_V = static_cast<Real>(0.8) * diam*diam;

Just remove the 0.8.

[efaulhaber]

Thanks a lot, that worked! Now I get the expected volume and density.