Reduce computing time ode system

2 Ansichten (letzte 30 Tage)
EldaEbrithil
EldaEbrithil am 21 Jun. 2020
Bearbeitet: EldaEbrithil am 22 Jun. 2020
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

Antworten (2)

Alan Stevens
Alan Stevens am 21 Jun. 2020
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 Kommentar
EldaEbrithil
EldaEbrithil am 21 Jun. 2020
Bearbeitet: EldaEbrithil am 21 Jun. 2020
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

Melden Sie sich an, um zu kommentieren.


Alan Stevens
Alan Stevens am 21 Jun. 2020
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 Kommentare
Alan Stevens
Alan Stevens am 22 Jun. 2020
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 am 22 Jun. 2020
Bearbeitet: EldaEbrithil am 22 Jun. 2020
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

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Programming finden Sie in Help Center und File Exchange

Produkte


Version

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by