Unable to classify the variable 'feq' in the body of the parfor-loop

Hi everybody,

I have a code that has the following error: Error: Unable to classify the variable 'feq' in the body of the parfor-loop. I don't understand the issue and how to solve it. Thank you very much for helping me.

% 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 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
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
% delete(gcp('nocreate')); % delete any existing parallel pool
% parpool; % start a parallel pool from the Parallel Computing Toolbox
% parpool('local', 12); % 12 is the number of workers (= number of cores)
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 = 1;
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.
% Main simulation loop
% parfor t=0;dt;tmax
% feq = zeros(1,q);
%     f = zeros(1,q);
%     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
parfor t=0:dt:tmax
    
    % 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)*feval(bc_func,t); % Option A
    f(1,:,2) = f(2,:,2) - 2*w(2)*f(1,:,2).*e(2,1)*feval(bc_func,t); % Option B
    f(1,:,2) = f(2,:,2) - 2*w(2)*f(1,:,2).*e(2,1)*feval(bc_func,t);
    f(1,:,5) = f(2,:,5) - 2*w(5)*f(1,:,5).*e(5,1)*feval(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)*feval(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)*feval(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)*feval(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)*feval(bc_func,t));
    % Boundary condition for right wall
    f(Lx,:,4) = f(Lx-1,:,4) - 2*w(4)*f(Lx,:,4).*e(4,1)*feval(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)*feval(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)*feval(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)*feval(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)*feval(bc_func,t));
    % Boundary condition for bottom wall
    f(:,1,3) = f(:,2,3) - 2*w(3)*f(:,1,3).*e(3,2)*feval(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)*feval(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)*feval(bc_func,t));
    f(:,1,2) = f(:,2,2) - 2*w(2)*f(:,1,2).*e(2,2)*feval(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)*feval(bc_func,t));
    % Boundary condition for top wall
    f(:,Ly,5) = f(:,Ly-1,5) - 2*w(5)*f(:,Ly,5).*e(5,2)*feval(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)*feval(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)*feval(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)*feval(bc_func,t));
    f(:,Ly,1) = f(:,Ly-1,1) - 2*w(1)*f(:,Ly,1).*e(1,2)*feval(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)*feval(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

Answer 1

[Thierry Rebetez]

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

[Walter Roberson]

paralell workers and background workers can never (directly) display graphics. They are not connected to the user's display. They can construct graphs, but unless the graphs are captured or the graphics objects are saved to file or copied to the managing thread, the graphics cannot be visible to the user.

... And if they could be displayed, you would have trouble with the fact that multiple workers might be trying to write to the same figures.

Answer 2

[Walter Roberson]

initializing feq before the parfor loop implies that you expect feq to be an output, that you expect to be able to access feq after the parfor. However you do not index feq by the loop variable t

The implication is that you want feq to be shared between all of the workers, so that if one worker happens to write to feq(:, :, 3) and a different worker with a different t happens to write to there as well before the first worker happens to read back from feq(:, :, 3) that you want the revised values to be read back. The implication is that you want the results to depend upon continual races between the workers to be the most recent worker to write to the array.

parfor does not support that kind of deliberate race between workers.

If you did not intend that, if you intend that feq should be local to each worker, then initialize feq on the worker rather than before the parfor.