MATLAB Answers

0

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

Asked by yassine amadane on 18 Sep 2017
Latest activity Commented on by Walter Roberson
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

  4 Comments

Show 1 older comment

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

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

Sign in to comment.

Tags

2 Answers

Answer by Chad Greene
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

  1 Comment

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.

Sign in to comment.


Answer by OCDER
on 18 Sep 2017
Edited by OCDER
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

  0 Comments

Sign in to comment.