Fortran ERROR

M_H215 · November 14, 2020, 05:03

Hello
I am trying to write a code for Lid-Driven-Cavity in FORTRAN (Very short Code). Actually, I am new with programming. The problem is, when I run the code Software apparently do the job so fast and the result created but there isn't any result...!!
I mean that the PLT format file created but when I want to open it with tecplot, it gives me an error..!!
when I check the result with Note-Pad , All of the data are Zero...!!

Any help will be appreciated...

the code is :
..............................

Program Lid

implicit none
Integer :: I,J,nx, ny, L, W, Iteration, Max_Iteration , Re, M, N, dt
Real :: Delta, dx, dy
Real, allocatable :: u(:,:-), v(:,:-), p(:,:-), u_old(:,:-), v_old(:,:-), p_old(:,:-), X(:-), Y(:.)
!************************************************* **!

PRINT *, "ENTER THE DESIRED POINTS ..."
PRINT *, "... IN X DIRECTION: SUGGESTED RANGE (20-200)"
READ*, M
PRINT *, "... IN Y DIRECTION: SUGGESTED RANGE (10-100)"
READ*, N

! Define Geometry
dt = 0.001
Delta = 2
Re = 100
L = 10
W = 10
dx = L /(M-1)
dy = W /(N-1)

ALLOCATE (X(M),Y(N),u(M,N),u_old(M,N),v(M,N),v_old(M,N),p(M ,N),p_old(M,N))

! Grid Generation

Do I = 1, M
x(I) = (I-1)* dx
End Do

Do J=1 , N
y(J) = (J-1) * dy
End Do

! Boundray Condition
Do I=1 , M
u(I,1) = 0
u(1,I) = 0
u(N,I) = 0
u(I,N) = 1 ! Lid Velocity
End Do

Do J=1, N
v(J,1) = 0
v(1,J) = 0
v(J,M) = 0
v(M,J) = 0
End Do

! Initialization

u = 0; v = 0 ; p = 0

! Solver
Do I=2, M-1
Do J=2, N-1
u_old(I,J) = u(I,J)
v_old(I,J) = v(I,J)
p_old(I,J) = p(I,J)

u(I,J) = - dt / 4* dx * (( u(I, J+1)+ u_old(I,J))**2 - (u_old(I, J)+u(I,J-1))**2) - dt / 4* dy * ((u_old(I,J)+ u(I-1,J)) &
* (v(I-1,J) + v(I-1, J+1)) - (u_old(I,J) + u(I+1,J)) * (v_old(I,J) + v(I,J+1))) - dt / dx *(p(I, J+1) - p_old(I,J)) &
+ dt / Re * ((u(I+1,J) - 2 * u_old(I,J) + u(I-1,J)) / dx**2 + (u(I,J+1) - 2 * u_old(I,J) + u(I,J+1)) / dy**2) + u_old(I,J)

v(I,J) = - dt / 4* dy * (( v(I-1, J)+ v(I-1,J+1))**2 - (v_old(I, J)+v(I,J+1))**2) - dt / 4* dx * ((u_old(I,J)+ u(I,J+1)) &
* (v(I,J+1) + v(I-1, J+1)) - (u_old(I,J) + u(I,J-1)) * (v_old(I,J) + v(I-1,J))) - dt / dy *(p(I, J+1) - p(I,J)) &
+ dt / Re * ((v(I+1,J) - 2 * v_old(I,J) + v(I-1,J)) / dx**2 + (v(I,J+1) - 2 * v_old(I,J) + v(I,J+1)) / dy**2) + v_old(I,J)

p(I,J) = - Delta * dt / 2 * ((u(I,J+1)+ u_old(I,J)) - (u_old(I,J) + u(I,J-1))) - Delta * dt / 2 &
* ((v(I-1,J)+ v(I-1,J+1)) - (v_old(I,J) + v(I,J+1)))

End Do

End Do

!-----------------------OUTPUTS GENERATION-----------------------------
OPEN (1,FILE='FIELD.PLT')
WRITE (1,*) 'VARIABLES=X,Y,u,v,p'
WRITE (1,*) 'ZONE I=',M,' J=',N
DO J=1,N
DO I=1,M
WRITE (1,*) X(I),Y(J),u(I,J),v(I,J),p(I,J)
END DO
END DO

End Program Lid

flotus1 · November 14, 2020, 05:21

There are a few questionable pieces of code in here.
But for the immediate problem at hand: the initialization over-writes the boundary condition.

Next time you want to post code here, try the code format. This way, your code will look like this:

Code:

REAL, DIMENSION(:,:), ALLOCATABLE :: u
do i=1, N
    u(i)=0.0
end do

instead of this
REAL, DIMENSION(:,

, ALLOCATABLE :: u
do i=1, N
u(i)=0.0
end do

M_H215 · November 14, 2020, 06:29

Thank you so much....

Actually I tried to send my code like something that you mentioned but I could not find a way to do that...

How I can send my code like something that you Mentioned ? (Like Code)

Also, Could you please tell me that how I can prevent this issue ? (initialization over-written my boundaries)

Thank you

M_H215 · November 14, 2020, 06:44

I changed the Initialization section (Add a Do Loop) in order to prevent the over-written but still that problem exist...Do you think still over-written happen? (Fortunately I found how I can send my post like a code)

Code:

Program Lid
    
    implicit none
    Integer :: I,J,nx, ny, L, W, Iteration, Max_Iteration , Re, M, N, dt
    Real :: Delta, dx, dy
    Real, allocatable :: u(:,:), v(:,:), p(:,:), u_old(:,:), v_old(:,:), p_old(:,:), X(:), Y(:)
    !***************************************************!

    PRINT *, "ENTER THE DESIRED POINTS ..." 
    PRINT *, "... IN X DIRECTION:  SUGGESTED RANGE (20-200)" 
    READ*, M 
    PRINT *, "... IN Y DIRECTION:  SUGGESTED RANGE (10-100)" 
    READ*, N 
    
    ! Define Geometry
    dt = 0.001
    Delta = 2
    Re = 100
    L = 10
    W = 10
    dx = L /(M-1)
    dy = W /(N-1)
    
    ALLOCATE (X(M),Y(N),u(M,N),u_old(M,N),v(M,N),v_old(M,N),p(M,N),p_old(M,N))
    
    ! Grid Generation
    
    Do I = 1, M
        x(I) = (I-1)* dx
    End Do
    
    Do J=1 , N
        y(J) = (J-1) * dy
    End Do
     
    ! Boundray Condition
    Do I=1 , M
        u(I,1) = 0
        u(1,I) = 0
        u(N,I) = 0
        u(I,N) = 1       ! Lid Velocity
    End Do
    
    Do J=1, N
        v(J,1) = 0
        v(1,J) = 0
        v(J,M) = 0
        v(M,J) = 0
    End Do
    
    ! Initialization
    Do I= 2, M-1
        Do J=2, N-1
            u(I,J) = 0
            v(I,J) = 0
            p(I,J) = 0
            
        End do
    End do
                
    ! Solver
        Do I=2, M-1
            Do J=2, N-1
                u_old(I,J) = u(I,J) 
                v_old(I,J) = v(I,J)
                p_old(I,J) = p(I,J)
            
                u(I,J) = - dt / 4* dx * (( u(I, J+1)+ u_old(I,J))**2 - (u_old(I, J)+u(I,J-1))**2) - dt / 4* dy * ((u_old(I,J)+ u(I-1,J)) &
                    * (v(I-1,J) + v(I-1, J+1)) - (u_old(I,J) + u(I+1,J)) * (v_old(I,J) + v(I,J+1))) - dt / dx *(p(I, J+1) - p_old(I,J)) &
                    + dt / Re * ((u(I+1,J) - 2 * u_old(I,J) + u(I-1,J)) / dx**2 + (u(I,J+1) - 2 * u_old(I,J) + u(I,J+1)) / dy**2) + u_old(I,J)
            
                v(I,J) = - dt / 4* dy * (( v(I-1, J)+ v(I-1,J+1))**2 - (v_old(I, J)+v(I,J+1))**2) - dt / 4* dx * ((u_old(I,J)+ u(I,J+1)) &
                    * (v(I,J+1) + v(I-1, J+1)) - (u_old(I,J) + u(I,J-1)) * (v_old(I,J) + v(I-1,J))) - dt / dy *(p(I, J+1) - p(I,J)) &
                    + dt / Re * ((v(I+1,J) - 2 * v_old(I,J) + v(I-1,J)) / dx**2 + (v(I,J+1) - 2 * v_old(I,J) + v(I,J+1)) / dy**2) + v_old(I,J)
            
                p(I,J) = - Delta * dt / 2 * ((u(I,J+1)+ u_old(I,J)) - (u_old(I,J) + u(I,J-1))) - Delta * dt / 2 &
                    * ((v(I-1,J)+ v(I-1,J+1)) - (v_old(I,J) + v(I,J+1)))
            
            End Do
        
        End Do
        
    
    !-----------------------OUTPUTS GENERATION----------------------------- 
OPEN (1,FILE='FIELD.PLT') 
WRITE (1,*) 'VARIABLES=X,Y,u,v,p' 
WRITE (1,*) 'ZONE I=',M,' J=',N 
DO J=1,N 
  DO I=1,M 
    WRITE (1,*) X(I),Y(J),u(I,J),v(I,J),p(I,J)
  END DO 
END DO 
    
End Program Lid

flotus1 · November 14, 2020, 06:55

There are many ways to fix your initialization issue. Two of them would be
1) put the boundary conditions after the initialization
or
2) restrict initialization to data that does not belong to the boundary.

For code formatting, either mark your code and click on the code formatting button. Assuming you have at least the standard editor enabled.
Or type code brackets manually.
(CODE)code goes here(/CODE)
Replace round brackets () with square ones [] for full effect

M_H215 · November 14, 2020, 07:07

Thank you....I changed the Initialization section (Add a Do Loop) in order to prevent the over-written issue but still that problem exist...Do you think still over-written issue exist?

Code:

Program Lid
    
    implicit none
    Integer :: I,J,nx, ny, L, W, Iteration, Max_Iteration , Re, M, N, dt
    Real :: Delta, dx, dy
    Real, allocatable :: u(:,:), v(:,:), p(:,:), u_old(:,:), v_old(:,:), p_old(:,:), X(:), Y(:)
    !***************************************************!

    PRINT *, "ENTER THE DESIRED POINTS ..." 
    PRINT *, "... IN X DIRECTION:  SUGGESTED RANGE (20-200)" 
    READ*, M 
    PRINT *, "... IN Y DIRECTION:  SUGGESTED RANGE (10-100)" 
    READ*, N 
    
    ! Define Geometry
    dt = 0.001
    Delta = 2
    Re = 100
    L = 10
    W = 10
    dx = L /(M-1)
    dy = W /(N-1)
    
    ALLOCATE (X(M),Y(N),u(M,N),u_old(M,N),v(M,N),v_old(M,N),p(M,N),p_old(M,N))
    
    ! Grid Generation
    
    Do I = 1, M
        x(I) = (I-1)* dx
    End Do
    
    Do J=1 , N
        y(J) = (J-1) * dy
    End Do
     
    ! Boundray Condition
    Do I=1 , M
        u(I,1) = 0
        u(1,I) = 0
        u(N,I) = 0
        u(I,N) = 1       ! Lid Velocity
    End Do
    
    Do J=1, N
        v(J,1) = 0
        v(1,J) = 0
        v(J,M) = 0
        v(M,J) = 0
    End Do
    
    ! Initialization
    Do I= 2, M-1
        Do J=2, N-1
            u(I,J) = 0
            v(I,J) = 0
            p(I,J) = 0
            
        End do
    End do
                
    ! Solver
        Do I=2, M-1
            Do J=2, N-1
                u_old(I,J) = u(I,J) 
                v_old(I,J) = v(I,J)
                p_old(I,J) = p(I,J)
            
                u(I,J) = - dt / 4* dx * (( u(I, J+1)+ u_old(I,J))**2 - (u_old(I, J)+u(I,J-1))**2) - dt / 4* dy * ((u_old(I,J)+ u(I-1,J)) &
                    * (v(I-1,J) + v(I-1, J+1)) - (u_old(I,J) + u(I+1,J)) * (v_old(I,J) + v(I,J+1))) - dt / dx *(p(I, J+1) - p_old(I,J)) &
                    + dt / Re * ((u(I+1,J) - 2 * u_old(I,J) + u(I-1,J)) / dx**2 + (u(I,J+1) - 2 * u_old(I,J) + u(I,J+1)) / dy**2) + u_old(I,J)
            
                v(I,J) = - dt / 4* dy * (( v(I-1, J)+ v(I-1,J+1))**2 - (v_old(I, J)+v(I,J+1))**2) - dt / 4* dx * ((u_old(I,J)+ u(I,J+1)) &
                    * (v(I,J+1) + v(I-1, J+1)) - (u_old(I,J) + u(I,J-1)) * (v_old(I,J) + v(I-1,J))) - dt / dy *(p(I, J+1) - p(I,J)) &
                    + dt / Re * ((v(I+1,J) - 2 * v_old(I,J) + v(I-1,J)) / dx**2 + (v(I,J+1) - 2 * v_old(I,J) + v(I,J+1)) / dy**2) + v_old(I,J)
            
                p(I,J) = - Delta * dt / 2 * ((u(I,J+1)+ u_old(I,J)) - (u_old(I,J) + u(I,J-1))) - Delta * dt / 2 &
                    * ((v(I-1,J)+ v(I-1,J+1)) - (v_old(I,J) + v(I,J+1)))
            
            End Do
        
        End Do
        
    
    !-----------------------OUTPUTS GENERATION----------------------------- 
OPEN (1,FILE='FIELD.PLT') 
WRITE (1,*) 'VARIABLES=X,Y,u,v,p' 
WRITE (1,*) 'ZONE I=',M,' J=',N 
DO J=1,N 
  DO I=1,M 
    WRITE (1,*) X(I),Y(J),u(I,J),v(I,J),p(I,J)
  END DO 
END DO 
    
End Program Lid

flotus1 · November 14, 2020, 07:21

Boundary values getting over-written is probably no longer an issue here.
But you should go and double-check -and add as a comment to your code- which boundary (top, bottom, left, right) each of the 4 lines in your boundary condition definition corresponds to. The loop variables "I" and "J" appearing both as the first and second index is a bug if N is not equal to M.

M_H215 · November 14, 2020, 07:39

Thank you so much....I add some comments for my Boundary Conditions. Also, I changed the name of M and N in each loop of Boundary Condition in order to prevent some bugs that may arise...But I don't know what is the problem in my code...

!!

Code:

Program Lid
    
    implicit none
    Integer :: I,J,nx, ny, L, W, Iteration, Max_Iteration , Re, M, N, dt
    Real :: Delta, dx, dy
    Real, allocatable :: u(:,:), v(:,:), p(:,:), u_old(:,:), v_old(:,:), p_old(:,:), X(:), Y(:)
!****************************************************************!

    PRINT *, "ENTER THE DESIRED POINTS ..." 
    PRINT *, "... IN X DIRECTION:  SUGGESTED RANGE (20-200)" 
    READ*, M 
    PRINT *, "... IN Y DIRECTION:  SUGGESTED RANGE (10-100)" 
    READ*, N 
!****************************************************************!
    ! Define Geometry
    dt = 0.001
    Delta = 2
    Re = 100
    L = 10
    W = 10
    dx = L /(M-1)
    dy = W /(N-1)
    
    ALLOCATE (X(M),Y(N),u(M,N),u_old(M,N),v(M,N),v_old(M,N),p(M,N),p_old(M,N))
    
!****************************************************************!
    ! Grid Generation
    
    Do I = 1, M
        x(I) = (I-1)* dx     ! Grid Generation in X Direction
    End Do
    
    Do J=1 , N
        y(J) = (J-1) * dy    ! Grid Generation in Y Direction
    End Do
    
!****************************************************************!
    ! Boundray Condition
    Do I=1 , M
        u(I,1) = 0       ! Bottom Wall
        u(1,I) = 0       ! Left Wall
        u(M,I) = 0       ! Right Wall
        u(I,M) = 1       ! Lid Velocity (Top Wall)
    End Do
    
    Do J=1, N
        v(J,1) = 0       ! Bottom Wall
        v(1,J) = 0       ! Left Wall
        v(J,N) = 0       ! Top Wall 
        v(N,J) = 0       ! Right Wall
    End Do
!****************************************************************!
    
    ! Initialization
    Do I= 2, M-1
        Do J=2, N-1
            u(I,J) = 0
            v(I,J) = 0
            p(I,J) = 0
            
        End do
    End do
    
!****************************************************************!
    
    ! Solver
        Do I=2, M-1
            Do J=2, N-1
                u_old(I,J) = u(I,J) 
                v_old(I,J) = v(I,J)
                p_old(I,J) = p(I,J)
            
                u(I,J) = - dt / 4* dx * (( u(I, J+1)+ u_old(I,J))**2 - (u_old(I, J)+u(I,J-1))**2) - dt / 4* dy * ((u_old(I,J)+ u(I-1,J)) &
                    * (v(I-1,J) + v(I-1, J+1)) - (u_old(I,J) + u(I+1,J)) * (v_old(I,J) + v(I,J+1))) - dt / dx *(p(I, J+1) - p_old(I,J)) &
                    + dt / Re * ((u(I+1,J) - 2 * u_old(I,J) + u(I-1,J)) / dx**2 + (u(I,J+1) - 2 * u_old(I,J) + u(I,J+1)) / dy**2) + u_old(I,J)
            
                v(I,J) = - dt / 4* dy * (( v(I-1, J)+ v(I-1,J+1))**2 - (v_old(I, J)+v(I,J+1))**2) - dt / 4* dx * ((u_old(I,J)+ u(I,J+1)) &
                    * (v(I,J+1) + v(I-1, J+1)) - (u_old(I,J) + u(I,J-1)) * (v_old(I,J) + v(I-1,J))) - dt / dy *(p(I, J+1) - p(I,J)) &
                    + dt / Re * ((v(I+1,J) - 2 * v_old(I,J) + v(I-1,J)) / dx**2 + (v(I,J+1) - 2 * v_old(I,J) + v(I,J+1)) / dy**2) + v_old(I,J)
            
                p(I,J) = - Delta * dt / 2 * ((u(I,J+1)+ u_old(I,J)) - (u_old(I,J) + u(I,J-1))) - Delta * dt / 2 &
                    * ((v(I-1,J)+ v(I-1,J+1)) - (v_old(I,J) + v(I,J+1)))
            
            End Do
        
        End Do
        
    
    !-----------------------OUTPUTS GENERATION----------------------------- 
OPEN (1,FILE='FIELD.PLT') 
WRITE (1,*) 'VARIABLES=X,Y,u,v,p' 
WRITE (1,*) 'ZONE I=',M,' J=',N 
DO J=1,N 
  DO I=1,M 
    WRITE (1,*) X(I),Y(J),u(I,J),v(I,J),p(I,J)
  END DO 
END DO 
    
End Program Lid

flotus1 · November 14, 2020, 07:56

Assuming you still have the problem that all variables are zero...
Your computation loop contains another bug. You are only copying values to _old one-by-one, also leaving out the boundary nodes. Which leads to accessing uninitialized portions of the _old variables.

It should instead look something like this:

Code:

do t=1, n_iter
    u_old(:,:) = u(:,:)
    do i=2, m-1
        do j=2, n-1
            ! computation of u based on u_old
        end do
    end do
end do

There are probably more issues. I won't search and point out each of them. Instead, I would like to encourage you to take a step back, and go over your code a few more times.

M_H215 · November 14, 2020, 08:20

Thank you.... You mean that I should separate the loops... For example, first one for u...next one for v .. last one for p...?

since in the loop example that you sent there is just u...there isn't any v and p....

Also, is there any differences between

Code:

u(:,:) = u_old(:,:)

and

Code:

u(I,J) = u_old(I,J)

?? since you wrote the first model

Code:

u(:,:)

Thank you

flotus1 · November 14, 2020, 08:35

The reason why I did not recreate your whole code with my implementation: I'm lazy. And the short example should be enough to convey the general idea.

The difference between the single lines

Code:

u_old(:,:) = u(:,:)

and

Code:

u_old(i,j) = u(i,j)

should be fairly obvious: the former copies the whole array, the latter only copies values for indices i and j.
In the context of this code, the placement of these two lines also matters a lot. I put mine outside of the computation loop, in order to copy the values before they are accessed. Which is the right thing to do here.
Go through your original code, and try to trace back which values of u_old etc. are accessed, and whether they always hold the intended values at that point. The don't.

M_H215 · November 15, 2020, 00:46

Thank you... But even by changing them the problem still exist...

Code:

u(I,J), v(I,J), and p(I,J)

to

Code:

u(:,:), v(:,:), p(:,:)

Code:

Program Lid
    
    implicit none
    Integer :: I,J,nx, ny, L, W, Iteration, Max_Iteration , Re, M, N, dt
    Real :: Delta, dx, dy
    Real, allocatable :: u(:,:), v(:,:), p(:,:), u_old(:,:), v_old(:,:), p_old(:,:), X(:), Y(:)
!****************************************************************!

    PRINT *, "ENTER THE DESIRED POINTS ..." 
    PRINT *, "... IN X DIRECTION:  SUGGESTED RANGE (20-200)" 
    READ*, M 
    PRINT *, "... IN Y DIRECTION:  SUGGESTED RANGE (10-100)" 
    READ*, N 
!****************************************************************!
    ! Define Geometry
    dt = 0.001
    Delta = 2
    Re = 100
    L = 10
    W = 10
    dx = L /(M-1)
    dy = W /(N-1)
    Max_Iteration = 2000
    
    ALLOCATE (X(M),Y(N),u(M,N),u_old(M,N),v(M,N),v_old(M,N),p(M,N),p_old(M,N))
    
!****************************************************************!
    ! Grid Generation
    
    Do I = 1, M
        x(I) = (I-1)* dx     ! Grid Generation in X Direction
    End Do
    
    Do J=1 , N
        y(J) = (J-1) * dy    ! Grid Generation in Y Direction
    End Do
    
!****************************************************************!
    ! Boundray Condition
    Do I=1 , M
        u(I,1) = 0       ! Bottom Wall
        u(1,I) = 0       ! Left Wall
        u(M,I) = 0       ! Right Wall
        u(I,M) = 1       ! Lid Velocity (Top Wall)
    End Do
    
    Do J=1, N
        v(J,1) = 0       ! Bottom Wall
        v(1,J) = 0       ! Left Wall
        v(J,N) = 0       ! Top Wall 
        v(N,J) = 0       ! Right Wall
    End Do
!****************************************************************!
    
    ! Initialization
    Do I= 2, M-1
        Do J=2, N-1
            u(I,J) = 0
            v(I,J) = 0
            p(I,J) = 0
            
        End do
    End do
    
!****************************************************************!
    
    ! Solver
    Do Iteration = 1, Max_Iteration
        
        u_old(:,:) = u(:,:) 
        v_old(:,:) = v(:,:)
        p_old(:,:) = p(:,:)
        
        Do I=2, M-1
            Do J=2, N-1
            
                u(I,J) = - dt / 4* dx * (( u(I, J+1)+ u_old(I,J))**2 - (u_old(I, J)+u(I,J-1))**2) - dt / 4* dy * ((u_old(I,J)+ u(I-1,J)) &
                    * (v(I-1,J) + v(I-1, J+1)) - (u_old(I,J) + u(I+1,J)) * (v_old(I,J) + v(I,J+1))) - dt / dx *(p(I, J+1) - p_old(I,J)) &
                    + dt / Re * ((u(I+1,J) - 2 * u_old(I,J) + u(I-1,J)) / dx**2 + (u(I,J+1) - 2 * u_old(I,J) + u(I,J+1)) / dy**2) + u_old(I,J)
            
                v(I,J) = - dt / 4* dy * (( v(I-1, J)+ v(I-1,J+1))**2 - (v_old(I, J)+v(I,J+1))**2) - dt / 4* dx * ((u_old(I,J)+ u(I,J+1)) &
                    * (v(I,J+1) + v(I-1, J+1)) - (u_old(I,J) + u(I,J-1)) * (v_old(I,J) + v(I-1,J))) - dt / dy *(p(I, J+1) - p(I,J)) &
                    + dt / Re * ((v(I+1,J) - 2 * v_old(I,J) + v(I-1,J)) / dx**2 + (v(I,J+1) - 2 * v_old(I,J) + v(I,J+1)) / dy**2) + v_old(I,J)
            
                p(I,J) = - Delta * dt / 2 * ((u(I,J+1)+ u_old(I,J)) - (u_old(I,J) + u(I,J-1))) - Delta * dt / 2 &
                    * ((v(I-1,J)+ v(I-1,J+1)) - (v_old(I,J) + v(I,J+1)))
            
            End Do
        
        End Do
        
    End Do
    
        
    
    !-----------------------OUTPUTS GENERATION----------------------------- 
OPEN (1,FILE='FIELD.PLT') 
WRITE (1,*) 'VARIABLES=X,Y,u,v,p' 
WRITE (1,*) 'ZONE I=',M,' J=',N 
DO J=1,N 
  DO I=1,M 
    WRITE (1,*) X(I),Y(J),u(I,J),v(I,J),p(I,J)
  END DO 
END DO 
    
End Program Lid

flotus1 · November 15, 2020, 08:47

Like I said, there are a lot more questionable (i.e. wrong) pieces of code present. I would have to be paid to find them all at this point, and you wouldn't learn much.

These are the pointers I can give you:
1) use debug flags when compiling your code. Which ones depends on the compiler you are using
2) be aware of integer arithmetic. 1/4 is not the same as 1.0/4.0 in fortran. And don't use integers for variables that really should not be integers. examples in your code: dt, dx...

M_H215 · November 15, 2020, 11:40

Thank you...Actually I am using the debug flags but since the Fortran can not find any error, so debug flags as well (does not find any error...!)

Definitely there is something wrong with my code but Fortran says that there is not any error...!!

Code runs so fast (less than a second..!!) without any error but there is zero result.
I mean that I have the output but all of the data are zero...!!!

flotus1 · November 15, 2020, 13:09

Two things worth noting here:
1) I compiled your code with gfortran and all the useful debug flags I know. I got one serious compile-time warning about dt. It is declared as an integer, and you try to assign a value of 0.001 to it. Then during run-time, I got a segfault pointing to your implementation of boundary conditions in cases where M is not equal to N. Something I had already pointed out earlier.
2) The fact that a compiler does not find any errors is no indication that a code is free of bugs. Again, accidental use of integer arithmetic is a big one here. Look at how you compute values dx and dy. This is using integer arithmetic. If you add write(*,*) dx, dy right after you compute them, you will see that they are both 0. Leading to division by 0 later in your compute loop.
There is no other way: you have to go over your whole code, and eliminate accidental use of integer arithmetic. A compiler can't do that for you.

November 14, 2020, 05:03	Fortran ERROR	#1
M_H215 New Member Mohammad Join Date: Apr 2015 Posts: 17 Rep Power: 11	Hello I am trying to write a code for Lid-Driven-Cavity in FORTRAN (Very short Code). Actually, I am new with programming. The problem is, when I run the code Software apparently do the job so fast and the result created but there isn't any result...!! I mean that the PLT format file created but when I want to open it with tecplot, it gives me an error..!! when I check the result with Note-Pad , All of the data are Zero...!! Any help will be appreciated... the code is : .............................. Program Lid implicit none Integer :: I,J,nx, ny, L, W, Iteration, Max_Iteration , Re, M, N, dt Real :: Delta, dx, dy Real, allocatable :: u(:,:-), v(:,:-), p(:,:-), u_old(:,:-), v_old(:,:-), p_old(:,:-), X(:-), Y(:.) !*********************************************** ! PRINT , "ENTER THE DESIRED POINTS ..." PRINT , "... IN X DIRECTION: SUGGESTED RANGE (20-200)" READ, M PRINT , "... IN Y DIRECTION: SUGGESTED RANGE (10-100)" READ, N ! Define Geometry dt = 0.001 Delta = 2 Re = 100 L = 10 W = 10 dx = L /(M-1) dy = W /(N-1) ALLOCATE (X(M),Y(N),u(M,N),u_old(M,N),v(M,N),v_old(M,N),p(M ,N),p_old(M,N)) ! Grid Generation Do I = 1, M x(I) = (I-1) dx End Do Do J=1 , N y(J) = (J-1) * dy End Do ! Boundray Condition Do I=1 , M u(I,1) = 0 u(1,I) = 0 u(N,I) = 0 u(I,N) = 1 ! Lid Velocity End Do Do J=1, N v(J,1) = 0 v(1,J) = 0 v(J,M) = 0 v(M,J) = 0 End Do ! Initialization u = 0; v = 0 ; p = 0 ! Solver Do I=2, M-1 Do J=2, N-1 u_old(I,J) = u(I,J) v_old(I,J) = v(I,J) p_old(I,J) = p(I,J) u(I,J) = - dt / 4* dx * (( u(I, J+1)+ u_old(I,J))2 - (u_old(I, J)+u(I,J-1))2) - dt / 4* dy * ((u_old(I,J)+ u(I-1,J)) & * (v(I-1,J) + v(I-1, J+1)) - (u_old(I,J) + u(I+1,J)) * (v_old(I,J) + v(I,J+1))) - dt / dx (p(I, J+1) - p_old(I,J)) & + dt / Re ((u(I+1,J) - 2 * u_old(I,J) + u(I-1,J)) / dx*2 + (u(I,J+1) - 2 u_old(I,J) + u(I,J+1)) / dy*2) + u_old(I,J) v(I,J) = - dt / 4 dy * (( v(I-1, J)+ v(I-1,J+1))2 - (v_old(I, J)+v(I,J+1))2) - dt / 4* dx * ((u_old(I,J)+ u(I,J+1)) & * (v(I,J+1) + v(I-1, J+1)) - (u_old(I,J) + u(I,J-1)) * (v_old(I,J) + v(I-1,J))) - dt / dy (p(I, J+1) - p(I,J)) & + dt / Re ((v(I+1,J) - 2 * v_old(I,J) + v(I-1,J)) / dx*2 + (v(I,J+1) - 2 v_old(I,J) + v(I,J+1)) / dy*2) + v_old(I,J) p(I,J) = - Delta dt / 2 * ((u(I,J+1)+ u_old(I,J)) - (u_old(I,J) + u(I,J-1))) - Delta * dt / 2 & * ((v(I-1,J)+ v(I-1,J+1)) - (v_old(I,J) + v(I,J+1))) End Do End Do !-----------------------OUTPUTS GENERATION----------------------------- OPEN (1,FILE='FIELD.PLT') WRITE (1,) 'VARIABLES=X,Y,u,v,p' WRITE (1,) 'ZONE I=',M,' J=',N DO J=1,N DO I=1,M WRITE (1,*) X(I),Y(J),u(I,J),v(I,J),p(I,J) END DO END DO End Program Lid

November 14, 2020, 05:21		#2
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,427 Rep Power: 49	There are a few questionable pieces of code in here. But for the immediate problem at hand: the initialization over-writes the boundary condition. Next time you want to post code here, try the code format. This way, your code will look like this: Code: REAL, DIMENSION(:,:), ALLOCATABLE :: u do i=1, N u(i)=0.0 end do instead of this REAL, DIMENSION(:,, ALLOCATABLE :: u do i=1, N u(i)=0.0 end do M_H215 likes this.

November 14, 2020, 06:55		#5
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,427 Rep Power: 49	There are many ways to fix your initialization issue. Two of them would be 1) put the boundary conditions after the initialization or 2) restrict initialization to data that does not belong to the boundary. For code formatting, either mark your code and click on the code formatting button. Assuming you have at least the standard editor enabled. Or type code brackets manually. (CODE)code goes here(/CODE) Replace round brackets () with square ones [] for full effect M_H215 and aero_head like this.

November 14, 2020, 07:21		#7
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,427 Rep Power: 49	Boundary values getting over-written is probably no longer an issue here. But you should go and double-check -and add as a comment to your code- which boundary (top, bottom, left, right) each of the 4 lines in your boundary condition definition corresponds to. The loop variables "I" and "J" appearing both as the first and second index is a bug if N is not equal to M. M_H215 likes this.

November 14, 2020, 07:56		#9
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,427 Rep Power: 49	Assuming you still have the problem that all variables are zero... Your computation loop contains another bug. You are only copying values to _old one-by-one, also leaving out the boundary nodes. Which leads to accessing uninitialized portions of the _old variables. It should instead look something like this: Code: do t=1, n_iter u_old(:,:) = u(:,:) do i=2, m-1 do j=2, n-1 ! computation of u based on u_old end do end do end do There are probably more issues. I won't search and point out each of them. Instead, I would like to encourage you to take a step back, and go over your code a few more times. M_H215 and aero_head like this.

November 14, 2020, 06:29		#3
M_H215 New Member Mohammad Join Date: Apr 2015 Posts: 17 Rep Power: 11	Thank you so much.... Actually I tried to send my code like something that you mentioned but I could not find a way to do that... How I can send my code like something that you Mentioned ? (Like Code) Also, Could you please tell me that how I can prevent this issue ? (initialization over-written my boundaries) Thank you

November 14, 2020, 08:20		#10
M_H215 New Member Mohammad Join Date: Apr 2015 Posts: 17 Rep Power: 11	Thank you.... You mean that I should separate the loops... For example, first one for u...next one for v .. last one for p...? since in the loop example that you sent there is just u...there isn't any v and p.... Also, is there any differences between Code: u(:,:) = u_old(:,:) and Code: u(I,J) = u_old(I,J) ?? since you wrote the first model Code: u(:,:) Thank you

November 14, 2020, 08:35		#11
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,427 Rep Power: 49	The reason why I did not recreate your whole code with my implementation: I'm lazy. And the short example should be enough to convey the general idea. The difference between the single lines Code: u_old(:,:) = u(:,:) and Code: u_old(i,j) = u(i,j) should be fairly obvious: the former copies the whole array, the latter only copies values for indices i and j. In the context of this code, the placement of these two lines also matters a lot. I put mine outside of the computation loop, in order to copy the values before they are accessed. Which is the right thing to do here. Go through your original code, and try to trace back which values of u_old etc. are accessed, and whether they always hold the intended values at that point. The don't. M_H215 likes this.

November 15, 2020, 08:47		#13
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,427 Rep Power: 49	Like I said, there are a lot more questionable (i.e. wrong) pieces of code present. I would have to be paid to find them all at this point, and you wouldn't learn much. These are the pointers I can give you: 1) use debug flags when compiling your code. Which ones depends on the compiler you are using 2) be aware of integer arithmetic. 1/4 is not the same as 1.0/4.0 in fortran. And don't use integers for variables that really should not be integers. examples in your code: dt, dx... M_H215 likes this.

November 15, 2020, 11:40		#14
M_H215 New Member Mohammad Join Date: Apr 2015 Posts: 17 Rep Power: 11	Thank you...Actually I am using the debug flags but since the Fortran can not find any error, so debug flags as well (does not find any error...!) Definitely there is something wrong with my code but Fortran says that there is not any error...!! Code runs so fast (less than a second..!!) without any error but there is zero result. I mean that I have the output but all of the data are zero...!!!

November 15, 2020, 13:09		#15
flotus1 Super Moderator Alex Join Date: Jun 2012 Location: Germany Posts: 3,427 Rep Power: 49	Two things worth noting here: 1) I compiled your code with gfortran and all the useful debug flags I know. I got one serious compile-time warning about dt. It is declared as an integer, and you try to assign a value of 0.001 to it. Then during run-time, I got a segfault pointing to your implementation of boundary conditions in cases where M is not equal to N. Something I had already pointed out earlier. 2) The fact that a compiler does not find any errors is no indication that a code is free of bugs. Again, accidental use of integer arithmetic is a big one here. Look at how you compute values dx and dy. This is using integer arithmetic. If you add write(,) dx, dy right after you compute them, you will see that they are both 0. Leading to division by 0 later in your compute loop. There is no other way: you have to go over your whole code, and eliminate accidental use of integer arithmetic. A compiler can't do that for you. M_H215 likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
DPM udf error	haghshenasfard	FLUENT	0	April 13, 2016 07:35
[OpenFOAM] Native ParaView Reader Bugs	tj22	ParaView	270	January 4, 2016 12:39
[OpenFOAM.org] Compile OF 2.3 on Mac OS X .... the patch	gschaider	OpenFOAM Installation	225	August 25, 2015 20:43
CGNS lib and Fortran compiler	manaliac	Main CFD Forum	2	November 29, 2010 07:25
How to get the max value of the whole field	waynezw0618	OpenFOAM Running, Solving & CFD	4	June 17, 2008 06:07