Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Steady state Heat Conduction

  1. Aug 5, 2014 #1
    Hi guys,I am new for programming. I made my first attempt to solve an heat equation by finite difference method and wrote a code for it in Matlab.I got a solution but i need help from you guys to (vectorize) increase the speed of my program.

    ((∂^2 T)/(∂x^2 )+(∂^2 T)/(∂y^2 )+(∂^2 T)/(∂z^2 ))(k)=-Q(x,y,z)

    clear;
    close all;
    clc;
    n = 5;% n grids and has n - 2 interior points per dimension

    x = linspace(0,1,n);% line in x direction of 0 to 1 and divided into n parts

    dx = x(2)-x(1); % distance between grids

    y = x;

    dy = dx;

    k=1; %thermal conductivity

    Q=300 ; % heat source
    m=(n-2)^2;

    %Right hand side matrix

    B=zeros(m,1);

    if (mod(m,2)==0)
    B(m/2,1)=Q;
    else
    B((m+1)/2,1)=Q;
    end

    %Main Matrix

    Tsol = zeros(n);
    %boundry conditions
    Tsol(1,1:n) = 0; %TOP
    Tsol(n,1:n) = 0; %BOTTOM
    Tsol(1:n,1) = 0; %LEFT
    Tsol(1:n,n) = 0; %RIGHT

    tic
    error = 1; z = 0;
    while error > 0.000001
    z = z+1;
    c=1;
    Torg = Tsol;
    for i = 2:n-1
    for j = 2:n-1
    Tsol(i,j) =abs(((dx^2*B(c,1))+((Torg(i+1,j)+Torg(i-1,j)+Torg(i,j+1)+Torg(i,j-1))*k))/(4*k));
    c=c+1;
    end
    end

    error = max(max(abs(Torg-Tsol)));
    end


    toc

    %% plotting

    figure, contour(x,y,Tsol),
    title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar



    figure,pcolor(x,y,Tsol),shading interp,
    title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar
     
  2. jcsd
  3. Aug 6, 2014 #2

    kreil

    User Avatar
    Gold Member

    You initialize Tsol as a matrix of zeros, yet the next 4 lines set different boundary values to 0? Seems unnecessary.

    Can't you just use del2() or gradient() here?
     
  4. Aug 6, 2014 #3

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    "Initializing the solution" and "setting the boundary conditions" are two different steps in the solution algorithm. For this particular equation they happen to have the same numerical effect (setting some terms to zero twice) but optimizing that away would have a negligible effect on the total execution time.

    To optimize the code, you first have to find out what is taking the time. For a small problem (n = 5) you only have about (n-2)2 = 9 calculations of the elements of Tsol, in each iteration. So if your program runs for a long time (e.g. more than 1 second) you must be doing a huge number of iterations. Try counting the number of times the "while error > 0.000001" loop runs, and think about how to reduce that to a more "reasonable" number (e.g. 100 or less).

    Converting the "Tsol" calculation to an array operation will help, but if you want to optimize a program effectively, you first need to find where the real bottlenecks are, not just optimize bits of code "because you can".
     
  5. Aug 6, 2014 #4
    I have already checked my program and found that my ittirations are going way beond than 100000,it is the reason why it is taking long time. I am trying a lot to vectorize my while loop but i am not able to do because of my lack of knowledge and experience in programming skills.Thats why i needed a helping hand,i will be glad to take your help in my first experience.

    error = 1; z = 0;
    while error > 0.000001
    z = z+1;
    c=1;
    Torg = Tsol;
    for i = 2:n-1
    for j = 2:n-1
    Tsol(i,j) =abs(((dx^2*B(c,1))+((Torg(i+1,j)+Torg(i-1,j)+Torg(i,j+1)+Torg(i,j-1))*k))/(4*k));
    c=c+1;
    end
    end

    error = max(max(abs(Torg-Tsol)));
    end
     
  6. Aug 6, 2014 #5

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    Print out value of "error" in the first terations (say the first 100). Is it converging monotonically, or oscillating, or jumping about "randomly"?
     
  7. Aug 6, 2014 #6
    the error is exponentially decreasing and some times it is raising up at certain values of n,like you can see in the images down below at n = 25.
     

    Attached Files:

  8. Aug 6, 2014 #7

    AlephZero

    User Avatar
    Science Advisor
    Homework Helper

    Very strange.

    I don't understand the first value of "error". On the first iteration everything in the calculation of Tsol should be 0, apart from one value in the B array which is 300. dx = 0.25. So the only non-zero value in tsol should be
    300 * (0.25)^2 / 4 = about 4.7. But on your plots it is about 0.67, not 4.7.

    Maybe I misunderstood what your code does, or maybe you are not running the same code that you posted here.

    Incidentally, why is the abs(...) function in the calculation of tsol?
     
  9. Aug 6, 2014 #8

    kreil

    User Avatar
    Gold Member

    OK, I invested some time and looked in to this.

    While it's easy to vectorize the for loops in a general sense, you were using the [itex]c[/itex] variable as a sort of delta function to pick out the value from [itex]B[/itex] for the heat source. This means that your loops were doing nothing for the first [itex]m/2[/itex] iterations, where the heat source is finally picked up and starts to propagate. So you can get rid of those [itex]m/2[/itex] unnecessary iterations by just initializing the heat source and adding it in differently. I did this by reshaping [itex]B[/itex] into the interior of a matrix, [itex]A[/itex]. Since [itex]A[/itex] is the same size as [itex]Tsol[/itex], you can use it to initialize [itex]Tsol[/itex] and then add in the heat source after each iteration. This means there are no wasted iterations; the very first one is productive.

    Using MATLAB R2014a with n=25, the time differences on my laptop were 0.101265 seconds for your code, and 0.038494 seconds for the code below, so the improvement is significant. I believe this could get even more dramatic for larger n.

    Hopefully this helps!

    Note: While I still think del2 should be useful here (since it calculates vectorized finite differences via bsxfun), it uses single sided differences on the edges of the matrix, which messes with your boundary conditions.

    Code (Text):

    n = 25;% n grids and has n - 2 interior points per dimension
    x = linspace(0,1,n);% line in x direction of 0 to 1 and divided into n parts
    dx = x(2)-x(1); % distance between grids
    y = x;
    dy = dx;
    k=1; %thermal conductivity
    Q=300 ; % heat source
    m=(n-2)^2;

    %Right hand side matrix
    B=zeros(m,1);
    if (mod(m,2)==0)
    B(m/2,1)=Q;
    else
    B((m+1)/2,1)=Q;
    end

    %% Main Matrix

    %Initialize a heat source matrix
    A = zeros(n);
    A(2:n-1,2:n-1) = reshape(B,n-2,n-2);
    A = (dx^2*A)/(4*k);

    Tsol = A;
    %boundary conditions
    Tsol(1,1:n) = 0; %TOP
    Tsol(n,1:n) = 0; %BOTTOM
    Tsol(1:n,1) = 0; %LEFT
    Tsol(1:n,n) = 0; %RIGHT

    tic
    error = 1;
    z = 0;
    i = 2:n-1;
    j = 2:n-1;
    while error > 0.000001
        z = z+1;
        Torg = Tsol;
        Tsol(i,j) = abs(((Torg(i+1,j)+Torg(i-1,j)+Torg(i,j+1)+Torg(i,j-1))*k))/(4*k);
        Tsol = Tsol + A;
        error = max(max(abs(Torg-Tsol)));
    end
    toc

    figure, contour(x,y,Tsol),
    title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar
    figure,pcolor(x,y,Tsol),shading interp,
    title('Temperature (Steady State)'),xlabel('x'),ylabel('y'),colorbar
     
     
    Last edited: Aug 6, 2014
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Steady state Heat Conduction
Loading...