How to use a sparse Jacobian for a stiff ODE solver to avoid the error "array exceeds maximum array size preference"

I have a model that I want to run for 1000 years. It represents an area of the sea floor (x by y), 61 size classes of organisms (B), and the food source ("detritus"; R). When you consider both B and R, the model has 62 state variables.

I give the ode15s solver an x by y by 62 matrix, which it seems to vectorise internally, which means in my internal ode function I have to reshape the inputs back into the x by y by 62 matrix to then do my matrix multiplication equations.

On a small x by y grid (i.e. when the length of x is 7, and length of y is 14), my code works fine. I can run it for 1000 years in like a minute, using the odeset function like this:

options = odeset('Refine', 2, 'abstol', 1e-10, 'MaxStep', 100, 'InitialStep', 0.01)

But, scaling up from 7x14 grid to a larger 360x180 grid, I get this error:

"Requested 4017600x4017600 (120260.6GB) array exceeds maximum array size preference (15.9GB). This might

cause MATLAB to become unresponsive.

Error in odenumjac (line 126)

ydel = y(:,ones(1,ny)) + diag(del);

Error in ode15s (line 348)

[dfdy,Joptions.fac,nF] = odenumjac(odeFcn, {t,y,odeArgs{:}}, f0, Joptions); %#ok<CCAT>

Error in dynorun_log_2D (line 178)

[t,y] = ode15s(@benthic_dydt_log2D,[tstart tfinal], y0,options);"

It's trying to create a (360x180x62)^2 matrix, which I guess is the Jacobian for the whole system. For now, there are no interactions between points, so this could be a sparse identity matrix.

I've gone back to the smaller 7x14 grid size while I work out setting the Jacobian, and have set it using:

options = odeset('JPattern',S);

I've tried making S a fully sparse matrix of size (7x14x62)^2, which caused the model to take >30 minutes to run for 1000 years, when before it took a minute. I've also tried making S a sparse identity matrix of size (7x14x62)^2, which hasn't yet finished running...

Neither option seem suitable for the bigger grid size. I thought setting the Jacobian to be sparse would make the solver faster? I'm not sure what to try next.

[Torsten]

Are you sure a nonstiff solver like ODE45, together with a given tspan vector with more than two elements, will take too long to solve your equations ?

And I guess there is some exchange of B and R between the (x/y) cells so that you have to solve the complete nx * ny * 62 system in one run ?

[Sophy]

Thanks Torsten. I went with stiff, because at the beginning of the model run we need to be able to take tiny time steps. In the first few seconds/minutes, the smallest organisms (one end of B) change rapidly while the biggest don't react (the other end of B). Then with time the smallest organisms get to steady state, while the bigger ones then change over the next days/years. Towards the end of the model run as we approach steady state, nothing really changes, so we can end up taking a timestep of around 300 years. But until I've run the model I don't know exactly what timesteps to take that will best resolve the system, so I wouldn't know what tspan vector to give it. Plus how dynamic the system is depends on the initial conditions etc so really I need to solver to internally work out what timestep to take next. Can a nonstiff solver like ODE45 handle this?

And yes, eventually I'd like some exchange of B and R between the (x/y) cells, hence why I'm not just looping through the cells. But there won't be exchange between all (x/y) cells, just between neighbouring cells, so the jacobian will be mostly sparse, with off diagonal components.

[Torsten]

But until I've run the model I don't know exactly what timesteps to take that will best resolve the system, so I wouldn't know what tspan vector to give it. Plus how dynamic the system is depends on the initial conditions etc so really I need to solver to internally work out what timestep to take next. Can a nonstiff solver like ODE45 handle this?

Setting "tspan" doesn't prescribe the solver to take the timesteps specified in this vector. The internal choice of the time steps for both solvers (ode15s and ode45) is independent of the tspan vector. The reason why I suggested to prescribe the output times is that I fear that ODE45 will take much smaller step sizes than ODE15S, and if all the intermediate results are saved in the solution vector, you will get the same RAM problem as you had with using ODE15S. The advantage of using ODE45 is that it doesn't require a Jacobian, the disadvantage is that if your system is really stiff, it will take forever for a computation. So I suggest you just test it for a small system before thinking about possible solution with ODE15S.

[Sophy]

Ah I see. I've tested it like you suggested. With a 14x7 grid, ode15s takes 1.5 minutes, and ode45 (with tspan set as the same timestep outputs as ode15s chose) takes 9.5 minutes. Bit too long to scale up using ode45, plus I would like some (x/y) exchange eventually.

[Torsten]

As said, you can choose tspan arbitrarily - you don't need to use the one from ode15s. Only more than start and end point have to be specified.

And exchange between cells is also no problem when you use ode45.

Why do you think that without interactions between points, the Jacobian is the identity matrix ?

[Sophy]

Yes you're right - I only chose the same output timestepping for consistency in the comparison, and outputting at fewer timesteps doesn't make it much quicker (still >9 minutes).

If I vectorise all the (x/y) grid points into a 98 length vector, the matrix holding the information of the transfer between each grid cell would be a 98 by 98 identity matrix, if there were no exchange. Should it be the mass matrix, not the jacobian, that should be the identity matrix?

[Torsten]

If there is no exchange between the cells, the Jacobian matrix of the whole system is a block-diagonal matrix with 98 blocks each of size 62x62. Thus there should be 98*62*62 elements that differ from 0 in a (98*62)^2 matrix.

sparsity = 1-98*62*62/(98*62)^2

sparsity = 0.9898

The sparsity of the matrix will become smaller if there is exchange between the cells.

[Sophy]

Thanks. I've set the sparcity pattern like this:

n1 = 98;
n2 = 62;
C = cell(1,n1);
for n = 1:n1
    C{n} = sparse(ones(n2));
end
S = blkdiag(C{:}); 
options = odeset('JPattern',S);

When using this for the 14x7 model run with ode45, it finished in around 10 minutes like before, but in ode15s it doesn't look like it will finish. I tried running for just 1 year instead of 1000, and in ode15s with this new sparcity pattern it looks like the timestep just keeps getting smaller and smaller as it gets towards the end time.

I don't mind continuing with ode45 if 10 minutes is the quickest I can get with this 14x7 grid, especially because when I scale up I won't encounter the "array exceeds maximum array size preference" error again. But with ode15s before specifying the jacobian sparcity it was as quick as 1 minute, which would be much better for when I move to a bigger grid, but I'm blocked by the "array exceeds maximum array size preference" error unless I can get the sparcity pattern working quicker.

Answer 1

[Torsten]

Here is an example code how to create the sparsity pattern for a given (quite complicated) ODE function.

Use the function "S_pattern" to generate your individual sparsity pattern matrix S.

Note that ode15s will come into trouble if the sparsity pattern changes during the integration. So if, e.g., some exchange terms between the cells appear at a time t > 0 during integration that were not present at time t = 0, ode15s will start to make small time steps because the Jacobian it is working with has become incomplete.

ode45 does not need a Jacobian and thus no sparsity pattern. So it's no surprise that the computation time when using ode45 did not change.

% Set physical parameters
L = [160 10 25 10 160]*1e-6; % [m] MEA layer thicknesses
kt = [1.6 0.27 0.3 0.27 1.6]; % [W/mK] Thermal conductivity
rcp = [1.58 1.56 1.9 1.56 1.58]*1e6;    % [J/m3K] Volumetric heat capacity
Qtot = [0 10 100 120 15]*1e3;     % [W] Heat Source
%Qtot = zeros(1,5);
% Number of mesh points in zones
nx = [11 11 11 11 11];
% Set boundary conditions
T_start = 343;
T_end = T_start;
% Set initial conditions
T0 = 0;
% Set time span of integration
tspan = 0:0.001:0.025;

% Create grid in domains
% Number of domains
Nd = numel(L); 
% Grid vectors
xstart = 0;
xend = L(1);
x{1} = linspace(xstart,xend,nx(1)).';
for i = 2:Nd
  xstart = xend;
  xend = xstart + L(i);
  x{i} = linspace(xstart,xend,nx(i)).';
end

% Number of ODEs to be solved
nodes = sum(nx) - (Nd-1); 

% Create initial solution vector 
y0 = zeros(nodes,1);
y0(1) = T_start;
y0(2:nodes-1) = T0;
y0(nodes) = T_end;

S = S_pattern(@(t,y)fun(t,y,Nd,nx,x,kt,rcp,Qtot),tspan(1),y0);
figure(1)
spy(S)

% Call solver
[T,Y] = ode15s(@(t,y)fun(t,y,Nd,nx,x,kt,rcp,Qtot),tspan,y0,odeset('JPattern',S));

% Plot results
xplot = [x{1};x{2}(2:end);x{3}(2:end);x{4}(2:end);x{5}(2:end)];
figure(2)
plot(xplot,Y.');

% Function
function dy = fun(t,y,Nd,nx,x,kt,rcp,Qtot)

  % Write y in local variables
  ileft = 1;
  iright = nx(1);
  T{1} = y(ileft:iright);
  for j = 2:Nd
    ileft = iright ;
    iright = ileft + nx(j) - 1;
    T{j} = y(ileft:iright);
  end

  % Compute values in ghost points
  for j = 1:Nd-1
    T1 = T{j};
    T2 = T{j+1};
    x1 = x{j};
    x2 = x{j+1};
    lambda1 = kt(j);
    lambda2 = kt(j+1);
    rhocp1 = rcp(j);
    rhocp2 = rcp(j+1);
    nx1 = nx(j);
    nx2 = nx(j+1);
    Q1 = Qtot(j);
    Q2 = Qtot(j+1);

    xg1 = x1(end)+(x1(end)-x1(end-1));
    xg2 = x2(1)-(x2(2)-x2(1));
    % Set up linear system of equations for [Tg1, Tg2]
    a11 = lambda1/(xg1-x1(nx1-1));
    a12 = lambda2/(x2(2)-xg2);
    b1 = lambda1*T1(nx1-1)/(xg1-x1(nx1-1))+lambda2*T2(2)/(x2(2)-xg2);
    a21 = (lambda1/(xg1-x1(nx1))/((x1(nx1)+xg1)/2 - (x1(nx1)+x1(nx1-1))/2))*rhocp2;
    a22 = (-lambda2/(x2(1)-xg2)/((x2(2)+x2(1))/2 - (x2(1)+xg2)/2))*rhocp1;
    b2 = ((lambda1*T1(nx1)/(xg1-x1(nx1))+lambda1*(T1(nx1)-T1(nx1-1))/(x1(nx1)-x1(nx1-1)))/...
          ((x1(nx1)+xg1)/2-(x1(nx1)+x1(nx1-1))/2)+Q1)*rhocp2 +...
         ((lambda2*(T2(2)-T2(1))/(x2(2)-x2(1))-lambda2*T2(1)/(x2(1)-xg2))/...
          ((x2(2)+x2(1))/2 - (x2(1)+xg2)/2)+Q2)*rhocp1;
    A = [a11,a12;a21,a22];
    b = [b1;b2];
    % Solve linear system for [Tg1, Tg2]
    sol = A\b;
    Tg1(j) = sol(1);
    Tg2(j) = sol(2);
    Xg1(j) = xg1;
    Xg2(j) = xg2;
  end

  % Compute time derivatives
  dTdt{1}(1) = 0.0;
  for j = 1:Nd
    if j < Nd
      Tj = [T{j};Tg1(j)];
      xj = [x{j};Xg1(j)];
    else
      Tj = T{j};
      xj = x{j};
    end
    lambdaj = kt(j);
    rhocpj = rcp(j);
    Qj = Qtot(j);
    for i=2:numel(xj)-1
      dTdt{j}(i) = ((lambdaj*(Tj(i+1)-Tj(i))/(xj(i+1)-xj(i)) -...
                     lambdaj*(Tj(i)-Tj(i-1))/(xj(i)-xj(i-1)))/...
                     ((xj(i+1)+xj(i))/2-(xj(i)+xj(i-1))/2)+Qj)/(rhocpj);
    end
  end
  dTdt{Nd}(end+1) = 0.0;

  % Return time derivatives
  dy = dTdt{1}(1:end);
  for j = 2:Nd
    dy = [dy,dTdt{j}(2:end)];
  end
  dy = dy.';
 end
 
 function S = S_pattern(f,t,u)
  y = f(t,u);
  h = 1e-6;
  ir = [];
  ic = [];
  for i = 1:numel(u)
    u(i) = u(i) + h;
    yh = f(t,u);
    %J(:,i) = (yh-y)/h;
    Jh = (yh-y)/h;
    Jhnz = find(abs(Jh)>0);
    if ~isempty(Jhnz)
      ir = [ir;Jhnz];
      ic = [ic;i*ones(size(Jhnz,1),1)];   
    end
    u(i) = u(i) - h;
  end
  S = sparse(ir,ic,ones(numel(ir),1),numel(u),numel(u));
end