Pdepe cylindrical coordinates with Neumann BC

I try to solve 1D diffusion eqn with Neumann(flux) BC ar r=R_c. I couldn't use variable Conc in the BC function. It gives 'Unrecognized function or variable' error. I appreciate for any ideas on how to resolve.

Example input: AdvDiff_Cyl(4E-11,4.2E-10,0.4,[0 6; 1*3600 5; 6*3600 3.75; 20*3600 1.85; 48*3600 0.7; 95*3600 0.3]);

function Conc = AdvDiff_Cyl(Deff,P,Kav,NP1)

% Define the parameters
R_C = 19.2*10^-6/2/pi; % radius of microvessel [m]
R_O = 10*R_C;          % radius of surrounding outer tissue [m]
tsim = 24 * 3600       % simulation time in [s]
r = [R_C:1E-6:R_O];    % create spatial and temporal solution mesh
t = [0:5:tsim];

% Solve the PDE numerically using pdepe which solves the pde using ode15s
m = 1   ; %1D cylindrical model
[sol] = pdepe(m, @advdiff_pde, @advdiff_initial, @advdiff_bc,r,t);
Conc = sol(:,:,1);

colors =  [    0      0.4470    0.7410;
            0.8500    0.3250    0.0980;
            0.9290    0.6940    0.1250;
            0.4940    0.1840    0.5560;
            0.4660    0.6740    0.1880;
            0.3010    0.7450    0.9330;
            0.6350    0.0780    0.1840];

figure()
i = 1;

for time = [5/3600 1 2 5 10 20]*3600/5
plot(r/R_C,Conc(time,:),'-', 'color', colors(i,:),'LineWidth',2)
hold on
i = i + 1;
end
xlabel('Distance [r/R_N]', 'Interpreter','Tex');
ylabel('Concentration [\mug/ml]','Interpreter','Tex');
title('Numerical Solution of ADE using pdepe','Interpreter','Tex');
legend('Initial Condition', '1 hr', '2 hr','5 hr','10 hr','20 hr','Location','northeast')
lgd.FontSize = 6;
legend('boxoff')


% Define the DE
function [c,f,s] = advdiff_pde(r,t,Conc,dcdr)
c = 1;
f = Deff*dcdr;
s = 0
end

% Define the initial condition
function c0 = advdiff_initial(r)
    c0 = 0;
end

% Define the boundary condition
function [pl,ql,pr,qr] = advdiff_bc(xl,cl,xr,cr,t)
    % BC1: 
    pl = P * (interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap') - (Conc/Kav));
    ql = Deff; %ignored by solver since m=1

    % BC2: 
    pr = 0;
    qr = 1;
end
end

Answer 1

[Torsten]

You will have to use

    pl = P * (interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap') - (cl/Kav));
    ql = 1; %ignored by solver since m=1    

instead of

    pl = P * (interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap') - (Conc/Kav));
    ql = Deff; %ignored by solver since m=1

[Ali]

Great! thank you.

Does cl refers to the value at the left boundary? It fixed the error, but there is one thing I am not sure. I attached the solution for 20 hrs. After some point since the BC value at the left decreases, I expect negative slope at the interface. I believe when the slope is negative, the BC value on the left end supposed to be smaller than the value on the right of the interface (referring to the green line at t=20hrs I would expect a smaller value for the horizontal line). Do you have any idea on that?

%Sample Input: AdvDiff_Cyl(4E-14,4.2E-8,0.4,[0 6; 1*3600 5; 6*3600  3.75; 20*3600  1.85; 48*3600  0.7; 95*3600  0.3]);

function Conc = AdvDiff_Cyl(Deff,P,Kav,NP1)

% Define the parameters
R_C = 19.2*10^-6/2/pi; % radius of microvessel [m]
R_O = 10*R_C;          % radius of surrounding outer tissue [m]
tsim = 24 * 3600       % simulation time in [s]
r = [R_C:1E-7:R_O];    % create spatial and temporal solution mesh
t = [0:1:tsim];

% Solve the PDE numerically using pdepe which solves the pde using ode15s
m = 1   ; %1D cylindrical model
[sol] = pdepe(m, @advdiff_pde, @advdiff_initial, @advdiff_bc,r,t);
Conc = sol(:,:,1);

colors =  [    0      0.4470    0.7410;
            0.8500    0.3250    0.0980;
            0.9290    0.6940    0.1250;
            0.4940    0.1840    0.5560;
            0.4660    0.6740    0.1880;
            0.3010    0.7450    0.9330;
            0.6350    0.0780    0.1840];

figure()
i = 1;
% time = [1 2 5 10 20]*3600/5
for time = [1 2 5 10 20]*3600
h1(i) = plot(r/R_C,Conc(time,:),'-', 'color', colors(i,:),'LineWidth',2);
hold on
Conc_vessel= interp1(NP1(:,1), NP1(:,2), time, 'spline', 'extrap');
plot([0 1], [Conc_vessel Conc_vessel],'-', 'color', colors(i,:),'LineWidth',2);
hold on
i = i + 1;
end
rectangle('Position',[0.99 0 0.01 5])
xlabel('Distance [r/R_N]', 'Interpreter','Tex');
ylabel('Concentration [\mug/ml]','Interpreter','Tex');
title('Numerical Solution of ADE using pdepe','Interpreter','Tex');
legend([h1(1), h1(2), h1(3), h1(4), h1(5)], '1 hr', '2 hr','5 hr','10 hr','20 hr','Location','northeast')
lgd.FontSize = 6;
legend('boxoff')


% Define the DE
function [c,f,s] = advdiff_pde(r,t,Conc,dcdr)
c = 1;
f = Deff*dcdr;
s = 0;
end

% Define the initial condition
function c0 = advdiff_initial(r)
    c0 = 0;
end

% Define the boundary condition
function [pl,ql,pr,qr] = advdiff_bc(xl,cl,xr,cr,t)
    % BC1: Transvascular Transport at R = R_c
    pl = P * (interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap') - (cl/Kav));
    ql = 1; %ignored by solver since m=1    

    % BC2: Concentration gradient is zero at the oter edge 
    pr = 0;
    qr = 1;
end
end

[Torsten]

Does cl refers to the value at the left boundary?

Yes.

You should plot

Kav*interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap')

, not

interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap')

at the left interface to compare with Conc.

[Ali]

@Torsten Thank you!

Sorry for bothering you but I have one more problem that I couldn't adress.

On the surface plot below the shorter surface is the BC at r=R_c and the longer surface is the values inside the domain. It seems the values inside the domain are starting to be greater than the BC values after some time, which normally should be imposible. Do you think I am formulating sth wrong? For some other input parameters this error/difference between BC and domain even greater.

AdvDiff_Cyl(1E-11,4.2E-8,1.0,0,[0 1; 1*3600 1; 6*3600  1; 20*3600  1; 48*3600  1; 95*3600  1]);

Index in position 1 exceeds array bounds. Index must not exceed 51841.

Error in solution>AdvDiff_Cyl (line 39)
h1(i) = plot(r/R_C,Conc(time,:),'-', 'color', colors(i,:),'LineWidth',2);

function Conc = AdvDiff_Cyl(Deff,P,Kav,Phi,NP1)

% Define the parameters
R_C = 19.2*10^-6/2/pi; 
R_O = 50*R_C;          
tsim = 24 * 3 * 3600;  % simulation time in [s]
r = [R_C:1E-7:R_O];    % create spatial and temporal solution mesh
t = [0:5:tsim];

% Solve the PDE numerically using pdepe which solves the pde using ode15s
m = 1   ; %2-D Cylindrical coordinates with azimuthal angular symmetry
[sol] = pdepe(m, @advdiff_pde, @advdiff_initial, @advdiff_bc,r,t);
Conc = sol(:,:,1);

colors =  [    0      0.4470    0.7410;
            0.8500    0.3250    0.0980;
            0.9290    0.6940    0.1250;
            0.4940    0.1840    0.5560;
            0.4660    0.6740    0.1880;
            0.3010    0.7450    0.9330;
            0.6350    0.0780    0.1840];
figure() %2D Surface Kymograph Plot
% ZZ = Kav*interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap'); %Constant or Specified BC
ZZ = Kav *(1+sin(2*pi/24/3600*(t+Phi*3600)))/2;   %Sin BC
surf(r(1:10:numel(r))/R_C,t(1:10:numel(t))/3600,Conc((1:10:numel(t)),(1:10:numel(r))),'edgecolor','none')
hold on
surf([0 1],t/3600, ZZ'.*[1 1],'edgecolor','none')
xlabel('Distance [r/R_N]'); xlim([0 25])
ylabel('Time [hr]'); ylim([0 tsim/3600]); yticks([0:24:tsim/3600])
colorbar; caxis([0, 1]);
view(2)


figure() %1D plot
i = 1;
% time = [1 2 5 10 20]*3600/5
for time = [1 2 4 8 12 24]*3600
h1(i) = plot(r/R_C,Conc(time,:),'-', 'color', colors(i,:),'LineWidth',2);
hold on
Conc_vessel = interp1(NP1(:,1), NP1(:,2), time, 'spline', 'extrap'); %Constant or Specified BC
% Conc_vessel = (1+sin(2*pi/24/3600*time))/2;   %Sin BC
plot([0 1], Kav*[Conc_vessel Conc_vessel],'-', 'color', colors(i,:),'LineWidth',2);
hold on
i = i + 1;
end
rectangle('Position',[0.8 0 0.4 5],'EdgeColor',"none",'FaceColor',[1 0 0 0.3])
ylim([0 1])
xlim([0 20])
xlabel('Distance [r/R_N]', 'Interpreter','Tex');
ylabel('Concentration [\mug/ml]','Interpreter','Tex');
title('Numerical Solution of ADE using pdepe','Interpreter','Tex');
legend([h1(1), h1(2), h1(3), h1(4), h1(5) h1(6)],'1 hr','2 hr','4 hr' ,'8 hr','12 hr','24 hr','Location','northeast')
lgd.FontSize = 6;
legend('boxoff')



% Define the DE
function [c,f,s] = advdiff_pde(r,t,Conc,dcdr)
c = 1;
f = Deff*dcdr;
s = 0;
end

% Define the initial condition
function c0 = advdiff_initial(r)
    c0 = 0.5;
end

% Define the boundary condition
function [pl,ql,pr,qr] = advdiff_bc(xl,cl,xr,cr,t)
    % BC1: Transvascular Transport at R = R_c
    pl = P * (1+sin(2*pi/24/3600*t)/2 - (cl/Kav)); %Sin BC
    ql = 1;   

    % BC2: Concentration gradient is zero at the outer edge 
    pr = 0;
    qr = 1;
end
end

[Torsten]

Please correct your code (see above).

[Ali]

Sorry for that, I deleted the part causing the error. The problem I adressed is shown in the surface generated it has values greater than 1 (the maximum BC value).

AdvDiff_Cyl(1E-11,4.2E-8,1.0,0,[0 1; 1*3600 1; 6*3600  1; 20*3600  1; 48*3600  1; 95*3600  1]);

function Conc = AdvDiff_Cyl(Deff,P,Kav,Phi,NP1)
% Define the parameters
R_C = 19.2*10^-6/2/pi; 
R_O = 50*R_C;          
tsim = 24 * 3 * 3600;  % simulation time in [s]
r = [R_C:1E-7:R_O];    % create spatial and temporal solution mesh
t = [0:5:tsim];
% Solve the PDE numerically using pdepe which solves the pde using ode15s
m = 1   ; %2-D Cylindrical coordinates with azimuthal angular symmetry
[sol] = pdepe(m, @advdiff_pde, @advdiff_initial, @advdiff_bc,r,t);
Conc = sol(:,:,1);
colors =  [    0      0.4470    0.7410;
            0.8500    0.3250    0.0980;
            0.9290    0.6940    0.1250;
            0.4940    0.1840    0.5560;
            0.4660    0.6740    0.1880;
            0.3010    0.7450    0.9330;
            0.6350    0.0780    0.1840];

figure() %2D Surface Kymograph Plot
% ZZ = Kav*interp1(NP1(:,1), NP1(:,2), t, 'spline', 'extrap'); %Constant or Specified BC
ZZ = Kav *(1+sin(2*pi/24/3600*(t+Phi*3600)))/2;   %Sin BC
surf(r(1:10:numel(r))/R_C,t(1:10:numel(t))/3600,Conc((1:10:numel(t)),(1:10:numel(r))),'edgecolor','none')
hold on
surf([0 1],t/3600, ZZ'.*[1 1],'edgecolor','none')
xlabel('Distance [r/R_N]'); xlim([0 25])
ylabel('Time [hr]'); ylim([0 tsim/3600]); yticks([0:24:tsim/3600])
colorbar; caxis([0, 1]);
view(2)

% Define the DE
function [c,f,s] = advdiff_pde(r,t,Conc,dcdr)
c = 1;
f = Deff*dcdr;
s = 0;
end
% Define the initial condition
function c0 = advdiff_initial(r)
    c0 = 0.5;
end
% Define the boundary condition
function [pl,ql,pr,qr] = advdiff_bc(xl,cl,xr,cr,t)
    % BC1: Transvascular Transport at R = R_c
    pl = P * (1+sin(2*pi/24/3600*t)/2 - (cl/Kav)); %Sin BC
    ql = 1;   
    % BC2: Concentration gradient is zero at the outer edge 
    pr = 0;
    qr = 1;
end
end

[Torsten]

If you add the following lines at the end of your code

figure(2)
for i = 1:10:numel(t)
  [~,dudx(i)] = pdeval(m,r,sol(i,:),r(1));
  bc_cond(i) = Deff*dudx(i) + P * (1+sin(2*pi/24/3600*t(i))/2 - sol(i,1)/Kav);
end
plot(t(1:10:numel(t)),bc_cond(1:10:numel(t)))

you will see that bc_cond is approximately 0.

Thus your boundary condition at r = R_C that you prescribe as

Deff*dc/dr + P * (1+sin(2*pi/24/3600*t)/2 - c/Kav);

is satisfied.

I don't know how you conclude from this boundary condition that c should not exceed 1.

Maybe you mean

pl = P * ((1+sin(2*pi/24/3600*t))/2 - (cl/Kav)); %Sin BC
ql = 1;      

instead of

pl = P * (1+sin(2*pi/24/3600*t)/2 - (cl/Kav)); %Sin BC
ql = 1;   

[Ali]

I should clarify my point furher. What I try to model is cyclic concentration BC (1+sin(2*pi/24/3600*t)/2 at r=R_c which is in between 0 and 1, diffusing into the domain r>R_c.

That's why physically I expect sinusoidal concentration BC in between 0 and 1 for (0<r<R_c) to be greater or equal to the concentration values inside the domain (r>R_c) after being diffused according to the pde in cylindrical coordinates.

Thank you for your help

[Torsten]

(1+sin(2*pi/24/3600*t)/2

is between 1/2 and 3/2.

That's why I corrected

pl = P * (1+sin(2*pi/24/3600*t)/2 - (cl/Kav)); %Sin BC
ql = 1;   

to

pl = P * ((1+sin(2*pi/24/3600*t))/2 - (cl/Kav)); %Sin BC
ql = 1;