parfor-loop leads to wrong results of PDE solver

Hello,

I am currently using the PDE toolbox to simulate a convection-diffusion problem. To simulate multiple cases I want to speed up the process by using a parfor-loop executing the "solvepde"-function. Unfortunately the results change dramatically when using the parfor-loop instead of the normal for-loop. You can try it by replacing the parfor-command with the for-command.

The following programm repeats the same calculation 6 times. Coefficient c is a diffusion coefficient, coefficient f contains a convection and a sink term. The geometry represents a 1D flow in a pipe with Dirichlet BC at the inlet and Neumann=0 BC at all other edges. The result should be a nice smooth decrease in concentration along the x-Axis, just as the for-loop calculates it. Instead, the parfor-loop results are nearly zero.

My guess is that the iterations arent fully independent, but even after reading the documentation pages about parfor and some forum posts I have no idea where this could come from. Do you have any ideas?

clearvars, close all

clc

% Create PDE model with 1 equation

model = createpde(1);

model.SolverOptions.ReportStatistics = 'on';

model.SolverOptions.ResidualTolerance = 2e-3;

L = 1; % length

H = 1e-1; % arbitrary height

D = 1e-5; % 1e-5 % Diffusion coefficient m^2/s

v = 3e-3; % 3e-3 % velocity m/s

Cfeed = 0.1; % 0.1 % Feed concentration mol/m^3

C0 = Cfeed * 2e-1; % Cfeed * 2e-1 % Initial condition mol/m^3;

kinParam = 2e-1; % kinetic parameter

hmax = H; % H % Size of FEM element

% define geometry for fluid flow

F_Rct = [3 4 0 L L 0 0 0 H H]'; %Description Matrix

gd = F_Rct;

sf = 'F1';

ns = ('F1')';

geo = decsg(gd,sf,ns);

% assign geometry to model

geometryFromEdges(model,geo);

% generate mesh

generateMesh(model,'Hmax',hmax);

% pdemesh(model);

% assign coefficients of the PDE

% for coefficient f pass the additional parameter v for fluid velocity

specifyCoefficients(model,'m',0,...

'd',0,...

'c',D,...

'a',0,...

'f',@(location,state)f_coeff (location,state,v,kinParam));

% define boundary conditions, Edge 4 is inlet

applyBoundaryCondition(model,'dirichlet','Edge',4,'u',Cfeed);

applyBoundaryCondition(model,'neumann','Edge',[1 2 3]);

% define initial conditions

setInitialConditions(model,C0);

% solve PDE

parpool(6)

parfor i = 1:6 %%% for correct results, replace parfor by for

results(i) = solvepde(model);

end

delete(gcp('nocreate'));

xq = linspace(0,L,100);

yq = H/2*ones(1,length(xq));

for i = 1:length(results)

CAinterp(i,:) = interpolateSolution(results(i),xq,yq);

end

figure(3)

for i = 1:length(results)

hold on

plot(xq,CAinterp(i,:))

end

% ylim([0 Cfeed])

hold off

%% Function defining coefficient f

function f = f_coeff(~,state,v,kinParam)

% convective term state.ux devided by arbitrary factor state.u

f = -v * state.ux ./ state.u - kinParam * state.u;

end

Answer 1

[Ravi Kumar]

Hi Anselm,

This is a bug when the model gets transferred to the worker in the process of parfor execution. I apologize for the inconvenience. Here is a workaround:

Chage the line:

applyBoundaryCondition(model,'dirichlet','Edge',4,'u',Cfeed);

to

applyBoundaryCondition(model,'dirichlet','Edge',4,'r',Cfeed);

This should bypass the bug. I obtained the following solution of the suggested change:

Regards,

Ravi

[Anselm Dreher]

Thank you, this works exactly as intended.

Since I need to use a mixed boundary condition in my actual code, I'd like to add, that they can be defined like its shown in the documentation using h,r,q and g. Just switching u for r doesnt work!

Assume I have N=10 equations that get a Dirchlet BC and a last one gets a Neumann BC, so there are 11 in total. The Dirichlet BCs are defined in “BCDirichlet”, for simplicity just the numbers 1 to 10. The matrix h is an identity matrix for all Dirichlet BC except for the last one with a Neumann BC for which it is set 0. Since the Neumann BC for Equation M should equal zero, Q and G are just filled with zeros.

N = 10;

M = N+1;

H = eye(M,M);

H(M,M) = 0;

BCDirichlet = (1:10);

BCDirichlet(M) = 0;

Q = zeros(M,M);

G = zeros(M,1);

applyBoundaryCondition(model,'mixed','Edge',4,'h',H,'r',BCDirichlet(1:M),'q',Q,'g',G);

This also works in a parfor-loop.