[mechie94]

Hello everyone,

I am developing a code in FORTRAN to solve 2d steady state heat conduction equation(Laplace equation) problem using Gauss Seidel iterative method.
Right now I am running the solution upto 1000 iterations and getting a solution.
However, I want the solution to stop iterating as soon as it converges.
I am using a three dimensional array to store the values of temparature. The main code snippet for gauss seidel iterations is given below.(k: iterations, i: x nodes, j: y nodes):

DO K = 2, 1000
  DO I = 2, 39
    DO J = 2, 19
       T(K, I, J) = 0.25*(T(K, I-1, J) + T(K-1, I+1, J) + T(K, I, J-1) + T(K-1, I, J+1))
    END DO
  END DO
END DO

Please suggest how can I stop this as the solution converges.

[FMDenaro]

You have to use a "while" construction, then check the residual of the equations at the k-th iteration and if it (in some norm) is less than your threshold you quite the loop.

[Simbelmynë]

Originally Posted by mechie94

Hello everyone,

I am developing a code in FORTRAN to solve 2d steady state heat conduction equation(Laplace equation) problem using Gauss Seidel iterative method.
Right now I am running the solution upto 1000 iterations and getting a solution.
However, I want the solution to stop iterating as soon as it converges.
I am using a three dimensional array to store the values of temparature. The main code snippet for gauss seidel iterations is given below.(k: iterations, i: x nodes, j: y nodes):

DO K = 2, 1000
  DO I = 2, 39
    DO J = 2, 19
       T(K, I, J) = 0.25*(T(K, I-1, J) + T(K-1, I+1, J) + T(K, I, J-1) + T(K-1, I, J+1))
    END DO
  END DO
END DO

Please suggest how can I stop this as the solution converges.

One option is to:

resid=0.d0 ! Residual (not in a strict sense though)
 DO K = 2, 1000
   DO I = 2, 39
     DO J = 2, 19
        TNew= 0.25*(T(K, I-1, J) + T(K-1, I+1, J) + T(K, I, J-1) + T(K-1, I, J+1))
        resid=max(resid,abs(TNew-T(k,i,j))
        T(k,i,j)=TNew
     END DO
   END DO
  END DO

 eps=resid/1000 ! Here you set the convergence criterion, eps, to reduce the initial residual by three orders of magnitude (this could be changed depending on accuracy).

resid=eps+1
do while (resid.gt.eps)
 resid=0.d0
 DO K = 2, 1000
   DO I = 2, 39
     DO J = 2, 19
        TNew = 0.25*(T(K, I-1, J) + T(K-1, I+1, J) + T(K, I, J-1) + T(K-1, I, J+1))
        resid=max(resid,abs(TNew-T(k,i,j))
        T(k,i,j)=TNew
     END DO
   END DO
  END DO
enddo

[FMDenaro]

It is required to evaluate the vector of the residuals that are homogenoeus at the same iteration level (rk=A.xk-q), otherwise the evaluation of the convergence can be degraded.

[Simbelmynë]

Originally Posted by FMDenaro

It is required to evaluate the vector of the residuals that are homogenoeus at the same iteration level (rk=A.xk-q), otherwise the evaluation of the convergence can be degraded.

Yes I now see that the OP had a 2D problem, where K is the iteration level (although it starts at K=2, which is really confusing). I suggest that the problem is trimmed to a 2d array of T, i.e. T(j,i) instead of maintaining information about the previous iterations.

If the previous iteration levels need to be kept, then perhaps a pre-allocation of the 3D array with size of K being equal to N² (N being the number of nodes in the system). This way it would scale somewhat similar to the Gauss-Seidel method.

[FMDenaro]

This is an m-file I wrote for my students for the GS method applied to the solution of the Laplace equation discretized in 2D with second order central formula. I used just the same 2D array.

____________________
clear;format long e
Ni=11;Lx=1.;dx=Lx/(Ni-1);x=(1:Ni)*dx-dx;dx2=dx*dx;
Nj=11;Ly=1.;dy=Ly/(Nj-1);y=(1:Nj)*dy-dy;dy2=dy*dy;
% coefficienti della matrice A
ac=-(2/dx2+2/dy2);acp=1/dx2;acm=acp;acpp=1/dy2;acmm=acpp;
Fovest=1;Fest=0;Fnord=0;Fsud=0;% Cond. contorno Dirichlet
for i=1:Ni phi(i,1)=Fsud; phi(i,Nj)=Fnord; end
for j=2:Nj-1 phi(1,j)=Fovest; phi(Ni,j)=Fest; end
% Vettore di tentativo iniziale x0,calcolo del residuo r0
iterazione=0;for i=2:Ni-1;for j=2:Nj-1;phi(i,j)=0;q(i,j)=0;
r(i,j)=ac*phi(i,j)+acp*phi(i+1,j)+acm*phi(i-1,j)+...
acpp*phi(i,j+1)+acmm*phi(i,j-1)-q(i,j);
end; end; rmax=norm(r,inf); eps_tolleranza=10^-7;
% Ciclo iterativo con metodo di Gauss-Seidel
while rmax > eps_tolleranza iterazione=iterazione+1
for i=2:Ni-1 for j=2:Nj-1 % Incognite interne
sum=acp*phi(i+1,j)+acm*phi(i-1,j)+acpp*phi(i,j+1)+...
acmm*phi(i,j-1);
phi(i,j)=(q(i,j)-sum)/ac;
end; end
for i=2:Ni-1 for j=2:Nj-1
r(i,j)=ac*phi(i,j)+acp*phi(i+1,j)+acm*phi(i-1,j)+...
acpp*phi(i,j+1)+acmm*phi(i,j-1)-q(i,j);
end; end; rmax=norm(r,inf);res(iterazione)=rmax;
end % fine ciclo while
semilogy(1:iterazione,res,'o'), xlabel('n. iterazione')...
,ylabel('|r|_i_n_f'): pause
contour(x,y,phi,16);shading interp;colorbar('vert');