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
