Asked by yassine amadane
on 18 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

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

OCDER
on 18 Sep 2017

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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 4 Comments

## OCDER (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/357252-hi-this-program-contains-no-errors-i-do-not-know-why-the-compilaton-is-long#comment_485466

## yassine amadane (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/357252-hi-this-program-contains-no-errors-i-do-not-know-why-the-compilaton-is-long#comment_485469

## OCDER (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/357252-hi-this-program-contains-no-errors-i-do-not-know-why-the-compilaton-is-long#comment_485472

## Walter Roberson (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/357252-hi-this-program-contains-no-errors-i-do-not-know-why-the-compilaton-is-long#comment_485527

Sign in to comment.