Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

converting a code from Matlab to Fortran, but Fortran diverged

Status
Not open for further replies.

eSAFARI

Technical User
Oct 17, 2010
3
0
0
IR
thread214-1651477

I am dealing with the same problem. I translated a code from Matlab to Fortran ( CVF). The Matlab code works pretty well but in iteration loop in Fortran the error is accumulated to a unacceptable number. I used double precision in Fortran. I don't understand why Matlab can calculate the result with just a 4-5 iteration but the Fortran diverges.
Is the anyone who can help me out with this?

Ehsan Safari,
Reservoir Engineer
ICOFC
 
Maybe you could forget to initialize some variable(s), or you didn't use proper numerical method for this task, ... etc, ... etc
To find the error, you need to debug your program.
IMO, nobody could help you, without seeing an example where it happens.
 
Thank you mikrom!
what do you mean by initialization of variables?
here is the part of the code that the divergence starts. The main code is too big to post it here.
the problems occurs in the do-while loop as in Matlab the results is obtained after 4-5 iteration but in Fortran it will never get out of the loop. I am sure all the subroutine inside the loop are working properly except the subroutine VapourPhase which its results have some discrepancy of 0.000000001%.

Code:
This Is The Matlab Code
===============================================================
function [ Xi , Yi , Ki_new , nV , nL  ] = CheckFlash(n,Zi,P,T,Ki,Pci,Tci,wi,Kij,OmegaA,OmegaB,Mixture)
 
if isempty(Ki)
    %  Estimating Equilibrium Ratios
    [ Ki ] = EstimateEquilibriumRatios(T,P,Pci,Tci,wi,Mixture);
end
% Initializing Equation Of State Mixture Dependents Parameters
[ ai , bi , ui , Wi ] = feval(@PR,T,Tci,Pci,wi,OmegaA,OmegaB);
% Calculating Liquid and Vapour Phase Compositions Solving Rachford-Rice Equation
MinnV   = 1/(1-max(Ki));
MaxnV   = 1/(1-min(Ki));
nV      = (MaxnV+MinnV)/2;
[ Xi , Yi , nV , nL ] = SolveRachfordRice(Zi,Ki,nV);
% Calculating Liquid and Vapour Phase Fugacitis and Fugacity Coefficients
[ phiL ] = LiquidPhase(n,Xi,P,T,ai,bi,ui,Wi,Kij);
[ phiV ] = VapourPhase(n,Yi,P,T,ai,bi,ui,Wi,Kij);
% Calculating New Equilibrium Ratios
Ki_old = Ki;
Ki_new = phiL./phiV;
% Calculating Error
Err = sum((1-Ki_old./Ki_new).^2);
% Final Loop
Itr  = 0;

while Err > 1e-12
    Itr = Itr+1;
    
    [ Xi , Yi , nV , nL ] = SolveRachfordRice(Zi,Ki_new,nV);
    [ phiL ]              = LiquidPhase(n,Xi,P,T,ai,bi,ui,Wi,Kij);
    [ phiV ]              = VapourPhase(n,Yi,P,T,ai,bi,ui,Wi,Kij);
    
    Ki_old = Ki_new;
    Ki_new = phiL./phiV;
    Err    = sum((1-Ki_old./Ki_new).^2);
    
    if Itr == 100
        break
    end
end
 
Report(T,P,Zi,Xi,Yi,Ki_new,nV,nL)
 
end

Code:
Here is the Fortran code:
===============================================================
Subroutine CheckFlash(nComp,n,Zi,P,T,Ki,Pci,Tci,wi,Kij,OmegaA,OmegaB,Xi , Yi , Ki_new , nV , nL  )
implicit none
Integer:: i,j,num,nComp,nargin,Itr,PWR2
Real (kind=8):: n,R,a,b ,u,w,Ti,RT ,AA ,BB ,UU,WW,C1,C2 ,C3 ,C4 ,Z,P,V,Roots(3),S,aprime,bprime,uprime,RTi,&
part1,part3, part4 ,fact1,part6 ,part7,fact5,eps,T, nV , nL,Err,Mixture_GLB(100), Zi_GLB(100), Pci_GLB(100), &
Tci_GLB(100), wi_GLB(100) ,Kij_GLB(100,100) ,OmegaA_GLB(100), OmegaB_GLB(100),MinnV , MaxnV,nV1,one,mult2,mult4,&
kvalue_test(21),zi_test(21),Eps1,Eps2
Real(KIND=8),Dimension(nComp)::Zi,ai,bi,ui,Wii,fi,phi,Yii,yi,part2,part5,fact2,fact3,dArdni,alphai,&
fact4,part8,Ki,OmegaA,OmegaB,Pci,Tci,wi,Mixture,mi,aci,Xi,Ki_out,phiL,phiV,Ki_old,Ki_new,Eps3
Real(KIND=8),Dimension(nComp,nComp)::Kij
COMPLEX (KIND=8):: Rts(3)
CHARACTER:: LiquidPhaseRPT*30,VapourPhaseRPT*50,FlashRPT*50,CWD*255
character::IniKvalPlusFracMeth*20
parameter(R   = 0.0083144621,one=1.0,mult2=2.0,PWR2=2,mult4=4.0,eps = 1.0e-12)                                                                                                                             

common/IniKvalPlusFracMeth_CMN/IniKvalPlusFracMeth
common/INPUTS/  Zi_GLB, Pci_GLB, Tci_GLB, wi_GLB ,Kij_GLB ,OmegaA_GLB, OmegaB_GLB 
COMMON/Directory/CWD
 Zi     (1:nComp) =Zi_GLB     (1:nComp) 
 Pci    (1:nComp) =Pci_GLB    (1:nComp) 
 Tci    (1:nComp) =Tci_GLB    (1:nComp) 
 wi     (1:nComp) =wi_GLB     (1:nComp) 
 Kij    (1:nComp,1:nComp) =Kij_GLB    (1:nComp,1:nComp) 
 OmegaA (1:nComp) =OmegaA_GLB (1:nComp) 
 OmegaB (1:nComp) =OmegaB_GLB (1:nComp) 

if (Ki(1) < 0) then !isempty(Ki)
    !  Estimating Equilibrium Ratios
    call EstimateEquilibriumRatios(nComp,T,P,Pci,Tci,wi,Ki) 
endif

! Initializing Equation Of State Mixture Dependents Parameters
call PR(nComp,T,Tci,Pci,wi,OmegaA,OmegaB,ai , bi , ui , Wii , mi , aci,alphai)       
! Calculating Liquid and Vapour Phase Compositions Solving Rachford-Rice Equation

MinnV   = 1.0/(1.0-maxval(Ki));
MaxnV   = 1.0/(1.0-minval(Ki));
nV1      = (MaxnV+MinnV)/2.0;

	 MinnV=sum(Zi*(Ki-1.0))
	 MaxnV=sum(Zi*(Ki-1.0)/Ki)
     nV1=abs(MinnV/(MinnV-MaxnV))  !initial Nv  

call SolveRachfordRice (nComp,Zi,Ki,nV1,Xi , Yi , nV , nL)          
! Calculating Liquid and Vapour Phase Fugacitis and Fugacity Coefficients

  call LiquidPhase(nComp,n,Xi,P,T,ai,bi,ui,Wii,Kij,phiL , fi , V , a , b , u , w , Z,LiquidPhaseRPT ) ![ phiL ] = LiquidPhase(n,Xi,P,T,ai,bi,ui,Wi,Kij);


 call VapourPhase(nComp,n,Yi,P,T,ai,bi,ui,Wii,Kij,phiV , fi , V , a , b , u , w , Z,VapourPhaseRPT )  ![ phiV ] = VapourPhase(n,Yi,P,T,ai,bi,ui,Wi,Kij);
           
! Calculating New Equilibrium Ratios
Ki_old = Ki;
Ki_new = phiL/phiV;

! Calculating Error
Err = sum((one-Ki_old/Ki_new)**PWR2);

! Final Loop
nV1=nV  
Itr  = 0;

do while (Err > eps  )   ! EBI: 1e-12
    Itr = Itr+1;
    call SolveRachfordRice (nComp, Zi,Ki_new,nV1,Xi , Yi , nV , nL)   
   call LiquidPhase(nComp,n,Xi,P,T,ai,bi,ui,Wii,Kij,phiL , fi , V , a , b , u , w , Z,LiquidPhaseRPT ) ![ phiL ] = LiquidPhase(n,Xi,P,T,ai,bi,ui,Wi,Kij);
          
   call VapourPhase(nComp,n,Yi,P,T,ai,bi,ui,Wii,Kij,phiV , fi , V , a , b , u , w , Z,VapourPhaseRPT )  ![ phiV ] = VapourPhase(n,Yi,P,T,ai,bi,ui,Wi,Kij);          
    
    Ki_old = Ki_new;
    Ki_new = phiL/phiV;
    Err    = sum((one-Ki_old/Ki_new)**PWR2);
	
	
    if (Itr == 100) then
       exit !break
    endif


enddo

call Report(nComp,T,P,Zi,Xi,Yi,Ki_new ,nV,nL) 

end Subroutine CheckFlash

Ehsan Safari,
Reservoir Engineer
ICOFC
 
IMO the error could be in the computation of the arrays phiL and/or phiV.
You could try to print the arrays after each step and compare them with the same arrays computed by Matlab.
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top