3D Waterwave Simulation using Python

I used Numpy  Matplotlib with Animation and 3d Plot module on my OS X Yosemite.

The code is still messy and clearly not efficient (there's slow loop here and there) but it works, :)
Here The Result
The Code

import numpy as np

n = 8;
g = 9.8;
dt = 0.02;
dx = 1.0;
dy = 1.0;

h = np.ones((n+2,n+2))
u = np.zeros((n+2,n+2))
v = np.zeros((n+2,n+2))

hx = np.zeros((n+1,n+1))
ux = np.zeros((n+1,n+1))
vx = np.zeros((n+1,n+1))

hy = np.zeros((n+1,n+1))
uy = np.zeros((n+1,n+1))
vy = np.zeros((n+1,n+1))

nsteps = 0
h[1,1] = .5;

def reflective():
    h[:,0] = h[:,1]
    h[:,n+1] = h[:,n]
    h[0,:] = h[1,:]
    h[n+1,:] = h[n,:]
    u[:,0] = u[:,1]
    u[:,n+1] = u[:,n]
    u[0,:] = -u[1,:]
    u[n+1,:] = -u[n,:]
    v[:,0] = -v[:,1]
    v[:,n+1] = -v[:,n]
    v[0,:] = v[1,:]
    v[n+1,:] = v[n,:]

def proses():
    #hx = (h[1:,:]+h[:-1,:])/2-dt/(2*dx)*(u[1:,:]-u[:-1,:])
    for i in range (n+1):
        for j in range(n):
            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])
            ux[i,j] = (u[i+1,j+1]+u[i,j+1])/2- dt/(2*dx)*((pow(u[i+1,j+1],2)/h[i+1,j+1]+ g/2*pow(h[i+1,j+1],2))- (pow(u[i,j+1],2)/h[i,j+1]+ g/2*pow(h[i,j+1],2)))
            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]))

    for i in range (n):
        for j in range(n+1):
            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])
            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]))
            vy[i,j] = (v[i+1,j+1]+v[i+1,j])/2 - dt/(2*dy)*((pow(v[i+1,j+1],2)/h[i+1,j+1] + g/2*pow(h[i+1,j+1],2)) - (pow(v[i+1,j],2)/h[i+1,j] + g/2*pow(h[i+1,j],2)))
    
    for i in range (1,n+1):
        for j in range(1,n+1):
            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])
            u[i,j] = u[i,j] - (dt/dx)*((pow(ux[i,j-1],2)/hx[i,j-1] + g/2*pow(hx[i,j-1],2)) - (pow(ux[i-1,j-1],2)/hx[i-1,j-1] + g/2*pow(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]))
            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)*((pow(vy[i-1,j],2)/hy[i-1,j] + g/2*pow(hy[i-1,j],2)) - (pow(vy[i-1,j-1],2)/hy[i-1,j-1] + g/2*pow(hy[i-1,j-1],2)))

#dh = dt/dt*(ux[1:,:]-ux[:-1,:])+ dt/dy*(vy[:,1:]-vy[:,:-1])
    reflective()
    return h,u,v
'''
for i in range (17):
    #print h
    proses(1)
'''

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from mpl_toolkits.mplot3d import Axes3D
a = n
x = np.arange(n+2)
y = np.arange(n+2)
x,y = np.meshgrid(x,y)

fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')

def plotset():
    ax.set_xlim3d(0, a)
    ax.set_ylim3d(0, a)
    ax.set_zlim3d(0.5,1.5)
    ax.set_autoscalez_on(False)
    ax.zaxis.set_major_locator(LinearLocator(10))
    ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
    cset = ax.contour(x, y, h, zdir='x', offset=0 , cmap=cm.coolwarm)
    cset = ax.contour(x, y, h, zdir='y', offset=n , cmap=cm.coolwarm)
    cset = ax.contour(x, y, h, zdir='z', offset=.5, cmap=cm.coolwarm)

plotset()

surf = ax.plot_surface(x, y, h,rstride=1, cstride=1,cmap=cm.coolwarm,linewidth=0,antialiased=False, alpha=0.7)

fig.colorbar(surf, shrink=0.5, aspect=5)


from matplotlib import animation


def data(k,h,surf):
    proses()
    ax.clear()
    plotset()
    surf = ax.plot_surface(x, y, h,rstride=1, cstride=1,cmap=cm.coolwarm,linewidth=0,antialiased=False, alpha=0.7)
    return surf,

ani = animation.FuncAnimation(fig, data, fargs=(h,surf), interval=10, blit=False)
#ani.save('air.mp4', bitrate=512)
plt.show()

and the snapshot

3D Surface Plot Animation using Matplotlib in Python

And here's the animation

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D

def data(i, z, line):
    z = np.sin(x+y+i)
    ax.clear()
    line = ax.plot_surface(x, y, z,color= 'b')
    return line,

n = 2.*np.pi
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

x = np.linspace(0,n,100)
y = np.linspace(0,n,100)
x,y = np.meshgrid(x,y)
z = np.sin(x+y)
line = ax.plot_surface(x, y, z,color= 'b')

ani = animation.FuncAnimation(fig, data, fargs=(z, line), interval=30, blit=False)

plt.show()

The result

The snapshot

3D Surface Plot using Matplotlib in Python

It's slightly modified from before

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D

n = 2.*np.pi
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

x = np.linspace(0,n,100)
y = np.linspace(0,n,100)
x,y = np.meshgrid(x,y)
z = np.sin(x+y)
line = ax.plot_surface(x, y, z,color= 'b')

plt.show()

the result


the snapshot

Matplotlib Animation in Python

Here is the update from before

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

def simData():
    t_max = n
    dt = 1./8
    k = 0.0
    t = np.linspace(0,t_max,100)
    while k < t_max:
        x = np.sin(np.pi*t+np.pi*k)
        k = k + dt
        yield x, t

def simPoints(simData):
    x, t = simData[0], simData[1]
    line.set_data(t, x)
    return line
n = 2.
fig = plt.figure()
ax = fig.add_subplot(111)
line, = ax.plot([], [], 'b')
ax.set_ylim(-1, 1)
ax.set_xlim(0, n)

ani = animation.FuncAnimation(fig, simPoints, simData, blit=False,\
                              interval=100, repeat=True)
plt.show()

and the result

Playing with Matplotlib Animation in Python

Coding like this

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import animation

fig = plt.figure()
n = 10
x = np.linspace(0,2*np.pi,100)



def init():
    pass
def animate(k):
    h = np.sin(x+np.pik)
    plt.plot(x,h)


ax = plt.axes(xlim=(0, 2*np.pi), ylim=(-1.1, 1.1))

anim = animation.FuncAnimation(fig, animate,init_func=init,frames=360,interval=20,blit=False)

plt.show()

The result

Playing (again) with 'Home Made' Vector in Delphi

Here it is. I create a vector as new type, which is in itself is three dimension array.

Then I declared u as vector with three dimension;
u (h,i,j)

where h = 0, 1, 2  as physical component (eg: height, velocity, momentum)
i , j = 0, 1, 2, ..., n as row n column


So if we read u[0,1,1], it means height value at coordinate (1,1); u[1,1,1] is the velocity value; [2,1,1] is the momentum value at the same coordinate.


Trying some of properties of it. I found out that we can initialize all component of vector-u with this one line code

u:=fu(h[i,j],i,j);

so the component u(h,i,j) will filled. Notice that the function has vector (or in this case array) return value.



The code below show how I fill the value of component u(0, i, j)

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;
const n=3;
type vector=array[0..2,0..n,0..n]of real;

type
  TForm1 = class(TForm)
    Button1: TButton;
    Memo1: TMemo;
    function fu(a:real;i,j:integer):vector;
    procedure FormCreate(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation{$R *.dfm}

function tform1.fu(a:real;i,j:integer):vector;
begin
  fu[0,i,j]:=a;
end;

procedure TForm1.FormCreate(Sender: TObject);
var i,j:integer;
h:array[0..n,0..n]of real;
u:vector;
begin
  for i:=0 to n do begin
    for j:=0 to n do begin
       h[i,j]:=1;
       u:=fu(h[i,j],i,j);
    end;
  end;
  memo1.Text:='';
  memo1.Lines.Append('h[1,1]='+floattostr(h[1,1]));
  memo1.Lines.Append('u[0,1,1]='+floattostr(u[0,1,1]));
  memo1.Lines.Append('u[0,2,1]='+floattostr(u[0,2,1]));
end;

end.

:)

Returning Function as Array in Delphi

Do you wonder how to do vector operation in Delphi? No, of course, :).

We could go like this.

function tform1.adv(a,b:real):real;
begin
  adv:=a+b;
end;

The problem is the return is real, which is single value only. We want a and b as vector. Wait...

How we define vector in Delphi? I don't know. I used to treat a vector in Delphi as array. So I coded it like this

var a,b:array[0..1]of real;

So far I had no problem. Lately, I am going crazy with overuse functions in Delphi, and trying operating vectors using function too.

But if I write the code like this

function tform1.adv(a,b:array[0..1]of real):real;
begin
  adv:=a[0]+b[0];  
  {a[1]+b[1]?}
end;

It will only return one value. So I improvised by modify it



like this

function tform1.adv(a,b:array[0..1]of real):array[0..1]ofreal;
begin
  adv[0]:=a[0]+b[0];  
  adv[1]:=a[1]+b[1];  
end;

But it won't compile. (it will give error message "identifier expected but ARRAY found"). So I try another approach

type
  vector=array[0..1] of real;

function tform1.adv(a,b:vector):vector;
begin
  adv[0]:=a[0]+b[0];
  adv[1]:=a[1]+b[1];
end;

It works, :).

Here my last night tinkering with "vector" in Delphi, :)

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, Math;
type
  vector=array[0..1] of real;

type
  TForm1 = class(TForm)
    procedure proses;
    function mux(a:real;b:vector):vector;
    function dot(a,b:vector):real;
    function norm(a,b:vector):vector;
    function adv(a,b:vector):vector;
    function suv(a,b:vector):vector;
    procedure FormCreate(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;
const
  n=3;
var
  Form1: TForm1;
  r,v:array[1..n] of vector;
implementation

{$R *.dfm}
function tform1.adv(a,b:vector):vector;
begin
  adv[0]:=a[0]+b[0];
  adv[1]:=a[1]+b[1];
end;
function tform1.suv(a,b:vector):vector;
begin
  suv[0]:=a[0]-b[0];
  suv[1]:=a[1]-b[1];
end;
function tform1.mux(a:real;b:vector):vector;
begin
  mux[0]:=a*b[0];
  mux[1]:=a*b[1];
end;

function tform1.dot(a,b:vector):real;
begin
  dot:=a[0]*b[0]+a[1]*b[1];
end;
function tform1.norm(a,b:vector):vector;
var mag,i,j:real;
begin
  i:=b[0]-a[0];
  j:=b[1]-a[1];
  mag:=sqrt(power(i,2)+power(j,2));
  if mag<>0 then begin
    norm[0]:=i/mag;
    norm[1]:=j/mag;
  end;
end;
procedure tform1.proses;
var direction:vector;
vi,vj,swap:real;
i,j:integer;
begin
  j:=2;i:=1;
  direction:=norm(r[j],r[i]);//call function
  vi:=dot(v[i],direction);
  vj:=dot(v[j],direction);
  swap:=vj-vi;
  v[i]:=adv(v[i],mux(swap,direction));
  v[j]:=suv(v[j],mux(swap,direction));
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
  proses;
end;

end.