Mouse over the particle field above. Link to source files and editor.

let xvels = [];
let yvels = [];
let SCALE = 0.09;
let DIFFUSE = 0.1;
let MOUSE_MULTIPLIER = 4;
let SPEED = 5;
let PIXELS = 25;
let N = 1200;
let PARTICLE_SPEED = 16;
let RADIUSz = 70;

let WIDTH;
let HEIGHT;

let particles = [];
let colors = [];

function setup() {
    w = (int(800/25)+1)*25;
    h = (int(800/25)+1)*25;
    createCanvas(w, h);
  
    noStroke();
    smooth();
    
    WIDTH = width/PIXELS;
    HEIGHT = height/PIXELS;
    xvels = [];
    yvels = [];
    for ( x=0; x<width; x++) {
      xvels[x] = [];
      yvels[x] = [];
        for ( y=0; y<height; y++) {
            xvels[x][y] = 0;//SPEED*(noise(x*SCALE, y*SCALE)-0.4);
            yvels[x][y] = 0;//SPEED*(noise(x*SCALE, y*SCALE)-0.4);
        }
    }
     
    for (i=0; i<N; i++) {
        particles[i] = new p5.Vector(random(width), random(height), 0);
        colors[i] = color(random(255), random(255), random(255));
    }
  
}

let old_mousex;
let old_mousey;

function draw() {
  background(0);
  let mousex = int((mouseX-1) / PIXELS);
  let mousey = int((mouseY-1) / PIXELS);
     
    if (mousex == old_mousex && mousey == old_mousey) {
        old_mousex = 0;
        old_mousey = 0;
    } else {
        if (old_mousex > 0 && old_mousey > 0 &&
            mousex < WIDTH && mousey < HEIGHT) {
          xvels[mousex][mousey] = MOUSE_MULTIPLIER * max(-20, min(20, mousex - old_mousex));
          yvels[mousex][mousey] = MOUSE_MULTIPLIER * max(-20, min(20, mousey - old_mousey));
        }
        old_mousex = mousex;
        old_mousey = mousey;
    }
     
    update(0.1);
     
     
    for (let i=0; i<N; i++) {
        particle = particles[i];
        velx = sample_surrounding(xvels, particle.x/PIXELS, particle.y/PIXELS);
        vely = sample_surrounding(yvels, particle.x/PIXELS, particle.y/PIXELS);
        particle.x = (particle.x + velx*PARTICLE_SPEED + width) % width;
        particle.y = (particle.y + vely*PARTICLE_SPEED + height) % height;
        particle.z *= 0.95;
        particle.z = max(particle.z, pow((abs(velx) + abs(vely)),0.75)*RADIUSz);
        fill(colors[i]);
        ellipse(particle.x, particle.y, particle.z, particle.z);
    }
  
}

function sample(arr, x, y) {
    let w = WIDTH;
    let h = HEIGHT;
    return arr[(x + w) % w][(y + h) % h];
}

function sample_surrounding(arr,x,y) {
    let w = WIDTH-1;
    let h = HEIGHT-1;
    x = (x + 10*w) % w;
    y = (y + 10*h) % h;
    let fx = floor(x);
    let fy = floor(y);
    let cx = fx+1;
    let cy = fy+1;
    let wx1 = cx - x;
    let wx2 = x - fx;
    let wy1 = cy - y;
    let wy2 = y - fy;
    return  arr[fx][fy] * wx1*wy1 +
            arr[fx][cy] * wx1*wy2 +
            arr[cx][fy] * wx2*wy1 +
            arr[cx][cy] * wx2*wy2;
}

function diffuse(arr,dt,drag) {
    let ret = [];
    let wsurrounding = dt*DIFFUSE;
    let wcenter = drag-(wsurrounding*4);
    for (let x=0; x<WIDTH; x++) {
        ret[x] = [];
        for (let y=0; y<HEIGHT; y++) {
            ret[x][y] = wcenter*sample(arr, x, y) +
                          wsurrounding*sample(arr, x, y-1) +
                          wsurrounding*sample(arr, x, y+1) +
                          wsurrounding*sample(arr, x-1, y) +
                          wsurrounding*sample(arr, x+1, y);
        }
    }
    return ret;
}

function project(xarr,yarr,dt) {
    let proj = [];
     
    for (let x=0; x<WIDTH; x++) {
      proj[x] = [];
        for (let y=0; y<HEIGHT; y++) {
            proj[x][y] = -0.5 * ((sample(xarr, x+1, y) - sample(xarr, x-1, y)) +
                                 (sample(yarr, x, y+1) - sample(yarr, x, y-1)));
        }
    }
     
    proj = diffuse(proj, dt, 1);
     
    for (let x=0; x<WIDTH; x++) {
        for (let y=0; y<HEIGHT; y++) {
            xarr[x][y] += -0.5 * (sample(proj, x+1, y) - sample(proj, x-1, y));
            yarr[x][y] += -0.5 * (sample(proj, x, y+1) - sample(proj, x, y-1));
        }
    }
}

function update(dt) {
    xvels1 = [];
    yvels1 = [];
     
    xvels = diffuse(xvels, dt, 0.98);
    yvels = diffuse(yvels, dt, 0.98);
     
    project(xvels, yvels, dt);
     
    // advect
    for (let x=0; x<WIDTH; x++) {
      xvels1[x] = [];
      yvels1[x] = [];
        for (y=0; y<HEIGHT; y++) {
            xvels1[x][y] = sample_surrounding(xvels, x-xvels[x][y]*dt*SPEED, y-yvels[x][y]*dt*SPEED);
            yvels1[x][y] = sample_surrounding(yvels, x-xvels[x][y]*dt*SPEED, y-yvels[x][y]*dt*SPEED);
        }
    }
     
    project(xvels1, yvels1, dt);
     
    xvels = xvels1;
    yvels = yvels1;
}