Thank you Walter Robertson for your answer. So I modified the code as follows: I put the initialization of feq, f, etc. after the parfor loop as follows:
parfor t=0:tmax
% Initialize fluid and floaters
f = zeros(Lx,Ly,q); % distribution function
feq = zeros(Lx,Ly,q); % equilibrium distribution function
ux = zeros(Lx,Ly); % fluid x-velocity
uy = zeros(Lx,Ly); % fluid y-velocity
x = x0; % floater x-positions
y = y0; % floater y-positions
vx = zeros(m,1); % floater x-velocities
vy = zeros(m,1); % floater y-velocities
% Collision step
for i=1:q
rho = ones(Lx,Ly)*rho0; % fluid density
feq(:,:,i) = w(i)*rho.*(1 + 3.*(e(i,1).*ux + e(i,2).*uy)./c + 9./2.*(e(i,1).*ux + ...
e(i,2).*uy).^2./c.^2 - 3/2.*(ux.^2 + uy.^2)./c.^2); % TODO check ^2
end
...
The code runs without an error but I don't get any plotting at the end of the code. Why?
Here again the updated code:
% Set simulation parameters
Lx = 50; % length of container in x-direction
Ly = 50; % length of container in y-direction
dx = 1; % grid spacing
dt = 0.1; % time step
c = dx/dt; % speed of sound (?)
tmax = 100; % total simulation time
omega = 1.8; % relaxation parameter
rho0 = 1; % fluid density
mu = 0.1; % fluid viscosity
nu = mu/rho0; % fluid kinematic viscosity
m = 100; % number of floaters
R = 1; % radius of floaters
x0 = rand(m,1)*Lx; % initial x-positions of floaters
y0 = rand(m,1)*Ly; % initial y-positions of floaters
% Set lattice parameters
q = 9; % number of lattice velocities
e = [0,0;0,1;1,0;0,-1;-1,0;1,1;1,-1;-1,-1;-1,1]; % lattice velocities % TODO: check order of arrows
w = [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]; % lattice weights
% Initialize fluid and floaters
% f = zeros(Lx,Ly,q); % distribution function
% feq = zeros(Lx,Ly,q); % equilibrium distribution function
rho = ones(Lx,Ly)*rho0; % fluid density
% ux = zeros(Lx,Ly); % fluid x-velocity
% uy = zeros(Lx,Ly); % fluid y-velocity
% x = x0; % floater x-positions
% y = y0; % floater y-positions
% vx = zeros(m,1); % floater x-velocities
% vy = zeros(m,1); % floater y-velocities
% Set up boundary conditions
A = 0.01; % amplitude of Faraday wave perturbation
fFarad = 5; % frequency of Faraday wave perturbation
bc_func = @(t) A*sin(2*pi*fFarad*t); % boundary function
% Set lattice parameters
q = 9; % number of lattice velocities
e = [0,0;0,1;1,0;0,-1;-1,0;1,1;1,-1;-1,-1;-1,1]; % lattice velocities % TODO: check order of arrows
w = [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]; % lattice weights
poolobj = gcp('nocreate'); % get the pool object, and do it avoiding creating a new one.
if isempty(poolobj) % check if there is not a pool.
poolsize = 12;
else
poolsize = poolobj.NumWorkers;
delete( gcp('nocreate')); % delete the current pool object.
end
parpool( poolsize, 'IdleTimeout',Inf); % create a new pool with the previous poolsize and new specs.
parfor t=0:tmax
% Initialize fluid and floaters
f = zeros(Lx,Ly,q); % distribution function
feq = zeros(Lx,Ly,q); % equilibrium distribution function
ux = zeros(Lx,Ly); % fluid x-velocity
uy = zeros(Lx,Ly); % fluid y-velocity
x = x0; % floater x-positions
y = y0; % floater y-positions
vx = zeros(m,1); % floater x-velocities
vy = zeros(m,1); % floater y-velocities
% Collision step
for i=1:q
rho = ones(Lx,Ly)*rho0; % fluid density
feq(:,:,i) = w(i)*rho.*(1 + 3.*(e(i,1).*ux + e(i,2).*uy)./c + 9./2.*(e(i,1).*ux + ...
e(i,2).*uy).^2./c.^2 - 3/2.*(ux.^2 + uy.^2)./c.^2); % TODO check ^2
end
for i=1:q
f(:,:,i) = (1-omega)*f(:,:,i) + omega*feq(:,:,i) + omega*w(i)*(3*(e(i,1)*ux + e(i,2)*uy)./c + ...
9/2*(e(i,1)*ux + e(i,2)*uy).^2./c^2 - 3/2*(ux.^2 + uy.^2)./c^2);
end
% External forcing
% bc = [1,2,3,4]; % boundary nodes
% ux(bc,:) = bc_func(t); % TODO: What is this supposed to do?
% Streaming step
for i=1:q
f(:,:,i) = circshift(f(:,:,i), e(i,:));
end
% Compute macroscopic quantities
rho = sum(f,3);
ux = sum(f.*reshape(e(:,1),[1 1 q]),3)./rho;
uy = sum(f.*reshape(e(:,2),[1 1 q]),3)./rho;
% Update floater positions and velocities
for i=1:m
ix = max(1, min(floor(x(i)/dx) + 1, Lx));
iy = max(1, min(floor(y(i)/dx) + 1, Ly));
vx(i) = (ux(ix,iy) - vx(i));
vy(i) = (uy(ix,iy) - vy(i));
% ix = floor(x(i)/dx) + 1;
% iy = floor(y(i)/dx) + 1;
% vx(i) = (ux(ix,iy) - vx(i));
% vy(i) = (uy(ix,iy) - vy(i));
x(i) = x(i) + vx(i)*dt;
y(i) = y(i) + vy(i)*dt;
% Check for collisions with walls
if x(i) < R
x(i) = R;
vx(i) = -vx(i);
elseif x(i) + R > Lx
x(i) = Lx - R;
vx(i) = -vx(i);
end
if y(i) < R
y(i) = R;
vy(i) = -vy(i);
elseif y(i) + R > Ly
y(i) = Ly - R;
vy(i) = -vy(i);
end
% Check for collisions with other floaters
for j=1:m
if i ~= j && norm([x(i) y(i)] - [x(j) y(j)]) < 2*R % TODO überprüfen
theta = atan2(y(i)-y(j),x(i)-x(j));
v1 = [vx(i)*cos(theta) + vy(i)*sin(theta) vx(i)*-sin(theta) + vy(i)*cos(theta)];
v2 = [vx(j)*cos(theta) + vy(j)*sin(theta) vx(j)*-sin(theta) + vy(j)*cos(theta)];
v1f = (v1*(R-R) + 2*R*v2)./(R+R);
v2f = (v2*(R-R) + 2*R*v1)./(R+R);
vx(i) = v1f(1)*cos(theta) + v1f(2)*-sin(theta);
vy(i) = v1f(1)*sin(theta) + v1f(2)*cos(theta);
vx(j) = v2f(1)*cos(theta) + v2f(2)*-sin(theta);
vy(j) = v2f(1)*sin(theta) + v2f(2)*cos(theta);
x(i) = x(i) + vx(i)*dt;
y(i) = y(i) + vy(i)*dt;
x(j) = x(j) + vx(j)*dt;
y(j) = y(j) + vy(j)*dt;
end
end
end
% Apply boundary conditions
% Boundary condition for left wall
%% f(1,:,2) = f(2,:,2) - 2*w(2)*rho(1,:).*e(2,1)*bc_func(t); % Option A
f(1,:,2) = f(2,:,2) - 2*w(2)*f(1,:,2).*e(2,1)*bc_func(t); % Option B
f(1,:,2) = f(2,:,2) - 2*w(2)*f(1,:,2).*e(2,1)*bc_func(t);
f(1,:,5) = f(2,:,5) - 2*w(5)*f(1,:,5).*e(5,1)*bc_func(t);
f(1,:,6) = f(2,:,6) - 2*w(6)*f(1,:,6).*(e(6,1)*feval(bc_func,t)+ e(6,2)*bc_func(t));
f(1,:,7) = f(2,:,7) - 2*w(7)*f(1,:,7).*(-e(7,1)*feval(bc_func,t)+ e(7,2)*bc_func(t));
f(1,:,8) = f(2,:,8) - 2*w(8)*f(1,:,8).*(-e(8,1)*feval(bc_func,t)- e(8,2)*bc_func(t));
f(1,:,9) = f(2,:,9) - 2*w(9)*f(1,:,9).*(e(9,1)*feval(bc_func,t)- e(9,2)*bc_func(t));
% Boundary condition for right wall
f(Lx,:,4) = f(Lx-1,:,4) - 2*w(4)*f(Lx,:,4).*e(4,1)*bc_func(t);
f(Lx,:,8) = f(Lx-1,:,8) - 2*w(8)*f(Lx,:,8).*(-e(8,1)*feval(bc_func,t)- e(8,2)*bc_func(t));
f(Lx,:,9) = f(Lx-1,:,9) - 2*w(9)*f(Lx,:,9).*(e(9,1)*feval(bc_func,t)- e(9,2)*bc_func(t));
f(Lx,:,6) = f(Lx-1,:,6) - 2*w(6)*f(Lx,:,6).*(e(6,1)*feval(bc_func,t)+ e(6,2)*bc_func(t));
f(Lx,:,3) = f(Lx-1,:,3) - 2*w(3)*f(Lx,:,3).*(e(3,1)*feval(bc_func,t)+ e(3,2)*bc_func(t));
% Boundary condition for bottom wall
f(:,1,3) = f(:,2,3) - 2*w(3)*f(:,1,3).*e(3,2)*bc_func(t);
f(:,1,6) = f(:,2,6) - 2*w(6)*f(:,1,6).*(e(6,1)*feval(bc_func,t)+ e(6,2)*bc_func(t));
f(:,1,7) = f(:,2,7) - 2*w(7)*f(:,1,7).*(-e(7,1)*feval(bc_func,t)+ e(7,2)*bc_func(t));
f(:,1,2) = f(:,2,2) - 2*w(2)*f(:,1,2).*e(2,2)*bc_func(t);
f(:,1,9) = f(:,2,9) - 2*w(9)*f(:,1,9).*(e(9,1)*feval(bc_func,t)- e(9,2)*bc_func(t));
% Boundary condition for top wall
f(:,Ly,5) = f(:,Ly-1,5) - 2*w(5)*f(:,Ly,5).*e(5,2)*bc_func(t);
f(:,Ly,8) = f(:,Ly-1,8) - 2*w(8)*f(:,Ly,8).*(-e(8,1)*feval(bc_func,t)- e(8,2)*bc_func(t));
f(:,Ly,6) = f(:,Ly-1,6) - 2*w(6)*f(:,Ly,6).*(e(6,1)*feval(bc_func,t)+ e(6,2)*bc_func(t));
f(:,Ly,7) = f(:,Ly-1,7) - 2*w(7)*f(:,Ly,7).*(-e(7,1)*feval(bc_func,t)+ e(7,2)*bc_func(t));
f(:,Ly,1) = f(:,Ly-1,1) - 2*w(1)*f(:,Ly,1).*e(1,2)*bc_func(t);
f(:,Ly,9) = f(:,Ly-1,9) - 2*w(9)*f(:,Ly,9).*(e(9,1)*feval(bc_func,t)- e(9,2)*bc_func(t));
% Plot fluid and floaters
if mod(t,10) == 0
figure(1)
clf
subplot(2,2,1)
imagesc(ux')
colorbar
title('Fluid x-velocity')
subplot(2,2,2)
imagesc(uy')
colorbar
title('Fluid y-velocity')
subplot(2,2,3)
scatter(x,y,10,'filled')
axis equal
xlim([0 Lx])
ylim([0 Ly])
title('Floaters')
subplot(2,2,4)
quiver(ux',uy')
axis equal
title('Fluid velocity field')
drawnow
end
end