Sparse matrix division results in "array exceeds maximum array size preference" in one case but not another

I've written a finite-difference frequency-domain program to simulate a silicon ring resonator. The basics of it is that I'm creating a matrix that performs the Laplacian operator on a 2D grid which has a geometry meshed across it.

I've got an absorbing boundary condition set up, and by default it's a PML on the edges of the grid, but for my project I'm testing how it works with a circular PML (there's a paper that claims it handles reflections better and thus solves faster). So far, everything else works perfectly. I can visualize the wave in the waveguide and how its transmission gets blocked or allowed by the ring resonator, it's very cool. But when I change the geometry of the absorbing boundary, Matlab won't even calculate it.

The line of interest is

A = DHX/URyy*DEX + DHY/URxx*DEY + ERzz;

These are all sparse matrices. The 2D geometric grid is full 1292 x 1129. I wrote an absorbing boundary function which adds a PML to the grid either on the edges, or creates a circular PML (both are ~10 units thick).

This results in sparse derivative (DHX, DEX, DHY, DEX) and dielectric (URyy, URxx, ERzz) matrices that are 1,458,668 squared. These are sparse, the derivative matrices are only populated on the 0 and 1 diagonal, and the dielectric matrices are only populated on the diagonal. Regardless of how the PML is added, it's transformed from the full 2D geometric grid to the "unrolled" sparse matrix with the following:

ERxx = diag(sparse(ERxx(:)));

For some reason, when I add the boundary condition at the edges, the A = DHX/URyy etc line runs no problem. But when I create the circular boundary condition, it results in the following error:

"Error using /

Requested 1458668x1458668 (15852.7GB) array exceeds maximum array size preference (15.8GB). This might cause MATLAB to become unresponsive."

Doing some debugging, for some reason Matlab is treating URyy and URxx as full matrices only when I set up the original 2D grid with a circular boundary condition.

I cannot find a good reason this is happening. In both cases, a thin region of the 2D grid is some complex value and 1 everywhere else, and then this 2D grid is turned into a large sparse diagonal matrix with only values on the main diagonal. I've thoroughly checked and there is no other difference. I can understand the final results being incorrect and my FDFD simulation giving odd results, but it literally won't do the calculation because it's treating URxx as a full "1458668x1458668 (15852.7GB) array" only in the case of interest, for no apparent reason.

The code where I implement the boundary condition is below. You can see where it implements either the rectangular boundary (no error) or circular boundary (works, but results in error when dividing). Also in case it matters, I wrote the bulk of it at home using 2024b, but wrote the circular PML on my work laptop which has 2022a where I'm getting the error. Any help would be very appreciated, let me know if any more info is needed. Thanks!!

function [eps_xx, eps_yy, eps_zz, mu_xx, mu_yy, mu_zz] = cs_upml(eps_map, mu_map, pml_thickness, circ)
% ADDUPML2D - Apply UPML (Uniaxial Perfectly Matched Layer) to a 2D Yee grid
%
% Inputs:
%   eps_map       - Relative permittivity on 2x Yee grid
%   mu_map        - Relative permeability on 2x Yee grid
%   pml_thickness - [x_lo x_hi y_lo y_hi] PML region thickness (1x grid units)
%   circ          - 1 for circular PML, 0 for rectangular
%
% Outputs:
%   eps_xx, eps_yy, eps_zz - PML-modified permittivity tensor components
%   mu_xx, mu_yy, mu_zz    - PML-modified permeability tensor components

%% --- PML configuration ---
stretch_max = 4;      % Maximum coordinate stretching
loss_max    = 1;      % Maximum conductivity scaling
poly_order  = 3;      % Polynomial order for grading

%% --- Grid setup ---
[nx, ny] = size(eps_map);

% Convert PML thickness from 1x to 2x grid
[x_lo, x_hi, y_lo, y_hi] = deal(2 * pml_thickness(1), 2 * pml_thickness(2), 2 * pml_thickness(3), 2 * pml_thickness(4));

% Initialize stretch tensors
s_x = ones(nx, ny);
s_y = ones(nx, ny);

%% -- Rectangular PML (circ = 0)
if circ == 0
    % PML - x
    for m = 1:x_lo
        scale = 1 + (stretch_max - 1) * (m / x_lo)^poly_order;
        loss  = loss_max * sin(0.5 * pi * m / x_lo)^2;
        s_x(x_lo - m + 1, :) = scale * (1 - 1i * 60 * loss);
    end

    for m = 1:x_hi
        scale = 1 + (stretch_max - 1) * (m / x_hi)^poly_order;
        loss  = loss_max * sin(0.5 * pi * m / x_hi)^2;
        s_x(nx - x_hi + m, :) = scale * (1 - 1i * 60 * loss);
    end

    % PML - y
    for n = 1:y_lo
        scale = 1 + (stretch_max - 1) * (n / y_lo)^poly_order;
        loss  = loss_max * sin(0.5 * pi * n / y_lo)^2;
        s_y(:, y_lo - n + 1) = scale * (1 - 1i * 60 * loss);
    end

    for n = 1:y_hi
        scale = 1 + (stretch_max - 1) * (n / y_hi)^poly_order;
        loss  = loss_max * sin(0.5 * pi * n / y_hi)^2;
        s_y(:, ny - y_hi + n) = scale * (1 - 1i * 60 * loss);
    end

end

%% -- Circular PML (circ = 1)
if circ == 1
    centerx = nx/2;
    centery = ny/2;
    %s_r = ones(nx,ny);

    for n = 1:y_hi-1
        scale = 1 + (stretch_max - 1) * (n / y_hi)^poly_order;
        loss  = loss_max * sin(0.5 * pi * n / y_hi)^2;
        radius = (ny - y_hi + n) - centery;
        for theta = 1:18000
            x1 = floor(centerx - radius*cos((theta/100-90)*pi/180));
            x2 = floor(centerx + radius*cos((theta/100-90)*pi/180));
            y = floor(centery + radius*sin((theta/100-90)*pi/180));

            s_x(x1,y) = scale * (1 - 1i * 60 * loss) * cos((theta/100-90)*pi/180);
            s_x(x2,y) = scale * (1 - 1i * 60 * loss) * cos((theta/100-90)*pi/180);

            s_y(x1,y) = scale * (1 - 1i * 60 * loss) * sin((theta/100-90)*pi/180);
            s_y(x2,y) = scale * (1 - 1i * 60 * loss) * sin((theta/100-90)*pi/180);
        end

    end

end

%% --- Apply coordinate stretching to tensors ---
eps_tensor_xx = eps_map ./ s_x .* s_y;
eps_tensor_yy = eps_map .* s_x ./ s_y;
eps_tensor_zz = eps_map .* s_x .* s_y;

mu_tensor_xx  = mu_map  ./ s_x .* s_y;
mu_tensor_yy  = mu_map  .* s_x ./ s_y;
mu_tensor_zz  = mu_map  .* s_x .* s_y;

%% --- Extract Yee grid-aligned tensor components ---
eps_xx = eps_tensor_xx(2:2:end, 1:2:end);  % E-field x-component
eps_yy = eps_tensor_yy(1:2:end, 2:2:end);  % E-field y-component
eps_zz = eps_tensor_zz(1:2:end, 1:2:end);  % E-field z-component

mu_xx  = mu_tensor_xx(1:2:end, 2:2:end);   % H-field x-component
mu_yy  = mu_tensor_yy(2:2:end, 1:2:end);   % H-field y-component
mu_zz  = mu_tensor_zz(2:2:end, 2:2:end);   % H-field z-component

end

[Matt J]

The information is insufficient for troubleshooting, since you give us neither the complete stack of error messages, including line numbers, nor the steps to run the code and reproduce the error. Regardless, you should be able trap the point where the error is occurring with,

>> dbstop if error

(or run with Pause On Errors set in the editor gui)

[Shridhar]

Thank you for your response.

I know where the error is, dbstop if error doesn't help. mrdivide and trying to multiply by inv(URyy) don't work either. inv(URyy) results in the same error, that's where the issue is happening.

Basically, I have a 2D grid that's 1 everywhere, but in one case it's some complex values at the borders, whereas in another it's complex values in a circle. I then "unroll" the 2D grid into one 1D array, and then create a sparse diagonal from it and use that in the calculation. In the first case, inv(diag(sparse(ERxx(:)))) is fine, in the other inv(diag(sparse(ERxx(:)))) gives the "Requested 1458668x1458668 (15852.7GB) array exceeds maximum array size" error as if Matlab is trying to treat it like a full array instead of sparse.

Full project is attached, MRR_Test.m is the main file. Run the first two sections. circ_pml = 1 gives the error, circ_pml = 0 results in no error. Full error is:

Error using inv

Requested 1458668x1458668 (15852.7GB) array exceeds maximum array size preference (15.8GB). This might cause MATLAB to become unresponsive.

Error in mrr_sim (line 162)

A = DHX/URyy*DEX + DHY/URxx*DEY + ERzz;

Error in MRR_Test (line 12)

[neff, R, T] = mrr_sim(lambda0, n_nom, rib_dim, ring_dim, graph, circ_pml);

Answer 1

[Matt J]

When circ_pml = 1, your URyy matrix conatins Inf values. This confuses inv,mrdivide, and other similar operators.

[Shridhar]

Yup this was it. A diagonal matrix is singular if and only if one of the primary elements is 0. I was sweeping the angle and creating a radial boundary condition and calculating its x and y components with cos and sin, which could be 0, making it uninvertible and also causing divide by 0 errors if I tried element operators.

Fixed for now by simply making the circle bigger than the area of interest so it's never at pi/2, and those areas just get a flat PML. Simulation works great now!

[Matt J]

I'm glad, but regardless, it would probably be more efficient to dispense with all the diag(sparse(___)) conversions. Diagonal matrix operations can usually be replaced with element-wise operators, like,

A = DHX./diag(URyy)'*DEX + DHY/diag(URxx)'*DEY + ERzz;