Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
#!/bin/bash
gfortran -o $2 -fno-underscoring $1 libhgraph.a libg2c.so -L/usr/X11R6/lib64 -lX11
: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
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(nvortex,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
function test(nvortex, nsteps, delta_t, xyg)
vxy_old = zeros(nvortex,2)
xy_time = zeros(nsteps, nvortex, 2)
for step=1:nsteps
vxy = find_vel(nvortex,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
end
@time nvortex, nsteps, delta_t, xyg = load_dat("vortex.dat")
@time test(nvortex, nsteps, delta_t, xyg )
@time test(nvortex, nsteps, delta_t, xyg )
@time test(nvortex, nsteps, delta_t, xyg )
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()
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.