How to convert MATLAB code to C++ code?

500 views (last 30 days)
How do I convert my MATLAB code to a C++ code?
  1 Comment
Ankit Singh Jaswal
Ankit Singh Jaswal on 25 Sep 2023
%MAC algorithm based code developed by Rohit Kanchi
%Contact: +91 91081 52106
%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

Sign in to comment.

Accepted Answer

Ameer Hamza
Ameer Hamza on 20 Apr 2018

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.

  4 Comments
PatrizioGraziosi
PatrizioGraziosi on 1 Oct 2020
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
Ameer Hamza on 1 Oct 2020
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.

Sign in to comment.

More Answers (3)

Ahmad Nadeem
Ahmad Nadeem on 5 Jul 2021
I want to Change MATLAB code to C++ can anyone tell me how can I convert this? I am working on phase field simulation
  3 Comments
Ankit Singh Jaswal
Ankit Singh Jaswal on 25 Sep 2023
%MAC algorithm based code developed by Rohit Kanchi
%Contact: +91 91081 52106
%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
Steven Lord on 25 Sep 2023
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.

Sign in to comment.


Raghu Boggavarapu
Raghu Boggavarapu on 19 Jan 2022
Use MATLAB Coder to convert your MATLAB Code into C/C++

Israa Alhuniti
Israa Alhuniti on 16 Mar 2022
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

Categories

Find more on Introduction to Installation and Licensing in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!