MATLAB Answers

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

2 views (last 30 days)
Adam Wasserman
Adam Wasserman on 19 Mar 2020
Answered: Edric Ellis on 20 Mar 2020
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?

  0 Comments

Sign in to comment.

Accepted Answer

Edric Ellis
Edric Ellis on 20 Mar 2020
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

  0 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by