Solution of the 2D-Laplace's Equation (Heat diffusion) with Finite Difference (Jacobi-IM, GS-IM, PSOR) FORTRAN CODE

 The steady-state temperature distribution over a 1m x 1m square plate for different boundary conditions are found out by solving the simple Laplace equation using the finite difference method.

 

In this study, our main objective is to solve the Laplace equation of heat diffusion for a 2D square solid domain. This is a very basic and benchmark problem in the field of computational heat transfer.
This computational problem is solved with the help of numeric simulation technique written in FORTRAN 95.

 

Source Code: (Fortran 95)

PROGRAM ELLIPTICSOLVER_NMN_CASE3

IMPLICIT NONE

REAL, DIMENSION(:,:), ALLOCATABLE  ::X,Y,TOLD,TNEW

REAL ::DX,DY,R

INTEGER ::NX,NY,METHOD

NX=40                       !INPUT

NY=40                       !INPUT

ALLOCATE(X(NX,NY),Y(NX,NY),TOLD(NX,NY),TNEW(NX,NY))

 

CALL GRID(X,Y,NX,NY,DX,DY)

R=(DX/DY)

TOLD=0.0

TNEW=0.0

CALL BOUNDARYCONDITION(TOLD,NX,NY,X)

TNEW=TOLD

METHOD=1                    !INPUT 1=GAUSSSEIDEL, 2=JACOBI, 3=PSOR

 

SELECT CASE(METHOD)

    CASE(1)

            CALL GAUSSSEIDEL(TOLD,TNEW,NX,NY,R,X)           

    CASE(2)

            CALL JACOBI(TOLD,TNEW,NX,NY,R,X)

    CASE(3)

            CALL PSOR(TOLD,TNEW,NX,NY,R,X)

    CASE DEFAULT

            CALL GAUSSSEIDEL(TOLD,TNEW,NX,NY,R,X)

END SELECT

 

CALL DATAWRITER(NX,NY,X,Y,TNEW)

 

END PROGRAM ELLIPTICSOLVER_NMN_CASE3

 

 

SUBROUTINE GRID(X,Y,NX,NY,DX,DY)

IMPLICIT NONE

INTEGER, INTENT(IN) ::NX,NY

REAL, DIMENSION(NX,NY) ::X,Y

REAL  :: LX,LY

REAL, INTENT(OUT) ::DX,DY

INTEGER I,J

LX=1

LY=1

DX=LX/REAL(NX-1)

DY=LY/REAL(NY-1)

    DO J=1,NY

        DO I=1,NX

            X(I,J)=(I-1)*DX

            Y(I,J)=(J-1)*DY

        END DO

    END DO

 

END SUBROUTINE GRID

 

 

SUBROUTINE BOUNDARYCONDITION(TOLD,NX,NY,X)

INTEGER ,INTENT(IN) ::NX,NY

REAL, DIMENSION(NX,NY)::TOLD,X

REAL :: TRW

INTEGER I,J

TRW=0                 

DO I=1,NX

    DO J=1,NY

        IF(J==1) THEN

        TOLD(I,J)=100*(1+X(I,J))

        ELSE IF(I==NX) THEN

        TOLD(I,J)=TRW

        ELSE

        TOLD(I,J)=0.0

        END IF

    END DO

END DO

 

END SUBROUTINE BOUNDARYCONDITION

 

 

SUBROUTINE BOUNDARYUPDATE(TNEW,NX,NY,X)

INTEGER, INTENT(IN) ::NX,NY

REAL, DIMENSION(NX,NY)::TNEW,X

REAL ::TRW

INTEGER I,J

!TTW=0 ADIABATIC

TRW=0

!TLW=0 ADIABATIC

!TBW=100

 

DO I=1,NX

    DO J=1,NY

        IF(I==1) THEN

        TNEW(I,J)=TNEW(I+1,J)

        ELSE IF(J==1) THEN

        TNEW(I,J)=100*(1+X(I,J))

        ELSE IF(I==NX) THEN

        TNEW(I,J)=TRW

        ELSE IF(J==NY) THEN

        TNEW(I,J)=TNEW(I,J-1)

        END IF

    END DO

END DO

!CORNER TREATMENT

TNEW(1,1)=(TNEW(1,2)+TNEW(2,1))/2  

TNEW(1,NY)=(TNEW(1,NY-1)+TNEW(2,NY))/2

TNEW(NX,1)=(TNEW(NX,2)+TNEW(NX-1,1))/2

TNEW(NX,NY)=(TNEW(NX,NY-1)+TNEW(NX-1,NY))/2 

 

END SUBROUTINE BOUNDARYUPDATE

 

 

SUBROUTINE DATAWRITER(NX,NY,X,Y,TNEW)

INTEGER,INTENT(IN) ::NX,NY

REAL, DIMENSION(NX,NY)::TNEW,X,Y

INTEGER I,J

OPEN(20,FILE="DATA.txt")

    DO J=1,NY

        DO I=1,NX

            WRITE(20,*)X(I,J),Y(I,J),TNEW(I,J)

        END DO

    END DO

CLOSE(20)

END SUBROUTINE DATAWRITER

 

!######################################## Gauss seidel ###########################

SUBROUTINE GAUSSSEIDEL(TOLD,TNEW,NX,NY,R,X)

INTEGER, INTENT(IN) ::NX,NY

REAL,DIMENSION(NX,NY) ::X

REAL,DIMENSION(NX,NY) ::TOLD,TNEW

REAL, INTENT(IN) ::R

REAL ::ERROR,ERRORMAX

INTEGER ::COUNTER,I,J

 

COUNTER=0

ERROR=0.0

ERRORMAX=0.00001

 

 

DO

    ERROR=0.0

        DO I=2,NX-1

            DO J=2,NY-1

           

            TNEW(I,J)=(TOLD(I+1,J)+TNEW(I-1,J)+ (R*R)*(TOLD(I,J+1)+TNEW(I,J-1)))/(2.0+2.0*(R*R))

            ERROR= ERROR + (TNEW(I,J)-TOLD(I,J))**2

           

            END DO

        END DO

   

    CALL BOUNDARYUPDATE(TNEW,NX,NY,X)   

    TOLD=TNEW

    COUNTER=COUNTER+1

    PRINT *,SQRT(ERROR)/(NX*NY)

    IF (SQRT(ERROR)/(NX*NY)<ERRORMAX) EXIT

END DO

 

END SUBROUTINE GAUSSSEIDEL

 

!######################################## Jacobi ############################

SUBROUTINE JACOBI(TOLD,TNEW,NX,NY,R,X)

INTEGER, INTENT(IN) ::NX,NY

REAL,DIMENSION(NX,NY) ::X

REAL,DIMENSION(NX,NY) ::TOLD,TNEW

REAL, INTENT(IN) ::R

REAL ::ERROR,ERRORMAX

INTEGER ::COUNTER,I,J

 

COUNTER=0

ERROR=0.0

ERRORMAX=0.00001

 

 

DO

    ERROR=0.0

        DO I=2,NX-1

            DO J=2,NY-1

           

            TNEW(I,J)=(TOLD(I+1,J)+TOLD(I-1,J)+ (R*R)*(TOLD(I,J+1)+TOLD(I,J-1)))/(2.0+2.0*(R*R))

            ERROR= ERROR + (TNEW(I,J)-TOLD(I,J))**2

           

            END DO

        END DO

   

    CALL BOUNDARYUPDATE(TNEW,NX,NY,X)   

    TOLD=TNEW

    COUNTER=COUNTER+1

    PRINT *,SQRT(ERROR)/(NX*NY)

    IF (SQRT(ERROR)/(NX*NY)<ERRORMAX) EXIT

END DO

 

END SUBROUTINE JACOBI

 

!######################################## POINT SOR ###############################

SUBROUTINE PSOR(TOLD,TNEW,NX,NY,R,X)

INTEGER, INTENT(IN) ::NX,NY

REAL,DIMENSION(NX,NY) ::X

REAL,DIMENSION(NX,NY) ::TOLD,TNEW

REAL, INTENT(IN) ::R

REAL ::ERROR,ERRORMAX,W

INTEGER ::COUNTER,I,J

 

COUNTER=0

ERROR=0.0

ERRORMAX=0.000000001

W=1.5

 

DO

    ERROR=0.0

        DO I=2,NX-1

            DO J=2,NY-1

           

            TNEW(I,J)=(1-W)*TOLD(I,J)+(W*(TOLD(I+1,J)+TNEW(I-1,J)+ (R*R)*(TOLD(I,J+1)+TNEW(I,J-1))))/(2.0+2.0*(R*R))

           

            ERROR= ERROR + (TNEW(I,J)-TOLD(I,J))**2

           

            END DO

        END DO

   

    CALL BOUNDARYUPDATE(TNEW,NX,NY,X)   

    TOLD=TNEW

    COUNTER=COUNTER+1

    PRINT *,SQRT(ERROR)/(NX*NY)

    IF (SQRT(ERROR)/(NX*NY)<ERRORMAX) EXIT

END DO

 

END SUBROUTINE PSOR

 

Results: