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