Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

FORTRAN: I have no idea what is happening

  1. Jun 18, 2011 #1
    I am freaking out man :smile: I have this program and in it I have the variable nPT declared as Integer. Somehow, after nPT is passed to subroutine, its value changes to a very large negative integer value. I cannot make sense of it. In the subroutine pasted at the bottom, the value of N (which is passed the value nPT in the subroutine) has the value of 452 immediately before the Do loop and a value of -1992004147 immediately after the Do loop.

    I have not done anything to the value of N in the Do loop within this subroutine, so I do not know how it is changing.

    Here is the pertinent part of the Main Program

    Code (Text):


        nPT = j        [B]! At this point nPT = 452[/B]
       
    c   ********************Calculate Thermodynamic Properties****************
            [B]! At this point nPT = 452 still[/B]

    c   Unburned Gas Temperature   
        Call T_unburned(nPT,T_0,P_0,P,a_mix,R_mix,T_u)
        [B]! At this point nPT = 452 still[/B]   

    c   Unburned Gas Specific Volume
        Call v_unburned(nPT,T_u,P,R_mix,v_u)
    [COLOR="Red"]       [B]! At this point nPT =  -1992004147[/B][/COLOR]
     

    The Subroutine v_unburned is as follows:

    Code (Text):

        Subroutine v_unburned(N,T,P,R,v)
    c////////////////////////////////////////////////////////////////////
    c   This subroutine calculates the specific volume (m**3/kg) of
    c   the nuburned gas using the ideal gas equation of state:
    c   pv = RT
    c
    c   INPUT:
    c
    c       P           = pressure array
    c       T           = temperature array
    c       R           = gas constant of mixture
    c       N           = number of elements
    c      
    c   OUTPUT:
    c
    c       v           = specific volume (m**3 / kg) array
    c
    c////////////////////////////////////////////////////////////////////
        Real*8  T(1000),P(1000),v(1000),R
        Integer i, N
       
        [B]! At this point N = 452 [/B] N is nPT in the subroutine   
        Do i = 1,N
            v(i)=R * T(i) / (P(i) * 101.325 / 14.7)
        End Do
           [COLOR="Red"]       [B]! At this point nPT =  -1992004147[/B][/COLOR]

        Return
       
        End
     
    Any thoughts :confused:
     
  2. jcsd
  3. Jun 18, 2011 #2

    Hurkyl

    User Avatar
    Staff Emeritus
    Science Advisor
    Gold Member

    Typically, such problems come from buffer overflows. Are the arrays T_u, P, and v_u actually 452 elements long?

    I'm not too familiar with FORTRAN syntax; don't you have to do something to indicate what subroutine arguments are inputs and which are outputs?
     
  4. Jun 18, 2011 #3

    gb7nash

    User Avatar
    Homework Helper

    Just out of curiousity, if you replace N in the do loop with some arbitrary number, is the value of N still 452 when you exit? Also, do these other values in the do loop rely on N in any way?
     
  5. Jun 18, 2011 #4

    Mark44

    Staff: Mentor

    Can you show us a complete listing of your code? I don't see anything in your code for v_unburned that would cause it to change the value of N from 452 to -1992004147. In fact, I don't see anything that changes N at all in this routine.

    The only thing that comes to mind is a mismatch in the types of the actual parameters and the formal parameters in the v_unburned routine. It would be helpful to see the declarations of nPT,T_u,P,R_mix, and v_u (the actual parameters) and compare them with the declarations of N,T,P,R, and v (the formal parameters), that are used in the routine.
     
  6. Jun 18, 2011 #5
    When I run in too a similar bug in my code for the first time, it took me hours to find out what the problem was. As hurkyl pointed out, such errors are caused when you assign values to an array beyond the allocated memory.
    Now, I'm not familiar with fortran, but it doesn't seem the case that you're exceeding the array length, unless you have declared the arrays outside the subroutine with a different size (1000 in this case) which would mean that memory is reallocated and written over the existing and substituted with that strange negative number.
     
  7. Jun 18, 2011 #6
    Here is the Main Program:

    Code (Text):
        Program Two Zone

    c////////////////////////////////////////////////////////////////////
    c   This program is used to calculate the ....<snip>
    c
    c       T_0       = initial Temp (K)       
    c       P_0       = initial pressure (psi) 
    c       PHI       = equivalence ratio      
    c       FAM       = total # moles per mole of fuel
    c       X(i)      = mole fraction of i-th species
    c       SA        = stoichiometric air
    c       a(i,j)    = i x j(=7) matrix containing reactant
    c                   NASA coefficients
    c       a_mix(i)  = mole-weighted NASA coefficients of mixture
    c       MM_mix    = mole-weighted molar mass of mixture (kg/kmol)
    c       R_mix     = gas constant of mixture (kJ/kg-K)
    c       temp_t    = temporary time array (for sorting)
    c       temp_P    = temporary pressure array (for sorting)
    c       P(i)      = pressure array (psi)
    c       t(i)      = time array (msec)
    c       nPT       = number of points
    c       T_u(i)    = unburned gas temperature array (K)
    c       v_u(i)    = unburned gas specific vol array (m**3/kg)
    c       u_u(i)    = unburned gas specific energy array (kJ/kg)
    c   ********************Initialize Problem************************
        Real*8  T_0, P_0, PHI, FAM, X_air, X(10), a(10,10), a_mix(10),
         &      MM_mix, R_mix, temp_P(5000), temp_t(5000),P(2000),
         &      t(2000),T_u(1000),TIME,CH7
         
        Integer i, j, nPT

        Character   Molecule*10


        T_0 = 298.0
        P_0 = 14.58131852
        PHI = 1.0
        SA  = 11.9

    c    Summarize Problem Statement

        Write(*,*) " --------------------------------- "
        Write(*,*) "   Initial Thermodynamic Condition "
        Write(*,*) " --------------------------------- "
        Write(*,*) "  "
        Write(*,*) "   Initial Temperature  = ", T_0," (K)"
        Write(*,*) "  "
        Write(*,*) "   Initial Pressure     =", P_0," (psi)"
        Write(*,*) "  "
        Write(*,*) "   Equivalence Ratio    =", PHI
        Write(*,*) "  "
        Write(*,*) " --------------------------------- "

        Write(*,*) "   The Fuel Molecule is C2H4F2"

    c   Calculate mole fractions
        FAM   = 1 + SA / PHI
        X_air = (FAM - 1) / FAM     !mole fraction air
       
        X(1)  = 0.21 * X_air            !mole fraction O2
        X(2)  = 0.79 * X_air            !mole fraction N2
        X(3)  = 1 - X_air               !mole fraction fuel

    c   Calculate NASA coefficients of Mixture
       
        Data a / 100 * 0.0 /        !initialize species coeff array to zero
        Data a_mix / 10 * 0.0 /     !initialize mixture coeff array to zero
       
        Open(1, File = 'reactants.dat')     !open file of coeffs
         
        i = 1
        Do While (.NOT.EOF(1))

            Read(1,*)Molecule               !Reads species name
            Read(1,*)(a(i,j), j = 1,7)      !Reads 7 coefficients
            i = i + 1
               
        End Do
        Close(1)                            !close file of coeffs


    c   Calculate mole-weighted mixture coefficients

        Do j = 1,7
            Do i = 1,3
                a_mix(j) = a_mix(j) + X(i)*a(i,j)
            End Do
        End Do
       
    c   Calculate molar mass of mixture (kg/kmol)
       
        MM_mix = X(1) * (2 * 15.9994)
         &          + X(2) * 2 * 14.0067
         &          + X(3) * 66.0496
       
    c   Calculate gas constant of mixture (kJ/kg-K)
       
        R_mix = 8.314 / MM_mix


    c   ************Read in pressure history from file************************
    c   ***********pressure.dat must have units of psi, msec******************

        Data temp_t / 5000 * 0.0 /      !initialize time array to zero
        Data temp_P / 5000 * 0.0 /      !initialize press array to zero
       
        Open(2, File = 'pressure.dat') 
       
        Read(2,*)TIME,CH7           !*********REMOVE THIS LINE WHEN DONE COMPARING

        i = 0
        Do While (.NOT.EOF(2))
           
            i = i + 1
            Read(2,*) temp_t(i), temp_P(i)
           
        End Do
        Close(2)     
       
        nPT = i                 !number of data points in file
       
    c   Reduce Number of data points (by 3)
        j = 0  
        Do i = 1, nPT, 3

            j    = j + 1
            t(j) = temp_t(i)
            P(j) = temp_P(i)

        End Do

        nPT = j
       
    c   ********************Calculate Thermodynamic Properties****************

    c   Unburned Gas Temperature   
        Call T_unburned(nPT,T_0,P_0,P,a_mix,R_mix,T_u)
             
    c   Unburned Gas Specific Volume
        Call v_unburned(nPT,T_u,P,R_mix,v_u)
       
    c   Unburned Gas Specific Internal Energy
        Call u_unburned(nPT,R_mix,a_mix,T_u,u_u)





        End
    Subroutine T_unburned:

    Code (Text):
        Subroutine T_unburned(N,T_i,P_i,P_sub,A_MIX,R_MIX,T_u_sub)
    c////////////////////////////////////////////////////////////////////
    c   This subroutine is used to calculate the temperature of the
    c   isentropically compressed unburned gas. For each time step
    c   an initial guess is made. The entropy change is calculated
    c   s(T_guess,P) - s(T_0,P_0). If the change is sufficiently close
    c   zero, T_guess is accepted, else iteration is employed.
    c
    c   INPUT:
    c
    c       N           = number of elements in pressure array
    c       T_i         = initial pressure
    c       P_i         = initial pressure
    c       P_sub       = pressure array in subroutine
    c       A_MIX       = mole-weighted NASA coefficients of mixture
    c       R_MIX       = gas constant of mixture
    c      
    c   OUTPUT:
    c
    c       T_u_sub     = temperature of unburned gas  
    c
    c////////////////////////////////////////////////////////////////////
        Real*8      A_MIX(10),R_MIX, T_u_sub(1000),P_sub(1000),T_i,
         &          P_i,delta_S,dsdT,T_temp
           
        Integer     i
       
        Do i = 1,N
            !Give initial values (arbitrary)
            delta_S = 2.0           !delta_S = s(T_guess,P) - s(T_0,P_0)
            dsdT    = 1.0           !dsdT    = ds/dT

            If (i.EQ.1) Then
               
                T_temp     = T_i + 0.1          !guess that T_2 = T_i + 0.1            
               
            Else

                T_temp = T_u_sub(i-1) + 0.1     !guess T_n = T_(n-1) + 0.1

            End If 

    c   Calculate entropy change to ensure isentropic compression

            Do While(Abs(delta_S).GT.0.001)

    c   Call Entropy subroutine to calculate entropy change and ds/dT  

                Call Entropy(T_i,P_i,P_sub(i),T_temp,A_MIX,R_MIX,
         &                  delta_S,dsdT)  

                If (dsdT.EQ.0.0) Then
                    Write(*,*) 'INVALID DERIVATIVE'
                    Stop   
                Else

                    T_temp = T_temp - (delta_S/dsdT)    !decrement T_guess

                End If
           
            End Do
           
            T_u_sub(i) = T_temp     !when delta_S = 0, accept T_guess
            T_u_sub(1) = T_i        !T(1) = T_initial

          End Do
       
        Return
       
        End
    Subroutine v_unburned

    Code (Text):
        Subroutine v_unburned(N,T,P,R,v)
    c////////////////////////////////////////////////////////////////////
    c   This subroutine calculates the specific volume (m**3/kg) of
    c   the unburned gas using the ideal gas equation of state:
    c   pv = RT
    c
    c   INPUT:
    c
    c       P           = pressure array
    c       T           = temperature array
    c       R           = gas constant of mixture
    c       N           = number of elements
    c      
    c   OUTPUT:
    c
    c       v           = specific volume (m**3 / kg) array
    c
    c////////////////////////////////////////////////////////////////////
        Real*8  T(1000),P(1000),v(1000),R
        Integer i, N
       
        Write(*,*)'N (Pre Do) = ', N
        Do i = 1,N
            v(i)=R * T(i) / (P(i) * 101.325 / 14.7)
        End Do
        Write(*,*)'N (Post Do) = ', N
        Return
       
        End
     
  8. Jun 18, 2011 #7
    I don't see v_u declared (or my idea?). Also I see you declare P with a size of 2000 and in the subroutine it has a size 1000.
     
  9. Jun 18, 2011 #8
    You guys are good! Somehow, it had to do with the fact that "v_u" in
    Call v_unburned(nPT,T_u,P,R_mix,v_u)
    was never declared in the Main code.

    I have know idea how that affected the value of N ... but it did.

    Also, I don't like how you can use a variable that is not declared and Fortran does not tell you, like newer languages. I guess I just have to be more careful.

    Thanks guys!!!
     
  10. Jun 18, 2011 #9

    jtbell

    User Avatar

    Staff: Mentor

    Use IMPLICIT NONE at the beginning of each program or subroutine. This forces you to declare all variables.
     
  11. Jun 18, 2011 #10
    Nice! I like it. :smile: Thanks jt
     
  12. Jun 19, 2011 #11

    Mark44

    Staff: Mentor

    Since v_u wasn't declared in your Main routine, its default type was real, and this isn't consistent with the declared type of the corresponding formal parameter in your v_unburned subroutine, which is an array of type real*8. jtbell's recommendation that you use IMPLICIT NONE would have prevented that error from happening.

    Single-stepping through the code with a debugger would show how nPT got overwritten. I'm not as familiar with the parameter passing mechanism in Fortran as I am in C/C++, but I suspect that the actual parameters are pushed onto the stack in some way. In order for the called routine to be able to pick out the parameters, it must know the types of the parameters that were placed on the stack. Since v_u wasn't declared, its type was different from the formal parameter v, and things predictably got screwed up.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: FORTRAN: I have no idea what is happening
Loading...