Libgdx port of Coding Challenge #132: Fluid Simulation : https://www.youtube.com/watch?v=alhpH6ECFvQ

Fluid.java

import com.badlogic.gdx.graphics.Color;
import com.badlogic.gdx.graphics.glutils.ShapeRenderer;
import com.badlogic.gdx.graphics.glutils.ShapeRenderer.ShapeType;
import com.badlogic.gdx.math.MathUtils;

public class Fluid {
	
	public static final int N = 128 / 1;
	public static final int SCALE = 1;
	public static final int iter = 4; //16 / 8;

	
	public static int IX(int x, int y) {
		  x = MathUtils.clamp(x, 0, N-1);
		  y = MathUtils.clamp(y, 0, N-1);
		  return x + (y * N);
		}
	
	  int size;
	  float dt;
	  float diff;
	  float visc;

	  float[] s;
	  float[] density;

	  float[] Vx;
	  float[] Vy;

	  float[] Vx0;
	  float[] Vy0;
	  
	  public Fluid(float dt, float diffusion, float viscosity) {

	    this.size = N;
	    this.dt = dt;
	    this.diff = diffusion;
	    this.visc = viscosity;

	    this.s = new float[N*N];
	    this.density = new float[N*N];

	    this.Vx = new float[N*N];
	    this.Vy = new float[N*N];

	    this.Vx0 = new float[N*N];
	    this.Vy0 = new float[N*N];
	  }

	  public void step() {
	    int N          = this.size;
	    float visc     = this.visc;
	    float diff     = this.diff;
	    float dt       = this.dt;
	    float[] Vx      = this.Vx;
	    float[] Vy      = this.Vy;
	    float[] Vx0     = this.Vx0;
	    float[] Vy0     = this.Vy0;
	    float[] s       = this.s;
	    float[] density = this.density;

	    diffuse(1, Vx0, Vx, visc, dt);
	    diffuse(2, Vy0, Vy, visc, dt);

	    project(Vx0, Vy0, Vx, Vy);

	    advect(1, Vx, Vx0, Vx0, Vy0, dt);
	    advect(2, Vy, Vy0, Vx0, Vy0, dt);

	    project(Vx, Vy, Vx0, Vy0);

	    diffuse(0, s, density, diff, dt);
	    advect(0, density, s, Vx, Vy, dt);
	  }

	  public void addDensity(int x, int y, float amount) {
	    int index = IX(x, y);
	    this.density[index] += amount;
	  }

	  public void addVelocity(int x, int y, float amountX, float amountY) {
	    int index = IX(x, y);
	    this.Vx[index] += amountX;
	    this.Vy[index] += amountY;
	  }

	  public void renderD(ShapeRenderer renderer) {
		  renderer.begin(ShapeType.Filled);

	    for (int i = 0; i < N; i++) {
	      for (int j = 0; j < N; j++) {
	        float x = i * SCALE;
	        float y = j * SCALE;
	        float d = this.density[IX(i, j)];
//	        renderer.setColor((((d * 8) % 256f) / 255f) * .3f,0,d/ 255f, 1f);
	        renderer.setColor((d/255f) ,(d/255f),(d/255f), 1f);
	        renderer.rect(x, y, SCALE, SCALE);
	      }
	    }
	    renderer.end();
	  }

	  public void renderV(ShapeRenderer renderer) {
		  renderer.setColor(Color.WHITE);
		  renderer.begin(ShapeType.Line);
		  
	    for (int i = 0; i < N; i++) {
	      for (int j = 0; j < N; j++) {
	        float x = i * SCALE;
	        float y = j * SCALE;
	        float vx = this.Vx[IX(i, j)];
	        float vy = this.Vy[IX(i, j)];

	        if (!(Math.abs(vx) < 0.1 && Math.abs(vy) <= 0.1)) {
	        	renderer.line(x, y, x+vx*SCALE, y+vy*SCALE );
	        }
	      }
	    }
	    renderer.end();
	  }

	  public void fadeD() {
	    for (int i = 0; i < this.density.length; i++) {
	      float d = density[i];
	      density[i] = MathUtils.clamp(d-0.02f, 0, 255);
	    }
	  }
	  
	  public static void diffuse (int b, float[] x, float[] x0, float diff, float dt) {
		  float a = dt * diff * (N - 2) * (N - 2);
		  lin_solve(b, x, x0, a, 1 + 6 * a);
		}

	public static void lin_solve(int b, float[] x, float[] x0, float a, float c) {
		  float cRecip = 1.0f / c;
		  for (int k = 0; k < iter; k++) {
		    for (int j = 1; j < N - 1; j++) {
		      for (int i = 1; i < N - 1; i++) {
		        x[IX(i, j)] =
		          (x0[IX(i, j)]
		          + a*(    x[IX(i+1, j)]
		          +x[IX(i-1, j)]
		          +x[IX(i, j+1)]
		          +x[IX(i, j-1)]
		          )) * cRecip;
		      }
		    }

		    set_bnd(b, x);
		  }
		}
		public static void project(float[] velocX, float[] velocY, float[] p, float[] div) {
		  for (int j = 1; j < N - 1; j++) {
		    for (int i = 1; i < N - 1; i++) {
		      div[IX(i, j)] = -0.5f*(
		        velocX[IX(i+1, j)]
		        -velocX[IX(i-1, j)]
		        +velocY[IX(i, j+1)]
		        -velocY[IX(i, j-1)]
		        )/N;
		      p[IX(i, j)] = 0;
		    }
		  }

		  set_bnd(0, div); 
		  set_bnd(0, p);
		  lin_solve(0, p, div, 1, 6);

		  for (int j = 1; j < N - 1; j++) {
		    for (int i = 1; i < N - 1; i++) {
		      velocX[IX(i, j)] -= 0.5f * (  p[IX(i+1, j)]
		        -p[IX(i-1, j)]) * N;
		      velocY[IX(i, j)] -= 0.5f * (  p[IX(i, j+1)]
		        -p[IX(i, j-1)]) * N;
		    }
		  }
		  set_bnd(1, velocX);
		  set_bnd(2, velocY);
		}


		public static void advect(int b, float[] d, float[] d0, float[] velocX, float[] velocY, float dt) {
		  float i0, i1, j0, j1;

		  float dtx = dt * (N - 2);
		  float dty = dt * (N - 2);

		  float s0, s1, t0, t1;
		  float tmp1, tmp2, x, y;

		  float Nfloat = N;
		  float ifloat, jfloat;
		  int i, j;

		  for (j = 1, jfloat = 1; j < N - 1; j++, jfloat++) { 
		    for (i = 1, ifloat = 1; i < N - 1; i++, ifloat++) {
		      tmp1 = dtx * velocX[IX(i, j)];
		      tmp2 = dty * velocY[IX(i, j)];
		      x    = ifloat - tmp1; 
		      y    = jfloat - tmp2;

		      if (x < 0.5f) x = 0.5f; 
		      if (x > Nfloat + 0.5f) x = Nfloat + 0.5f; 
		      i0 = MathUtils.floor(x); 
		      i1 = i0 + 1.0f;
		      if (y < 0.5f) y = 0.5f; 
		      if (y > Nfloat + 0.5f) y = Nfloat + 0.5f; 
		      j0 = MathUtils.floor(y);
		      j1 = j0 + 1.0f; 

		      s1 = x - i0; 
		      s0 = 1.0f - s1; 
		      t1 = y - j0; 
		      t0 = 1.0f - t1;

		      int i0i = (int)i0;
		      int i1i = (int)i1;
		      int j0i = (int)j0;
		      int j1i = (int)j1;

		      // DOUBLE CHECK THIS!!!
		      d[IX(i, j)] = 
		        s0 * (t0 * d0[IX(i0i, j0i)] + t1 * d0[IX(i0i, j1i)]) +
		        s1 * (t0 * d0[IX(i1i, j0i)] + t1 * d0[IX(i1i, j1i)]);
		    }
		  }

		  set_bnd(b, d);
		}



		public static void set_bnd(int b, float[] x) {
		  for (int i = 1; i < N - 1; i++) {
		    x[IX(i, 0  )] = b == 2 ? -x[IX(i, 1  )] : x[IX(i, 1 )];
		    x[IX(i, N-1)] = b == 2 ? -x[IX(i, N-2)] : x[IX(i, N-2)];
		  }
		  for (int j = 1; j < N - 1; j++) {
		    x[IX(0, j)] = b == 1 ? -x[IX(1, j)] : x[IX(1, j)];
		    x[IX(N-1, j)] = b == 1 ? -x[IX(N-2, j)] : x[IX(N-2, j)];
		  }

		  x[IX(0, 0)] = 0.5f * (x[IX(1, 0)] + x[IX(0, 1)]);
		  x[IX(0, N-1)] = 0.5f * (x[IX(1, N-1)] + x[IX(0, N-2)]);
		  x[IX(N-1, 0)] = 0.5f * (x[IX(N-2, 0)] + x[IX(N-1, 1)]);
		  x[IX(N-1, N-1)] = 0.5f * (x[IX(N-2,N-1)] + x[IX(N-1, N-2)]);
		}
}

FluidGame.java

import com.badlogic.gdx.Game;
import com.badlogic.gdx.Gdx;
import com.badlogic.gdx.graphics.GL20;
import com.badlogic.gdx.graphics.glutils.ShapeRenderer;
import com.badlogic.gdx.math.MathUtils;
import com.badlogic.gdx.math.Vector2;
import com.badlogic.gdx.utils.viewport.FitViewport;
import com.badlogic.gdx.utils.viewport.Viewport;

public class FluidGame extends Game
{
	float t = 0;
	
	public static final int SIZE = Fluid.N * Fluid.SCALE;
	public static final int width = SIZE;
	public static final int height = SIZE;

	private Viewport viewport;
	private ShapeRenderer renderer;
	Fluid fluid;

	@Override
	public void create() {
		fluid = new Fluid(1, .5f, 1e-4f); // 1e-7f);
		viewport = new FitViewport(SIZE, SIZE);
		renderer = new ShapeRenderer();
		
		float s = .1f;
		for(int y=0 ; y<Fluid.N ; y++){
			for(int x=0 ; x<Fluid.N ; x++){
				fluid.addDensity(x, y, 255f * ((MathUtils.random()*2f-1f) * .1f + .5f));
				fluid.addVelocity(x, y, MathUtils.random(-s,s), MathUtils.random(-s,s) );
			}
		}
	}
	
	@Override
	public void resize(int width, int height) {
		viewport.update(width, height, true);
	}
	
	@Override
	public void render() {
		fluid.diff = 1e-8f;
		fluid.visc = 1e-9f;
		
		Gdx.gl.glClear(GL20.GL_COLOR_BUFFER_BIT);
		
		float s = .01f;
		for(int y=1 ; y<Fluid.N-1 ; y++){
			for(int x=2 ; x<3 ; x++){
				//fluid.addVelocity(x, y, MathUtils.random(s/2, s), MathUtils.random(-s/3,s/3) );
			}
		}
	  int cx = (int)(0.5f*width/Fluid.SCALE);
	  int cy = (int)(0.5*height/Fluid.SCALE);
	  Vector2 p = viewport.unproject(new Vector2(Gdx.input.getX(), Gdx.input.getY())).scl(1f / Fluid.SCALE);
	  cx = (int)p.x;
	  cy = (int)p.y;
	  int kernel = 1;
	  for (int i = -kernel; i <= kernel; i++) {
	    for (int j = -kernel; j <= kernel; j++) {
	      fluid.addDensity(cx+i, cy+j, 112f); //MathUtils.random(50f, 150f));
	    }
	  }
	  for (int i = 0; i < 2; i++) {
	    float angle = t * MathUtils.PI2 * 2;
	    Vector2 v = new Vector2(1,0); //.setAngleRad(angle);
	    v.scl(0.2f);
	    t += Gdx.graphics.getDeltaTime() * .01f;
	    fluid.addVelocity(cx, cy, v.x, v.y );
	  }

	  fluid.step();

	  renderer.setProjectionMatrix(viewport.getCamera().combined);
	  fluid.renderD(renderer);
	 // fluid.renderV(renderer);
	 fluid.fadeD();
	}

}