Communicate Between Workers in SPMD Computations

Problem Overview

The 2D heat equation models transient heat conduction and is given by:

where:

You can solve the 2D heat equation by discretizing the $$x$$ and $$y$$ spatial domain using the finite difference method and applying Dirichlet boundary condition, where the temperature at the boundary is kept fixed.

To parallelize the computation, you partition the spatial domain among the workers along the $x$-direction. Each worker computes updates for its subdomain and workers exchange the boundary columns with their immediate neighbors at each time step. This communication allows each worker to update its boundary points using the most recent information from adjacent subdomains.

Define Parameters

Define the physical and numerical parameters for the simulation.

L = 1;           % Length of the plate (m)
H = 1;           % Height of the plate (m)
tmax = 2000;        % Total simulation time (s)
alpha = 1.11e-4;   % Thermal diffusivity for copper (m^2/s)
dx = 0.0025;       % Grid spacing in x-direction (m)
dy = 0.0025;       % Grid spacing in y-direction (m)
dt = 0.01;         % Time step size (s)
nt = tmax/dt;

Check that the values satisfy the stability condition required for the explicit method, otherwise, adjust dt, dx, or dy to maintain stability

% Stability parameters
rX = alpha*dt/dx^2;
rY = alpha*dt/dy^2;

% Assert stability condition for explicit scheme
assert(rX + rY < 0.5, ...
    "Stability condition violated: decrease dt or increase dx/dy.");

Set Up Grid and Boundary Conditions

In this example, you model a copper plate of size 1m by 1m. The first step is to discretize it using a 2-D grid. Define the number of grid points in each dimension, X and Y, from the predefined grid spacing.

nx = round(L/dx)+1;
ny = round(H/dy)+1;

Discretize the spatial dimensions with a 2-D grid by using the meshgrid function. Divide each dimension uniformly according to the number of points by using linspace.

[X,Y] = meshgrid(linspace(0,L,nx),linspace(0,H,ny));

Define the initial interior and boundary temperatures. Suppose the copper plate is resting on a heat source. The interior of the plate starts at 10K. The left and right boundaries have a constant temperature of 25K, the top boundary has a constant temperature of 0K, and the bottom boundary has a constant temperature of 50K. These boundary conditions specify the constant values that a solution must take along the boundary of the domain. This type of boundary condition is known as the Dirichlet boundary condition.

u0 = ones(size(X))*10;
u0(:,1) = 25;
u0(:,end) = 25;
u0(1,:) = 50;
u0(end,:) = 0;

Plot the initial temperature distribution of the plate when t is 0.

figure;
plotTemperature(X,Y,u0,0);

Partition Grid on Workers

Start a parallel pool of workers using a remote cluster profile. To try this example, replace the cluster profile name with your own cluster profile or use the local Process profile instead.

pool = parpool("myCluster",10);

Starting parallel pool (parpool) using the 'myCluster' profile ...
Connected to parallel pool with 10 workers.

Use an spmd block to manually partition the temperature variable u across the memory of the workers in your cluster. Each worker get a contiguous block of columns. You can also use codistributed arrays to distribute the temperature variable to the workers.

spmd
    colsPerWkr = floor(ny/spmdSize);
    extra = mod(ny,spmdSize);
    if spmdIndex <= extra
        localCols = colsPerWkr+1;
        startCol = (spmdIndex-1)*localCols + 1;
    else
        localCols = colsPerWkr;
        startCol = (extra*(colsPerWkr+1)) ...
            + (spmdIndex-extra-1)*colsPerWkr + 1;
    end
    endCol = startCol + localCols - 1;

    uLocal = u0(:,startCol:endCol);
end

Exchange Initial Worker Boundaries and Set Up Ghost Columns

Before starting the time-stepping loop, each worker must update its local temperature array, uLocal, to include ghost columns from its immediate neighbors. These ghost columns are essential because the finite difference method requires values from adjacent subdomains to correctly update boundary points. In this section, each worker exchanges its boundary columns with its left and right neighbors using spmdSendReceive. The workers then add the received columns to uLocal as ghost columns on the appropriate side.

When exchanging boundary columns between neighbors, using a symmetric communication pattern such as spmdSendReceive simplifies the code and distributes communication more evenly across workers. Define the index of the left and right neighbor workers for each worker. On each worker, the spmdIndex-1 worker is the left worker, and the spmdIndex+1 worker is the right worker. When you use modulo division, worker 1 receives data from the worker whose index is equal to spmdSize. Edge workers can simply not add the ghost columns, allowing all workers to follow the same logic without special‑case branching.

spmd
    offset = 1;
    leftWkr = mod((spmdIndex-offset)-1,spmdSize) + 1;
    rightWkr = mod((spmdIndex+offset)-1,spmdSize) + 1;
end

spmd

Extract local boundary columns.

    leftSend  = uLocal(:,1);
    rightSend = uLocal(:,end);

Send left boundary and receive right ghost. For the rightmost worker, do not append the right ghost column.

    rightRecv = spmdSendReceive(leftWkr,rightWkr,leftSend);
    if spmdIndex < spmdSize
        uLocal = [uLocal,rightRecv];
    end

Send right boundary and receive left ghost. For the leftmost worker, do not append the left ghost column.

    leftRecv = spmdSendReceive(rightWkr,leftWkr,rightSend);
    if spmdIndex > 1
        uLocal = [leftRecv,uLocal];
    end

end

Perform Time-Stepping Loop with Boundary Exchange

Determine the size of the local temperature array. This value is different on workers with additional ghost boundaries and you use this value when you update the local temperature values.

spmd
    [lx,ly] = size(uLocal);

Additionally, each worker records snapshots of its local temperature array at the quarter, half, three-quarters, and final time steps. Determine which time steps to collect snapshots and preallocate a cell array to store the temperature data.

    snapSteps = round([0.25,0.5,0.75,1]*nt);
    numSnaps = numel(snapSteps);
    uSnapshots = cell(1,numSnaps);
end

Perform the main time-stepping loop for the heat equation. At each time step, each worker updates its local temperature values using the finite difference method. After updating, workers exchange their boundary columns with their neighbors to keep ghost columns up to date for the next iteration. If the time step is a snapshot time step, each worker stores the current local temperature values in a cell array.

spmd
    snapIdx = 1;
    for t = 1:nt
        uNewLocal = uLocal;
        % Vectorized finite-difference update on the interior
        uNewLocal(2:lx-1,2:ly-1) = uLocal(2:lx-1,2:ly-1) ...
            + rX*(uLocal(2:lx-1,3:ly)   ...       % right
            - 2*uLocal(2:lx-1,2:ly-1) ...
            + uLocal(2:lx-1,1:ly-2)) ...          % left
            + rY*(uLocal(3:lx,2:ly-1) ...         % down
            - 2*uLocal(2:lx-1,2:ly-1) ...
            + uLocal(1:lx-2,2:ly-1));           % up

        uLocal = uNewLocal;

        leftSend  = uLocal(:,2);       % inner-left updated boundary
        rightSend = uLocal(:,end-1);   % inner-right updated boundary

        % Send updated left column and receive right ghost
        rightRecv = spmdSendReceive(leftWkr,rightWkr,leftSend);
        if spmdIndex < spmdSize
            uLocal(:,end) = rightRecv;
        end

        leftRecv = spmdSendReceive(rightWkr,leftWkr,rightSend);
        if spmdIndex > 1
            uLocal(:,1) = leftRecv;
        end

        % Collect snapshot if this is a designated time step
        if snapIdx <= numSnaps && t == snapSteps(snapIdx)
            uSnapshots{snapIdx} = uLocal;
            snapIdx = snapIdx+1;
        end
    end
end

Reconstruct Global Solution

After the time-stepping loop, each worker holds a portion of the temperature snapshot data that includes extra ghost columns at the subdomain boundaries. To reconstruct the full global solution, remove the ghost columns from each worker.

spmd
    if spmdIndex > 1
        for idx = 1:numSnaps
            uSnapshots{idx} = uSnapshots{idx}(:,2:end);
        end
    end

    if spmdIndex < spmdSize
        for idx = 1:numSnaps
            uSnapshots{idx} = uSnapshots{idx}(:,1:end-1);
        end
    end
end

Concatenate the remaining snapshot data across the workers in the correct column order using the spmdCat function and store the resulting data on worker 1.

spmd
    globalSnapshots = cell(1,numSnaps);
    for s = 1:numSnaps
        globalSnapshots{s} = spmdCat(uSnapshots{s},2,1);
    end
end

Retrieve the global snapshot results from worker 1.

snapshotTimes = snapSteps{1}*dt;
finalSnapshots = globalSnapshots{1};

Visualize Results

Plot the temperature distributions at the quarter, half, three-quarters, and final time steps.

figure;
tiledlayout(2,2);
for idx = 1:4
    nexttile
    plotTemperature(X,Y,finalSnapshots{idx},snapshotTimes(idx));
end

function plotTemperature(X,Y,u,t)
contourf(X,Y,u,20,LineColor="none");
title("Time = "+num2str(t)+"s");
meanT = mean(u,"all");
subtitle(compose("AvgTemp = %.2fK",meanT))
xlabel("x (m)");
ylabel("y (m)");
axis equal tight;
c = colorbar;
c.Label.String = "Temperature (K)";
end