INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Tek-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

#### Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

# converting a code from Matlab to Fortran, but Fortran diverged

## converting a code from Matlab to Fortran, but Fortran diverged

(OP)

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
www.icofc.org

### RE: converting a code from Matlab to Fortran, but Fortran diverged

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.

### RE: converting a code from Matlab to Fortran, but Fortran diverged

(OP)
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 --> matlab

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 --> Fortran

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
www.icofc.org

### RE: converting a code from Matlab to Fortran, but Fortran diverged

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.

#### Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

#### Red Flag Submitted

Thank you for helping keep Tek-Tips Forums free from inappropriate posts.
The Tek-Tips staff will check this out and take appropriate action.

Close Box

# Join Tek-Tips® Today!

Join your peers on the Internet's largest technical computer professional community.
It's easy to join and it's free.

Here's Why Members Love Tek-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!