Arrays have incompatible sizes for this operation. Error in (line 34), can you please tell me how to solve it?

% Define problem parameters

L = 1; % size of domain

Ul = 1; % velocity of upper lid

Re = 1000; % Reynolds number

nu = Ul*L/Re; % viscosity

nx = 40; % number of cells in x direction

ny = 40; % number of cells in y direction

dx = L/nx; % cell size in x direction

dy = L/ny; % cell size in y direction

dt = 0.001; % time step

nt = 100; % number of time steps

% Create grid

x = linspace(0, L, nx+1);

y = linspace(0, L, ny+1);

[X, Y] = meshgrid(x, y);

% Set initial conditions

u = zeros(nx+1, ny+2);

v = zeros(nx+2, ny+1);

p = zeros(nx+2, ny+2);

% Define coefficients for finite volume discretization

ae = repmat(1/nu*dt/dx^2, nx+2, ny+2);

aw = repmat(1/nu*dt/dx^2, nx+2, ny+2);

an = repmat(1/nu*dt/dy^2, nx+2, ny+2);

as = repmat(1/nu*dt/dy^2, nx+2, ny+2);

ap = ae+aw+an+as+repmat(1/dt, nx+2, ny+2);

b = zeros(nx+2, ny+2);

% Main loop for time integration

for n = 1:nt

% Calculate intermediate velocities

ue = u(2:end-1, 2:end-1) + dt*(-(p(2:end-1, 3:end)-p(2:end-1, 2:end-1))/dy + nu*(u(3:end, 2:end-1)-2*u(2:end-1, 2:end-1)+u(1:end-2, 2:end-1))/dx^2);

uw = u(1:end-2, 2:end-1) + dt*(-(p(2:end-1, 2:end-1)-p(2:end-1, 1:end-2))/dy + nu*(u(2:end-1, 2:end-1)-2*u(1:end-2, 2:end-1)+u(1:end-3, 2:end-1))/dx^2);

vn = v(2:end-1, 2:end-1) + dt*(-(p(3:end, 2:end-1)-p(2:end-1, 2:end-1))/dx + nu*(v(2:end-1, 3:end)-2*v(2:end-1, 2:end-1)+v(2:end-1, 1:end-2))/dy^2);

vs = v(2:end-1, 1:end-2) + dt*(-(p(2:end-1, 2:end-1)-p(1:end-2, 2:end-1))/dx + nu*(v(2:end-1, 2:end-1)-2*v(2:end-1, 1:end-2)+v(2:end-1, 1:end-3))/dy^2);

end

% Impose boundary conditions

% Velocity boundary conditions

u(1,:) = Ul; % Top wall

v(1,:) = 0;

u(end,:) = 0; % Bottom wall

v(end,:) = 0;

u(:,1) = 0; % Left wall

v(:,1) = 0;

u(:,end) = 0; % Right wall

v(:,end) = 0;

% Pressure boundary conditions

% Top wall (Neumann boundary condition)

for i = 2:Nx-1

Ap =zeros (1, i);

Ap(i,end) = 0; % diagonal element

Aw =zeros (1, i);

Aw(i,end) = 0; % west coefficient

Ae = zeros (1, i);

Ae(i,end) = 0; % east coefficient

An = zeros (1, i);

An(i,end) = 1/dy; % north coefficient

As = zeros (1, i);

As(i,end) = -1/dy; % south coefficient

b(i,end) = 0; % right-hand side

end

% Bottom wall (Neumann boundary condition)

for i = 2:Nx-1

Ap(i,1) = 0; % diagonal element

Aw(i,1) = 0; % west coefficient

Ae(i,1) = 0; % east coefficient

An(i,1) = -1/dy; % north coefficient

As(i,1) = 1/dy; % south coefficient

b(i,1) = 0; % right-hand side

end

% Left wall (Neumann boundary condition)

for j = 2:Ny-1

Ap(1,j) = 1; % diagonal element

Aw(1,j) = 0; % west coefficient

Ae(1,j) = -1/dx; % east coefficient

An(1,j) = 0; % north coefficient

As(1,j) = 0; % south coefficient

b(1,j) = 0; % right-hand side

end

% Right wall (Dirichlet boundary condition)

for j = 2:Ny-1

Ap(end,j) = 1; % diagonal element

Aw(end,j) = -1/dx; % west coefficient

Ae(end,j) = 0; % east coefficient

An(end,j) = 0; % north coefficient

As(end,j) = 0; % south coefficient

b(end,j) = 0; % right-hand side

end

[Dyuman Joshi]

As p has different size than u and v, using 2:end-1 or 3:end results in a larger array for p as compared to u and v.

I suggest you to modify the indices accordingly.

L = 1;      % size of domain
Ul = 1;     % velocity of upper lid
Re = 1000;  % Reynolds number
nu = Ul*L/Re;   % viscosity
nx = 40;    % number of cells in x direction
ny = 40;    % number of cells in y direction
dx = L/nx;  % cell size in x direction
dy = L/ny;  % cell size in y direction
nt = 100;   % number of time steps
% Set initial conditions
u = zeros(nx+1, ny+2);
v = zeros(nx+2, ny+1);
p = zeros(nx+2, ny+2);

%for n=1:nt
%Breaking down the formula of intermediate velocity
%ue
A1=-(p(2:end-1, 3:end)-p(2:end-1, 2:end-1))/dy

A1 = 40×40

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

A2=nu*(u(3:end, 2:end-1)-2*u(2:end-1, 2:end-1)+u(1:end-2, 2:end-1))/dx^2

A2 = 39×40

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

%end

Answer 1

[Ashu]

Hi Nusrat

I understand that you are facing error with incompatible array sizes.

Few corrections that you need to consider are -

• Rename nx, ny to Nx, Ny respectively, because in the further code you are refering them with the later variable names.

Nx = 40;    % number of cells in x direction
Ny = 40;    % number of cells in y direction

• Check your array operations and their sizes in the 'Main operation loop'

Example analysis -

ue = u(2:end-1, 2:end-1) + dt*(-(p(2:end-1, 3:end) - p(2:end-1, 2:end-1))/dy + nu*(u(3:end, 2:end-1)-2*u(2:end-1, 2:end-1)+u(1:end-2, 2:end-1))/dx^2);
    % size = 39 x 40                size = 40 x 40    size = 40 x 40         

It's clear here that the size for first operand and the second operand don't match.

A simple correction can be like this

    ue = u(2:end-1, 2:end-1) + dt*(-(p(2:end-2, 3:end) - p(2:end-2, 2:end-1))/dy + nu*(u(3:end, 2:end-1)-2*u(2:end-1, 2:end-1)+u(1:end-2, 2:end-1))/dx^2);
                                     %------changes in the index----------%

These changes will depend upon your use case, so you can change accordingly.