Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Vortex Performance comparison if Julia (Julia-lang), Fortran and Python
#!/bin/bash
gfortran -o $2 -fno-underscoring $1 libhgraph.a libg2c.so -L/usr/X11R6/lib64 -lX11
50,000 time steps:
:vortex || time ./compile_fortran.sh vortex.f vortexf;time ./vortexf
real 0m0.060s
user 0m0.039s
sys 0m0.022s
real 0m0.051s
user 0m0.051s
sys 0m0.000s
:vortex || time julia vortex.jl
Time taken to load dat file: 0.5606400966644287secs
Time taken to compute 50000 time steps with 11 vortices: 2.2556450366973877secs
Press enter to finish:
real 0m8.414s
user 0m8.294s
sys 0m0.127s
:vortex || time python vortex.py
Time taken to load dat file: 5.91278076172e-05 secs
Time taken to compute 50000 time steps with 11 vortices: 30.8461220264 secs
real 0m31.218s
user 0m31.152s
sys 0m0.016s
:vortex.git || time julia vortex_fast.jl
elapsed time: 0.209670419 seconds (101222184 bytes allocated)
Time taken to load dat file: 0.5140218734741211secs
Time taken to compute 50000 time steps with 11 vortices: 0.7584049701690674secs
Press enter to finish:
real 0m6.947s
user 0m6.830s
sys 0m0.124s
500,000 time steps:
:vortex.git || time ./compile_fortran.sh vortex.f vortexf;time ./vortexf
real 0m0.049s
user 0m0.035s
sys 0m0.014s
real 0m0.495s
user 0m0.496s
sys 0m0.000s
:vortex.git || time julia vortex_fast.jl
elapsed time: 1.547171999 seconds (994022184 bytes allocated)
Time taken to load dat file: 0.5128679275512695secs
Time taken to compute 500000 time steps with 11 vortices: 2.0945730209350586secs
Press enter to finish:
real 0m8.521s
user 0m8.389s
sys 0m0.139s
11 50000 No of vortices and No of steps
0.001 -10.0 Length of time step in seconds, source strength
0.5 0.7 1.0 X,Y coordinates and strength of vortices
1.0 1.0 -1.0 Repeat for each NVORTEX vortices
-1.0 1.0 1.0
-1.0 -1.0 -1.0
1.0 -1.0 1.0
0.5 0.5 2.0
1.0 0.0 -2.0
0.0 -1.0 -0.5
-0.5 0.5 0.5
.1 .1 -2.
-.1 -.1 2.0
C
C PROGRAM TO COMPUTE THE MOTION OF A CLOUD OF VORTICES
C
COMMON NVORTEX,GAMMA(100),X(100),Y(100),VX(100),VY(100),
& VX_OLD(100),VY_OLD(100),DELTA_T,NSTEPS
C
INTEGER STEP
C
OPEN(UNIT=1,FILE='vortex.dat')
C
CALL INPUT
C
CLOSE(1)
C
C OPEN THE PLOTTING UNIT USING GRAPHICS LIBRARY ROUTINE 'SELPLT'.
C
C Switch off plotting for performance test:
! CALL SELPLT(8)
C
C FIND THE SIZE OF THE DOMAIN CONTAINING THE VORTICES
C
XMIN = 0.0
XMAX = 0.0
YMIN = 0.0
YMAX = 0.0
DO 10 N = 1,NVORTEX
IF(X(N).GT.XMAX) XMAX = X(N)
IF(X(N).LT.XMIN) XMIN = X(N)
IF(Y(N).GT.YMAX) YMAX = Y(N)
IF(Y(N).LT.YMIN) YMIN = Y(N)
10 CONTINUE
C
C SET THE PLOTTING REGION USING GRAPHICS LIBRARY ROUTINE 'SETBOX'.
C
C Switch off plotting for performance test:
! CALL SETBOX(2*XMIN,2*XMAX,2*YMIN,2*YMAX,0,1)
C
DO 1000 STEP =1,NSTEPS
C
C OBTAIN THE CURRENT VELOCITY OF THE VORTICES FROM SUBROUTINE
C FIND_VEL
C
CALL FIND_VEL
C
C SET THE OLD VALUES OF VELOCITY EQUAL TO THE INITIAL
C VELOCITIES FOR THE FIRST STEP ONLY.
C
IF(STEP.EQ.1) THEN
DO 50 N=1,NVORTEX
VX_OLD(N) = VX(N)
VY_OLD(N) = VY(N)
50 CONTINUE
ENDIF
C
C FIND THE NEW POSITIONS OF THE VORTICES
C
DO 100 N=1,NVORTEX
X(N) = X(N) + 0.5*(VX(N) + VX_OLD(N))*DELTA_T
Y(N) = Y(N) + 0.5*(VY(N) + VY_OLD(N))*DELTA_T
VX_OLD(N) = VX(N)
VY_OLD(N) = VY(N)
100 CONTINUE
C
C NOW PLOT THE VORTEX POSITIONS USING GRAPHICS LIBRARY ROUTINE 'PLTSYM'.
C
C Switch off plotting for performance test:
! DO N=1,NVORTEX
! CALL PLTSYM(X(N),Y(N),'C')
! END DO
C
1000 CONTINUE
C
C CLOSE THE PLOTTING UNIT USING GRAPHICS LIBRARY ROUTINE 'BRKPLT'.
C
C Switch off plotting for performance test:
! CALL BRKPLT
C
STOP
END
C
C***************************************************************************
SUBROUTINE INPUT
C
COMMON NVORTEX,GAMMA(100),X(100),Y(100),VX(100),VY(100),
& VX_OLD(100),VY_OLD(100),DELTA_T,NSTEPS
C
READ(1,*) NVORTEX,NSTEPS
READ(1,*) DELTA_T
DO N = 1,NVORTEX
READ(1,*) X(N),Y(N),GAMMA(N)
END DO
RETURN
END
C
C**************************************************************************
SUBROUTINE FIND_VEL
C
COMMON NVORTEX,GAMMA(100),X(100),Y(100),VX(100),VY(100),
& VX_OLD(100),VY_OLD(100),DELTA_T,NSTEPS
C
TWOPI = 2*3.1415926
DO 100 N=1,NVORTEX
XVEL = 0.0
YVEL = 0.0
DO 90 M=1,NVORTEX
IF(M.EQ.N) GO TO 90
DX = X(N)-X(M)
DY = Y(N)-Y(M)
DISTSQ = DX*DX + DY*DY
XVEL = XVEL + GAMMA(M)/TWOPI*DY/DISTSQ
YVEL = YVEL - GAMMA(M)/TWOPI*DX/DISTSQ
90 CONTINUE
VX(N) = XVEL
VY(N) = YVEL
100 CONTINUE
RETURN
END
function load_dat(fname::String)
lines = String[]
open(fname,"r") do f
for line in eachline(f)
push!(lines, line)
end
end
line_split = split(lines[1])
nvortex = int(line_split[1])
nsteps = int(line_split[2])
line_split = split(lines[2])
delta_t = float(line_split[1])
strength = float(line_split[2])
xyg=zeros(nvortex,3)
for v = 1:nvortex
line_split = split(lines[v + 2])
for j = 1:3
xyg[v, j] = float(line_split[j])
end
end
return nvortex, nsteps, delta_t, xyg
end
function find_vel(xyg::Array{Float64,2})
vxy = zeros(nvortex,2)
twopi = 2*pi
for v1 = 1:nvortex
for v2 = 1:nvortex
if v1 == v2
continue
end
dx = xyg[v1,1] - xyg[v2, 1]
dy = xyg[v1,2] - xyg[v2, 2]
distsq = dx^2 + dy^2
vxy[v1, 1] += xyg[v2,3]/twopi*dy/distsq
vxy[v1, 2] -= xyg[v2,3]/twopi*dx/distsq
end
end
return vxy
end
t0 = time()
nvortex, nsteps, delta_t, xyg = load_dat("vortex.dat")
t1 = time()
println("Time taken to load dat file: ", t1-t0, "secs")
vxy_old = zeros(nvortex,2)
xy_time = zeros(nsteps, nvortex, 2)
for step=1:nsteps
vxy = find_vel(xyg)
if step == 1
vxy_old = vxy
end
for dim = 1:2
xyg[:, dim] += 0.5*(vxy[:, dim] + vxy_old[:, dim])*delta_t
xy_time[step,:,dim] = xyg[:,dim]
end
vxy_old = vxy
end
t1 = time()
println("Time taken to compute ", nsteps, " time steps with ", nvortex, " vortices: ", t1-t0, "secs")
using PyPlot
hold(true)
for v = 1:nvortex
plot(xy_time[:,v,1], xy_time[:,v,2])
end
#println("Press enter to finish: ")
#readline(STDIN)
import numpy as np
from time import time
def load_dat(fname):
lines = open(fname, 'r').read().split('\n')
line_split = [i for i in lines[0].split() if i != '']
nvortex = int(line_split[0])
nsteps = int(line_split[1])
line_split = [i for i in lines[1].split() if i != '']
delta_t = float(line_split[0])
strength = float(line_split[1])
xyg=np.zeros([nvortex,3])
for v in range(nvortex):
line_split = [i for i in lines[v + 2].split() if i != '']
for j in range(3):
xyg[v, j] = float(line_split[j])
return nvortex, nsteps, delta_t, xyg
def find_vel(xyg, nvortex):
vxy = np.zeros([nvortex,2])
for v1 in range(nvortex):
for v2 in range(nvortex):
if v1 == v2:
continue
dx = xyg[v1, 0] - xyg[v2, 0]
dy = xyg[v1, 1] - xyg[v2, 1]
distsq = dx*dx + dy*dy
vxy[v1, 0] += xyg[v2,2]/(2*np.pi)*dy/distsq
vxy[v1, 1] -= xyg[v2,2]/(2*np.pi)*dx/distsq
return vxy
t0 = time()
nvortex, nsteps, delta_t, xyg = load_dat("vortex.dat")
t1 = time()
print "Time taken to load dat file:", t1-t0, "secs"
vxy_old = np.zeros([nvortex,2])
xy_time = np.zeros([nsteps, 2, nvortex])
for step in range(nsteps):
vxy = find_vel(xyg, nvortex)
if step == 1:
vxy_old = np.copy(vxy)
for v in range(nvortex):
for dim in range(2):
xyg[v, dim] += 0.5*(vxy[v, dim] + vxy_old[v, dim])*delta_t
xy_time[step,dim,v] = xyg[v, dim]
vxy_old = np.copy(vxy)
t1 = time()
print "Time taken to compute", nsteps, "time steps with", nvortex, "vortices:", t1-t0, "secs"
import matplotlib.pyplot as plt
plt.hold(True)
for v in range(nvortex):
plt.plot(xy_time[:,0,v], xy_time[:,1,v])
#plt.show()
function load_dat(fname::String)
lines = String[]
open(fname,"r") do f
for line in eachline(f)
push!(lines, line)
end
end
line_split = split(lines[1])
nvortex = int(line_split[1])
nsteps = int(line_split[2])
line_split = split(lines[2])
delta_t = float64(line_split[1])
# strength = float(line_split[2])
x=zeros(Float64, nvortex)
y=zeros(Float64, nvortex)
g=zeros(Float64, nvortex)
for v = 1:nvortex
line_split = split(lines[v + 2])
x[v] = float(line_split[1])
y[v] = float(line_split[2])
g[v] = float(line_split[3])
end
return nvortex, nsteps, delta_t, x, y, g
end
function find_vel(nvortex::Int64, x::Array{Float64,1}, y::Array{Float64,1}, g::Array{Float64,1})
vx = zeros(Float64, nvortex)
vy = zeros(Float64, nvortex)
twopi = 2*pi
for v1 = 1:nvortex
for v2 = 1:nvortex
if v1 == v2
continue
end
dx = x[v1] - x[v2]
dy = y[v1] - y[v2]
distsq = dx*dx + dy*dy
vx[v1] += g[v2]/twopi*dy/distsq
vy[v1] -= g[v2]/twopi*dx/distsq
end
end
return vx, vy
end
function perform_calc(nvortex::Int64, nsteps::Int64, delta_t::Float64, x::Array{Float64,1}, y::Array{Float64,1}, g::Array{Float64,1})
vx_old = zeros(Float64, nvortex)
vy_old = zeros(Float64, nvortex)
x_time = zeros(Float64, nsteps, nvortex)
y_time = zeros(Float64, nsteps, nvortex)
for step=1:nsteps
vx, vy = find_vel(nvortex, x, y, g)
if step == 1
vx_old = vx
vy_old = vy
end
x += 0.5*(vx + vx_old)*delta_t
y += 0.5*(vy + vy_old)*delta_t
x_time[step,:] = x
y_time[step,:] = y
vx_old = vx
vy_old = vy
end
return x_time, y_time
end
t0 = time()
nvortex, nsteps, delta_t, x, y, g = load_dat("vortex.dat")
t1 = time()
@time x_time, y_time = perform_calc(nvortex, nsteps, delta_t, x, y, g)
t2 = time()
println("Time taken to load dat file: ", t1-t0, "secs")
println("Time taken to compute ", nsteps, " time steps with ", nvortex, " vortices: ", t2-t0, "secs")
using PyPlot
hold(true)
for v = 1:nvortex
plot(x_time[:,v], y_time[:,v])
end
println("Press enter to finish: ")
# readline(STDIN)
@ViralBShah

This comment has been minimized.

Copy link

commented Feb 15, 2014

Could you try rewriting the julia version devectorized?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.