## Hi.. this program contains no errors, I do not know why the compilaton is long !!!

on 18 Sep 2017
Latest activity Commented on by Walter Roberson

### Walter Roberson (view profile)

on 19 Sep 2017
clc
clear
format short
% Define the parameters
eps = 0.2; % Perturbation Parameter
% Define the number of grid points in x and y direction
% nx = number of grid points in x direction
% ny = number of grid points in y direction
nx = 101;
ny = 65;
t = (2.0 * cos(pi/(nx * ny)))^2;
% Calculate a optimum value of SOR parameter
omega1 = (16.0+sqrt((256.0-(64.0 * t))))/(2.0 * t);
omega2 = (16.0-sqrt((256.0-(64.0 * t))))/(2.0 * t);
oopt = min(omega1,omega2);
if ( (oopt <= 1.0) || (oopt >= 2.0 ))
oopt = 1.0;
end
omega = oopt;
% Define the dimension of the domain
Lx = 1.6; % Length in x of the computation region
Ly = 1.0; % Length in y of the computation region
% Calculate the mesh size
hx = Lx/(nx-1); % Grid spacing in x
hy = Ly/(ny-1); % Grid spacing in y
% Generate the mesh/grid
x=zeros(nx,1);
for i = 2:nx
x(i)=x(i-1)+hx;
end
y=zeros(ny,1);
for j = 2:ny
y(j)=y(j-1)+hy;
end
% Solve the Temperature equation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize the temperature field as zero
T = zeros(nx,ny);
% Max-Norm on Error {L-inf Error) Initialized
Linf = 1.0;
iteration = 0;
while (Linf > 1.0e-5)
iteration = iteration+1;
% Store the old values of T in Told
Told = T;
% Apply the boundary conditions
for i = 1:nx
% BC for the bottom boundary
T(i,1)=T(i,2)+hy;
% BC for the top boundary
if ( x(i) <= 0 )
T(i,ny) = 1;
else
T(i,ny) = 0;
end
end
for j = 1:ny
% BC for the Left boundary
T(1,j)=T(2,j);
% BC for the Left boundary
T(nx,j)=T(nx-1,j);
end
% Now compute the interior domain using 2nd order finite difference
for i = 2:nx-1
for j = 2:ny-1
term1 = ((eps/hx)^2) * (T(i-1,j)+T(i+1,j));
term2 = ((1.0/hy)^2) / (T(i,j-1)+T(i,j+1));
num = term1+term2;
den = 2.0 * (((eps/hx)^2)+((1.0/hy)^2));
Tgs = num/den;
T(i,j) = (omega * Tgs)+((1.0-omega) * Told(i,j));
end
end
% Calculate the error
Terr = abs(T-Told);
Linf = norm(Terr,2);
% Print the convergence history every 50 iterations
ccheck = round(iteration/100)-(iteration/100);
if ( (iteration == 1) || ( ccheck == 0) )
fprintf('%d \t %e \n', iteration, Linf)
end
end
% Plot the solutions
figure
[X,Y] = meshgrid(x,y);
clevel = [-1 -0.005 0 0.053 0.152 0.252 0.352 0.451 0.551 0.65 0.75 0.949 1.049 1.148 1.248 1.348 1.447 1.547 1.646 1.746 1.848 1.945 2.0 3.0];
contourf(X',Y',T,clevel)
colorbar
axis([-1.2 1.2 -0.8 1.8])
xlabel('x')
ylabel('y')
zlabel('T(x,y)')
hold on
yi = (0:0.01:1);
xi = zeros(length(yi));
plot(xi,yi,'black')
print -djpeg figures\chap4\temperature4

Show 1 older comment

on 18 Sep 2017
I mean running ...
OCDER

### OCDER (view profile)

on 18 Sep 2017
Oh ok. One last thing - highlight your question text, and push the {} Code button located above the text box you're using, so that we can see your code
Like this
for k = 1:100
end
Walter Roberson

### Walter Roberson (view profile)

on 19 Sep 2017
Are you asking why your code is taking a long time to converge? Or are you asking for hints on vectorizing your code?

### Tags

on 18 Sep 2017

I formatted the code to make it easier for folks to read.
Using the profiler for a few iterations indicates most of the time is spent inside the nested loops. I stopped the code on its 443rd time through the while loop, but by then, term1, term2, etc had already been calculated over 2.7 million times. I suggest rewriting the code ( vectorize) to remove these nested loops:
for i = 2:nx-1
for j = 2:ny-1
...
end
end

OCDER

### OCDER (view profile)

on 18 Sep 2017
Thanks Chad for formatting the code! Yassine, the profiler that Chad is suggesting is really helpful for optimizing code, so please do give the document a read. I use it a lot to find bottlenecks quickly, and then the tic/toc or timeit method for a quick check to see if a new code is faster/slower.

on 18 Sep 2017
Edited by OCDER

### OCDER (view profile)

on 18 Sep 2017

I'm getting that norm is the slowest step, follow closely by that nested for loop that Chad is pointing out. There's a sqrt in norm that is slowing things down. You can try rewriting the while loop to check Linf BEFORE the sqrt is taken, for instance:
while Linf2 > (1.0e-5)^2
....
Linf2 = sum(sum(T-Told).^2)); %Calc the norm^2 yourself, but careful about precision as it's not 100% the same as norm(T-Told, 2).
end
This sped things up 26 times for me, but then the nested for loop becomes the bottleneck. To be even faster, this article could be helpful: https://www.mathworks.com/company/newsletters/articles/programming-patterns-maximizing-code-performance-by-optimizing-memory-access.html