parfor-loop leads to wrong results of PDE solver

4 Ansichten (letzte 30 Tage)
Anselm Dreher
Anselm Dreher am 7 Jan. 2022
Kommentiert: Anselm Dreher am 12 Jan. 2022
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
% 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
% generate mesh
% pdemesh(model);
% assign coefficients of the PDE
% for coefficient f pass the additional parameter v for fluid velocity
'f',@(location,state)f_coeff (location,state,v,kinParam));
% define boundary conditions, Edge 4 is inlet
applyBoundaryCondition(model,'neumann','Edge',[1 2 3]);
% define initial conditions
% solve PDE
parfor i = 1:6 %%% for correct results, replace parfor by for
results(i) = solvepde(model);
xq = linspace(0,L,100);
yq = H/2*ones(1,length(xq));
for i = 1:length(results)
CAinterp(i,:) = interpolateSolution(results(i),xq,yq);
for i = 1:length(results)
hold on
% 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;

Akzeptierte Antwort

Ravi Kumar
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:
This should bypass the bug. I obtained the following solution of the suggested change:
  1 Kommentar
Anselm Dreher
Anselm Dreher am 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);
This also works in a parfor-loop.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)




Community Treasure Hunt

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

Start Hunting!

Translated by