![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/858990/image.png)
parfor-loop leads to wrong results of PDE solver
17 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Anselm Dreher
am 7 Jan. 2022
Kommentiert: Anselm Dreher
am 12 Jan. 2022
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 Kommentare
Akzeptierte Antwort
Ravi Kumar
am 10 Jan. 2022
Bearbeitet: Edric Ellis
am 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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/858990/image.png)
Regards,
Ravi
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Boundary Conditions finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!