Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
28794 Discussions

What is the solution for the analysis result being NaN? (Fortran)

alize
Beginner
971 Views

        Hi,

  • I try to solve hydrodynamic flow equation in Fortran. The codes include two steps. In the first step, hydraulic parameters solved by direct step method. In second step, I try to use this hydraulic parameters to find flow depth and velocity in unsteady flow analysis.
  • The code works well without any error. For the first part gives hydraulic table as it should be(first image) Hydraulic analysis Table
  • But for the second output for unsteady flow analysis it gives flow depth and flow velocity as NaN format.(second image) unsteady flow analysis
  • When I check the code I don't see any reason why the code result as NaN should be given, such as division by 0 or a result that is not an integer.
  • Please support me to solve this issue. I am sharing my code as part 1 and part 2 . And sharing my code as in the file.
  •  

  • WRITE(3,111)LTIME
    111 FORMAT (//,2X,'TIME=',I6,/)
    WRITE(3,110) Q
    110 FORMAT (2X,'DISCHARGE=',F10.4,/)
    WRITE(3,105)
    105 FORMAT (2X,'LENGTH',3X,'DEPTH',3X,'VELOCTIY',/)
    WRITE(3,106)(DX(K),Y(K),U(K),K=1,I)
    106 FORMAT (2X,F12.4,3X,F10.4,3X,F10.4)
    ENDIF
    ! REPLACING
    DO 50 K=1,I
    YO(K)=Y(K)
    UO(K)=U(K)
    50 CONTINUE
    20 CONTINUE
    RETURN
    END
  • ! Part 2 END ----------------------------
  •  
  •  
0 Kudos
1 Solution
andrew_4619
Honored Contributor III
528 Views

OK I copy/pasted  that source into visual studio and set the settings to not check things I would formally check (standards etc) to allow it to compile.

I tried running with the input you suggest but I think that is not correct:

(input: 1,1, 1,3.6,7,15,7,0,0.0015, 0.02,100,1,100,5, 20 and 7200, 5 and 10000, 12000, 100) 

The 11th item asked is IMETHOD and the data above suggests a value 100  but actually the only valid imethod value is 1 so I think the input is listed incorrectly or I am reading it wrong. Either way you ask too much for a help forum I should not have to guess or to try to understand the code to be able to help you, this is a Fortran forum not a hydrodynamics forum.

 

You should  run the program and copy/past the cmd window showing question and answer that would be easier to follow. I am sure that running this program in the debugger will show the problem very easily. 

 

 

 

View solution in original post

0 Kudos
16 Replies
JohnNichols
Valued Contributor III
905 Views
    dimension YO(400), TERSY(400)
    REAL N,L
    character*14 FILEOU
    AREA(B,YY,Z)=(B+Z*YY)*YY
    PERI(B,YY,Z)=B+2*YY*(1+Z**2)**.5
    WRITE(*,8)
8   format(/,19X,'PLEASE ENTER THE OUTPUT FILE NAME',//)
    READ(*,8008)FILEOU
8008 FORMAT(A14)
    OPEN(3,FILE=FILEOU)
    READ(*,*)QWE
    WRITE(*,*)'DOWNSTREAM BOUNDARY CONDITION:'
    WRITE(*,*)'SPILLWAY (1), BOTTOM OUTLET (2)'
    READ(*,*)ICHECK
    IF(ICHECK.EQ.1)THEN
        WRITE(*,*)'ENTER SPILLWAY HEIGHT, P(m)'
        read(*,*)P
        WRITE(*,*)'ENTER SPILLWAY CREST LENGTH, L(m)'
        READ(*,*)L
    ENDIF
    IF(ICHECK.EQ.2) THEN
        WRITE(*,*)'ENTER DISCHARGE COEFFICIENT'
        READ(*,*)C
        WRITE(*,*)'ENTER DIAMETER OF BOTTOM OUTLET(m)'
        READ(*,*)D
        WRITE(*,*)'ENTER THE NUMBER OF BOTTOM OUTLETS'
        READ(*,*)NNN
    ENDIF
    WRITE(*,*)'ENTER STEADY - STATE DISCHARGE (m^3/s)'
    READ(*,*)QQ
    WRITE(*,*)'ENTER BOTTOM WIDTH(m)'
    READ(*,*)BB
    WRITE(*,*)'ENTER SIDE INCLINATION'
    READ(*,*)ZZ
    WRITE(*,*)'ENTER BED SLOPE'
    READ(*,*)S0
    WRITE(*,*)'ENTER MANNING ROUGHNESS'
    READ(*,*)N
    WRITE(*,*)'ENTER Dx(m)'
    READ(*,*)DXX

    WRITE(*,*)'CHOOSE THE SOLUTION TECHNIQUE:'
    WRITE(*,*)'COMPLETE DYNAMIC WAVE(1)'

    READ(*,*)IMETHOD
    WRITE(*,*)'ENTER TIME INTERVAL Dt, (sec)'
    READ(*,*)DT
    WRITE(*,*)'ENTER THE HYDROGRAPH PARAMETERS :'
    WRITE(*,*)'Qo(m^3/S)'
    READ(*,*)QQ0
    WRITE(*,*)'Qp (m^3/s) AND Tp(sec)'
    READ(*,*)QQP
    READ(*,*)Tp
    WRITE(*,*)'Qf(m^3/s) AND Tn (sec)'
    READ(*,*)QQF
    READ(*,*)Tn
    WRITE(*,*)'Tmax(sec)'
    READ(*,*) Tmax


    YNN=0.01
99  ETOTAL=(QQ*N)/S0**.5
    ECHECK=AREA(BB,YNN,ZZ)**(5./3.)/PERI(BB,YNN,ZZ)**(2./3.)
    IF(ETOTAL.LT.ECHECK) GOTO 100
    YNN=YNN+.01
    GOTO 99
100 W=QQ/(3.132*BB**2.5)
    IF(ZZ.EQ.0) THEN
        YCOB=W**.672717*1.022785
    ELSE IF(ZZ.EQ..5) THEN
        YCOB=W**.602684*0.748173
    ELSE IF(ZZ.EQ.1) THEN
        YCOB=W**.574105*0.657397
    ELSE IF(ZZ.EQ.1.5) THEN
        YCOB=W**.564345*0.625201
    ELSE IF(ZZ.EQ.2) THEN
        YCOB=W**.552851*0.584043
    ELSE IF(ZZ.EQ.2.5) THEN
        YCOB=W**.545379*0.552333
    ELSE IF(ZZ.EQ.3) THEN
        YCOB=W**.537298*0.523298
    ELSE IF(ZZ.EQ.4) THEN
        YCOB=W**.524667*0.487638
    END IF
    YC=YCOB*BB
    IF(YNN.GT.YC)THEN
        WRITE(3,551)
551     FORMAT('HYDRAULIC PARAMETERS RESULT IN MILD SLOPE'   ,/)
    ELSE
        WRITE(3,552)
552     FORMAT('THIS PROGRAM IS NOT VALID FOR STEEP SLOPE CASE'   ,/)

    ENDIF
    IF(ICHECK.EQ.1) THEN
        H0=0.1
98      X=P/H0
        IF(X.GE.2.8) THEN
            CO=2.175
        ELSE
            C0=.018*X**7-.218*X**6+1.084*X**5-2.872*X**4+4.37*X**3-3.847*X**&
                2+1.889*X+1.72
        ENDIF
        QCHECK=C0*L*H0**1.5
        IF(QQ.LT.QCHECK) GOTO 101
        H0=H0+.05
        GOTO 98
101     YMAX = P+H0
    ELSE
        YMAX=0.0826*QQ**2/(NNN**2*C**2*D**4)
    ENDIF
    KZ=ZZ
    WRITE(3,15)QQ,YMAX,YNN,BB,KZ,S0,N,DXX,YC
15  FORMAT(10X,'Q(m3/s =',F6.1,/,10x,'Ymax(m)=',F6.3,/,10X,'Yn(&
        m)=',F6.3,/,/,10X,'B(m)=',f6.1,/,10x,'Z=1V:',I1,'H',/,10X,'So&
        =',F6.5,/,10X,'n=',F6.3,/,10X,'DX(m)=',F6.3,/,10X,'Yc(m)=&
        ',F6.3,///)
    WRITE(3,50)
50  FORMAT(1X,'I',5X,'Z',6X,'Y',5X,'A',5X,'P',5X,'R',6X,'U',7X,'Sf',6X&
        ,'H',7X,'HE',5X,'LENGTH',/,6X,'(m)'4X,'(m)',2X'(m2)',3X,'(m)',3X&
        ,'(m)',3X,'(m/s)',11X,'(m)',5X,'(m)',/,78('--'),/)
    T=0.
    Z1=0.
    Y1=YMAX
    I=1
    A1=AREA(BB,Y1,ZZ)
    P1=PERI(BB,Y1,ZZ)
    R1=A1/P1
    U1=QQ/A1
    S1=(U1*N/R1**(2./3.))**2
    H1=Z1+Y1+U1**2/19.62
    WRITE(3,51)I,Z1,Y1,A1,P1,R1,U1,S1,H1
51  FORMAT(1X,I3,1X,F6.3,1X,F6.3,1X,2(F7.1,1X),2(F6.3,1X),F7.6,1X,F6.3&
        ,1X,F6.2,3X,F7.1)
30  I=I+1
    T=T+DXX
    Z2=Z1+S0*DXX
    Y2=0.
10  Y2=Y2+0.01
    A2=AREA(BB,Y2,ZZ)
    P2=PERI(BB,Y2,ZZ)
    R2=A2/P2
    U2=QQ/A2
    S2=(U2*N/R2**(2./3.))**2
    H2=Z2+Y2+U2**2/19.62
    H22=H1+.5*DXX*(S1+S2)
    HE=H2-H22
    HHE=HE*100
    IF(ABS(HE).LT.0.0074) GOTO 20
    GOTO 10
20  FRU2 = U2/(9.81*Y2)**.5
    DELY=HE/(1-FRU2**2+(3*S2*DXX)/(2*R2))
    Y2=Y2-DELY
    WRITE(3,51)I,Z2,Y2,A2,P2,R2,U2,S2,H2,HHE,T
    IF(Y2.LE.YNN) GOTO 40
    WRITE(*,*)I,Y2
    TERSY(I)=Y2
    Y1=Y2
    A1=A2
    P1=P2
    R1=R2
    U1=U2
    S1=S2
    Z1=Z2
    H1=H2
    S1=S2
    GOTO 30
40  WRITE(3,75) T
75  FORMAT(//,10X,'TOTAL LENGTH IS',F7.1,'m',//)
    YO(1) = YNN
    YO(I) = YMAX
    WRITE(*,*)T,I
    DO 597 K =2,I-1
        YO(K)=TERSY(I-K+1)
597 CONTINUE
    WRITE(*,9)
9   FORMAT(//,15('!'),3X,'PLEASE WAIT COMPUTATION CONTUNIES',3X,15('!'&
        ))


    IF(IMETHOD.EQ.1)CALL COMPW(QQ,BB,S0,N,DXX,YNN,I,ZZ,YO,ICHECK,D,&
        C,NNN,QQ0,QQP,QQF,TP,TN,TMAX,DT,P,L)

553 STOP
    END
    SUBROUTINE&
        COMPW(QQ,BB,S0,N,DXX,YNN,I,ZZ,YO,ICHECK,D,C,NNN,QQ0,QQP&
        ,QQF,TP,TN,TMAX,DT,P,L)
    REAL N,L
    DIMENSION DX(400),UO(400),YO(400),Y(400),U(400)
    AREA(B,YY,Z)=(B+Z*YY)*YY
    PERI(B,YY,Z)=B+2*YY*(1+Z**2)**.5
    WRITE(3,*)'RESULTS OF COMPLETE DYNAMIC WAVE'
    QIN=QQ
    WID=BB
    S=S0
    CON=N
    DL=DXX
    DX(1)=0
    DO 15 K=1,I
        UO(K)=QIN/AREA(WID,YO(K),ZZ)
        DX(K+1)=DX(K)+DL

15  CONTINUE
    WRITE(3,199)
199 FORMAT (2X,'TIME=0',/)
    WRITE(3,101) QIN
101 FORMAT (2X, 'INITIAL DISCHARGE=',F9.4,/)
    WRITE(3,102)
102 FORMAT (2X, 'LENGTH',3X,'DEPTH ',3X,'VELOCITY',/)
    WRITE(3,103) (DX(K),YO(K),UO(K),K=1,I)
103 FORMAT (2X,F12.4,3X,F10.4,3X,F10.4)
    WRITE(*,*)'ENTER THE PRINT OUT TIME INTERVAL, (sec)'
    READ(*,*)J
    DO 20 LTIME=1,int(TMAX)
        WRITE(*,*)LTIME

        ! EQUATIONS OF FLOOD HYDROGRAPH
        DIS=QQ0+(QQP-QQ0)/TP*LTIME
        IF(LTIME.GT.TP) DIS=QQP-(QQP-QQF)/(TN-TP)*(LTIME-TP)
        IF(LTIME.GT.TN) DIS=QQF

        ! DETERMINATION OF VELOCTIES AND WATER DEPTHS AT TIME LEVEL NNN
        Q=DIS
        Y(1)=0.01
121     ETOTAL=(Q*CON)/S**.5
        ECHECK=(AREA(WID,Y(1),ZZ))**(5./3.)/(PERI(WID,Y(1),ZZ))**(2./3.)
        IF(ETOTAL.LT.ECHECK) GOTO 122
        Y(1) = Y(1)+.01
        GOTO 121
122     U(1)=Q/AREA(WID,Y(1),ZZ)
        DX(1)=0
        DO 30 K=2,I-1

            DX(K+1)=DX(K)+DL
            Y1=UO(K)*(YO(K-1)-YO(K+1))
            Y2=YO(K)*(UO(K-1)-UO(K+1))
            Y(K)=YO(K)+DT/2./DL*(Y1+Y2)
            GAM=Y(K)**(4./3)/9.81/DT/CON**2
            BET1=UO(K+1)-UO(K-1)
            BET2=9.81*DT/2./DL*(YO(K+1)-YO(K-1))
            BET=GAM*(UO(K)*DT/2./DL*BET1+BET2-9.81*DT*S-UO(K))
            U(K)=(-GAM+(GAM**2.-4*BET)**0.5)/2.

30      CONTINUE
        DO 40 K=2,I
            Y1=YO(K)*(UO(K-1)-UO(K))
            Y2=UO(K)*(YO(K-1)-YO(K))
            Y(K)=YO(K)+DT/DL*(Y1+Y2)




40      CONTINUE
        IF(J.EQ.LTIME)THEN
            J=J+1800
            WRITE(3,111)LTIME
111         FORMAT (//,2X,'TIME=',I6,/)
            WRITE(3,110) Q
110         FORMAT (2X,'DISCHARGE=',F10.4,/)
            WRITE(3,105)
105         FORMAT (2X,'LENGTH',3X,'DEPTH',3X,'VELOCTIY',/)
            WRITE(3,106)(DX(K),Y(K),U(K),K=1,I)
106         FORMAT (2X,F12.4,3X,F10.4,3X,F10.4)
        ENDIF
        ! REPLACING
        DO 50 K=1,I
            YO(K)=Y(K)
            UO(K)=U(K)
50      CONTINUE
20  CONTINUE
    RETURN
    END

There is a button to insert text so it is properly set out and colour coded. 

What is QWE?  You ask for it as the second input but there is no indication of what it should be? People can not read minds, it is important to remember this when writing code. 

You should start with a program unit and you should declare all your variables, it is easy to make a mistake.  

Provide a sample input and output so we can replicate same. 

 

0 Kudos
JohnNichols
Valued Contributor III
904 Views
AREA(B,YY,Z)=(B+Z*YY)*YY
PERI(B,YY,Z)=B+2*YY*(1+Z**2)**.5

 

The complier does not throw an error at this line,  but area is undefined in the code and the watch window will tell you that. 

 

include 'mkl_vsl.f90'
include "errcheck.inc"

!     ******************************************************************
!
PROGRAM ULARC
!
!     ******************************************************************
!
!     SAMPLE ELASTO-PLASTIC ANALYSIS OF PLANE FRAMES
!
!     PROGRAMMED BY A. SUDHAKAR, G. H. POWELL, G. ORR, R. WHEATON
!		    UNIVERSITY OF CALIFORNIA, BERKELEY, CALIFORNIA
!		    1972
!
!     ******************************************************************
!

use Base
use Scotia
use ACADDraw
Implicit None
character*60 ZFILE1,ZFILE2

INTEGER IBUST           !   Failure Code
INTEGER ICLOSE          !   Limi
INTEGER NCYC
INTEGER INEL
INTEGER JNEL
INTEGER NLC
INTEGER NELS
REAL XFAC
INTEGER NEQ
INTEGER NJTS
INTEGER NSJTS,NTERM
REAL DMAX
REAL FAX
INTEGER NLCS,MBAND, XX


REAL X(nodes),Y(nodes),RRTOT(nodes,3),XA(ij,mn_10),YA(ij,mn_10),LIMIT(nstates), mass(nstates)

REAL YMP(nstates),YMN(nstates),ROTCP(nstates),ROTCN(nstates)

If you are going to ask humans for help, at least make it to a standard that a masters student in 1972 got to, see sample for how to start a program in Fortran from 1972. Implicit none is critical.  

 

0 Kudos
alize
Beginner
882 Views

My friend I didn't write this code just I tried to convert to Matlab so  Please keep your advice to yourself and  be logical

0 Kudos
JohnNichols
Valued Contributor III
880 Views
  • I try to solve hydrodynamic flow equation in Fortran. The codes include two steps. In the first step, hydraulic parameters solved by direct step method. In second step, I try to use this hydraulic parameters to find flow depth and velocity in unsteady flow analysis.
  • The code works well without any error. For the first part gives hydraulic table as it should be(first image) Hydraulic analysis Table
  • But for the second output for unsteady flow analysis it gives flow depth and flow velocity as NaN format.(second image) unsteady flow analysis
  • When I check the code I don't see any reason why the code result as NaN should be given, such as division by 0 or a result that is not an integer.
  • Please support me to solve this issue. I am sharing my code as part 1 and part 2 . And sharing my code as in the file.

 

Where do you say you are looking at MATLAB?

I try to solve hydrodynamic flow equation in Fortran == in English this says you are doing it in Fortran 

I took the time to download the code, compile it and then tell you YOUR CODE WILL NOT RUN and if you ask for help then provide the information I need to run the program.  

You are using a board that has a facility to show the code in a readable format - learn to use it if you want EXPERT help.  

What is QWE as the input? A simple question and your code will not execute as shown with correct answers.  

Do you expect people to look at the code and tell you the error?  If so you are not really understanding how to code. 

 

You asked for advice, I am now finished with your problem - have a nice day. 

 

 

0 Kudos
alize
Beginner
878 Views
0 Kudos
JohnNichols
Valued Contributor III
866 Views

There is a button to insert text so it is properly set out and colour coded.  This is courtesy to your reader.  The fact that you did not use it means you are not commonly on the board or you do not care about someone reading it in legible fashion as I did for you and the others who might help. 

What is QWE?  You ask for it as the second input but there is no indication of what it should be? People can not read minds, it is important to remember this when writing code. 

You should start with a program unit and you should declare all your variables, it is easy to make a mistake.  

Provide a sample input and output so we can replicate same. 

--------------------------------------------------------------------------------------------------------------------------

Answer these questions and you might get some help - otherwise you are wasting your time as you cannot run the program if we do not know what the answers to the questions are?  It is called a sample input and output.

This cannot be the complete code as I noted variables are not defined.  

---------------------------------------------------------------------------------------------------------------------------

This board has very few readers on the weekend.  

----------------------------------------------------------------------------------------------------------------------------

If you want to argue go ahead.  This is my last post until I get QWE and a sample or if you wish to continue to be an interesting person. 

0 Kudos
andrew_4619
Honored Contributor III
509 Views

John,

AREA and PERI are obsolete "statement functions" and between  declarations and the first executable statement is the correct place to have them. Old code, old constructs....

0 Kudos
alize
Beginner
890 Views

......

0 Kudos
alize
Beginner
856 Views

.......

0 Kudos
JohnNichols
Valued Contributor III
837 Views

Screenshot 2024-10-26 123544.png

If I run your code it stops after I give it an output file name and it WAITS forever. 

If I look at the code it is stopped on a read 

Screenshot 2024-10-26 123608.png

It is stopped at the read(*,*)QWE

I asked in the first post what value do I put in for QWE so I can run our code and look for errors.  

I have read all your posts since you joined a fortnight ago, you are not getting answers to your questions, because they are not questions people can answer without the code.   The one answer was declare variables. 

Good luck on writing the MATLAB code.  

If you are going to proceed asking questions, provide the whole code in the standard presentation window and provide a sample input.  

Otherwise you are wasting your time.  Me I just enjoy watching someone fumble. 

 

0 Kudos
alize
Beginner
823 Views

My friend, you can put any number you want there. Afterwards, the code works without any problems and does not give any errors, as I already stated in the first message. The only problem with the table is that some numbers appear as "NaN".
(input: 1,1, 1,3.6,7,15,7,0,0.0015, 0.02,100,1,100,5, 20 and 7200, 5 and 10000, 12000, 100)

I uploaded also output as a zip file.

0 Kudos
JohnNichols
Valued Contributor III
820 Views

My dear friend:

1. Thank you for the input data.

2. It is not readily observed in the first post. 

I will try it later today. 

JMN

0 Kudos
andrew_4619
Honored Contributor III
529 Views

OK I copy/pasted  that source into visual studio and set the settings to not check things I would formally check (standards etc) to allow it to compile.

I tried running with the input you suggest but I think that is not correct:

(input: 1,1, 1,3.6,7,15,7,0,0.0015, 0.02,100,1,100,5, 20 and 7200, 5 and 10000, 12000, 100) 

The 11th item asked is IMETHOD and the data above suggests a value 100  but actually the only valid imethod value is 1 so I think the input is listed incorrectly or I am reading it wrong. Either way you ask too much for a help forum I should not have to guess or to try to understand the code to be able to help you, this is a Fortran forum not a hydrodynamics forum.

 

You should  run the program and copy/past the cmd window showing question and answer that would be easier to follow. I am sure that running this program in the debugger will show the problem very easily. 

 

 

 

0 Kudos
JohnNichols
Valued Contributor III
454 Views

Andrew:

I would tend to agree with you, but I actually did this type of analysis for a Fortran program in 1988 and I find the problem interesting.  But the code needs a bit of fixing to get to the problems identified easily with VS.  

But like all things, ULARC calls mightily in the night to solve real bridge problems, so the person will have to be patient, unless a real expert like Jim can look at the code and say ah....   I have known several of those people - interesting 4 or 5 depending on the way the North wind is blowing are on this forum.  

In 1948, in Boston, the ASCE published some sewer flow charts that were still used when I was designing sewers in the middle of a recession.  I coded the Manning's equation in Fortran and it gave the exact answers, once I guessed the roughness for the very low open flows. The moral of the story, no matter how long it takes to code in Fortran it is always faster than by hand and not matter how long it takes to code in MATLAB it is always faster by hand.  

Have a nice day.  

JMN

0 Kudos
alize
Beginner
521 Views

Okey guys, thank you for taking time. Hava e a nice day.

0 Kudos
Reply