Any way to speed up code?

3 views (last 30 days)
Aslan
Aslan on 13 Oct 2014
Commented: Stephen23 on 13 Oct 2014
Hey guys, I've generated a code for calculating the stream function and vorticity in a driven cavity flow, but the code I am running takes forever to iterate (should take roughly 20,000 iterations) is there any way to speed this code up? Cheers
clc
clear all
close all
nx = 196;
ny = 196;
dx = 1.0/nx;
dy = 1.0/ny;
dx2 = dx^2;
dy2 = dy^2;
relax1 = 0.2;
relax2 = 1.2;
re = 3200;
anu = 1/re;
rex = 2.0/dx2 + 2.0/dy2;
for i = 1:nx+1
for j = 1:ny+1
psi(i,j) = 0.0;
w(i,j) = 0.0;
x(i,j) = (i-1)*dx;
y(i,j) = (j-1)*dy;
end;
end;
for k = 1:50000
adu = 0.0;
adv = 0.0;
psimin = 1000000.0;
psimax = -1000000.0;
for i = 2:nx %sweep through all internal points on the grid
for j = 2:ny %updating until convergence.
wimj0 = w(i-1,j); %for the vorticity equations these values
wipj0 = w(i+1,j); %may have to be modified at the points
wi0jm = w(i,j-1); %adjacent to the walls
wi0jp = w(i,j+1);
%Inclusion of Boundary Conditions
if(i == 2) %Side Wall
wimj0 = -2.0*(psi(i,j))/dx2;
end;
if(i == nx)%Side Wall
wipj0 = -2.0*(psi(i,j))/dx2;
end;
if(j == 2)%Bottom Wall
wi0jm = -2.0*(psi(i,j))/dy2;
end;
if(j == ny)%Top Wall
wi0jp = -2.0*(psi(i,j)+dy)/dy2;
end;
wfunc1 = ((wipj0+wimj0)/dx2+(wi0jp+wi0jm)/dy2)*anu;
wfunc2 = (psi(i,j+1)-psi(i,j-1))*(wipj0-wimj0)/(4*dx*dy);
wfunc3 = (psi(i+1,j)-psi(i-1,j))*(wi0jp-wi0jm)/(4*dx*dy);
wij = (wfunc1-wfunc2+wfunc3)/(anu*rex); %Second Order FDE for Vorticity
dw = wij - w(i,j); %compute the difference between new and old value.
w(i,j) = w(i,j) + relax1*dw; %update using relaxation
Psifunc1 = (psi(i+1,j)+psi(i-1,j))/dx2+(psi(i,j+1)+psi(i,j-1))/dy2;
psiij = (Psifunc1+wij)/rex; %Second Order FDE for Streamfunction
dpsi = psiij - psi(i,j);
psi(i,j) = psi(i,j) + relax2*dpsi;
if(psi(i,j)>psimax) %determine maximum and minimum streamfunction
psimax = psi(i,j);
end;
if(psi(i,j)<psimin)
psimin = psi(i,j);
end;
ddu = abs(dpsi); %compute maximum change in solution over full sweep.
ddv = abs(dw);
if(ddu>adu)
adu = ddu;
end;
if(ddv>adv)
adv = ddv;
end;
end;
end;
k
adu
adv
if(adu<10.0e-7 && adv<10.0e-7) %exit loop if converged
break;
end;
end;
for i = 1:10
v(i) = -(i-1)*0.01;
end;
for i = 11:20
v(i) = (i-10)*psimax/10;
end;
contourf(x,y,psi,v);

Accepted Answer

Stephen23
Stephen23 on 13 Oct 2014
Edited: Stephen23 on 13 Oct 2014
Vectorization is the key. In MATLAB it is faster to perform operations and functions on whole arrays, rather than using loops. These pages will give you more information:
For checking how your code executes (eg. to find the slow parts), you can use the profiler . You will find examples and other useful info here:
  2 Comments
Aslan
Aslan on 13 Oct 2014
Thanks mate! This looks like what I'm after
Stephen23
Stephen23 on 13 Oct 2014
My pleasure. Vectorization is pretty straightforward, the trick is to think of it as maths, and not as code.
PS: You can also accept any answer that solves your problem ;-)

Sign in to comment.

More Answers (0)

Categories

Find more on Develop Apps Using App Designer in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!