INTELLIGENT WORK FORUMS
FOR COMPUTER PROFESSIONALS

Log In

Come Join Us!

Are you a
Computer / IT professional?
Join Tek-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • 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.

Jobs

converting a code from Matlab to Fortran, but Fortran diverged

converting a code from Matlab to Fortran, but Fortran diverged

converting a code from Matlab to Fortran, but Fortran diverged

(OP)
thread214-1651477: Fortran's precision vs. MATLAB's

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.

Reply To This Thread

Posting in the Tek-Tips forums is a member-only feature.

Click Here to join Tek-Tips and talk with other members!

Resources

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:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close