How to convert MATLAB code to C++ code?

260 次查看（过去 30 天）
goutham gumm
goutham gumm 2018-4-20
评论： Harshal Suresh 2023-12-18
How do I convert my MATLAB code to a C++ code?
  1 个评论
Ankit Singh Jaswal
Ankit Singh Jaswal 2023-9-25
%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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Ameer Hamza
Ameer Hamza 2018-4-20

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 个评论
显示 2更早的评论
PatrizioGraziosi
PatrizioGraziosi 2020-10-1
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 2020-10-1
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.

请先登录，再进行评论。

更多回答（3 个）

Raghu Boggavarapu
Raghu Boggavarapu 2022-1-19
Use MATLAB Coder to convert your MATLAB Code into C/C++
Ahmad Nadeem
Ahmad Nadeem 2021-7-5
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 个评论
显示 1更早的评论
Ankit Singh Jaswal
Ankit Singh Jaswal 2023-9-25
%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 2023-9-25
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.

请先登录，再进行评论。

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

类别

在 Help Center 和 File Exchange 中查找有关 Software Development Tools 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by