Reduce computing time ode system

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

Answer 1

[Alan Stevens]

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

[EldaEbrithil]

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

Answer 2

[Alan Stevens]

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

[EldaEbrithil]

Done

always too much time

[Alan Stevens]

You call function chambsinglebobb(...), but your listing above only shows function tubosinglebobb(...).

[EldaEbrithil]

sorry, wrong function name...

now it is correct

[EldaEbrithil]

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

I'll try to explain a little more how the above code works.

at each iteration the domain changes increasing size. At each iteration the triad of initial conditions change, and together with them change Ch and Fc. So ode23 solves my system for different domains and with a different set of initial conditions /Ch / Fc

[Alan Stevens]

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]

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