Fortran Successive Over Relaxation Method in FORTRAN

[Aun Muhammad]

Hey everyone.
I've attached a code I have written in FORTRAN implementing the SOR Method for a 2 D Laplace Equation.
Subroutines have been created, grid has been created and initial and BC's applied appropriately.
While computing the error of each iteration, I am unable to achieve the desired graph. In the graphs I obtained, the error after a few iterations, instantly dives to a very low value and remains there. However, the graph must display a lowest value of error for an 'optimum value' of omega.
I'll be grateful if somebody with FORTRAN experience could look into this and kindly point out where the error is being made. I've tried it a hundred times now and can't find where the error lies. Thanks.

 

[DrClaude]

Next time, please post your code directly, not as an attachment, and use code tags:

PROGRAM ELLIPTIC_SOLVER
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE

!
!*************            VARIABLES AND PARAMETERS            *************!
!

INTEGER(I4B)::i,j,k,l
INTEGER(I4B),PARAMETER::N=21,MAXITRS=20,maxw=501
REAL(DP)::DX,DY,w,ERROR
REAL(DP), PARAMETER::deltaw=0.002
REAL(DP),DIMENSION(:,:),ALLOCATABLE::X,Y,S,U0,UOLD,UNEW,R
ALLOCATE(X(N,N),Y(N,N),S(N,N),U0(N,N),UOLD(N,N),UNEW(N,N),R(N,N))

!
!*************                   GRID SIZING                  *************!
!

CALL GRID(X,Y,DX,DY,N)
    WRITE(*,*)'The Matrix X(i,j) is given as:'
    DO i=LBOUND(X,1),UBOUND(X,1)
        WRITE(*,*) (X(i,j), j=LBOUND(X,2),UBOUND(X,2))
    END DO
   
    WRITE(*,*)'The Matrix Y(i,j) is given as:'
    DO i=LBOUND(Y,1),UBOUND(Y,1)
        WRITE(*,*) (Y(i,j), j=LBOUND(Y,2),UBOUND(Y,2))
    END DO
   
!
!*************               ANALYTICAL SOLUTION              *************!
!

DO j=1,N
    DO i=1,N
        U0(i,j)=(X(i,j))**2-(Y(i,j))**2
    END DO
END DO

    WRITE(*,*)'The Analytical Solution of U0(i,j) is given as:'
    DO i=LBOUND(U0,1),UBOUND(U0,1)
        WRITE(*,*)(U0(i,j), j=LBOUND(U0,2),UBOUND(U0,2))
    END DO
   
!
!*************         INITIAL AND BOUNDARY CONDITIONS        *************!
!

DO i=1,N
    S(i,1)=U0(i,1)              !Boundary Condition from Analytical Solution
    S(i,N)=U0(i,N)              !Boundary Condition from Analytical Solution
END DO

DO j=1,N
    S(1,j)=U0(1,j)              !Boundary Condition from Analytical Solution
    S(N,j)=U0(N,j)              !Boundary Condition from Analytical Solution
END DO

DO i=2,N-1
    DO j=2,N-1
        S(i,j)=0.0                               !Initial-Value from Problem
    END DO
END DO

DO i=1,N
    DO j=1,N
        UOLD(i,j)=S(i,j)           !Equating Values from S(i,j) to UOLD(i,j)
    END DO
END DO

    WRITE(*,*)'The Initial Value Matrix UOLD(i,j) is given as:'
    DO i=LBOUND(UOLD,1),UBOUND(UOLD,1)
      WRITE(*,*) (UOLD(i,j), j=LBOUND(UOLD,2),UBOUND(UOLD,2))
    END DO
   
!
!*************               NUMERICAL SOLUTION               *************!
!

DO i=1,N                                !Boundary Conditions Remain Constant
    UNEW(i,1) = UOLD(i,1)
    UNEW(i,N) = UOLD(i,N)
END DO

DO j=1,N       
    UNEW(1,j) = UOLD(1,j)
    UNEW(N,j) = UOLD(N,j)
END DO

UNEW(i,j)=UOLD(i,j)
ERROR=0.0
w=1.0
DO l=1,maxw
    DO k=1,MAXITRS
        DO j=2,N-1
            DO i=2,N-1
                CALL PSOR(UOLD,UNEW)
                    CALL ERROR_CALCULATION(UOLD,UNEW,ERROR)
                        CALL RR(UNEW,U0,R)
                          
                UOLD(i,j)=UNEW(i,j)
            END DO
        END DO
         WRITE(*,*)'The Initial Value Matrix R(i,j) is given as:'
                            DO i=LBOUND(R,1),UBOUND(R,1)
                            WRITE(*,*) (R(i,j), j=LBOUND(R,2),UBOUND(R,2))
                            END DO
    END DO
    WRITE(*,*)w,ERROR
    ERROR=0.0
    w=w+deltaw
END DO

CONTAINS

!
!*************                   SUBROUTINES                  *************!
!

SUBROUTINE GRID(X,Y,DX,DY,N)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
    INTEGER(I4B)::i,j
    REAL(DP), DIMENSION(N,N)::X,Y
    INTEGER(I4B),INTENT(IN)::N
    REAL(DP),INTENT(OUT)::DX,DY
   
    DX=1.0/(FLOAT(N)-1.0)
    DY=1.0/(FLOAT(N)-1.0)
   
    DO j=1,N
        DO i=1,N
            X(i,j)=(i-1)*DX
            Y(i,j)=(j-1)*DY
        END DO
    END DO
END SUBROUTINE GRID

!

SUBROUTINE PSOR(UOLD,UNEW)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
    REAL(DP)::SUM1
    REAL(DP),DIMENSION(N,N),INTENT(IN)::UOLD
    REAL(DP),DIMENSION(N,N),INTENT(OUT)::UNEW

SUM1=0.25*(UOLD(i+1,j)+UNEW(i-1,j)+UOLD(i,j+1)+UNEW(i,j-1)-(4.0*UOLD(i,j)))
UNEW(i,j)=UOLD(i,j)+(w*SUM1)
END SUBROUTINE PSOR

!

SUBROUTINE ERROR_CALCULATION(UOLD,UNEW,ERROR)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
    REAL(DP),DIMENSION(N,N),INTENT(IN)::UOLD,UNEW
    REAL(DP),INTENT(OUT)::ERROR

    ERROR=ERROR+ABS(UNEW(i,j)-UOLD(i,j))
END SUBROUTINE ERROR_CALCULATION

SUBROUTINE RR(UNEW,U0,R)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
    REAL(DP),DIMENSION(N,N),INTENT(IN)::U0,UNEW
    REAL(DP),DIMENSION(N,N),INTENT(OUT)::R
R(i,j)=U0(i,j)-UNEW(i,j)
END SUBROUTINE RR

END PROGRAM ELLIPTIC_SOLVER

 

[DrClaude]

Not related to your problem, but the following code is extremely inefficient:

DO i=2,N-1
    DO j=2,N-1
        S(i,j)=0.0                               !Initial-Value from Problem
    END DO
END DO

DO i=1,N
    DO j=1,N
        UOLD(i,j)=S(i,j)           !Equating Values from S(i,j) to UOLD(i,j)
    END DO
END DO

2D arrays are stored in column-major order in Fortran, so the leftmost index such change first. And since you are using Fortran 90, you should simply write

S = 0.0 !Initial-Value from Problem

DO i=1,N
    S(i,1)=U0(i,1)              !Boundary Condition from Analytical Solution
    S(i,N)=U0(i,N)              !Boundary Condition from Analytical Solution
END DO

DO j=1,N
    S(1,j)=U0(1,j)              !Boundary Condition from Analytical Solution
    S(N,j)=U0(N,j)              !Boundary Condition from Analytical Solution
END DO

UOLD = S           !Equating Values from S(i,j) to UOLD(i,j)

 

[DrClaude]

There is a problem with this here:

SUBROUTINE ERROR_CALCULATION(UOLD,UNEW,ERROR)
USE NRTYPE; USE NRUTIL, ONLY:NRERROR
IMPLICIT NONE
    REAL(DP),DIMENSION(N,N),INTENT(IN)::UOLD,UNEW
    REAL(DP),INTENT(OUT)::ERROR

    ERROR=ERROR+ABS(UNEW(i,j)-UOLD(i,j))
END SUBROUTINE ERROR_CALCULATION

ERROR is declared with an intent of out, which means that its value on entry is discarded. Therefore it is executable statement there is equivalent to

    ERROR = ABS(UNEW(i,j)-UOLD(i,j))

 

[Aun Muhammad]

I've made the necessary changes.
However, the results are still the same.
I'm plotting 'omega' versus 'logError'. I'm unable to find wopt(optimum) from the graph.
The blue-lined graph below shows my results.
The red-lined graph is the desired output.

 

[DrClaude]

I can't comment much because I am not familiar with the method you are using. But looking at your graph, the values are in the range 10⁴ to 10^-12. That's 16 orders of magnitude, which is basically equal to the machine epsilon for a double precision real.

This might be a case where you can't do any better because you are running out of decimal places.

 

[Aun Muhammad]

What you are suggesting seems plausible.
The results however must display a general trend similar to the desired result. The Error must decrease, reach a minima and then begin to increase again. Sadly that is not true in my case.

 

[eys_physics]

Hey everyone.
I've attached a code I have written in FORTRAN implementing the SOR Method for a 2 D Laplace Equation.
Subroutines have been created, grid has been created and initial and BC's applied appropriately.
While computing the error of each iteration, I am unable to achieve the desired graph. In the graphs I obtained, the error after a few iterations, instantly dives to a very low value and remains there. However, the graph must display a lowest value of error for an 'optimum value' of omega.
I'll be grateful if somebody with FORTRAN experience could look into this and kindly point out where the error is being made. I've tried it a hundred times now and can't find where the error lies. Thanks.

Hey, I see the following the problems:
1) Your code isn't computing the error of each iterations, as you state. Instead you are summing up the errors of all iterations. That is, ERROR = Error after one iteration + Error after two iterations + ... . In order to solve this you need to put "Error =0.0d0" after the line DO k=1,MAXITRS . Also don't use "Error = 0.0", since "0.0" is in single-precision, so a lot of precision is lost.
2) According to your formula you should divide by (N-2)*(N-2), but I cannot see that it is done in the code.

 

[Aun Muhammad]

Thank you for your reply.
I've changed the value from single to double precision as you stated. But could you please reiterate what the formula attached below actually means?
In my understanding,
Omega, or w, is being varied from 1 to 2 with a step-size of 0.002. This makes 501 iterations.
For every value of w, I'm carrying out 20 iterations for the calculation of UNEW.
Within these 20 iterations, i and j are varied from 2 to N-1 for the whole grid.
I'm trying to calculate the value of error for every value of w. Therefore, the sum of errors for the 20 iterations has to be summed up right?

 

[eys_physics]

Hey, sorry for the delay of my reply. No, in my understanding the formula isn't saying that.
Written out it is $\mathrm{Error} = \frac{\sum_{ij}|U_{i,j}-S_{i,j}|}{(N-2)*(N-2))}$, where $U$ is the exact solution and $V$ is the approximate one.
In other words, for a fixed number of iterations, the error is computed as the average error over the considered mesh in $x$ and $y$.
Typically, you should therefore compute the error after in your case 20 iterations. If you want to check convergence with respect to the number of iterations, you keep #w# fixed and compare the errors for let's say 19 and 20 iterations. However, in none of the cases you sum up the the errors.
If you do so, the total error will be dominated by the error coming from the first iterations (where errors are large).