trying to simulate fluids

gistfile1.txt

var defaultCellSize = 48 * window.devicePixelRatio;

function createParticles(canvas, width, height, cellSize = defaultCellSize)
{
	canvas.width  = width;
	canvas.height = height;

	// defaults

	canvas.showCells      = false;
	canvas.showVelocity   = false;
	canvas.showDensity    = false;
	canvas.showForces     = true;
                          
	canvas.cellSize       = cellSize;
	canvas.xCells         = Math.ceil(canvas.width  / canvas.cellSize);
	canvas.yCells         = Math.ceil(canvas.height / canvas.cellSize);
	
	canvas.scale          = 10;

    //canvas.particleBounce = 1;
	canvas.wallBounce     = 0.2;
    canvas.damping        = 0;

	canvas.random         = 0;
	canvas.gravity        = 9.8;

	canvas.viscositySigma = 0;     // linear viscosity
	canvas.viscosityBeta  = 0.001; // quadratic viscosity
	
	canvas.restDensity    = 10;

	canvas.stiffness      = 1.5;
	canvas.nearStiffness  = 3;

	canvas.time           = performance.now();
	canvas.dt             = 1;

	canvas.mouseDown      = false;
	canvas.mouseX         = 0;
	canvas.mouseY         = 0;

    // cells

	canvas.cells = new Array();
	
	for (var y = 0; y < canvas.yCells; y++)
	{
		for (var x = 0; x < canvas.xCells; x++)
		{
			canvas.cells.push({ 
				particles: new Array(), 
				cx:        x,
				cy:        y,
				left:      x * canvas.cellSize, 
				top:       y * canvas.cellSize, 
				right:     (x+1) * canvas.cellSize, 
				bottom:    (y+1) * canvas.cellSize,
				color:     'rgb(' + Math.floor(64+192*Math.random()) + ', '
							      + Math.floor(64+192*Math.random()) + ', '
							      + Math.floor(64+192*Math.random()) + ')'
			});
		}
	}


	function saveCells()
	{
	    //canvas.oldCells = canvas.cells;
	    canvas.oldCells = new Array();

	    for (var i = 0; i < canvas.cells.length; i++)
	    {
	        var cell = canvas.cells[i];

	        var oldCell = { 
	            particles: new Array(), 
	            cx:        cell.cx,
	            cy:        cell.cy,
	            left:      cell.left, 
	            top:       cell.top, 
	            right:     cell.right, 
	            bottom:    cell.bottom,
	            color:     cell.color
	        };

	        for (var j = 0; j < cell.particles.length; j++)
	        {
	            var p = cell.particles[j];

	            oldCell.particles.push({
	                cell:         oldCell,
	                px:           p.px, 
	                py:           p.py,
	                ppx:          p.ppx,
	                ppy:          p.ppy,
	                fx:           p.fx,
	                fy:           p.fx,
	                vx:           p.vx,
	                vy:           p.vx,
	                color:        p.color,
	                density:      p.density,
	                nearDensity:  p.nearDensity
	            });
	        }

	        canvas.oldCells.push(oldCell);
	    }
	}

	
    function updateCells()
	{
	    for (var i = 0; i < canvas.cells.length; i++)
	    {
	        var cell = canvas.cells[i];

	        for (var j = cell.particles.length - 1; j >= 0; j--)
	        {
	            var p = cell.particles[j];

	            if (   p.px >= cell.left
					&& p.py >= cell.top
					&& p.px <  cell.right
					&& p.py <  cell.bottom
                    || p.px < 0 || p.py < 0)
	                continue;

	            cell.particles.splice(j, 1);
	            p.cell = null;

	            addParticle(p);
	        }
	    }
	};

    // particles

	function addParticle(p)
	{
		var cx = Math.floor(p.px / canvas.cellSize);
		var cy = Math.floor(p.py / canvas.cellSize);
 
		if (   cx < 0 
			|| cy < 0 
			|| cx >= canvas.xCells 
			|| cy >= canvas.yCells)
			return;

		var cell = canvas.cells[cy * canvas.xCells + cx];
			
		p.cell = cell;
		cell.particles.push(p);
	}

	var radius = 2;

    // update

	canvas.update = function()
	{
	    var time    = performance.now();
	    canvas.dt   = (time - canvas.time) / 1000 * canvas.scale; // in seconds
	    canvas.time = time;
	    
	    canvas.updateParticles();
		canvas.paint();
	};

	canvas.updateParticles = function() 
	{
	    canvas.applyRandom();
	    canvas.applyDamping();
	    canvas.applyGravity();

	    canvas.applyViscosity();

	    canvas.saveParticlePositions();

	    canvas.moveParticles();
	    //updateCells();

	    canvas.relaxParticleDensity();
	    //updateCells();

	    canvas.bounceOffWalls();
	    updateCells(); // housekeeping

	    canvas.updateParticleVelocities();
	};

    // move

	canvas.saveParticlePositions = function()
	{
		for (var i = 0; i < canvas.cells.length; i++)
		{
			var cell = canvas.cells[i];

			for (var j = 0; j < cell.particles.length; j++)
			{
				var p = cell.particles[j];

				p.ppx = p.px;
				p.ppy = p.py;
			}
		}
	};



	canvas.moveParticles = function()
	{
	    for (var i = 0; i < canvas.cells.length; i++)
	    {
	        var cell = canvas.cells[i];

	        for (var j = 0; j < cell.particles.length; j++)
	        {
	            var p = cell.particles[j];

                p.px += p.vx * canvas.dt,
                p.py += p.vy * canvas.dt;
	        }
	    }
	};


    // forces

	canvas.applyRandom = function()
	{
		for (var i = 0; i < canvas.cells.length; i++)
		{
			var cell = canvas.cells[i];

			for (var j = 0; j < cell.particles.length; j++)
			{
				var p = cell.particles[j];

				p.vx += (-1 + Math.random()*2) * canvas.random;
				p.vy += (-1 + Math.random()*2) * canvas.random;
			}
		}	
	};



	canvas.applyDamping = function()
	{
		for (var i = 0; i < canvas.cells.length; i++)
		{
			var cell = canvas.cells[i];

			for (var j = 0; j < cell.particles.length; j++)
			{
				var p = cell.particles[j];

				p.vx *= 1 - canvas.damping;
				p.vy *= 1 - canvas.damping;
			}
		}	
	};

	canvas.applyGravity = function()
	{
		for (var i = 0; i < canvas.cells.length; i++)
		{
			var cell = canvas.cells[i];

			for (var j = 0; j < cell.particles.length; j++)
				cell.particles[j].vy += canvas.gravity * canvas.dt;
		}	
	};

    // viscosity

	canvas.applyViscosity = function()
	{
	    var h = canvas.cellSize;

	    for (var i = 0; i < canvas.cells.length; i++)
	    {
	        var cell = canvas.cells[i];

	        var cyMin = Math.max(0, cell.cy - 1);
	        var cyMax = Math.min(cell.cy + 1, canvas.yCells - 1);
	        var cxMin = Math.max(0, cell.cx - 1);
	        var cxMax = Math.min(cell.cx + 1, canvas.xCells - 1);

	        for (var j = 0; j < cell.particles.length; j++)
	        {
	            var p = cell.particles[j];

	            for (var cy = cyMin; cy <= cyMax; cy++)
	            {
	                for (var cx = cxMin; cx <= cxMax; cx++)
	                {
	                    var ci = cy * canvas.xCells + cx;
	                    var checkCell = canvas.cells[ci];

	                    for (var k = 0; k < checkCell.particles.length; k++)
	                    {
	                        if (i == ci && j == k) // check if the test is against itself
	                            continue;

	                        var p2 = checkCell.particles[k];

	                        var dx = p2.px - p.px;
	                        var dy = p2.py - p.py;

	                        var dist2 = dx*dx + dy*dy;

	                        if (dist2 >= h*h) // check if particle is too far
	                            continue;

	                        var dist = Math.sqrt(dist2);
	                        var q = 1 - dist/h;

	                        var dvx = p.vx - p2.vx;
	                        var dvy = p.vy - p2.vy;

	                        var x1 = dx / nozero(dist);
	                        var y1 = dy / nozero(dist);

	                        var u = dvx*x1 + dvy*y1; // dot product

	                        if (u > 0) 
	                        {
	                            var I = 
                                      canvas.dt 
                                    * q
                                    * (  canvas.viscositySigma * u 
                                       + canvas.viscosityBeta  * u*u);

	                            var Ix = I * x1;
	                            var Iy = I * y1;

	                            p.vx -= Ix/2;
	                            p.vy -= Iy/2;
                                
	                            p2.vx += Ix/2;
	                            p2.vy += Iy/2;
	                        }
	                    }
	                }
	            }
	        }
	    }
	};

    // double density relaxation

	canvas.relaxParticleDensity = function()
	{
	    saveCells();

	    //

	    var h = canvas.cellSize;

	    for (var i = 0; i < canvas.cells.length; i++)
	    {
	        var oldCell = canvas.oldCells[i];
	        var cell    = canvas.cells[i];

	        var cyMin = Math.max(0, cell.cy - 1);
	        var cyMax = Math.min(cell.cy + 1, canvas.yCells - 1);
	        var cxMin = Math.max(0, cell.cx - 1);
	        var cxMax = Math.min(cell.cx + 1, canvas.xCells - 1);

	        for (var j = 0; j < cell.particles.length; j++)
	        {
	            var p = cell.particles[j];

	            p.density     = 0;
	            p.nearDensity = 0;

	            for (var cy = cyMin; cy <= cyMax; cy++)
	            {
	                for (var cx = cxMin; cx <= cxMax; cx++)
	                {
	                    var ci = cy * canvas.xCells + cx;

	                    var checkCell = canvas.cells[ci];

	                    for (var k = 0; k < checkCell.particles.length; k++)
	                    {
	                        if (i == ci && j == k) // check if the test is against itself
	                            continue;

	                        var p2 = checkCell.particles[k];

	                        var dx = p2.px - p.px;
	                        var dy = p2.py - p.py;

	                        var dist2 = dx*dx + dy*dy;

	                        if (dist2 >= h*h) // check if particle is too far
	                            continue;

	                        var dist = Math.sqrt(dist2);
	                        var q = 1 - dist/h;

	                        p.density     += Math.sqr (q);
	                        p.nearDensity += Math.cube(q);
	                    }
	                }
	            }

	            var pressure     = canvas.stiffness     * (p.density - canvas.restDensity);
	            var nearPressure = canvas.nearStiffness * p.nearDensity;

	            p.fx = 0;
	            p.fy = 0;

	            var pOld = oldCell.particles[j];

	            for (var cy = cyMin; cy <= cyMax; cy++)
	            {
	                for (var cx = cxMin; cx <= cxMax; cx++)
	                {
	                    var ci = cy * canvas.xCells + cx;

	                    var checkOldCell = canvas.oldCells[ci];
	                    var checkCell    = canvas.cells[ci];

	                    for (var k = 0; k < checkOldCell.particles.length; k++)
	                    {
	                        if (i == ci && j == k) // check if the test is against itself
	                            continue;

	                        var p2    = checkCell.particles[k];
	                        var p2old = checkOldCell.particles[k];

	                        var dx = p2old.px - pOld.px;
	                        var dy = p2old.py - pOld.py;

	                        var dist2 = dx*dx + dy*dy;

	                        if (dist2 >= h*h) // check if particle is too far
	                            continue;

	                        var dist = Math.sqrt(dist2);
	                        var q = 1 - dist/h;

	                        var D = 
	                              Math.sqr(canvas.dt)
	                            * (  pressure     * q
	                               + nearPressure * q*q);

	                        p2old.fx = D/2 * dx / nozero(dist); // multiplied by unit vector
                            p2old.fy = D/2 * dy / nozero(dist);

	                        p2.px += p2old.fx;
	                        p2.py += p2old.fy;

	                        p.fx -= p2old.fx;
	                        p.fy -= p2old.fy;
                        }
                    }
                }

                p.px += p.fx;
                p.py += p.fy;
            }
	    }
	};

    // collisions

	function dist(x1, y1, x2, y2)
	{
	    var dx = x2 - x1;
	    var dy = y2 - y1;

	    return Math.sqrt(dx*dx + dy*dy);
	}

	function dist2(x1, y1, x2, y2)
	{
	    var dx = x2 - x1;
	    var dy = y2 - y1;

	    return dx*dx + dy*dy;
	}

	canvas.bounceOffWalls = function()
	{
	    var left   = 0;
	    var top    = 0;
	    var right  = canvas.width;
	    var bottom = canvas.height;

	    for (var i = 0; i < canvas.cells.length; i++)
	    {
	        var cell = canvas.cells[i];

	        for (var j = 0; j < cell.particles.length; j++)
	        {
	            var p = cell.particles[j];

	            if (p.px < left)
	            {
	                var iy = p.ppy + (p.py - p.ppy) * (left - p.ppx) / nozero(p.px - p.ppx);
	                var ix = left;

	                var d = dist(ix, iy, p.px, p.py);
	                var v = dist(p.ppx, p.ppy, p.px, p.py);

	                var f = canvas.wallBounce * d / nozero(v);

	                p.px  = 2 * left - p.px;
	                p.ppx = 2 * left - p.ppx;

	                p.px = ix + (p.px - ix) * f;
	                p.py = iy + (p.py - iy) * f;

	                p.ppx = ix + (p.ppx - ix) * f;
	                p.ppy = iy + (p.ppy - iy) * f;
                }
	            else if (p.px > right) 
	            {
	                var iy = p.ppy + (p.py - p.ppy) * (right - p.ppx) / nozero(p.px - p.ppx);
	                var ix = right;

	                var d = dist(ix, iy, p.px, p.py);
	                var v = dist(p.ppx, p.ppy, p.px, p.py);

	                var f = canvas.wallBounce * d / nozero(v);

	                p.px  = 2 * right - p.px;
	                p.ppx = 2 * right - p.ppx;

	                p.px = ix + (p.px - ix) * f;
	                p.py = iy + (p.py - iy) * f;

	                p.ppx = ix + (p.ppx - ix) * f;
	                p.ppy = iy + (p.ppy - iy) * f;
                }

	            if (p.py < top)
	            {
	                var ix = p.ppx + (p.px - p.ppx) * (top - p.ppy) / nozero(p.py - p.ppy);
	                var iy = top;

	                var d = dist(ix, iy, p.px, p.py);
	                var v = dist(p.ppx, p.ppy, p.px, p.py);

	                var f = canvas.wallBounce * d / nozero(v);

	                p.py  = 2 * top - p.py;
	                p.ppy = 2 * top - p.ppy;

	                p.px = ix + (p.px - ix) * f;
	                p.py = iy + (p.py - iy) * f;

	                p.ppx = ix + (p.ppx - ix) * f;
	                p.ppy = iy + (p.ppy - iy) * f;
                }
	            else if (p.py > bottom)
	            {
	                var ix = p.ppx + (p.px - p.ppx) * (bottom - p.ppy) / nozero(p.py - p.ppy);
	                var iy = bottom;

	                var d = dist(ix, iy, p.px, p.py);
	                var v = dist(p.ppx, p.ppy, p.px, p.py);

	                var f = canvas.wallBounce * d / nozero(v);

	                p.py  = 2 * bottom - p.py;
	                p.ppy = 2 * bottom - p.ppy;

	                p.px = ix + (p.px - ix) * f;
	                p.py = iy + (p.py - iy) * f;

	                p.ppx = ix + (p.ppx - ix) * f;
	                p.ppy = iy + (p.ppy - iy) * f;
	            }
	        }
	    }	
    };


	canvas.updateParticleVelocities = function()
	{
	    for (var i = 0; i < canvas.cells.length; i++)
	    {
	        var cell = canvas.cells[i];

	        for (var j = 0; j < cell.particles.length; j++)
	        {
	            var p = cell.particles[j];

	            p.vx = (p.px - p.ppx) / nozero(canvas.dt);
	            p.vy = (p.py - p.ppy) / nozero(canvas.dt);
	        }
	    }
	};
    
    // paint

	canvas.paint = function()
	{
		c = canvas.getContext('2d');

		// background

		c.fillStyle = '#000411';
		c.fillRect(0, 0, canvas.width, canvas.height);

		// particles

		for (var i = 0; i < canvas.cells.length; i++)
		{
			var cell = canvas.cells[i];

			for (var j = 0; j < cell.particles.length; j++)
			{
				var p = cell.particles[j];

				// body 

				if (   !canvas.showVelocity
                    && !canvas.showDensity
                    && !canvas.showForces)
				{
				    c.beginPath();
				    c.arc(p.px, p.py, radius, 0, Tau, false);
				    c.fillStyle = canvas.showCells ? cell.color : p.color;

				    //if (p.colliding && canvas.showCollisions) 
				    //    c.fillStyle = '#ff4';

				    c.fill();
				}

				if (canvas.showVelocity)
				{
					var angle = getAngle(
						p.px, 
						p.py, 
						p.px + p.vx, 
						p.py + p.vy);

					var l = Math.sqrt(p.vx*p.vx + p.vy*p.vy);

					var df = 60;
					l = Math.pow(l / df, 0.75) * df;

					c.beginPath();
					c.moveTo(
                        p.px + l * Math.cos(angle) / 2, 
                        p.py + l * Math.sin(angle) / 2);
					c.lineTo(
						p.px - l * Math.cos(angle) / 2,
						p.py - l * Math.sin(angle) / 2);
					c.lineWidth = 1;
					c.strokeStyle = '#0f0';
					c.stroke();
				}

				else if (canvas.showDensity)
				{
				    var h = canvas.cellSize;

				    c.beginPath();
				    c.arc(p.px, p.py, h/20 * p.density / canvas.restDensity, 0, Tau, false);
				    c.lineWidth = 1;//Math.cube(p.density / canvas.restDensity) * 2;
			        c.strokeStyle = '#f44';
				    c.stroke();
				}

				else if (canvas.showForces)
				{
				    var angle = getAngle(
						p.px, 
						p.py, 
						p.px + p.fx, 
						p.py + p.fy);

				    var l = dist(0, 0, p.fx, p.fy);

				    //var df = 60;
				    //l = Math.pow(l / df, 0.75) * df;
				    l *= 2 * canvas.scale;

				    c.beginPath();
				    c.moveTo(
                        p.px + l * Math.cos(angle), 
                        p.py + l * Math.sin(angle));
				    c.lineTo(
						p.px - l * Math.cos(angle),
						p.py - l * Math.sin(angle));
				    c.lineWidth = 1;

				    var scale = 10;
				    var pressure = p.density - canvas.restDensity;

				    //c.strokeStyle = 'rgb('
                    //    + Math.min(Math.round(pressure / canvas.restDensity * 255 * scale), 255).toString()
                    //    + ', 0, '
                    //    + Math.min(Math.round(Math.abs(-pressure / canvas.restDensity) * 255 * scale), 255).toString()
                    //    + ')';

				    c.strokeStyle = '#ff0';
				    c.stroke();
				}
            }
		}	

		// cells

		if (canvas.showCells)
		{
			c.lineWidth    = 1;

			//c.font         = '16px Arial bold, sans-serif';
			//c.textAlign    = 'left';
			//c.textBaseline = 'top';

			for (var i = 0; i < canvas.cells.length; i++)
			{
				var cell = canvas.cells[i];

				if (cell.particles.length > 0)
				{
				    c.strokeStyle = cell.color;
				    c.fillStyle   = cell.color;

				    c.strokeRect(
					    0.5 + cell.left,
					    0.5 + cell.top, 
					    cell.right  - cell.left - 1,
					    cell.bottom - cell.top  - 1);

				    //c.fillText(
                    //    cell.particles.length,
                    //    cell.left + 2,
                    //    cell.top + 2);
                }
			}	
		}
	};
}