FORTRAN: I have no idea what is happening

[Saladsamurai]

I am freaking out man 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

	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:

	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

 

[Hurkyl]

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?

 

[gb7nash]

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?

 

[Mark44]

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.

 

[Grieverheart]

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.

 

[Saladsamurai]

Here is the Main Program:

	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:

	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

	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

 

[Grieverheart]

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.

 

[Saladsamurai]

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!

 

[jtbell]

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.

Use IMPLICIT NONE at the beginning of each program or subroutine. This forces you to declare all variables.

 

[Saladsamurai]

Use IMPLICIT NONE at the beginning of each program or subroutine. This forces you to declare all variables.

Nice! I like it. Thanks jt

 

[Mark44]

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.

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.

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!