compile_fortran.sh

#!/bin/bash
gfortran -o $2 -fno-underscoring $1 libhgraph.a libg2c.so -L/usr/X11R6/lib64 -lX11

results

: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.dat

   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
  

vortex.f

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

vortex.jl

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)

vortex.py

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()