Error in ode15s: incompatible sizes arrays, while no errors in function

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);

Re difference too big 

Arrays have incompatible sizes for this operation.

Error in ode15s (line 553)
rhs = hinvGak*ode(tnew,ynew) - Mtnew*(psi+difkp1);

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!

Answer 1

[Torsten]

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.

[Jort Puiman]

Thank you so much! I completely overlooked that return would have this effect.

For the other comment, I am aware that my 'distribution' method is not working properly. I have a post-ode fix, but this is not how it should be. I will look at how to fix this.