Parfor indexing for looping over edges and assigning to adjacent cell indices

I am working on a finite volume CFD code that operates on a triangular mesh. It loops over all edges between cells and adds or subtracts things from the values in the two adjacent cells to assemble the "flux" vector.

I am having a lot of trouble implementing this loop as a parfor loop because the loop variable is the edge number, but the assignment indices are the numbers of the neighboring cells. I've commented on the lines that don't work below.

flux = zeros(5,Nelem); % preallocate flux matrix
sdl = zeros(1,Nelem); % preallocate wavespeed matrix

% Loop over interior edges
parfor idxEdge = 1:NedgeInterior
    
    % Find L and R neighboring elements, their corresponding states, and the edge normal vector
    % and length.
    elemL = elemInterior(idxEdge,1); % cell number of left neighboring cell
    UL = U(:,elemL);
    elemR = elemInterior(idxEdge,2); % cell number of right neighboring cell
    UR = U(:,elemR);
    nEdge = normInterior(idxEdge,:);
    lEdge = lengthInterior(idxEdge);
    
    % Calculate flux over the edge
    [F,smag] = calcFlux(UL,UR,nEdge);
    
    % Add edge flux contribution to element-wise net flux matrix
    flux(:,elemL) = flux(:,elemL) + F.*lEdge; % doesn't work, Error: The variable flux in a parfor cannot be classified.
    flux(:,elemR) = flux(:,elemR) - F.*lEdge; % doesn't work, Error: The variable flux in a parfor cannot be classified.
    
    % Add edge wave speed contribution to element-wise wave speed vector
    sdl(elemL) = sdl(elemL) + smag.*lEdge; % doesn't work, Error: The variable sdl in a parfor cannot be classified.
    sdl(elemR) = sdl(elemR) + smag.*lEdge; % doesn't work, Error: The variable sdl in a parfor cannot be classified.
    
end

I know parfor doesn't like assignment indices that aren't loop indices, but I see no other way to structure the code. The loop index must be the edge number, but the assignment incides must be the neighboring cell numbers. Is there any workaround I could try to make this work?

Answer 1

[Edric Ellis]

I think this can be made to work, but it might be a little tricky. parfor is complaining because your output variables flux and sdl are not being correctly "sliced" - simply put, this means that it must be obvious to the parfor machinery that you've got independent loop iterations, and the only way that works is if the indexed assignments have the raw loop index in one position. More in the doc.

All is not lost though - because you're only ever incrementing flux and sdl, it might work out to make these be "reduction" variables. It's not a totally trivial change to make, and it's not clear that it will result in a parfor loop with good performance (depends how expensive calcFlux is, and how large flux is...).

The basic idea is this: on each loop iteration, you compute an increment to flux that is the same size as flux, and then say flux = flux + fluxIncrement. Something a bit like this (note this is completely untested since your original code isn't executable...)

parfor ...
    % ....        

    % Add edge flux contribution to element-wise net flux matrix
    fluxIncrement = zeros(5,Nelem); % fluxIncrement is a parfor "temporary" variable
    fluxIncrement(:,elemL) = F.*lEdge;
    fluxIncrement(:,elemR) = - F.*lEdge;
    flux = flux + fluxIncrement; % flux is now a parfor "reduction" variable.
    
    % ... and the equivalent for 'sdl'.
end

As an alternative, you could do something more like this where you compute all the increments in the parfor loop, store them in a cell array, and then apply them back on the client

fluxIncrement = cell(Nelem, 2);
parfor ...
    % store the increments to the left and right
    fluxIncrement(idx, :) = {F.*lEdge, -F.*lEdge};
end
% and then
for idx = 1:NedgeInterior
    % ... compute elemL and elemR
    flux(:, elemL) = flux(:, elemL) + fluxIncrement{idx, 1};
    flux(:, elemR) = flux(:, elemR) + fluxIncrement{idx, 2};
end