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

`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

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

You mean running it? or compiling it into a binary file? Compiling could be slow slow as matlab has to check for dependencies, convert, etc.

Also, can you do us big favor and hit that

{} Codebutton, and explain how long it takes to compile?## 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

I mean running ...

## 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

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

## 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

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.