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

Status
Not open for further replies.

#### eSAFARI

##### Technical User

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.

Replies
2
Views
150
Replies
7
Views
302
Replies
28
Views
465
Replies
11
Views
180
Replies
25
Views
308