Height field fluids simulation and shallow water equation

So in this project I am about to simulate a water surface with simple interaction.

Source code can be found at source.

There are basically two ways of doing it : the height field fluid approach and shallow water equation approach. We would talk about the first approach in this post.

A good tutorial on height field fluids can be found here: http://www.matthiasmueller.info/talks/gdc2008.pdf

Rendering is based on the code found following this link: http://www.lousodrome.net/opengl/

The base code uses simple sine waves to simulate water flows. This post focuses on how we can create more complicated wave simulations based on simple physics law.

Height field fluid considers the interaction force as proportional to the curvature of the height function. Thus according to newton's law:

f(i,j) = curvature(i,j)
v(i,j) = v_old(i,j) + delta_t * f(i,j)
h(i,j) = h_old(i,j) + delta_t * v(i,j)

This is a simple update function. Implementation can be found in main_height_fluid.c
Code for update the height field:

void    update_height(void)
{
int i,j;
for(j = 0; j < RESOLUTION; j++)
for(i = 0; i <= RESOLUTION; i++)
{
float height_east = height[MIN(i+1,RESOLUTION)][j];
float height_west = height[MAX(i-1,0)][j];
float height_south = height[i][MAX(j-1,0)];
float height_north = height[i][MIN(j+1, RESOLUTION-1)];
height_v[i][j] += ((height_west+ height_east + height_south + height_north)/4 - height[i][j]);
height_v[i][j] *= damping;
height[i][j] += height_v[i][j];
}
}

Desktop recording available at:

Shallow water equation is more complicated. Basically it considers the conservation of mass and momentum in two directions.  It also takes the velocity of wave in two directions into account.

Code for update the height field using shallow water equation. The code is in Matlab which comes from the reference I posted above. In the source code you would find my C++ /matlab version of this equation solver. :

% Reflective boundary conditions
H(:,1) = H(:,2);      U(:,1) = U(:,2);       V(:,1) = -V(:,2);
H(:,n+2) = H(:,n+1);  U(:,n+2) = U(:,n+1);   V(:,n+2) = -V(:,n+1);
H(1,:) = H(2,:);      U(1,:) = -U(2,:);      V(1,:) = V(2,:);
H(n+2,:) = H(n+1,:);  U(n+2,:) = -U(n+1,:);  V(n+2,:) = V(n+1,:);

% First half step

% x direction
i = 1:n+1;
j = 1:n;

% height
Hx(i,j) = (H(i+1,j+1)+H(i,j+1))/2 - dt/(2*dx)*(U(i+1,j+1)-U(i,j+1));

% x momentum
Ux(i,j) = (U(i+1,j+1)+U(i,j+1))/2 -  ...
dt/(2*dx)*((U(i+1,j+1).^2./H(i+1,j+1) + g/2*H(i+1,j+1).^2) - ...
(U(i,j+1).^2./H(i,j+1) + g/2*H(i,j+1).^2));

% y momentum
Vx(i,j) = (V(i+1,j+1)+V(i,j+1))/2 - ...
dt/(2*dx)*((U(i+1,j+1).*V(i+1,j+1)./H(i+1,j+1)) - ...
(U(i,j+1).*V(i,j+1)./H(i,j+1)));

% y direction
i = 1:n;
j = 1:n+1;

% height
Hy(i,j) = (H(i+1,j+1)+H(i+1,j))/2 - dt/(2*dy)*(V(i+1,j+1)-V(i+1,j));

% x momentum
Uy(i,j) = (U(i+1,j+1)+U(i+1,j))/2 - ...
dt/(2*dy)*((V(i+1,j+1).*U(i+1,j+1)./H(i+1,j+1)) - ...
(V(i+1,j).*U(i+1,j)./H(i+1,j)));
% y momentum
Vy(i,j) = (V(i+1,j+1)+V(i+1,j))/2 - ...
dt/(2*dy)*((V(i+1,j+1).^2./H(i+1,j+1) + g/2*H(i+1,j+1).^2) - ...
(V(i+1,j).^2./H(i+1,j) + g/2*H(i+1,j).^2));

% Second half step
i = 2:n+1;
j = 2:n+1;

% height
H(i,j) = H(i,j) - (dt/dx)*(Ux(i,j-1)-Ux(i-1,j-1)) - ...
(dt/dy)*(Vy(i-1,j)-Vy(i-1,j-1));
% x momentum
U(i,j) = U(i,j) - (dt/dx)*((Ux(i,j-1).^2./Hx(i,j-1) + g/2*Hx(i,j-1).^2) - ...
(Ux(i-1,j-1).^2./Hx(i-1,j-1) + g/2*Hx(i-1,j-1).^2)) ...
- (dt/dy)*((Vy(i-1,j).*Uy(i-1,j)./Hy(i-1,j)) - ...
(Vy(i-1,j-1).*Uy(i-1,j-1)./Hy(i-1,j-1)));
% y momentum
V(i,j) = V(i,j) - (dt/dx)*((Ux(i,j-1).*Vx(i,j-1)./Hx(i,j-1)) - ...
(Ux(i-1,j-1).*Vx(i-1,j-1)./Hx(i-1,j-1))) ...
- (dt/dy)*((Vy(i-1,j).^2./Hy(i-1,j) + g/2*Hy(i-1,j).^2) - ...
(Vy(i-1,j-1).^2./Hy(i-1,j-1) + g/2*Hy(i-1,j-1).^2)); 