CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   Main CFD Forum (https://www.cfd-online.com/Forums/main/)
-   -   Convergence of Gauss Seidel Iterative Method (https://www.cfd-online.com/Forums/main/193227-convergence-gauss-seidel-iterative-method.html)

mechie94 September 19, 2017 15:49
Convergence of Gauss Seidel Iterative Method
 
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):

Code:
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 September 19, 2017 15:55
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ė September 20, 2017 03:02
Quote:
Originally Posted by mechie94 (Post 664942)
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):

Code:
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:

Code:


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 September 20, 2017 07:59
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ė September 20, 2017 09:42
Quote:
Originally Posted by FMDenaro (Post 665026)
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 September 20, 2017 11:35
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');


All times are GMT -4. The time now is 09:27.