Error in ode15s: incompatible sizes arrays, while no errors in function
3 views (last 30 days)
Show older comments
Hi everyone,
I have script where I try to calculate the position of multiple particles over time. This is not a problem by itself, but the velocity of the particles is dependent on the sum of the particles that are present at each position. After a while, I managed to make a function that at least runs, but now I have an error in ode15s itself. I get similar errors when using other solvers (23s, 45).
The function file is a bit more than 100 lines, but the most important part is at the bottom, where I have my mass balances. Since I want to calculate everything for multiple particle sizes, I give MatLab a 100*4 matrix. MatLab makes it a 400*1 matrix, which I redistribute to the original 100*4 size for every calculation. At the end, I change it back to a 400*1 matrix.
clear
clc
method = 'distribution';
%% parameters
% time
tstep =100;
tend = tstep*100;
t = 0:tstep:tend;
% feed properties
rhow = 1000; % Density water
rhod = 1200; % Preliminary guess density particles
rhop = 1100; % Preliminary guess density protein
etaL = 0.001; % Preliminary guess viscosity continuous phase
g = 9.81; % Gravity
cmax = 0.6; % Maximal volumetric fraction
cini = 0.2; % Initial volumetric fraction
cprot = 0.55; % Preliminary guess protein fraction in algae
rho_in = cini*rhod + (1-cini)*rhow;
% centrifuge properties
r = 0.076; % Centrifugation radius
z = 0.5; % Centrifugation height (m)
Nr = 100; % Parts over radius
Nz = 100; % Parts over height
Ny = 100; % Parts over perpendicular to radius
dr = r./Nr; % Size of radius part (m)
dz = z./Nz; % Size of height part (m)
Vt= pi*r^2*z; % Total volume (m^3)
dy = zeros(length(linspace(0,360,Ny)),1);
Vcell = zeros(Ny,1);
mprot_ini = Vt*(1-cini)*cprot*rhop;
for i = 1:length(dy) % Calculation for size of parts perpendicular to radius
if i>1
dy(i) = pi*0.5*dr*i/Ny;
Vcell(i) = pi*(i*dr)^2*dz/Ny-pi*((i-1)*dr)^2*dz/Ny;
else
Vcell(i) = pi*dr^2*dz/Ny;
end
end
% run parameters
f = 40000/60; % Rotation frequency (1/s)
gc = (2*pi*f)^2*r; % Centrifugal force (m/s^2)
% energetic parameters
Volt = 230;
Ampere = 15;
% economical values
Cp = 5; % Price of protein (€/kg)
Ce = 0.39; % price of energy (€/kWh)
%% Choose method for model
switch method
case 'one size'
dp = 0.5e-6; % Particle size
c = zeros(Nr,1);
c(:) = cini; % Initial condition
case 'distribution'
dp1 = 0.15*10^-6; % Particle size distribution
dp2 = 0.25*10^-6;
dp3 = 0.35*10^-6;
dp4 = 0.45*10^-6;
dp = [dp1 dp2 dp3 dp4];
distribution = [0.1 0.3 0.4 0.2]; % Distribution percentage
c = ones(Nr,1);
cini = cini .* distribution;
c = c.*cini; % Initial conditions
result = zeros(length(t), length(c), length(dp));
end
%% Solver
tic % Measure time of calculation
tolerance = 1e-6;
options = odeset('RelTol',tolerance,'AbsTol',tolerance,'NonNegative',1:100);
% ODE
[t,C] = ode15s((@(t,c)odebatch31(t,c, Nr,dr,cmax,rhod,rhow,dp,g,f,method)),t,c,options);
for i = 0:length(dp)-1
result(:,:,i+1) = C(:,Nr*i+1:Nr*i+Nr);
end
toc % Measure time of calculation
sum_result = sum(result, 3); % Add distributions together to get local concentration
% Mass distribution loop
for j = size(sum_result,2):-1:1
for i = 1:size(sum_result,1)
if sum_result(i,j)>cmax
excess(i,j) = sum_result(i,j)-cmax;
sum_result(i,j) = cmax;
sum_result(i,j-1) = sum_result(i,j-1)+excess(i,j);
end
end
end
cmax_reached = zeros(length(t),1);
for i = 1:length(t)
[cmax, max_col] = max(sum_result(i,:));
cmax_reached(i) = max_col;
if cmax_reached(i) == 1
cmax_reached(i) = 100;
end
end
%% Control on c
Mbegin = sum(sum_result(1,:))*rhod + (1-sum(sum_result(1,:)))*rhow;
Mend = sum(sum_result(end,:))*rhod + (1-sum(sum_result(end,:)))*rhow;
Mbal = sum( Mend) - sum( Mbegin );
if Mbegin-Mend>1e-3 || Mbegin-Mend < -1e-3
fprintf('Unequal mass at beginning and end \n')
fprintf('Mass balance = %f kg/m3 \n',Mbal)
end
%% figure for check
figure(1)
xaxis = linspace(0,r*100,Nr);
yaxis = t/60;
surf(xaxis,yaxis,sum_result)
xlabel('r [cm]')
ylabel('time [min]')
zlabel('volumetric fraction [-]')
colorbar
%% Protein concentration and figures for checking
c_water = 1-sum_result;
c_prot = c_water*cprot;
pellet_thickness = (Nr-cmax_reached)*dr;
polynomal = polyfit(t, pellet_thickness, 3);
polyval = polyval(polynomal,t);
pos = 0.2*Nr; % choose percentage of time
tpos = pos*tstep;
figure(2)
plot(xaxis,c_prot(pos,:), xaxis,sum_result(pos,:))
title(['concentration through centrifuge at ', num2str(tpos/60,'%.1f'),' min'])
legend('protein','cell debris')
xlabel('r [cm]')
ylabel('fraction[-]')
% figure(3)
% plot(t, pellet_thickness,t,polyval, 'or')
% legend('pellet size', 'polynomal')
% xlabel('time [s])')
% ylabel('pellet size [m]')
%% calculation for harvestable proteins
% Let's say that proteins can be harvested when the concentration of cell
% cebris is lower than 0.1. Now, we can calculate the protein in
% each cell, translate it to volumetric cells and therefore to a mass, and
% lastly calculate the total mass that can be harvested.
c_prot_harvest = c_water*cprot;
c_prot_harvest(sum_result>0.1)=0;
m_prot_harvest = c_prot_harvest.*Vcell'*rhop*Ny*Nz;
figure(4)
yyaxis left
plot(xaxis,m_prot_harvest(pos,:))
ylabel('mass [kg]')
yyaxis right
plot(xaxis,c_prot_harvest(pos,:))
title(['concentration through centrifuge at ', num2str(tpos/60,'%.1f'),' min'])
legend('protein mass','protein fraction')
xlabel('r [cm]')
ylabel('fraction [-]')
harvest = zeros(length(t),1); slope = zeros(length(t)-1,1);
for i = 1:length(t)
harvest(i) = sum(m_prot_harvest(i,:));
end
flag = false;
for i = 20:length(t)
slope(i) = (harvest(i)-harvest(i-1))/tstep;
if all(harvest(i)-harvest(i-1)<0.00005) && ~flag
disp(['Optimal centrifugal time is ', num2str((i*tstep)/60, '%.1f'), ' min'])
disp(['Protein harvest at this time is ', num2str(harvest(i), '%.2f'),' kg'])
disp(['The yield on protein is ', num2str(((harvest(i)/mprot_ini)*100), '%.2f'),'%'])
flag = true;
revenue(i) = harvest(i)*Cp;
energy = (rho_in*(2*pi*f*r)^2*Vt./t); % J
costs = energy/1000*Ce;
profit(i) = revenue(i)-costs;
disp(['The profit is €',num2str(profit(i), '%.2f')])
end
end
revenue = harvest*Cp;
energy = rho_in*(2*pi*f*r)^2*Vt./t; % J
costs = energy/1000*Ce;
profit = revenue-costs;
figure(5)
plot(t/60,(harvest/mprot_ini))
xlabel('time [min]')
ylabel('mass [kg]')
title('Protein harvest over time')
yyaxis right
plot(t/60, profit)
ylabel('profit [€]')
function [dcdt] = odebatch31(~,c, Nr,dr,cmax,rhod,rhow,dp,g,f,method)
%% zeros matrices
v = zeros(Nr,length(dp));
excess = zeros(Nr,1);
dcdt = zeros(Nr,length(dp));
Re = zeros(Nr,length(dp));
eta = zeros(Nr,1);
vgs = zeros(Nr,length(dp));
Ar = zeros(Nr,length(dp));
n_check = zeros(Nr,length(dp));
n = zeros(Nr,length(dp));
v_check = zeros(Nr,length(dp));
Re_check = zeros(Nr,length(dp));
csum = zeros(Nr,1);
c = reshape(c, Nr, length(dp));
r = zeros(Nr,1);
gc = zeros(Nr,1);
%% Mass distribution loop
switch method
case 'one size'
for i=Nr:-1:2
if c(i) >= cmax
excess(i) = c(i)-cmax;
c(i) = cmax;
c(i-1) = c(i-1)+excess(i);
end
end
% Move location of wall condition
saturated = zeros(Nr, 1);
for i = 1:Nr
if c(i) == cmax
saturated(i) = 1;
end
end
wall_coords = find(saturated == 1);
if isempty(wall_coords)
wall_condition = Nr;
else
wall_condition = wall_coords(1)-1;
end
case 'distribution'
for i=Nr:-1:2
clocal(i,:) = c(i,:);
clocal_sum(i) = sum(clocal(i));
if clocal_sum(i) >= cmax
distribution(i,:) = clocal(i,:)./clocal(i)*cmax;
c(i,:) = distribution(i,:);
excess(i,:) = clocal(i,:)-distribution(i,:);
c(i-1,:) = c(i-1,:)+excess(i,:);
end
end
% Move location of wall condition
saturated = zeros(Nr, 1);
for i = 1:Nr
if c(i) == cmax
saturated(i) = 1;
end
end
wall_coords = find(saturated == 1);
if isempty(wall_coords)
wall_condition = Nr;
else
wall_condition = wall_coords(1)-1;
end
end
%% Main loop
for i=1:wall_condition
%% find values for velocity and Reynolds check
for j = 1:length(dp)
dparticle(j) = dp(1,j);
eta(i) = exp(2.5*(sum(c(i,j)))./(1-sum(c(i,j))/cmax));
r(i) = r(i)+dr;
gc(i) = (2*pi*f)^2*r(i);
eta(i) = 0.001;
vgs(i,j) = ((rhod-rhow).*dparticle(j).^2.*gc(i))./(18*eta(i));
Ar(i,j) = (dparticle(j).^3*rhow*(rhod-rhow)*g)./eta(i).^3;
n_check(i,j) = (4.8+2.4*0.015.*Ar(i,j).^0.5)/(0.015.*Ar(i,j).^0.5+1);
v_check(i,j) = (1-c(i,j)).^n_check(i,j).*vgs(i,j);
Re_check(i,j) = (rhod*v_check(i,j).*dparticle(j))./eta(i);
n(i,j) = (4.8-2.4*Re_check(i,j).^(3/4)*0.015)/(1+Re_check(i,j).^(3/4)*0.015);
v(i,j) = (1-sum(c(i,j))).^n(i,j).*vgs(i,j);
Re(i,j) = (rhod.*v(i,j).*dparticle(j))./eta(i);
if Re(i,j)/Re_check(i,j)<0.95 || Re(i,j)/Re_check(i,j)>1.05
fprintf('Re difference too big \n')
return
end
if Re(i,j)>1
fprintf('Re too big for stokes \n')
return
end
%% ODE calculation
if i > 1
dcdt(i,j) = (1/dr).*((c(i-1,j).*v(i-1,j))-c(i,j).*v(i,j)); % general ode
else
dcdt(i,j)=(1/dr).*-(c(i,j).*v(i,j)); % initial condition
end
if i == wall_condition % wall condition
dcdt(i,j) = (1/dr).*(c(i-1,j).*v(i-1,j));
end
csum(i) = sum(c(i,:));
end
end
dcdt = dcdt(:);
end
When running the code, I do not get any errors, only after everything is run through the function and the ode is calculated, I get the following error:
Arrays have incompatible sizes for this operation.
Error in ode15s (line 553)
rhs = hinvGak*ode(tnew,ynew) - Mtnew*(psi+difkp1);
Error in centrifugation31 (line 83)
[t,C] = ode15s((@(t,c)odebatch31(t,c, Nr,dr,cmax,rhod,rhow,dp,g,f,method)),t,c,options);
I attached my script and function, so that you can maybe understand it a bit better.
I am looking forward to a reply, and a big thanks to anyone that wants to take a look!
0 Comments
Accepted Answer
Torsten
on 2 Oct 2023
Edited: Torsten
on 2 Oct 2023
if Re(i,j)/Re_check(i,j)<0.95 || Re(i,j)/Re_check(i,j)>1.05
fprintf('Re difference too big \n')
return
end
if Re(i,j)>1
fprintf('Re too big for stokes \n')
return
end
You cannot simply return from odebatch31 ; the solver will crash in this case because it expects at least any output for dcdt of size 400x1.
When I comment these lines of your code, I get an error in this if-statement:
if clocal_sum(i) >= cmax
distribution(i,:) = clocal(i,:)./clocal(i)*cmax;
c(i,:) = distribution(i,:);
excess(i,:) = clocal(i,:)-distribution(i,:);
c(i-1,:) = c(i-1,:)+excess(i,:);
end
The line
distribution(i,:) = clocal(i,:)./clocal(i)*cmax;
looks incorrect.
More Answers (0)
See Also
Categories
Find more on Mathematics and Optimization 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!