![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/858990/image.png)
parfor-loop leads to wrong results of PDE solver
17 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
采纳的回答
Ravi Kumar
2022-1-10
编辑:Edric Ellis
2022-1-11
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/858990/image.png)
Regards,
Ravi
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Boundary Conditions 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!