Error in assignment in parfor

5 views (last 30 days)
Flavien Saenz Molina
Flavien Saenz Molina on 27 Nov 2015
Commented: Walter Roberson on 27 Nov 2015
Hello everyone !
I am struggling with a piece of code.
x = sym('x', [1,H]);
eqm=gamma*p_isa.*x.^2==(WoS/C_D0)*[(ToW_sl*0.76*(0.907+0.262*(abs(x-0.5)).^1.5).*(rho_isa/rho_ssl).^0.7)-sqrt((ToW_sl*0.76*(0.907+0.262*(abs(x-0.5)).^1.5).*(rho_isa/rho_ssl).^0.7).^2-4*K*C_D0)];
eqM=gamma*p_isa.*x.^2==(WoS/C_D0)*[(ToW_sl*0.76*(0.907+0.262*(abs(x-0.5)).^1.5).*(rho_isa/rho_ssl).^0.7)+sqrt((ToW_sl*0.76*(0.907+0.262*(abs(x-0.5)).^1.5).*(rho_isa/rho_ssl).^0.7).^2-4*K*C_D0)];
parfor i=1:H
M_rmin(i)=double(vpasolve(eqm(i),x(i),M_min(i)));
M_rmax(i)=double(vpasolve(eqM(i),x(i),M_max(i)));
V_rmax(i)=M_rmax(i)*c(i);
V_rmin(i)=M_rmin(i)*c(i);
end
I defined 2 vector equations, each row corresponding to a given altitude. All the coefficients have been computed before and the only thing I change is the altitude step, hence H which is the length of my altitude vector. The code works fine for high steps (100m is ok) but i need to make it work with a smaller step, so I tried 10m. The thing is I get this error at the line of parfor : In an assignment A(:) = B, the number of elements in A and B must be the same.
The loop stopped at some point. I checked all the dimensions, they are all matching. I solved the equation outside the loop for a rank higher than the one at which it had stopped and I get a result, so this should work.
Could you help me ?
Thanks !

Answers (2)

Walter Roberson
Walter Roberson on 27 Nov 2015
You did not happen to provide any sample values, and it is not clear from your posting which variable is the "step".
Generally speaking, you could run into the problem you did if vpasolve() were to return something other than 1 result. That might occur if it returned no results (no solution found); for equations that it determines to be polynomials it can also return multiple (more than 1) value.

Flavien Saenz Molina
Flavien Saenz Molina on 27 Nov 2015
Edited: Flavien Saenz Molina on 27 Nov 2015
Thank you for your answer !
I thought about that, that's why I did check the next iteration manually and I found out that there is only one solution to it ! By step, I meant iteration, sorry for the misunderstanding. What I don't get is why my code works for an altitude step of 50m (like in the one I gave you here) but not for 10m or 20m ... By the way, I can't figure out why but the for loop exits with an error always for h=10000.
For instance, if you run my code below with an altitude step of 50, it works but if you use 20 instead, it doesn't, there is an error in the parfor at step 501, which exactly corresponds to h=10000. What I don't understand is now if outside the loop, I solve the equation for i=502, I have only one result, so there should not be any dimension issue ! Besides, if I change the step to 10, i have an error in the parfor after another number of iterations which also corresponds to h=10000.
Here is the all code I am using :
clear
format long g
%%%%%%%1 to compute %%%%%%%
plot_V=0; %V_max, V_min, V_stall
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% INDATA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C_D0 = 0.022;
K = 0.064;
WoS = 800; %N/m^2
ToW_sl = 0.30;
C_Lmax = 1.50;
rho_ssl = 1.225; %kg/m^3
gamma=1.40;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ISA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[h, z, T_isa, p_isa, rho_isa, c, mu]=isatm(50,15000);
s=size(h);
H=s(2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simple thrust model %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
parfor i=1:H
ToW_alt(i)=(rho_isa(i)/rho_ssl)*ToW_sl;
V_min(i)=sqrt(WoS/(rho_isa(i)*C_D0)*(ToW_alt(i)-sqrt(ToW_alt(i)^2-4*K*C_D0)));
V_max(i)=sqrt(WoS/(rho_isa(i)*C_D0)*(ToW_alt(i)+sqrt(ToW_alt(i)^2-4*K*C_D0)));
V_stall(i)=sqrt(2*WoS/(rho_isa(i)*C_Lmax));
M_min(i)=V_min(i)/c(i);
M_max(i)=V_max(i)/c(i);
M_stall(i)=V_stall(i)/c(i);
end
L=find(imag(V_max)~=0);
i_abs=min(L)-1;
h_abs=h(i_abs)
[int,i_int]=min(abs(V_min-V_stall));
h_int=h(i_int)
if plot_V==1
figure
subplot(1,1,1)
plot(V_max(1:i_abs),h(1:i_abs),'r')
hold on
plot(V_min(1:i_abs),h(1:i_abs),'b')
hold on
plot(V_stall(1:i_abs),h(1:i_abs),'g')
legend('V_{max}', 'V_{min}', 'V_{stall}')
xlabel('V (m/s)')
ylabel('h (m)')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Refined thrust model %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = sym('x', [1,H])
eqm=gamma*p_isa.*x.^2==(WoS/C_D0)*[(ToW_sl*0.76*(0.907+0.262*(abs(x-0.5)).^1.5).*(rho_isa/rho_ssl).^0.7)-sqrt((ToW_sl*0.76*(0.907+0.262*(abs(x-0.5)).^1.5).*(rho_isa/rho_ssl).^0.7).^2-4*K*C_D0)];
eqM=gamma*p_isa.*x.^2==(WoS/C_D0)*[(ToW_sl*0.76*(0.907+0.262*(abs(x-0.5)).^1.5).*(rho_isa/rho_ssl).^0.7)+sqrt((ToW_sl*0.76*(0.907+0.262*(abs(x-0.5)).^1.5).*(rho_isa/rho_ssl).^0.7).^2-4*K*C_D0)];
parfor i=1:H
M_rmin(i)=double(vpasolve(eqm(i),x(i),M_min(i)));
M_rmax(i)=double(vpasolve(eqM(i),x(i),M_max(i)));
V_rmax(i)=M_rmax(i)*c(i);
V_rmin(i)=M_rmin(i)*c(i);
end
and the isatm function :
function [h, z, T_isa, p_isa, rho_isa, c, mu]=isatm(p,h_max)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ISA constants %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g0 = 9.80665; %m/s^2
p_ssl = 101.325*1E3; %Pa
T_ssl = 288.15; %K
rho_ssl = 1.225; %kg/m^3
r = 6.356766E6; %m
R=p_ssl/(rho_ssl*T_ssl); %J/(kg*K)
gamma=1.40;
mu0=1.7894E-2; %g/(m*s)
T0=288.15; %K
S=110; %K
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Altitude %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=0:p:h_max;
s=size(h);
H=s(2);
for i=1:H
z(i)=r/(r-h(i))*h(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ISA %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if h_max>11000
I=find(h==11000);
% h = 0 ... 11000
a=-6.5E-3;
for i=1:I
T_isa(i)=T_ssl+a*h(i);
p_isa(i)=p_ssl*(T_isa(i)/T_ssl)^(-g0/(a*R));
rho_isa(i)=rho_ssl*(T_isa(i)/T_ssl)^(-((g0/(a*R))+1));
c(i)=sqrt(gamma*R*T_isa(i));
mu(i)=mu0*(T_isa(i)/T0)^(3/2)*((T0+S)/(T_isa(i)+S));
end
if h_max>20000
J=find(h==20000);
% h = 11000 ... 20000
a=0;
for i=I+1:J
T_isa(i)=T_isa(I)+a*(h(i)-h(I));
p_isa(i)=p_isa(I)*exp(-(g0/(R*T_isa(I)))*(h(i)-h(I)));
rho_isa(i)=rho_isa(I)*exp(-(g0/(R*T_isa(I)))*(h(i)-h(I)));
c(i)=sqrt(gamma*R*T_isa(i));
mu(i)=mu0*(T_isa(i)/T0)^(3/2)*((T0+S)/(T_isa(i)+S));
end
% h = 20000 ... h_max
a=1.0E-3;
for i=J+1:H
T_isa(i)=T_isa(J)+a*(h(i)-h(J));
p_isa(i)=p_isa(J)*(T_isa(i)/T_isa(J))^(-g0/(a*R));
rho_isa(i)=rho_isa(J)*(T_isa(i)/T_isa(J))^(-((g0/(a*R))+1));
c(i)=sqrt(gamma*R*T_isa(i));
mu(i)=mu0*(T_isa(i)/T0)^(3/2)*((T0+S)/(T_isa(i)+S));
end
else
% h = 11000 ... h_max
a=0;
for i=I+1:H
T_isa(i)=T_isa(I)+a*(h(i)-h(I));
p_isa(i)=p_isa(I)*exp(-(g0/(R*T_isa(I)))*(h(i)-h(I)));
rho_isa(i)=rho_isa(I)*exp(-(g0/(R*T_isa(I)))*(h(i)-h(I)));
c(i)=sqrt(gamma*R*T_isa(i));
mu(i)=mu0*(T_isa(i)/T0)^(3/2)*((T0+S)/(T_isa(i)+S));
end
end
else
% h = 0 ... h_max
a=-6.5E-3;
for i=1:H
T_isa(i)=T_ssl+a*h(i);
p_isa(i)=p_ssl*(T_isa(i)/T_ssl)^(-g0/(a*R));
rho_isa(i)=rho_ssl*(T_isa(i)/T_ssl)^(-((g0/(a*R))+1));
c(i)=sqrt(gamma*R*T_isa(i));
mu(i)=mu0*(T_isa(i)/T0)^(3/2)*((T0+S)/(T_isa(i)+S));
end
end
end
  1 Comment
Walter Roberson
Walter Roberson on 27 Nov 2015
Do you get the same error if you use "for" instead of "parfor" ?
I recommend,
parfor i=1:H
t_M_rmin = double(vpasolve(eqm(i),x(i),M_min(i)));
t_M_rmax = double(vpasolve(eqM(i),x(i),M_max(i)));
if length(t_M_rmin) ~= 1 || length(t_M) ~= 1
fprintf('iteration %d failed, received %d results for rmin and %d results for rmax\n', i, length(t_M_rmin), length(t_M_rmax));
M_rmin(i) = nan;
M_rmax(i) = nan;
V_rmax(i) = nan;
V_rmin(i) = nan;
else
M_rmin(i) = t_M_rmin;
M_rmax(i) = t_M_rmax;
V_rmax(i) = t_M_rmax * c(i);
V_rmin(i) = t_M_rmin(i) * c(i);
end
end

Sign in to comment.

Categories

Find more on Parallel for-Loops (parfor) in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!