MATLAB Answers

Converting working Fortran code to MATLAB

4 views (last 30 days)
Joe Check
Joe Check on 10 Jun 2020
Answered: Ben Barrowes on 11 Jun 2020
Hello,
I am trying to convert a Fortran code into a Matlab code. Most of the fortran code consists of matrices / arrays such as T(i,k), u(i,k), x(i), dpdx(x).
I'm wondering what is the best way to replicate those arrays in matlab. For instance the matrix T(i,k) in the for loop, shouldn't I use the zeros command instead of writing T(i,k) = 0 which is part of my initial condtions for my problem??
I'm a programming newbie so some code modifications are not as obvious to me.
One piece of the Fortran code is as following:
"
imin = 0
imax = 100
kmin = 0
kmax = 100
%%%
do i = imin, imax, 1
x(i) = xmin + real (i - imin, kr) / real (imax - imin, kr) * (xmax - xmin)
h(i) = hm + x(i) ** 2.0_kr / (2.0_kr * Rad)
dpdx(i) = 6.0_kr * visc / h(i) ** 3.0_kr * ((U_ls + U_us) * h(i) + cons1)
do k = kmin, kmax, 1
z(k) = zmin + real (k - kmin, kr) / real (kmax - kmin, kr) * (h(i) - zmin)
T(i,k) = 0.0_kr
A1 = dpdx(i) / (2.0_kr * visc)
A2 = (- dpdx(i) / (2.0_kr * visc) * h(i) + (U_us - U_ls) / h(i))
u(i,k) = A1 * z(k) ** 2.0_kr + A2 * z(k) + U_ls
dudz(i,k) = 2.0_kr * A1 * z(k) + A2
enddo
enddo
%%%
My matlab equivalent currently looks like this:
for i = 1:imax
x(i)= xmin + (i - imin) / (imax - imin) * (xmax - xmin);
h(i) = hmin + x(i).^2 / (2 * Rad);
dpdx = 6 * visc /(h(i).^3) * ((U_ls * U_us) * h(i) + cons1);
for k = kmin:kmax
z(i) = zmin + (k - kmin) / (kmax - kmin) * (h(i) - zmin);
T(i,k) = 0
A1 = dpdx / (2 * visc);
A2 = (-dpdx / (2 * visc)) * h(i) + (U_us - U_ls) / h(i);
u(i,k) = A1 * z(i).^2 + A2 * z(i) + U_ls;
dudz(i,k) = 2 * A1 * z(i) + A2;
end
end
I'm not sure about the i and k arrays in the for loop which in the original code run from 0 to 100 but in matlab they go from 1 to 100, how do I correct this?
Please ignore the other variables I listed like hmin, visc, Rad etc. I would really appreciate it if someone could tell me that the structure of my matlab code is correct and would yield me the same output as the original Fortran code.

  1 Comment

KSSV
KSSV on 10 Jun 2020
In MATLAB the indices should be posititve integers. You can run the loop from i = 1 to 101.

Sign in to comment.

Accepted Answer

Ben Barrowes
Ben Barrowes on 11 Jun 2020
Here is your fortran code cleaned up and compilable:
program loops
integer, parameter :: kr = 8
real x(0:100), h(0:100), dpdx(0:100), z(0:100)
real T(0:100,0:100),u(0:100,0:100),dudz(0:100,0:100)
imin = 0
imax = 100
kmin = 0
kmax = 100
do i = imin, imax, 1
x(i) = xmin + real (i - imin, kr) / real (imax - imin, kr) * (xmax - xmin)
h(i) = hm + x(i) ** 2.0_kr / (2.0_kr * Rad)
dpdx(i) = 6.0_kr * visc / h(i) ** 3.0_kr * ((U_ls + U_us) * h(i) + cons1)
do k = kmin, kmax, 1
z(k) = zmin + real (k - kmin, kr) / real (kmax - kmin, kr) * (h(i) - zmin)
T(i,k) = 0.0_kr
A1 = dpdx(i) / (2.0_kr * visc)
A2 = (- dpdx(i) / (2.0_kr * visc) * h(i) + (U_us - U_ls) / h(i))
u(i,k) = A1 * z(k) ** 2.0_kr + A2 * z(k) + U_ls
dudz(i,k) = 2.0_kr * A1 * z(k) + A2
enddo
enddo
end program loops
And here is the result when I run it through f2matlab:
function loops(varargin)
persistent a1 a2 cons1 dpdx dudz firstCall h hm i imax imin k kmax kmin kr rad t u u_ls u_us visc x xmax xmin z zmin
;
if isempty(firstCall),firstCall=1;end;
if firstCall;
kr = 8;
a1=0;
a2=0;
cons1=0;
dpdx=zeros(1,100+1);
dudz=zeros(100+1,100+1);
h=zeros(1,100+1);
hm=0;
i=0;
imax=0;
imin=0;
k=0;
kmax=0;
kmin=0;
rad=0;
t=zeros(100+1,100+1);
u=zeros(100+1,100+1);
u_ls=0;
u_us=0;
visc=0;
x=zeros(1,100+1);
xmax=0;
xmin=0;
z=zeros(1,100+1);
zmin=0;
end
firstCall=0;
imin = 0;
imax = 100;
kmin = 0;
kmax = 100;
for i = imin: 1: imax
x(i+1) = xmin + real(i-imin)./real(imax-imin).*(xmax-xmin);
h(i+1) = hm + x(i+1).^2.0./(2.0.*rad);
dpdx(i+1) = 6.0.*visc./h(i+1).^3.0.*((u_ls+u_us).*h(i+1)+cons1);
for k = kmin: 1: kmax
z(k+1) = zmin + real(k-kmin)./real(kmax-kmin).*(h(i+1)-zmin);
t(i+1,k+1) = 0.0;
a1 = dpdx(i+1)./(2.0.*visc);
a2 =(-dpdx(i+1)./(2.0.*visc).*h(i+1)+(u_us-u_ls)./h(i+1));
u(i+1,k+1) = a1.*z(k+1).^2.0 + a2.*z(k+1) + u_ls;
dudz(i+1,k+1) = 2.0.*a1.*z(k+1) + a2;
end
k = fix(kmax+ 1);
end
i = fix(imax+ 1);
end %program loops
This could be further optimized and cleaned up in matlab, but this gets you working code you can profile and tinker with.
PM me if you have more code or further troubles. Posting larger codes here can get unwieldy.

  0 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by