parfor-loop leads to wrong results of PDE solver

3 views (last 30 days)
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

Accepted Answer

Ravi Kumar
Ravi Kumar on 10 Jan 2022
Edited: Edric Ellis on 11 Jan 2022
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
  1 Comment
Anselm Dreher
Anselm Dreher on 12 Jan 2022
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.

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