How to convert MATLAB code to C++ code?

How do I convert my MATLAB code to a C++ code?

[Ankit Singh Jaswal]

%MAC algorithm based code developed by Rohit Kanchi

%Contact: [phone]

%All units are in SI System kg,m,s,K,Pa

clc

clear

close all

savepath='D:\CFD_MEE4006\Project\2D MAC Algorithm FDM Code\Simulation_Files\';

fprintf("This is an FDM MAC algorithm based code to predict 2-D unsteady incompressible viscous flows. \n\n");

fprintf("Please make sure to enter the input parameters correctly");

nx=100;

ny=40;

deltax=1;

deltay=1;

a=deltax;

b=deltay;

deltaT=0.001;

timesteps=1000;

mu=10^(-3);

rho=1000;

%Grid Point Generator

[Pgrid,Ugrid,Vgrid]=gridpointgenerator(nx,ny,deltax,deltay);

%Boundary Conditions:

%Top: Slip Wall [Axis of rotation]

%Bottom: No Slip Wall, Adiabatic

%Left: Velocity Inlet

Inletvel=0.1;

%Right: Pressure Boundary Condition

ExitP=101325;

%Grid Points Initialization at t=0s

%Pressure Grid Points:

l=length(Pgrid(:,1));

for i=1:l

Pgrid(i,3)=ExitP;

end

%u vel grid points:

l=length(Ugrid(:,1));

for i=1:l

if Ugrid(i,1)==(-(deltax/2))

Ugrid(i,3)=Inletvel;

else

Ugrid(i,3)=0;

end

end

%v vel grid points:

l=length(Vgrid(:,1));

for i=1:l

Vgrid(i,3)=0;

end

%%%% START MAC ALGORITHM

utild=Ugrid+0; %Initializing Provisional Velocity matrices

vtild=Vgrid+0;

Pprime=Pgrid+0;

Ugridtemp=Ugrid+0;

Vgridtemp=Vgrid+0;

for timestep=1:timesteps

flag=0; %This var, if zero, means residual are not sufficiently low.

%When the residual comes down in the while loop below, a break-

%-statement is encountered and the soln proceeds to next time-

%-step after updating values

iter=1;

while flag==0

%Calculate provisional velocities:

%Calculate U~

for i=1:length(Ugrid(:,1))

if Ugrid(i,1)==(-(deltax/2)) %Left Boundary

continue

elseif Ugrid(i,1)==((nx*deltax)-(deltax/2)) %Right Boundary

continue

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif Ugrid(i,2)==0 && Ugrid(i,1)~=((nx*deltax)-(deltax/2)) %Bottom Boundary

u1=Ugrid(i,3); %i,j

u2=Ugrid(i+1,3); %i+1,j

u3=Ugrid(i-1,3); %i-1,j

%i,j+1

for k=i:length(Ugrid(:,1))

if(Ugrid(k,1))==(Ugrid(i,1))

u4=Ugrid(k,3);

break

end

end

%i,j-1

u5=0-u1;

%v Velocity points for averaging

for k=1:length(Vgrid(:,1))

if Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg1=Vgrid(k,3); %i,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg2=Vgrid(k,3); %i+1,j

end

end

vforavg3=0;

vforavg4=0;

%Calculate Pressures:

for k=1:length(Pgrid(:,1))

if Pgrid(k,1)==Ugrid(i,1)+(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P1up=Pgrid(k,3); %P i+1,j

end

if Pgrid(k,1)==Ugrid(i,1)-(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P2up=Pgrid(k,3); %P i,j

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif Ugrid(i,2)==((ny*deltay)-deltay) && Ugrid(i,1)~=((nx*deltax)-(deltax/2)) %Top Boundary

u1=Ugrid(i,3); %i,j

u2=Ugrid(i+1,3); %i+1,j

u3=Ugrid(i-1,3); %i-1,j

%i,j+1

u4=u1+0;

%i,j-1

for k=i:-1:1

if(Ugrid(k,1))==(Ugrid(i,1))

u5=Ugrid(k,3);

break

end

end

%v Velocity points for averaging

for k=1:length(Vgrid(:,1))

if Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg1=Vgrid(k,3); %i,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg2=Vgrid(k,3); %i+1,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)

vforavg3=Vgrid(k,3); %i+1,j-1

elseif Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)

vforavg4=Vgrid(k,3); %i-1,j-1

end

end

%Calculate Pressures:

for k=1:length(Pgrid(:,1))

if Pgrid(k,1)==Ugrid(i,1)+(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P1up=Pgrid(k,3); %P i+1,j

end

if Pgrid(k,1)==Ugrid(i,1)-(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P2up=Pgrid(k,3); %P i,j

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

else %Interior Grid Points

u1=Ugrid(i,3); %i,j

u2=Ugrid(i+1,3); %i+1,j

u3=Ugrid(i-1,3); %i-1,j

%i,j+1

for k=i:length(Ugrid(:,1))

if(Ugrid(k,1))==(Ugrid(i,1))

u4=Ugrid(k,3);

break

end

end

%i,j-1

for k=i:-1:1

if(Ugrid(k,1))==(Ugrid(i,1))

u5=Ugrid(k,3);

break

end

end

%v Velocity points for averaging

for k=1:length(Vgrid(:,1))

if Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg1=Vgrid(k,3); %i,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg2=Vgrid(k,3); %i+1,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)

vforavg3=Vgrid(k,3); %i+1,j-1

elseif Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)

vforavg4=Vgrid(k,3); %i-1,j-1

end

end

%Calculate Pressures:

for k=1:length(Pgrid(:,1))

if Pgrid(k,1)==Ugrid(i,1)+(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P1up=Pgrid(k,3); %P i+1,j

end

if Pgrid(k,1)==Ugrid(i,1)-(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P2up=Pgrid(k,3); %P i,j

end

end

end

delcU= rho*((u1*((u2-u3)/(2*a))) + (((vforavg1+vforavg2+vforavg3+vforavg4)/4)*((u4-u5)/(2*b))));

deldU= mu*( ((u2-(2*u1)+u3)/(a^2)) + ((u4-(2*u1)+u5)/(b^2)));

%Calc provisional Up:

utild(i,3)= u1 + ((deltaT/rho)* (deldU-delcU + ((P2up-P1up)/a)));

end

%Calculate V~

for i=1:length(Vgrid(:,1))

if Vgrid(i,1)==0 || Vgrid(i,1)==((nx-1)*deltax) || Vgrid(i,2)==(0-(deltay/2)) || Vgrid(i,2)==((ny*deltay)-(deltay/2)) %left, right, bottom and top boundaries

continue

else %Interior grid points

v1=Vgrid(1,3); %i,j

for k=1:length(Vgrid(:,1)) %i,j+1

if Vgrid(k,1)==Vgrid(i,1) && Vgrid(k,2)==(Vgrid(i,2)+deltay)

v2=Vgrid(k,3);

break;

end

end

for k=1:length(Vgrid(:,1)) %i,j-1

if Vgrid(k,1)==Vgrid(i,1) && Vgrid(k,2)==(Vgrid(i,2)-deltay)

v3=Vgrid(k,3);

break;

end

end

v4=Vgrid(i+1,3); %i+1,j

v5=Vgrid(i-1,3); %i-1,j

%calc u vel avg

for k=1:length(Ugrid(:,1))

if Ugrid(k,1)==Vgrid(i,1)-(deltax/2) && Ugrid(k,2)==Vgrid(i,2)+(deltay/2)

uforavg1=Ugrid(k,3); %i,j

elseif Ugrid(k,1)==Vgrid(i,1)+(deltax/2) && Ugrid(k,2)==Vgrid(i,2)+(deltay/2)

uforavg2=Ugrid(k,3); %i+1,j

elseif Ugrid(k,1)==Vgrid(i,1)+(deltax/2) && Ugrid(k,2)==Vgrid(i,2)-(deltay/2)

uforavg3=Ugrid(k,3); %i+1,j-1

elseif Ugrid(k,1)==Vgrid(i,1)-(deltax/2) && Ugrid(k,2)==Vgrid(i,2)-(deltay/2)

uforavg4=Ugrid(k,3); %i-1,j-1

end

end

%Calc Pressures

for k=1:length(Pgrid(:,1))

if Pgrid(k,1)==Vgrid(i,1) && Pgrid(k,2)==(Vgrid(i,2)-(deltay/2))

P1vp=Pgrid(k,3); %P i,j

end

if Pgrid(k,1)==Vgrid(i,1) && Pgrid(k,2)==(Vgrid(i,2)+(deltay/2))

P2vp=Pgrid(k,3); %P i,j+1

end

end

end

delcV=rho*( (v1*((v2-v3)/(2*b))) + (((uforavg1+uforavg2+uforavg3+uforavg4)/4)*((v4-v5)/(2*a))) );

deldV=mu*( ((v4-(2*v1)+v5)/(a^2)) + ((v2-(2*v1)+v3)/(b^2)) );

vtild(i,3)=v1 + ((deltaT/rho)* (deldV-delcV + ((P1vp-P2vp)/b)));

end

%Calculate P' field

for i=1:length(Pgrid(:,1))

px1=Pgrid(i,1)+(deltax/2);

px2=Pgrid(i,1)-(deltax/2);

py1=Pgrid(i,2)+(deltay/2);

py2=Pgrid(i,2)-(deltay/2);

for k=1:length(utild(:,1))

if utild(k,1)==px1 && utild(k,2)==Pgrid(i,2)

up=utild(k,3);

end

if utild(k,1)==px2 && utild(k,2)==Pgrid(i,2)

uw=utild(k,3);

end

end

for k=1:length(vtild(:,1))

if vtild(k,1)==Pgrid(i,1) && vtild(k,2)==py1

vp=vtild(k,3);

end

if vtild(k,1)==Pgrid(i,1) && vtild(k,2)==py2

vs=vtild(k,3);

end

end

Pprime(i,3)= -(( ((up-uw)/a) +((vp-vs)/b) )/( ((2*deltaT)/rho)*((1/(a^2)) + (1/(b^2))) ));

end

%Updating utild and vtild to n+1 temporarily to check for-

%-convergence

for i=1:length(Ugrid(:,1))

if Ugrid(i,1)==(-(deltax/2)) || Ugrid(i,1)==((nx*deltax)-(deltax/2)) %Left Boundary and right boundary

continue

end

for k=1:length(Pprime(:,1))

if Pprime(k,1)==Ugrid(i,1)-((deltax)/2) && Pprime(k,2)==Ugrid(i,2)

Ugridtemp(i,3)=utild(i,3) - ( (deltaT/(rho*a))*(Pprime(k+1,3)-Pprime(k,3)));

end

end

end

%For right boundary we need special treatment for uvel:

for i=1:length(Ugrid(:,1))

if Ugrid(i,1)==((nx*deltax)-(deltax/2))

Ugridtemp(i,3)=Ugridtemp(i-1,3);

end

end

for i=1:length(Vgrid(:,1))

if Vgrid(i,1)==0 || Vgrid(i,1)==((nx-1)*deltax) || Vgrid(i,2)==(0-(deltay/2)) || Vgrid(i,2)==((ny*deltay)-(deltay/2)) %left, right, bottom and top boundaries

continue

end

for k=1:length(Pprime(:,1))

if Pprime(k,1)==Vgrid(i,1) && Pprime(k,2)==(Vgrid(i,2)-(deltay/2))

paa=Pprime(k,3);

elseif Pprime(k,1)==Vgrid(i,1) && Pprime(k,2)==(Vgrid(i,2)+(deltay/2))

pbb=Pprime(k,3);

end

end

Vgridtemp(i,3)=vtild(i,3) - ( (deltaT/(rho*b))*(pbb-paa));

end

%For right and top boundary we need special treatment for vvel:

for i=1:length(Vgrid(:,1))

if Vgrid(i,1)==((nx-1)*deltax) %right

Vgridtemp(i,3)=Vgrid(i-1,3);

elseif Vgrid(i,1)==((ny*deltay)-(deltay/2)) %top

for k=1:length(Vgrid(:,1))

if Vgrid(k,1)==Vgrid(i,1) && Vgrid(k,2)==(Vgrid(i,2)-deltay)

Vgridtemp(i,3)=Vgrid(k,3);

end

end

end

end

%Pressure correction

for i=1:length(Pgrid(:,1))

if Pgrid(i,1)==((nx-1)*deltax)

continue

else

Pgrid(i,3)=Pgrid(i,3)+Pprime(i,3);

end

end

%To ensure dp/dy at boundary is zero:

%Check for convergence using temporary values:

%[Convergence Check: L1 norm]

sum=0;

for i=1:length(Pgrid(:,1))

px1=Pgrid(i,1)+(deltax/2);

px2=Pgrid(i,1)-(deltax/2);

py1=Pgrid(i,2)+(deltay/2);

py2=Pgrid(i,2)-(deltay/2);

for k=1:length(Ugridtemp(:,1))

if Ugridtemp(k,1)==px1 && Ugridtemp(k,2)==Pgrid(i,2)

up=Ugridtemp(k,3);

end

if Ugridtemp(k,1)==px2 && Ugridtemp(k,2)==Pgrid(i,2)

uw=Ugridtemp(k,3);

end

end

for k=1:length(Vgridtemp(:,1))

if Vgridtemp(k,1)==Pgrid(i,1) && Vgridtemp(k,2)==py1

vp=Vgridtemp(k,3);

end

if Vgridtemp(k,1)==Pgrid(i,1) && Vgridtemp(k,2)==py2

vs=Vgridtemp(k,3);

end

end

residual=( ((up-uw)/a) +((vp-vs)/b) );

sum=sum+abs(residual);

end

m=sum/length(Pgrid(:,1));

disp(m);

%disp(iter);

if m<=0.001 %Converged

flag=1;

Vgrid=Vgridtemp; %Proceeding to n+1

Ugrid=Ugridtemp;

else %Not converged

iter=iter+1;

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Proceeding to next time step n -> n+1:

%The points in the u and v grids are first updated.

%While updation, points such as the inlet grid points for uvel are-

%-not touched because they are constant vel points and is a -

%-dirichlet bc.

%Similarly the Pgrid is also updated and the exit P points are not-

%-touched because it is a constant pressure Dirichlet BC

%Save files at n+1 for the current time-step

if mod((timestep*deltaT),0.01)==0

file_name = sprintf('U_vel_timestep_%d.txt',timestep);

id=fopen([savepath file_name],'w');

for i=1:length(Ugrid(:,1))

fprintf(id,'%f %f %f\n',Ugrid(i,1),Ugrid(i,2),Ugrid(i,3));

end

fclose(id);

file_name = sprintf('V_vel_timestep_%d.txt',timestep);

id=fopen([savepath file_name],'w');

for i=1:length(Vgrid(:,1))

fprintf(id,'%f %f %f\n',Vgrid(i,1),Vgrid(i,2),Vgrid(i,3));

end

fclose(id);

file_name = sprintf('P_timestep_%d.txt',timestep);

id=fopen([savepath file_name],'w');

for i=1:length(Pgrid(:,1))

fprintf(id,'%f %f %f\n',Pgrid(i,1),Pgrid(i,2),Pgrid(i,3));

end

fclose(id);

%Saving complete and proceed to next time-step

end

end

function [Pgrid,Ugrid,Vgrid] = gridpointgenerator(nx,ny,deltax,deltay)

% Grid point generator

%Row by row, upward

%Pressure collocated grid points

Pgridx(1)=0;

Pgridy(1)=0;

for i=2:nx

Pgridx(i)=Pgridx(i-1)+deltax;

end

for i=2:ny

Pgridy(i)=Pgridy(i-1)+deltay;

end

count=1;

for i=1:ny

for j=1:nx

Pgridxp(count)=Pgridx(j);

Pgridyp(count)=Pgridy(i);

count=count+1;

end

end

%x vel collocated grid points

nx2=nx+1;

ugridx(1)=-(deltax/2);

ugridy(1)=0;

for i=2:nx2

ugridx(i)=ugridx(i-1)+deltax;

end

for i=2:ny

ugridy(i)=ugridy(i-1)+deltay;

end

count=1;

for i=1:ny

for j=1:nx2

ugridxp(count)=ugridx(j);

ugridyp(count)=ugridy(i);

count=count+1;

end

end

%y vel collocated grid points

ny2=ny+1;

vgridy(1)=-(deltay/2);

vgridx(1)=0;

for i=2:nx

vgridx(i)=vgridx(i-1)+deltax;

end

for i=2:ny2

vgridy(i)=vgridy(i-1)+deltay;

end

count=1;

for i=1:ny2

for j=1:nx

vgridxp(count)=vgridx(j);

vgridyp(count)=vgridy(i);

count=count+1;

end

end

Pgridxp=Pgridxp';

Pgridyp=Pgridyp';

ugridxp=ugridxp';

ugridyp=ugridyp';

vgridxp=vgridxp';

vgridyp=vgridyp';

figure

plot(Pgridxp,Pgridyp,'o','MarkerSize',5,'Markeredgecolor','k','markerfacecolor','k');

hold on

plot(ugridxp,ugridyp,'x','MarkerSize',5,'Markeredgecolor','b','markerfacecolor','b');

hold on

plot(vgridxp,vgridyp,'*','MarkerSize',5,'Markeredgecolor','g','markerfacecolor','g');

title('FDM Grid Points')

legend('Pressure Points','u vel points','v vel points');

Pgrid=[Pgridxp,Pgridyp];

Ugrid=[ugridxp,ugridyp];

Vgrid=[vgridxp,vgridyp];

end

Answer 1

[Ameer Hamza]

You will need MATLAB Coder to perform the conversion. If it is already installed, you can use MATLAB coder app under Apps tab in MATLAB. It will guide you through steps to generate C or C++ code.

[PatrizioGraziosi]

Hi Ameer,

thanks for the answer.

There are mainly two reasons, as far as I see, to port a code from MATLAB to C++:

1) to speed it up by using it inside MATLAB environment,

2) to use it outside MATLAB, without having MATLAB installed.

MATLAB Coder should work for the first, does it work also for the second?

I don' t have any experience in this task but I'm going to face this challenge with a big code (tens of thousands of MATLAB code lines in many files) in the forthcoming months.

Thanks in advance

Patrizio

[Ameer Hamza]

Yes, you can use the generated C/C++ code in external projects without the need to have MATLAB installed. See here: https://www.mathworks.com/help/coder/deployment.html?s_tid=CRUX_lftnav, also check this example: https://www.mathworks.com/help/coder/ug/generate-and-modify-an-example-cc-main-function.html.

[PatrizioGraziosi]

Hi Ameer,

thank you!

the C code generated in this way can be used in any platform, while the mex files called in matlab are sensitive to the operating system, isn't it?

Thank

Patrizio

[Ameer Hamza]

Yes, mex files are platform specific. I think that the generated C code should work on other platforms too, but I am not entirely sure.

Answer 2

[Ahmad Nadeem]

I want to Change MATLAB code to C++ can anyone tell me how can I convert this? I am working on phase field simulation

[Raghu Boggavarapu]

As mentioned above, please use MATLAB Coder - MATLAB (mathworks.com) for the purpose.

[Ankit Singh Jaswal]

%MAC algorithm based code developed by Rohit Kanchi

%Contact: [phone]

%All units are in SI System kg,m,s,K,Pa

clc

clear

close all

savepath='D:\CFD_MEE4006\Project\2D MAC Algorithm FDM Code\Simulation_Files\';

fprintf("This is an FDM MAC algorithm based code to predict 2-D unsteady incompressible viscous flows. \n\n");

fprintf("Please make sure to enter the input parameters correctly");

nx=100;

ny=40;

deltax=1;

deltay=1;

a=deltax;

b=deltay;

deltaT=0.001;

timesteps=1000;

mu=10^(-3);

rho=1000;

%Grid Point Generator

[Pgrid,Ugrid,Vgrid]=gridpointgenerator(nx,ny,deltax,deltay);

%Boundary Conditions:

%Top: Slip Wall [Axis of rotation]

%Bottom: No Slip Wall, Adiabatic

%Left: Velocity Inlet

Inletvel=0.1;

%Right: Pressure Boundary Condition

ExitP=101325;

%Grid Points Initialization at t=0s

%Pressure Grid Points:

l=length(Pgrid(:,1));

for i=1:l

Pgrid(i,3)=ExitP;

end

%u vel grid points:

l=length(Ugrid(:,1));

for i=1:l

if Ugrid(i,1)==(-(deltax/2))

Ugrid(i,3)=Inletvel;

else

Ugrid(i,3)=0;

end

end

%v vel grid points:

l=length(Vgrid(:,1));

for i=1:l

Vgrid(i,3)=0;

end

%%%% START MAC ALGORITHM

utild=Ugrid+0; %Initializing Provisional Velocity matrices

vtild=Vgrid+0;

Pprime=Pgrid+0;

Ugridtemp=Ugrid+0;

Vgridtemp=Vgrid+0;

for timestep=1:timesteps

flag=0; %This var, if zero, means residual are not sufficiently low.

%When the residual comes down in the while loop below, a break-

%-statement is encountered and the soln proceeds to next time-

%-step after updating values

iter=1;

while flag==0

%Calculate provisional velocities:

%Calculate U~

for i=1:length(Ugrid(:,1))

if Ugrid(i,1)==(-(deltax/2)) %Left Boundary

continue

elseif Ugrid(i,1)==((nx*deltax)-(deltax/2)) %Right Boundary

continue

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif Ugrid(i,2)==0 && Ugrid(i,1)~=((nx*deltax)-(deltax/2)) %Bottom Boundary

u1=Ugrid(i,3); %i,j

u2=Ugrid(i+1,3); %i+1,j

u3=Ugrid(i-1,3); %i-1,j

%i,j+1

for k=i:length(Ugrid(:,1))

if(Ugrid(k,1))==(Ugrid(i,1))

u4=Ugrid(k,3);

break

end

end

%i,j-1

u5=0-u1;

%v Velocity points for averaging

for k=1:length(Vgrid(:,1))

if Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg1=Vgrid(k,3); %i,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg2=Vgrid(k,3); %i+1,j

end

end

vforavg3=0;

vforavg4=0;

%Calculate Pressures:

for k=1:length(Pgrid(:,1))

if Pgrid(k,1)==Ugrid(i,1)+(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P1up=Pgrid(k,3); %P i+1,j

end

if Pgrid(k,1)==Ugrid(i,1)-(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P2up=Pgrid(k,3); %P i,j

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif Ugrid(i,2)==((ny*deltay)-deltay) && Ugrid(i,1)~=((nx*deltax)-(deltax/2)) %Top Boundary

u1=Ugrid(i,3); %i,j

u2=Ugrid(i+1,3); %i+1,j

u3=Ugrid(i-1,3); %i-1,j

%i,j+1

u4=u1+0;

%i,j-1

for k=i:-1:1

if(Ugrid(k,1))==(Ugrid(i,1))

u5=Ugrid(k,3);

break

end

end

%v Velocity points for averaging

for k=1:length(Vgrid(:,1))

if Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg1=Vgrid(k,3); %i,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg2=Vgrid(k,3); %i+1,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)

vforavg3=Vgrid(k,3); %i+1,j-1

elseif Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)

vforavg4=Vgrid(k,3); %i-1,j-1

end

end

%Calculate Pressures:

for k=1:length(Pgrid(:,1))

if Pgrid(k,1)==Ugrid(i,1)+(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P1up=Pgrid(k,3); %P i+1,j

end

if Pgrid(k,1)==Ugrid(i,1)-(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P2up=Pgrid(k,3); %P i,j

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

else %Interior Grid Points

u1=Ugrid(i,3); %i,j

u2=Ugrid(i+1,3); %i+1,j

u3=Ugrid(i-1,3); %i-1,j

%i,j+1

for k=i:length(Ugrid(:,1))

if(Ugrid(k,1))==(Ugrid(i,1))

u4=Ugrid(k,3);

break

end

end

%i,j-1

for k=i:-1:1

if(Ugrid(k,1))==(Ugrid(i,1))

u5=Ugrid(k,3);

break

end

end

%v Velocity points for averaging

for k=1:length(Vgrid(:,1))

if Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg1=Vgrid(k,3); %i,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)+(deltay/2)

vforavg2=Vgrid(k,3); %i+1,j

elseif Vgrid(k,1)==Ugrid(i,1)+(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)

vforavg3=Vgrid(k,3); %i+1,j-1

elseif Vgrid(k,1)==Ugrid(i,1)-(deltax/2) && Vgrid(k,2)==Ugrid(i,2)-(deltay/2)

vforavg4=Vgrid(k,3); %i-1,j-1

end

end

%Calculate Pressures:

for k=1:length(Pgrid(:,1))

if Pgrid(k,1)==Ugrid(i,1)+(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P1up=Pgrid(k,3); %P i+1,j

end

if Pgrid(k,1)==Ugrid(i,1)-(deltax/2) && Pgrid(k,2)==Ugrid(i,2)

P2up=Pgrid(k,3); %P i,j

end

end

end

delcU= rho*((u1*((u2-u3)/(2*a))) + (((vforavg1+vforavg2+vforavg3+vforavg4)/4)*((u4-u5)/(2*b))));

deldU= mu*( ((u2-(2*u1)+u3)/(a^2)) + ((u4-(2*u1)+u5)/(b^2)));

%Calc provisional Up:

utild(i,3)= u1 + ((deltaT/rho)* (deldU-delcU + ((P2up-P1up)/a)));

end

%Calculate V~

for i=1:length(Vgrid(:,1))

if Vgrid(i,1)==0 || Vgrid(i,1)==((nx-1)*deltax) || Vgrid(i,2)==(0-(deltay/2)) || Vgrid(i,2)==((ny*deltay)-(deltay/2)) %left, right, bottom and top boundaries

continue

else %Interior grid points

v1=Vgrid(1,3); %i,j

for k=1:length(Vgrid(:,1)) %i,j+1

if Vgrid(k,1)==Vgrid(i,1) && Vgrid(k,2)==(Vgrid(i,2)+deltay)

v2=Vgrid(k,3);

break;

end

end

for k=1:length(Vgrid(:,1)) %i,j-1

if Vgrid(k,1)==Vgrid(i,1) && Vgrid(k,2)==(Vgrid(i,2)-deltay)

v3=Vgrid(k,3);

break;

end

end

v4=Vgrid(i+1,3); %i+1,j

v5=Vgrid(i-1,3); %i-1,j

%calc u vel avg

for k=1:length(Ugrid(:,1))

if Ugrid(k,1)==Vgrid(i,1)-(deltax/2) && Ugrid(k,2)==Vgrid(i,2)+(deltay/2)

uforavg1=Ugrid(k,3); %i,j

elseif Ugrid(k,1)==Vgrid(i,1)+(deltax/2) && Ugrid(k,2)==Vgrid(i,2)+(deltay/2)

uforavg2=Ugrid(k,3); %i+1,j

elseif Ugrid(k,1)==Vgrid(i,1)+(deltax/2) && Ugrid(k,2)==Vgrid(i,2)-(deltay/2)

uforavg3=Ugrid(k,3); %i+1,j-1

elseif Ugrid(k,1)==Vgrid(i,1)-(deltax/2) && Ugrid(k,2)==Vgrid(i,2)-(deltay/2)

uforavg4=Ugrid(k,3); %i-1,j-1

end

end

%Calc Pressures

for k=1:length(Pgrid(:,1))

if Pgrid(k,1)==Vgrid(i,1) && Pgrid(k,2)==(Vgrid(i,2)-(deltay/2))

P1vp=Pgrid(k,3); %P i,j

end

if Pgrid(k,1)==Vgrid(i,1) && Pgrid(k,2)==(Vgrid(i,2)+(deltay/2))

P2vp=Pgrid(k,3); %P i,j+1

end

end

end

delcV=rho*( (v1*((v2-v3)/(2*b))) + (((uforavg1+uforavg2+uforavg3+uforavg4)/4)*((v4-v5)/(2*a))) );

deldV=mu*( ((v4-(2*v1)+v5)/(a^2)) + ((v2-(2*v1)+v3)/(b^2)) );

vtild(i,3)=v1 + ((deltaT/rho)* (deldV-delcV + ((P1vp-P2vp)/b)));

end

%Calculate P' field

for i=1:length(Pgrid(:,1))

px1=Pgrid(i,1)+(deltax/2);

px2=Pgrid(i,1)-(deltax/2);

py1=Pgrid(i,2)+(deltay/2);

py2=Pgrid(i,2)-(deltay/2);

for k=1:length(utild(:,1))

if utild(k,1)==px1 && utild(k,2)==Pgrid(i,2)

up=utild(k,3);

end

if utild(k,1)==px2 && utild(k,2)==Pgrid(i,2)

uw=utild(k,3);

end

end

for k=1:length(vtild(:,1))

if vtild(k,1)==Pgrid(i,1) && vtild(k,2)==py1

vp=vtild(k,3);

end

if vtild(k,1)==Pgrid(i,1) && vtild(k,2)==py2

vs=vtild(k,3);

end

end

Pprime(i,3)= -(( ((up-uw)/a) +((vp-vs)/b) )/( ((2*deltaT)/rho)*((1/(a^2)) + (1/(b^2))) ));

end

%Updating utild and vtild to n+1 temporarily to check for-

%-convergence

for i=1:length(Ugrid(:,1))

if Ugrid(i,1)==(-(deltax/2)) || Ugrid(i,1)==((nx*deltax)-(deltax/2)) %Left Boundary and right boundary

continue

end

for k=1:length(Pprime(:,1))

if Pprime(k,1)==Ugrid(i,1)-((deltax)/2) && Pprime(k,2)==Ugrid(i,2)

Ugridtemp(i,3)=utild(i,3) - ( (deltaT/(rho*a))*(Pprime(k+1,3)-Pprime(k,3)));

end

end

end

%For right boundary we need special treatment for uvel:

for i=1:length(Ugrid(:,1))

if Ugrid(i,1)==((nx*deltax)-(deltax/2))

Ugridtemp(i,3)=Ugridtemp(i-1,3);

end

end

for i=1:length(Vgrid(:,1))

if Vgrid(i,1)==0 || Vgrid(i,1)==((nx-1)*deltax) || Vgrid(i,2)==(0-(deltay/2)) || Vgrid(i,2)==((ny*deltay)-(deltay/2)) %left, right, bottom and top boundaries

continue

end

for k=1:length(Pprime(:,1))

if Pprime(k,1)==Vgrid(i,1) && Pprime(k,2)==(Vgrid(i,2)-(deltay/2))

paa=Pprime(k,3);

elseif Pprime(k,1)==Vgrid(i,1) && Pprime(k,2)==(Vgrid(i,2)+(deltay/2))

pbb=Pprime(k,3);

end

end

Vgridtemp(i,3)=vtild(i,3) - ( (deltaT/(rho*b))*(pbb-paa));

end

%For right and top boundary we need special treatment for vvel:

for i=1:length(Vgrid(:,1))

if Vgrid(i,1)==((nx-1)*deltax) %right

Vgridtemp(i,3)=Vgrid(i-1,3);

elseif Vgrid(i,1)==((ny*deltay)-(deltay/2)) %top

for k=1:length(Vgrid(:,1))

if Vgrid(k,1)==Vgrid(i,1) && Vgrid(k,2)==(Vgrid(i,2)-deltay)

Vgridtemp(i,3)=Vgrid(k,3);

end

end

end

end

%Pressure correction

for i=1:length(Pgrid(:,1))

if Pgrid(i,1)==((nx-1)*deltax)

continue

else

Pgrid(i,3)=Pgrid(i,3)+Pprime(i,3);

end

end

%To ensure dp/dy at boundary is zero:

%Check for convergence using temporary values:

%[Convergence Check: L1 norm]

sum=0;

for i=1:length(Pgrid(:,1))

px1=Pgrid(i,1)+(deltax/2);

px2=Pgrid(i,1)-(deltax/2);

py1=Pgrid(i,2)+(deltay/2);

py2=Pgrid(i,2)-(deltay/2);

for k=1:length(Ugridtemp(:,1))

if Ugridtemp(k,1)==px1 && Ugridtemp(k,2)==Pgrid(i,2)

up=Ugridtemp(k,3);

end

if Ugridtemp(k,1)==px2 && Ugridtemp(k,2)==Pgrid(i,2)

uw=Ugridtemp(k,3);

end

end

for k=1:length(Vgridtemp(:,1))

if Vgridtemp(k,1)==Pgrid(i,1) && Vgridtemp(k,2)==py1

vp=Vgridtemp(k,3);

end

if Vgridtemp(k,1)==Pgrid(i,1) && Vgridtemp(k,2)==py2

vs=Vgridtemp(k,3);

end

end

residual=( ((up-uw)/a) +((vp-vs)/b) );

sum=sum+abs(residual);

end

m=sum/length(Pgrid(:,1));

disp(m);

%disp(iter);

if m<=0.001 %Converged

flag=1;

Vgrid=Vgridtemp; %Proceeding to n+1

Ugrid=Ugridtemp;

else %Not converged

iter=iter+1;

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Proceeding to next time step n -> n+1:

%The points in the u and v grids are first updated.

%While updation, points such as the inlet grid points for uvel are-

%-not touched because they are constant vel points and is a -

%-dirichlet bc.

%Similarly the Pgrid is also updated and the exit P points are not-

%-touched because it is a constant pressure Dirichlet BC

%Save files at n+1 for the current time-step

if mod((timestep*deltaT),0.01)==0

file_name = sprintf('U_vel_timestep_%d.txt',timestep);

id=fopen([savepath file_name],'w');

for i=1:length(Ugrid(:,1))

fprintf(id,'%f %f %f\n',Ugrid(i,1),Ugrid(i,2),Ugrid(i,3));

end

fclose(id);

file_name = sprintf('V_vel_timestep_%d.txt',timestep);

id=fopen([savepath file_name],'w');

for i=1:length(Vgrid(:,1))

fprintf(id,'%f %f %f\n',Vgrid(i,1),Vgrid(i,2),Vgrid(i,3));

end

fclose(id);

file_name = sprintf('P_timestep_%d.txt',timestep);

id=fopen([savepath file_name],'w');

for i=1:length(Pgrid(:,1))

fprintf(id,'%f %f %f\n',Pgrid(i,1),Pgrid(i,2),Pgrid(i,3));

end

fclose(id);

%Saving complete and proceed to next time-step

end

end

function [Pgrid,Ugrid,Vgrid] = gridpointgenerator(nx,ny,deltax,deltay)

% Grid point generator

%Row by row, upward

%Pressure collocated grid points

Pgridx(1)=0;

Pgridy(1)=0;

for i=2:nx

Pgridx(i)=Pgridx(i-1)+deltax;

end

for i=2:ny

Pgridy(i)=Pgridy(i-1)+deltay;

end

count=1;

for i=1:ny

for j=1:nx

Pgridxp(count)=Pgridx(j);

Pgridyp(count)=Pgridy(i);

count=count+1;

end

end

%x vel collocated grid points

nx2=nx+1;

ugridx(1)=-(deltax/2);

ugridy(1)=0;

for i=2:nx2

ugridx(i)=ugridx(i-1)+deltax;

end

for i=2:ny

ugridy(i)=ugridy(i-1)+deltay;

end

count=1;

for i=1:ny

for j=1:nx2

ugridxp(count)=ugridx(j);

ugridyp(count)=ugridy(i);

count=count+1;

end

end

%y vel collocated grid points

ny2=ny+1;

vgridy(1)=-(deltay/2);

vgridx(1)=0;

for i=2:nx

vgridx(i)=vgridx(i-1)+deltax;

end

for i=2:ny2

vgridy(i)=vgridy(i-1)+deltay;

end

count=1;

for i=1:ny2

for j=1:nx

vgridxp(count)=vgridx(j);

vgridyp(count)=vgridy(i);

count=count+1;

end

end

Pgridxp=Pgridxp';

Pgridyp=Pgridyp';

ugridxp=ugridxp';

ugridyp=ugridyp';

vgridxp=vgridxp';

vgridyp=vgridyp';

figure

plot(Pgridxp,Pgridyp,'o','MarkerSize',5,'Markeredgecolor','k','markerfacecolor','k');

hold on

plot(ugridxp,ugridyp,'x','MarkerSize',5,'Markeredgecolor','b','markerfacecolor','b');

hold on

plot(vgridxp,vgridyp,'*','MarkerSize',5,'Markeredgecolor','g','markerfacecolor','g');

title('FDM Grid Points')

legend('Pressure Points','u vel points','v vel points');

Pgrid=[Pgridxp,Pgridyp];

Ugrid=[ugridxp,ugridyp];

Vgrid=[vgridxp,vgridyp];

end

[Steven Lord]

If you're trying to ask how to convert that large section of MATLAB code to C code, use MATLAB Coder as @Raghu Boggavarapu said in their answer. If while doing so you encounter difficulties, post as a new question the specific error and/or warning messages you receive and ask for help in that new question.

Answer 3

[Raghu Boggavarapu]

Use MATLAB Coder to convert your MATLAB Code into C/C++

[Harshal Suresh]

Does it come free with the standard student license?

Answer 4

[Israa Alhuniti]

function [Top_predator_fit,Top_predator_pos,Convergence_curve]=MPA(SearchAgents_no,Max_iter,lb,ub,dim,fobj)

Top_predator_pos=zeros(1,dim);

Top_predator_fit=inf;

Convergence_curve=zeros(1,Max_iter);

stepsize=zeros(SearchAgents_no,dim);

fitness=inf(SearchAgents_no,1);

Prey=initialization(SearchAgents_no,dim,ub,lb);

Xmin=repmat(ones(1,dim).*lb,SearchAgents_no,1);

Xmax=repmat(ones(1,dim).*ub,SearchAgents_no,1);

Iter=0;

FADs=0.2;

P=0.5;

while Iter<Max_iter

%------------------- Detecting top predator -----------------

for i=1:size(Prey,1)

Flag4ub=Prey(i,:)>ub;

Flag4lb=Prey(i,:)<lb;

Prey(i,:)=(Prey(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;

fitness(i,1)=fobj(Prey(i,:));

if fitness(i,1)<Top_predator_fit

Top_predator_fit=fitness(i,1);

Top_predator_pos=Prey(i,:);

end

end

%------------------- Marine Memory saving -------------------

if Iter==0

fit_old=fitness; Prey_old=Prey;

end

Inx=(fit_old<fitness);

Indx=repmat(Inx,1,dim);

Prey=Indx.*Prey_old+~Indx.*Prey;

fitness=Inx.*fit_old+~Inx.*fitness;

fit_old=fitness; Prey_old=Prey;

%------------------------------------------------------------

Elite=repmat(Top_predator_pos,SearchAgents_no,1); %(Eq. 10)

CF=(1-Iter/Max_iter)^(2*Iter/Max_iter);

RL=0.05*levy(SearchAgents_no,dim,1.5); %Levy random number vector

RB=randn(SearchAgents_no,dim); %Brownian random number vector

for i=1:size(Prey,1)

for j=1:size(Prey,2)

R=rand();

%------------------ Phase 1 (Eq.12) -------------------

if Iter<Max_iter/3

stepsize(i,j)=RB(i,j)*(Elite(i,j)-RB(i,j)*Prey(i,j));

Prey(i,j)=Prey(i,j)+P*R*stepsize(i,j);

%--------------- Phase 2 (Eqs. 13 & 14)----------------

elseif Iter>Max_iter/3 && Iter<2*Max_iter/3

if i>size(Prey,1)/2

stepsize(i,j)=RB(i,j)*(RB(i,j)*Elite(i,j)-Prey(i,j));

Prey(i,j)=Elite(i,j)+P*CF*stepsize(i,j);

else

stepsize(i,j)=RL(i,j)*(Elite(i,j)-RL(i,j)*Prey(i,j));

Prey(i,j)=Prey(i,j)+P*R*stepsize(i,j);

end

else

stepsize(i,j)=RL(i,j)*(RL(i,j)*Elite(i,j)-Prey(i,j));

Prey(i,j)=Elite(i,j)+P*CF*stepsize(i,j);

end

end

end

for i=1:size(Prey,1)

Flag4ub=Prey(i,:)>ub;

Flag4lb=Prey(i,:)<lb;

Prey(i,:)=(Prey(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;

fitness(i,1)=fobj(Prey(i,:));

if fitness(i,1)<Top_predator_fit

Top_predator_fit=fitness(i,1);

Top_predator_pos=Prey(i,:);

end

end

if Iter==0

fit_old=fitness; Prey_old=Prey;

end

Inx=(fit_old<fitness);

Indx=repmat(Inx,1,dim);

Prey=Indx.*Prey_old+~Indx.*Prey;

fitness=Inx.*fit_old+~Inx.*fitness;

fit_old=fitness; Prey_old=Prey;

if rand()<FADs

U=rand(SearchAgents_no,dim)<FADs;

Prey=Prey+CF*((Xmin+rand(SearchAgents_no,dim).*(Xmax-Xmin)).*U);

else

r=rand(); Rs=size(Prey,1);

stepsize=(FADs*(1-r)+r)*(Prey(randperm(Rs),:)-Prey(randperm(Rs),:));

Prey=Prey+stepsize;

end

Iter=Iter+1;

Convergence_curve(Iter)=Top_predator_fit;

end