diffusion from point source boundary conditions

I am using MATLAB's pdepe partial diffeq solver to solve diffusion of a solute. I have a point source at r=0 where there should be intake of solute (so a sink, or negative source). if i use the "s" variable for the solver, does this input it at r=0, or for every value of r? if the latter, then should i incorporate the source term in my BC somehow?

any help is appreciated. thanks.

Answer 1

[Torsten]

if i use the "s" variable for the solver, does this input it at r=0, or for every value of r?

For every value of r (at least if you don't restrict s to special r-values by an if-statement).

if the latter, then should i incorporate the source term in my BC somehow?

The boundary condition at r=0 must be dc/dr = 0 - all other boundary conditions will be ignored.

In my opinion, you have two options:

Define a sphere with a hole in the middle of a certain radius and set the sink term in the boundary conditions part of the code as a flux over the inner surface (mol/(m^2*s)) or

Define a volume source in the s-term (mol/m^3*s) and restrict it to a specified inner radius:

rmin = ...;
if r <= rmin
  s = ...
end 

In either case: These are difficult error-prone and numerically challenging settings and you should verify your results by making a global molar or mass balance.

[Juliana]

hi, thanks again for the help, Torsten.

I am viewing the spatial variable as distance from surface of the cell, so i just restricted the source term to r=0. i think that does the job!

for this system though its more accurate if the source term could update for every value of t, since the source term is dependent irl on surrounding concentration.

do you think this is possible to do with a for loop?

[Torsten]

I am viewing the spatial variable as distance from surface of the cell, so i just restricted the source term to r=0. i think that does the job!

I think to get the source/sink in mol/s, you have to multiply the value you prescribed by the volume of the cell from r=0 to r=r(2). But check it from the results you get.

for this system though its more accurate if the source term could update for every value of t, since the source term is dependent irl on surrounding concentration.

do you think this is possible to do with a for loop?

Not with a for loop. One input parameter to pdefun is the time t so that you can define the source term directly as a function of t.

But again my advice: make a global mass balance.

[Juliana]

okay i see. but the source term is not dependent on time, its dependent on concentration at r=0. like, the cell uptakes more nutrients if there are more available nutrients, as dictated by michaelis-mentin kinetics

[Juliana]

this is what i get without defining the source at radius=0

vs with defining the following:

if r==0

s=-1;

else

s=0;

im not sure the second one makes much sense intuitively, since it should show decreasing concentrations by r=0

also if i change the s=-1 to s=1, the figure does not change

[Torsten]

If you change rmin to 1/37, nothing happens in the code below.

The function "heatsphere" is not called with the boundary values for r (thus r=0 and r=1).

The minimum r value with which it is called is 1/36. I don't know how MATLAB comes up with this value. It's neither part of the grid vector r nor a cell center. Maybe because m=2 instead of m=0, MATLAB includes a certain propagation factor in the r-values.

Since you also don't know the volume of the cells the source term is connected with, it is nearly impossible to dose the source according to what you want to set.

I suggest to include the source as a boundary source at a value r=rmin > 0.

Then you have the area A = 4*pi*rmin^2 over which the substance diffuses and you can set the flux accordingly as D*dn/dr = ndot/A.

r = linspace(0,1,25);
t = linspace(0,1000,25);
m = 2;
sol = pdepe(m,@heatsphere,@heatic,@heatbc,r,t);
u = sol(:,:,1);
surf(r,t,u)
xlabel('r')
ylabel('t')
zlabel('u(r,t)')
view([150 25])

plot(t,sol(:,1))
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at center of sphere')

function [c,f,s] = heatsphere(r,t,u,dudr)
D = 1e-4;
c = 1;
f = D*dudr;
s = 0;
rmin = 1/36;
if r <= rmin
  s = -12;
end
end
%----------------------------------------------
function u0 = heatic(r)
u0 = 300;
end
%----------------------------------------------
function [pl,ql,pr,qr] = heatbc(rl,ul,rr,ur,t)
pl = 0.0;
ql = 1.0;
pr = ur - 300;
qr = 0;
end
%----------------------------------------------

[Juliana]

Thanks again, the value of 1/36 worked well.

i am attempting to update the source value for every few time values, since the pdepe needs at least 3 values of time to eval for

i am having trouble with the syntax of it. not sure where the loop should go, because right now i get an error that the variable s1 is not recognized, but if i were to put the loop within the function, then i cant do multiple time step vectors. my code is pasted below.

%prompt = "What is concentration of nutrient in medium? ";    %user inputs
medium = 100; %prompt;
%prompt = "What is diffusivity of nutrient in medium? ";
D = .0001; %prompt;
%prompt = "What is the maximum radial distance for the system? ";
rinput = .5; %prompt;
%prompt = "What is maximum time span for the system? ";
tinput = 100; %prompt;
r = linspace(0,rinput,1000);        %space vector
u=linspace(100,0,100);
tn=1;
tm=4;
     
 for tm=4:tinput
    t0=linspace(tn,tm,4);   %make vector with small t values to do steps of diffusion

    rmin=1/36;              %smallest value that matlab computes
    if r<=rmin              %specifies source at approxr=0
    s1=-2*u(tm,rmin);     %rbitraty function for source that depends on u
    else                    %no sink anywhere else
    s1=0;
    end
Unew=diffusio(medium,D,r,t0);      %new diffusion calculation each time
% u= find some way to store new values of matrix with previous iterations
     tn=tn+4;               %add 4 to time vector so we get next terms
     tm=tm+4;

    end     
    



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure (1)                          %surfplot
surf(r,t,u,'EdgeColor', 'none')     %so grid doesnt make image all black
colorbar
clim ([0 100])
xlabel('r')
ylabel('t')
zlabel('u(r,t)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure (2)
colormap cool           %aesthetic choice
imagesc(r,t,u)
colorbar 
clim ([0 100])           %fixed colorbar range
xlabel('Radial Distance r','interpreter','latex')
ylabel('Time t','interpreter','latex')
title('Diffusion','interpreter','latex')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure (3)
plot(t,u(:,1))      %not sure if this is showing at surface of cell?
xlabel('Time')
ylabel('Concentration')
title('Concentration change at center of sphere')

%%%%%%%%%%%%%%%%%% FUNCTIONS %%%%%%%%%%%%%%%%%%%%

function [u]=diffusio(medium,D,r,t)

%medium=100;               %for if we dont wanna input each time
%D=.0001;
%rinput=.5;
%tinput=500;

m = 2;                          %indicates spherical sym

sol = pdepe(m,@diffpde,@icfun,@bcfun,r,t);
u = sol(:,:,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [c,f,s] = diffpde(r,t,u,dudr)     %pde definition
c = 1;
f = D*dudr;                 %diffusivity times gradient
s=s1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u0 = icfun(r)     %initial condition
  u0 = medium;              %at t=0, all r, conc=medium
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [pl,ql,pr,qr] = bcfun(rl,ul,rr,ur,t)   %boundary conditions
pl = 0;                 %ignored bc spherical   
ql = 0;                 %ignored bc spherical   
pr = 0;                 %indicates no flux on boundary
qr = 1;                 %indicates no flux on boundary
end

end

[Torsten]

I don't understand why you make things that complicated.

You can directly set the sink term in the function itself depending on u.

You also have access to the time t in the function if you want to set something time-dependent.

medium = 100;               %for if we dont wanna input each time
D = .0001;
rend = .5;
nr = 1000;
r = linspace(0,rend,nr);
rmin = 1/36;

tend = 500;
nt = 100;
t = linspace(0,tend,100);

u = diffusio(medium,D,rmin,r,t);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure (1)                          %surfplot
surf(r,t,u,'EdgeColor', 'none')     %so grid doesnt make image all black
colorbar
clim ([0 100])
xlabel('r')
ylabel('t')
zlabel('u(r,t)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure (2)
colormap cool           %aesthetic choice
imagesc(r,t,u)
colorbar 
clim ([0 100])           %fixed colorbar range
xlabel('Radial Distance r','interpreter','latex')
ylabel('Time t','interpreter','latex')
title('Diffusion','interpreter','latex')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure (3)
plot(t,u(:,1))      %not sure if this is showing at surface of cell?
xlabel('Time')
ylabel('Concentration')
title('Concentration change at center of sphere')

function u = diffusio(medium,D,rmin,r,t);
  m = 2;                          %indicates spherical sym
  sol = pdepe(m,@(r,t,u,dudr)diffpde(r,t,u,dudr,D,rmin),@(r)icfun(r,medium),@bcfun,r,t);
  u = sol(:,:,1);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
function [c,f,s] = diffpde(r,t,u,dudr,D,rmin)     %pde definition
  c = 1;
  f = D*dudr;                 %diffusivity times gradient
  s = 0.0;
  if r <=rmin
    s = -2e-1*u(1);
  end   
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u0 = icfun(r,medium)     %initial condition
  u0 = medium;                      %at t=0, all r, conc=medium
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [pl,ql,pr,qr] = bcfun(rl,ul,rr,ur,t)   %boundary conditions
  pl = 0;                 %ignored bc spherical   
  ql = 1;                 %ignored bc spherical   
  pr = 0;                 %indicates no flux on boundary
  qr = 1;                 %indicates no flux on boundary
end

[Juliana]

let me explain a bit better:

i am monitoring diffusion of, say, glucose into a cell.

the cell's rate of intake (source term of pde) depends on the surrounding concentration (so, concentration at r=rmin)

therefore, for each new value of concentration at each time t, the source term will change dependent on the new concentration.

i am trying to find a way to solve the pde multiple times over and over, updating the source term for each new concentration value.

so say i want to solve the pde from t=linspace(0,10,1000) so delta t = 0.01 s

so for t=0, i will get concentration=c1 at rmin

so for t=.01, i want to take c1 and input it into my source term function, lets say s=-20*(concentration at r=rmin). so now s=-20*c1. i will get concentration=c2

then for t=0.02, i want the source to be s=-20*c2

and so on.

but i cannot do this in a number of ways for a number of reasons. I am first trying to solve the pde at each time t, but the pde needs a time span of at least three values.

so now i will try and solve the pde in the script with a for loop, creating a time span of 3 values per my delta t of 0.01. and run the pde solver over and over.

but this does not work because the syntax is wrong. i cannot update the source term outside the function that defines it. and i cannot update the source term within the dunction because it does not recognize the output matrix, obviously.

how can I solve this?

[Torsten]

But what you try to do is fully equivalent to setting the source term to s = -20*u(1) for 0<=r<=rmin. u(1) is the (varying) concentration, and the source term is changing with changing concentration according to your above explained law.

Setting s to -20*u is like setting a first-order reaction.

[Juliana]

i updated the source to be -20*u(:,1) and I seem to get a result but I'm not convinced this isnt just taking u(last value,1). it just doesnt look right to me

[Torsten]

Setting s = -20*u(1) means that a 1st order reaction takes places that consumes u at a rate -20*u.

I don't understand what you mean by "u(last value)". The "last value" is always the "actual value", and this actual value is taken when you set s = -20*u(1). If the source term depends on a concentration in the past (u(t-tau) for some tau > 0), you have a delay PDE. Such a construct cannot be solved with pdepe.

What doesn't look right to you ? The concentration profile ?

[Juliana]

i see, youre right its the same as u(1), but yeah i would expect the concentration profile to look much different, but i guess im wrong idk