Translation of Fortran code to related Matlab code in Internal Fluid Dynamics

Hi,

I am a kind of new to Matlab. My supervisor gave me a task in Fortran. I am supposed to transfer the code into a Matlab code. The matlab code that I have written is:

clear all
close all
clc

dt=0.001;
t = 0:dt:3;
delta_x=0.01; %step size
x=0:delta_x:100;

rmass=0.5;
rmomi=1.0;
% dx=1/(ii-1);
pi=4.*atan(1);
fik=0.4;
H1=0.5;
H1D=0;
A=0;
AD=0.1;
ADinit=AD;
c1b=0.5;
c2b=1-c1b;

ub = cell(2, 1);
ini_cond = [0,0];
for i = 1:2;
 ub{i} = zeros(1,length(x));
 ub{i}(:, length(x)) = ini_cond(i) * rand(1, 1);
end    

for i=1:length (x);
%      x=i*delta_x ;
      fikness=fik*sin(pi.*x);
      ub{1}=(c1b-H1D*(x-0.5)+AD/2.*(x-0.5).^2)./(H1-0.5*fikness-A*(x-0.5));
      ub{2}=(c2b+H1D*(x-0.5)-AD/2.*(x-0.5).^2)./(1.-H1+0.5*fikness+A*(x-0.5));

end

c1=c1b;
c2=1-c1;

u = cell(2, 1);
ini_cond = [0,0];
for i = 1:2;
 u{i} = zeros(1,length(x));
 u{i}(:, length(x)) = ini_cond(i) * rand(1, 1);
end    

for i=1:length (x);
%      x=i*delta_x ;
      fikness=fik*sin(pi.*x);
      u{1}=(c1-H1D*(x-0.5)+AD/2.*(x-0.5).^2)./(H1-0.5*fikness-A*(x-0.5));
      u{2}=(c2+H1D*(x-0.5)-AD/2.*(x-0.5).^2)./(1.-H1+0.5*fikness+A*(x-0.5));

end

% p = cell(2, 1);
% q = cell(2, 1);
% % ini_cond = [0,0];
% 
% for i = 1:2;
%    p{i} = zeros(1,length(t));
%    p{i}(:, length(t)) = ini_cond(i) * rand(1, 1);
%    q{i} = zeros(1,length(t));
%    q{i}(:, length(t)) = ini_cond(i) * rand(1, 1);
% end    

p{1}(1)=0.5*(1.-u{1}(1).^2);
q{1}(1)=0;

for i=2:length(x)

q{1}(i)=q{1}(i-1)-delta_x*(u{1}(i-1)-ub{1}(i-1))./dt;
p{1}(i)=0.5*(1.-u{1}(i).^2)+q{1}(i);
end
p{2}(1)=0.5*(1.-u{2}(1).^2);
q{2}(1)=0;

for i=2:length(x)

q{2}(i)=q{2}(i-1)-delta_x*(u{2}(i-1)-ub{2}(i-1))./dt;
p{2}(i)=0.5*(1.-u{2}(i).^2)+q{2}(i);
end

sum = cell(2, 1);
st = cell(2, 1);
ini_cond = [0,0];

for i = 1:2;
 sum{i} = zeros(1,length(t));
 sum{i}(:, length(t)) = ini_cond(i) * rand(1, 1);
 st{i} = zeros(1,length(t));
 st{i}(:, length(t)) = ini_cond(i) * rand(1, 1);
end 

% for j=1:length (t);
% st(j)=p{1}(j)-p{2}(j);
% end

dc=0.001;
c1=c1+dc;
c1=(c1*st(1)-(c1-dc)*st(2))/(st(1)-st(2));

sum{1}(1)=0.5*(p{2}(1)-p{1}(1));
sum{2}(1)=0.5*(p{2}(1)-p{1}(1)).*(-1/2);
ADDOT(1)=sum{2}(1)*delta_x./rmomi;
H1DDOT(1)=-sum{1}(1).*delta_x./rmass;
H1(i)=H1(i-1)+dt*H1D(i-1);
ADDOT(i)=sum{2}(i)*delta_x./rmomi;
AD(i)=AD(i-1)+dt*ADDOT(i-1);
A(i)=A(i-1)+dt*AD(i-1);

for i=2:length(t)
% x=(i-1)*dx;
sum{1}(i)=sum{1}(i-1)+(p{2}(i)-p{1}(i));
sum{2}(i)=sum{2}(i-1)+(p{2}(i)-p{1}(i)).*(x-1/2);
H1DDOT(i)=-sum{1}(i).*delta_x./rmass;
H1D(i)=H1D(i-1)+dt*H1DDOT;
H1(i)=H1(i-1)+dt*H1D(i-1);
ADDOT(i)=sum{2}(i)*delta_x./rmomi;
AD(i)=AD(i-1)+dt*ADDOT(i-1);
A(i)=A(i-1)+dt*AD(i-1);
end

H1L=H1+A*1/2;
H1R=H1-A*1/2;
H2=1.-H1;
rat1=AD/ADinit;
rat2=0;
rat2=ADDOT/AD;

for i=1:length(x)
ub{1}(i)=u{1}(i);
ub{2}(i)=u{2}(i);
c1b=c1;
end

figure (1);
subplot(1,2,1);
plot(t,H1D,'.:',t,H1DDOT,'.:',t,c2,'*');
axis([0 1 0 2.5]);
xlabel('time');
ylabel('body leading and trailing edges');
hold on
subplot(1,2,2);
plot(t,A,t,AD,t,ADDOT);
xlabel('time');
ylabel('the angles');

figure (2);
subplot(1,2,1);
plot(x,u{1},x,u{2},':', 'LineWidth',1.5);
axis([0 10 0 2.5]);
xlabel('x lateral position');
ylabel('velocity profiles u1, u2 in the gaps');
hold on
subplot(1,2,2);
plot(x,p{1},':',x, p{2},':m', 'LineWidth', 2);
axis([0 10 -3 1]);
xlabel('x lateral position');
ylabel('pressure p1, p2 in the gaps');
hold on
suptitle('1 grain in a fluid (N=2)');

But 'sum' and 'st' parts (Newton-Raphson Method) do not work.

The fortran code is this:

program grab3
!!from grair2.f90 on 09.01.2009                 
!!1 grain in air between 2 walls: inviscid
!!iterates more than grair1
!!Here A = - theta

implicit double precision (a-h,o-z)
common u1(1001),u2(1001),u1b(1001),u2b(1001)&
,p1(1001),q1(1001),p2(1001),q2(1001),st(9)

ii=101
b=1.
dx=1./(ii-1.)
dt=0.0001
tmax=3.000000001
tmax=2.3999999999
!!tmax=4.500000001
!!tmax=dt*1.000001
nt=0
ntp=1000
ntp1=200
tstart=1.8999999999999
!!tstart=4.2999999999999
rmass=1.
rmomi=2.
!!rmomi=1.
!!rmass=10.
!!rmomi=20.

pi=4.*atan(1.)
!!thickness-effect parameter  --actually camber in this prog
fik=0.4
fik=-0.4

t=0.
H1=0.5
H1D=0.
A=0.
AD=0.1
!!AD=0.001
ADinit=AD
c1b=0.5
c2b=1.-c1b
do 1 i=1,ii
x=(i-1.)*dx
fikness=fik*sin(pi*x)
u1b(i)=(c1b-H1D*(x-b/2.)+AD/2.*(x-b/2.)**2)/(H1-0.5*fikness-A*(x-b/2.))
u2b(i)=(c2b+H1D*(x-b/2.)-AD/2.*(x-b/2.)**2)/(1.-H1+0.5*fikness+A*(x-b/2.))
!!note sign above
1 continue
dc=0.001

88 t=t+dt
nt=nt+1
mmm=1
c1=c1b
666 m=1
77 c2=1.-c1

do 3 i=1,ii
x=(i-1.)*dx
fikness=fik*sin(pi*x)
u1(i)=(c1-H1D*(x-b/2.)+AD/2.*(x-b/2.)**2)/(H1-0.5*fikness-A*(x-b/2.))
u2(i)=(c2+H1D*(x-b/2.)-AD/2.*(x-b/2.)**2)/(1.-H1+0.5*fikness+A*(x-b/2.))
3 continue

p1(1)=0.5*(1.-u1(1)**2)
q1(1)=0.
do 4 i=2,ii
q1(i)=q1(i-1)-dx*(u1(i-1)-u1b(i-1))/dt
p1(i)=0.5*(1.-u1(i)**2)+q1(i)
4 continue
p2(1)=0.5*(1.-u2(1)**2)
q2(1)=0.
do 44 i=2,ii
q2(i)=q2(i-1)-dx*(u2(i-1)-u2b(i-1))/dt
p2(i)=0.5*(1.-u2(i)**2)+q2(i)
44 continue
!!write(6,*)m,m,m,m
st(m)=p1(ii)-p2(ii)
!!write(6,*)m,m
m=m+1
if(m.eq.3)go to 51
if(m.eq.4)go to 52
c1=c1+dc
go to 77
51 c1=(c1*st(1)-(c1-dc)*st(2))/(st(1)-st(2))
go to 77
52 continue 
mmm=mmm+1
if(mmm.eq.2)go to 666
!!write(6,*)m,m,m,m
!!really need to iterate more
!!write(6,*)c1,st(3)  

sum1=0.5*(p2(1)-p1(1))
sum2=0.5*(p2(1)-p1(1))*(-b/2.)
!!write(6,*)m,m,m,m
do 93 i=2,ii-1
x=(i-1.)*dx
sum1=sum1+(p2(i)-p1(i))
93 sum2=sum2+(p2(i)-p1(i))*(x-b/2.)
!!write(6,*)m,m,m,m
H1DDOT=-sum1*dx/rmass
H1D=H1D+dt*H1DDOT
H1=H1+dt*H1D
ADDOT=sum2*dx/rmomi
AD=AD+dt*ADDOT
A=A+dt*AD
!!write(6,*)m,m,m,m,m
!!if(((nt/ntp)*ntp).ne.nt)go to 5499
if(((nt/ntp)*ntp).ne.nt)go to 7021
write(6,*)t,H1,A,st(3) 

H1L=H1+A*b/2.
H1R=H1-A*b/2.
H2=1.-H1
rat1=AD/ADinit
rat2=0.
if(abs(AD).lt.0.0000001)go to 9753
rat2=ADDOT/AD
9753 continue
write(5,*)t,H1L,H1R,c1
close(5)
open(5,file='g3bD.dat',status='unknown',position='append')
write(5,*)t,A,AD,ADDOT
close(5)

7021 continue
!!if(ii.gt.0)go to 5499
if(t.lt.tstart)go to 7023
if(((nt/ntp1)*ntp1).ne.nt)go to 7023
do 7022 i=1,ii
x=(i-1.)*dx
fikness=fik*sin(pi*x)
gap1=(H1-0.5*fikness-A*(x-b/2.))
gap2=(1.-H1+0.5*fikness+A*(x-b/2.))
open(5,file='g3cD.dat',status='unknown',position='append')
write(5,*)x,gap1,gap2
close(5)
7022 continue

7023 continue
if(ii.gt.0)go to 5499
if(t.lt.tstart)go to 5499
if(((nt/ntp1)*ntp1).ne.nt)go to 5499
do 702 i=1,ii
x=(i-1.)*dx
open(5,file='gibaS.dat',status='unknown',position='append')
write(5,*)x,u1(i),u2(i)
close(5)
open(5,file='gibbS.dat',status='unknown',position='append')
write(5,*)x,p1(i),p2(i)
close(5)
702 continue

5499 if(t.gt.tmax)go to 549
do 801 i=1,ii
u1b(i)=u1(i)
u2b(i)=u2(i)
801 continue
c1b=c1
go to 88

549 continue
stop
end

Could you please please please help me..

Answer 1

[Walter Roberson]

naming a variable "sum" almost always leads to trouble. Avoid naming variables the same as commonly-used MATLAB library routines.

"clear all" should only be used by experts.

You might want to try http://sourceforge.net/projects/f2matlab/ to do the conversion.

[Meva]

Thank you very much. Can I apply the newton_raphson numerical algorith inside the integral. Actually, in the code ADDOT and H1DDOT represent the integrals and sum1 and sum2 are inside of those integrals.