converting a code from Matlab to Fortran, but Fortran diverged

[eSAFARI]

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

http://www.icofc.org

 

[mikrom]

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.

 

[eSAFARI]

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%.

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

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

http://www.icofc.org

 

[mikrom]

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.