parfor-loop leads to wrong results of PDE solver

5 visualizaciones (últimos 30 días)
Anselm Dreher
Anselm Dreher el 7 de En. de 2022
Comentada: Anselm Dreher el 12 de En. de 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

Respuesta aceptada

Ravi Kumar
Ravi Kumar el 10 de En. de 2022
Editada: Edric Ellis el 11 de En. de 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 comentario
Anselm Dreher
Anselm Dreher el 12 de En. de 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.

Iniciar sesión para comentar.

Más respuestas (0)

Productos


Versión

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by