Reduce computing time ode system

1 次查看（过去 30 天）
EldaEbrithil
EldaEbrithil 2020-6-21
编辑： EldaEbrithil 2020-6-22
Hi all
i would like to reduce the computing time of my code in such a way as to increase the dimensions of integration domain (L), hradialchmaber vetor, and M01 vector
hradialchamber=(1.5:0.5:2);%mm
for i=1:length(hradialchamber)
Hradialchamber(i)=hradialchamber(i)/1000;%m
end
rextbobbin=9.525:0.3:17.5;%
Rextbobbin=rextbobbin/1000;%m
Atomicweight=131.29;%
gamma=1.667;%
R=8314/Atomicweight;%
Cp=(gamma/(gamma-1))*R;%e
dconduttore=(0.43);%mm
Dconduttore=dconduttore/1000;%m
throatRadius=0.100:0.0005:0.5;%mm
throatradius=throatRadius/1000;
exitRadius=0.600:0.001:3;%mm
exitradius=exitRadius/1000;
exitarea=pi*exitradius.^2;%m^2
Mexit=5.88;%INPUT
Tcoldgas=300;%K
for i=1:length(throatRadius)
for j=1:length(exitRadius)
radiusparameter=0.5;
Mcheck(i,j)=((0.50*exitradius(j)/(throatradius(i)))^2)-...
((1/Mexit)*((2/(gamma+1))*...
(1+(((gamma-1)/2)*Mexit^2)))^((gamma+1)/(2*(gamma-1))));%
diffM(i,j)=Mcheck(i,j)-Mexit;
end
end
DiffminM2=min(diffM(diffM>0));
diffminM2=max(diffM(diffM<0));%
if DiffminM2-abs(diffminM2)>0%
[righ2,colh2]=find(diffM==diffminM2);
for i=1:length(righ2)
Mcheck1(i)=Mcheck(righ2(i),colh2(i));
end
exitradius1=exitradius(colh2);
throatradius1=throatradius(righ2);
else DiffminM2-abs(diffminM2)<0;
[RigH,ColH]=find(diffM==DiffminM2);
for i=1:length(RigH)
Mcheck1(i)=Mcheck(RigH(i),ColH(i));
end
exitradius1=exitradius(ColH);
throatradius1=throatradius(RigH);
end
exitradiusfinal=max(exitradius1);
throatradiusfinal=max(throatradius1);
throatarea=pi*throatradiusfinal^2;%m^2
exitareafinal=pi*exitradiusfinal^2;%m^2
pressure=2.2;%bar
Pressure=pressure*10^5;%Pa
theater=803;%K
emittancetantalum=0.16;
dexttantalum=1.5;%mm
Dexttantalum=dexttantalum/1000;%m
for i=1:length(hradialchamber)
Aeffettiva(i)=(Dexttantalum*(Hradialchamber(i)))-(pi*(Dexttantalum^2)/4);
A2(i)=Aeffettiva(i);
M2=0.001:0.0001:0.8;
for j=1:length(M2)
M2check(j,i)=(A2(i)/(throatarea))-...
((1/M2(j))*((2/(gamma+1))*...
(1+(((gamma-1)/2)*M2(j)^2)))^((gamma+1)/(2*(gamma-1))));
end
end
Diff_minM2=min(M2check(M2check>0));
diff_minM2=max(M2check(M2check<0));
if Diff_minM2-abs(diff_minM2)>0%
[righ_2,colh_2]=find(M2check==diff_minM2);
M2final=M2(righ_2);
else Diff_minM2-abs(diff_minM2)<0
[RigH_2,ColH_2]=find(M2check==Diff_minM2);
M2final=M2(RigH_2);
end
L=0.2:0.1:0.3;
M01=(0.1:0.1:M2final);
dynamicviscosity=54.89*10^-6;
dynamicviscositycold=23.23*10^-6;
dynamicviscosityaverage=(dynamicviscosity+dynamicviscositycold)/2;
thermalconductivityxenonhot=13.07*10^-3;
thermalconductivityxenoncold=5.555*10^-3;
thermalconductivityxenonaverage=(thermalconductivityxenonhot+...
thermalconductivityxenoncold)/2;
Pr=dynamicviscosityaverage*Cp/thermalconductivityxenonaverage;
ro1=Pressure/(R*Tcoldgas);
for y=1:length(hradialchamber)
Aeffettiva(y)=(Dexttantalum*(Hradialchamber(y)))-(pi*(Dexttantalum^2)/4);
Lq(y)=(Aeffettiva(y)^0.5);
Perimetro(y)=Hradialchamber(y)*2+Dexttantalum*2+(2*pi*(Dexttantalum/2));
Dh(y)=(4*(Aeffettiva(y))/Perimetro(y));
sqrtAeffettiva(y)=Aeffettiva(y)^0.5;
for k=1:length(M01)
u1(k)=M01(k)*(gamma*R*Tcoldgas)^0.5;
mdotstart(k,y)=ro1*u1(k)*Aeffettiva(y);
Re(k,y)=(4*mdotstart(k,y))/(dynamicviscosityaverage*pi*Dh(y));
Deannumber(k,y)=Re(k,y)*(thermalconductivityxenonaverage)^0.5;
for i=1:length(rextbobbin)
delta(i)= Dexttantalum/(Rextbobbin(i));
Recritic(i)=2100*(1+12*(delta(i)^0.5));
Rediffer(i,k,y)=Re(k,y)-Recritic(i);
%%%%From here to end the critical part%%%
if Re(k,y)>Recritic(i)
RE(k,y)=Re(k,y);
[RErig,REcol]=find(RE>0);
A(y)=Aeffettiva(y);
Mo1(k)=M01(k);
NuT(k,y)=Re(k,y)*(Pr^(1/3))*(1/(2*(3.6*log10(6.115/Re(k,y)))^2));
fT(k,y)=0.0767/(Re(k,y)^0.25);
hT(k,y)=NuT(k,y)*thermalconductivityxenonaverage/(sqrtAeffettiva(y));
ChT(k,y)=NuT(k,y)/(Re(k,y)*Pr);
end
end
end
end
Nc=0.001;
P0=2.200000e+05;
T0=300;
Tt0=T0;
Pt0=P0;
for y=1:length(hradialchamber)
Aeffettiva_cell(y)={Aeffettiva(y)};
Perimetro_cell(y)={Perimetro(y)};
for j=1:length(M01)
Y0(j)={[Tt0,Pt0,Mo1(j)]};
ChT_cell(j,y)={ChT(j,y)};
ft_cell(j,y)={fT(j,y)};
end
end
for i=1:length(L)
Dall(i) ={[0:Nc:L(i)]};
end
for iDom = 1:numel(Dall)
xRange = Dall{iDom};
for iInitial = 1:numel(Y0)
for iAeff=1:numel(Aeffettiva_cell)
Ch=ChT_cell{iInitial,iAeff};%it doesn't depend by iDom
Fc=ft_cell{iInitial,iAeff};
Aeff=Aeffettiva_cell{iAeff};
Perim=Perimetro_cell{iAeff};
[xSol{iDom,iInitial,iAeff},YSol{iDom,iInitial,iAeff}]=ode23(@(x,Y) ...
chambsinglebobb(x,Y,Ch,Aeff,Perim,Fc) ,xRange,Y0{iInitial});
end
end
end
function dYdx=chambsinglebobb(x,Y,Ch,Aeff,Perim,Fc)
gamma=1.667;
Dexttantalum=0.001500000000000;
Tw=803;
Tt=Y(1);
Pt=Y(2);
M=Y(3);
dTtdx=Ch*(Tw-Tt)*(Perim/Aeff);
dPtdx=Pt*((-gamma*((M^2)/2)*Fc*(Perim/Aeff))-...
(gamma*((M^2)/2)*dTtdx*(1/Tt)));
dMdx=M*(((1+((gamma-1)/2)*M^2)/(1-M^2))*((gamma*(M^2)*Fc*...
Perim)/(2*Aeff))+...
((0.5*((gamma*M^2)+1))*(1+((gamma-1)/2)*M^2)/(1-M^2))*...
(Ch*(Tw-Tt)*Perim)/(Aeff*Tt));
dYdx=[dTtdx;dPtdx;dMdx];
end
In the order, the operations of the code are the calculation of:
  • throatradiusfinal
  • exitradiusfinal
  • throatarea
  • M2final
  • NuT
  • ChT
  • fT
  • YSol
The variable parameters are: hardialchamber, rextbobbin, M01, L. The problem related to high computational time is connected (i think) with the the calculation of the ode system solution. In particular i have written the code in such a way the ODE system is solved numerous time namely for different IC (showbelow), different value of Ch in the differential equations, different value of Fc in the differential equations and of course different integration domain (L):
Y0(j)={[Tt0,Pt0,Mo1(j)]}
Ch=ChT_cell{iInitial,iAeff};
Fc=ft_cell{iInitial,iAeff};
xRange = Dall{iDom};
The problem (i think) is connected with ChT that it is physically independent of L, this leads to increase very much the number of combinations allowed. I have written the same code but with the only difference of ChT dependent on L, this code is much faster than the privious one BUT of course it is not physically acceptable. So i ask you if you know some shrewdness for increasing the calculation time for this problem.
Regards

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

回答（2 个）

Alan Stevens
Alan Stevens 2020-6-21
Here's a start, looking at your triply nested y, k and i loops.
% Remove expressions that don't depend on y, k or i from the loops
dynamicviscosity=54.89*10^-6;
dynamicviscositycold=23.23*10^-6;
dynamicviscosityaverage=(dynamicviscosity+dynamicviscositycold)/2;
thermalconductivityxenonhot=13.07*10^-3;
thermalconductivityxenoncold=5.555*10^-3;
thermalconductivityxenonaverage=(thermalconductivityxenonhot+...
thermalconductivityxenoncold)/2;
Pr=dynamicviscosityaverage*Cp/thermalconductivityxenonaverage;
ro1=Pressure/(R*Tcoldgas);
for y=1:length(hradialchamber)
% Put expressions that depend only on y here
Aeffettiva(y)=(Dexttantalum*(Hradialchamber(y)))-(pi*(Dexttantalum^2)/4);
Lq(y)=(Aeffettiva(y)^0.5);
Perimetro(y)=Hradialchamber(y)*2+Dexttantalum*2+(2*pi*(Dexttantalum/2));
Dh(y)=(4*(Aeffettiva(y))/Perimetro(y));
sqrtAeffettiva(y)=Aeffettiva(y)^0.5;
for k=1:length(M01)
% Put expressions that depend on k or k and y here
u1(k)=M01(k)*(gamma*R*Tcoldgas)^0.5;
mdotstart(k,y)=ro1*u1(k)*Aeffettiva(y);
Re(k,y)=(4*mdotstart(k,y))/(dynamicviscosityaverage*pi*Dh(y));
Deannumber(k,y)=Re(k,y)*(thermalconductivityxenonaverage)^0.5;
for i=1:length(rextbobbin)
% Put expressions that depend on i, i and y, i and k or i, y and k here
delta(i)= Dexttantalum/(Rextbobbin(i));
Recritic(i)=2100*(1+12*(delta(i)^0.5));
Rediffer(i,k,y)=Re(k,y)-Recritic(i);
%%%%From here to end the critical part%%%
if Re(k,y)>Recritic(i)
RE(k,y)=Re(k,y);
[RErig,REcol]=find(RE>0);
A(y)=Aeffettiva(y);
Mo1(k)=M01(k);
NuT(k,y)=Re(k,y)*(Pr^(1/3))*(1/(2*(3.6*log10(6.115/Re(k,y)))^2));
fT(k,y)=0.0767/(Re(k,y)^0.25);
hT(k,y)=NuT(k,y)*thermalconductivityxenonaverage/(sqrtAeffettiva(y));
ChT(k,y)=NuT(k,y)/(Re(k,y)*Pr);
end
end
end
end
  1 个评论
EldaEbrithil
EldaEbrithil 2020-6-21
编辑：EldaEbrithil 2020-6-21
Thank you for reply me
ok i have done it, but the computational time is still too high
if you have any question please let me know

请先登录，再进行评论。

Alan Stevens
Alan Stevens 2020-6-21
Look at the other triple loop. Similar improvements can be made:
for y=1:length(hradialchamber)
Aeffettiva_cell(y)={Aeffettiva(y)};
Perimetro_cell(y)={Perimetro(y)};
for j=1:length(M01)
Y0(j)={[Tt0,Pt0,Mo1(j)]};
ChT_cell(j,y)={ChT(j,y)};
ft_cell(j,y)={fT(j,y)};
end
end
for i=1:length(L)
Dall(i) ={[0:Nc:L(i)]};
end
  6 个评论
显示 4更早的评论
Alan Stevens
Alan Stevens 2020-6-22
You should pre-allocate space for some of your larger matrices. However, this will only help a little. Look at the steps you require for your throat radii. Do you really need that number? If so, I'm afraid you need to put up with long running times! You could try compiling the code, but I suspect run-times will still be long.
EldaEbrithil
EldaEbrithil 2020-6-22
编辑：EldaEbrithil 2020-6-22
ok i have found the problem: the ode23 solve the system which gives some values of
YSol{iDom,iInitial}(iDom,3)
bigger or equal to 1. This is not allowed by the equations of the system, this create a lot of problems in the solving process and therefore the long computational times are explained. Is there a way to prevent ode23 to solve the system as soon as
YSol{iDom,iInitial}(iDom,3)
values are near 1 or maybe bigger than one?
something like that but i don't know:
if YSol{i,j}(i,3)>1
break;
end

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by