Clear Filters
Clear Filters

I have an indices error in my code

3 views (last 30 days)
Thierry
Thierry on 2 Jul 2023
Answered: Star Strider on 2 Jul 2023
Dear Matlab community,
I have written a code that presents the following error at line 145: index in position 1 is invalid. Array indices must be positive integers or logical values. I don't know how to get rid of this issue. Thank you for your help. Regards, Thierry.
set(0, 'DefaultFigurePosition', [100, 100, 800, 800]); % Adjust the figure size as needed
% BLOCK 1 *********************************************************
% Set up the external forcing
g = 9.81; % Gravity
A = 0.1*g; % amplitude of Faraday wave perturbation
fFaradX = 4; % frequency of Faraday wave perturbation in the x direction (odd values only!)
fFaradY = 11; % frequency of Faraday wave perturbation in the y direction (odd values only!)
% Set simulation parameters
radius = 10; % radius of circular container
center = [0, 0]; % center coordinates of the circular container
c_s = 1./(sqrt(3)); % lattice "speed of sound" used for the equilibrium state computation feq
tmax = 50; % total simulation time
dt = 1; % time step
dx = 1; % grid spacing
omega = 1.9; % relaxation parameter (=1/(relaxation time), omega=1/tau)
rho0 = 1; % fluid density
mu = 0.1; % fluid dynamic viscosity
nu = mu/rho0; % fluid kinematic viscosity
m = 10000; % number of floaters
R = 0.1; % radius of floaters
theta0 = rand(m,1)*2*pi; % initial angles of floaters
% Convert circular container to a rectangular grid
Lx = ceil(radius) + 2; % length of container in x-direction
Ly = ceil(radius) + 2; % length of container in y-direction
% Calculate grid coordinates
[xGrid, yGrid] = meshgrid(1:Lx, 1:Ly);
xGrid = xGrid - ceil(radius) - 1;
yGrid = yGrid - ceil(radius) - 1;
% Calculate distance from the center of the circular container
dist = sqrt(xGrid.^2 + yGrid.^2);
% Create a mask to represent the circular container
containerMask = dist <= radius;
% BLOCK 2 *****************************************************
% Set lattice parameters
q = 9; % number of lattice velocities // number of f(i) discrete populations
e = [0,0; 1,0; 0,1; -1,0; 0,-1; 1,1; -1,1; -1, -1; 1, -1]; % Discrete lattice Boltzmann velocity vectors
w = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; % lattice weights
% BLOCK 3 ******************************************************
% Initialize fluid and floaters
rho = ones(Lx,Ly); % fluid density
f = zeros(Lx,Ly,q); % distribution function
feq = zeros(Lx,Ly,q); % equilibrium distribution function
ux = zeros(Lx,Ly); % fluid x-velocity (macroscopic velocity)
uy = zeros(Lx,Ly); % fluid y-velocity (macroscopic velocity)
x = rand(m,1)*Lx; % floater x-positions
y = rand(m,1)*Ly; % floater y-positions
vx = zeros(m,1); % floater x-velocity
vy = zeros(m,1); % floater y-velocity
% Compute the center coordinates of the circular container
cx = Lx/2;
cy = Ly/2;
% Compute the radius of the circular container
radius = min(Lx, Ly) / 2;
% Compute the distance from each lattice point to the center of the container
distance = sqrt((xGrid - cx).^2 + (yGrid - cy).^2);
% Identify the lattice points that are inside the circular container
insideContainer = distance <= radius;
% Set the fluid density and velocity to zero outside the circular container
rho(~insideContainer) = 0;
ux(~insideContainer) = 0;
uy(~insideContainer) = 0;
% Initialize the distribution function and equilibrium distribution function
for ix = 1:Lx
for iy = 1:Ly
if insideContainer(ix, iy)
ux(ix, iy) = A * (cos(fFaradX * pi * (ix / Lx))^2);
uy(ix, iy) = A * (cos(fFaradY * pi * (iy / Ly))^2);
for j = 1:q
feq(ix, iy, j) = w(j) * rho(ix, iy) * (1 + (e(j, 1) * ux(ix, iy) + e(j, 2) * uy(ix, iy)) / c_s^2 + ...
(e(j, 1) * ux(ix, iy) + e(j, 2) * uy(ix, iy))^2 / (2 * c_s^4) - ...
(ux(ix, iy)^2 + uy(ix, iy)^2) / (2 * c_s^2));
f(ix, iy, j) = feq(ix, iy, j);
end
end
end
end
% Main simulation loop
for t = 1:tmax
for ix = 1:Lx
for iy = 1:Ly
if insideContainer(ix, iy)
% Compute equilibrium distribution
rho(ix, iy) = sum(f(ix, iy, :));
ux(ix, iy) = (sum(f(ix, iy, [2, 5, 8])) - sum(f(ix, iy, [4, 7, 9]))) / rho(ix, iy);
uy(ix, iy) = (sum(f(ix, iy, [3, 6, 7])) - sum(f(ix, iy, [1, 5, 8]))) / rho(ix, iy);
for j = 1:q
% Streaming step
ix_n = mod(ix + e(j, 1), Lx);
iy_n = mod(iy + e(j, 2), Ly);
f(ix_n, iy_n, j) = f(ix, iy, j);
end
end
end
end
% Update macroscopic variables
for ix = 1:Lx
for iy = 1:Ly
if insideContainer(ix, iy)
rho(ix, iy) = sum(f(ix, iy, :));
ux(ix, iy) = (sum(f(ix, iy, [2, 5, 8])) - sum(f(ix, iy, [4, 7, 9]))) / rho(ix, iy);
uy(ix, iy) = (sum(f(ix, iy, [3, 6, 7])) - sum(f(ix, iy, [1, 5, 8]))) / rho(ix, iy);
end
end
end
% Collision of distribution functions (BGK type)
for j = 1:q
for ix = 1:Lx
for iy = 1:Ly
if insideContainer(ix, iy)
feq(ix, iy, j) = w(j) * rho(ix, iy) * (1 + (e(j, 1) * ux(ix, iy) + e(j, 2) * uy(ix, iy)) / c_s^2 + ...
(e(j, 1) * ux(ix, iy) + e(j, 2) * uy(ix, iy))^2 / (2 * c_s^4) - ...
(ux(ix, iy)^2 + uy(ix, iy)^2) / (2 * c_s^2));
f(ix, iy, j) = (1 - omega) * f(ix, iy, j) + omega * feq(ix, iy, j);
end
end
end
end
% Update floater positions
for i = 1:m
% Update the floater velocity based on the fluid velocity at its position
ix = mod(round(x(i)), Lx) + 1; % Apply modulo and add 1 to ensure indices are within the valid range
iy = mod(round(y(i)), Ly) + 1; % Apply modulo and add 1 to ensure indices are within the valid range
vx(i) = (sum(f(ix, iy, [2, 5, 8])) - sum(f(ix, iy, [4, 7, 9]))) / rho(ix, iy);
vy(i) = (sum(f(ix, iy, [3, 6, 7])) - sum(f(ix, iy, [1, 5, 8]))) / rho(ix, iy);
% Update the floater position based on its velocity
x(i) = x(i) + vx(i) * dt;
y(i) = y(i) + vy(i) * dt;
% Apply periodic boundary conditions for floater positions
x(i) = mod(x(i), Lx);
y(i) = mod(y(i), Ly);
end
% Plot the fluid and floater positions (optional)
if mod(t, 2) == 0
% Plot fluid velocity field
quiver(1:Lx, 1:Ly, ux', uy');
hold on;
% Plot floater positions
scatter(x, y, 'filled', 'r');
hold off;
axis([0 Lx 0 Ly]);
title(sprintf('Fluid and Floater Positions at Time Step %d', t));
drawnow;
end
end

Answers (2)

Mihir
Mihir on 2 Jul 2023
The error you mentioned, "index in position 1 is invalid," is likely caused by the indexing in the line where you update the floater positions. To fix this issue, you need to modify the indexing for the `ix` and `iy` variables to ensure they are within the valid range. Here's the corrected code:
% Update floater positions
for i = 1:m
% Update the floater velocity based on the fluid velocity at its position
ix = mod(round(x(i)), Lx) + 1; % Apply modulo and add 1 to ensure indices are within the valid range
iy = mod(round(y(i)), Ly) + 1; % Apply modulo and add 1 to ensure indices are within the valid range
vx(i) = (sum(f(ix, iy, [2, 5, 8])) - sum(f(ix, iy, [4, 7, 9]))) / rho(ix, iy);
vy(i) = (sum(f(ix, iy, [3, 6, 7])) - sum(f(ix, iy, [1, 5, 8]))) / rho(ix, iy);
% Update the floater position based on its velocity
x(i) = x(i) + vx(i) * dt;
y(i) = y(i) + vy(i) * dt;
% Apply periodic boundary conditions for floater positions
x(i) = mod(x(i), Lx);
y(i) = mod(y(i), Ly);
end
Unrecognized function or variable 'm'.
In the code above, I added `+ 1` to the `ix` and `iy` indexing lines to shift the indices from the range `[0, Lx-1]` and `[0, Ly-1]` to `[1, Lx]` and `[1, Ly]`, respectively. This adjustment ensures that the indices are valid for MATLAB array indexing.

Star Strider
Star Strider on 2 Jul 2023
The problem is that in the ‘Update floater positions’ loop, both ‘ix’ and ‘iy’ are NaN at some point, and that throws the indexing error.
Adding this if, elseif, else block avoids the indexing error, however there are still problems because nothing plots afterward, because ‘f’ and ‘rho’ are zero everywhere, so the calculations involving them are all NaN, and only finite (and real) results plot.
% Update floater positions
for i = 1:m
% Update the floater velocity based on the fluid velocity at its position
ix = mod(round(x(i)), Lx) + 1; % Apply modulo and add 1 to ensure indices are within the valid range
iy = mod(round(y(i)), Ly) + 1; % Apply modulo and add 1 to ensure indices are within the valid range
if ~all(isnan([ix iy]))
vx(i) = (sum(f(ix, iy, [2, 5, 8])) - sum(f(ix, iy, [4, 7, 9]))) / rho(ix, iy);
vy(i) = (sum(f(ix, iy, [3, 6, 7])) - sum(f(ix, iy, [1, 5, 8]))) / rho(ix, iy);
% Update the floater position based on its velocity
x(i) = x(i) + vx(i) * dt;
y(i) = y(i) + vy(i) * dt;
% Apply periodic boundary conditions for floater positions
x(i) = mod(x(i), Lx);
y(i) = mod(y(i), Ly);
else
break
end
end
So this is just a termporary fix for the indexing problem. You will need to correct the ‘ix’ and ‘iy’ as well as the other calculations to provide a defninitive solution.
Also, in the set call, it is preferable to use groot rather than 0.
.

Categories

Find more on Fluid Dynamics in Help Center and File Exchange

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!