Discussion:
pdepe solver - Time dependent boundary condition
(too old to reply)
Kamel
2010-04-21 18:25:04 UTC
Permalink
Hello,

I try to impose a boundary condition to resolve a parabolic equation using the pdepe solver. Do you know how to impose a discrete time boundary conition for the right BC like this (series of ramp displacements) :

for t<500s, ur = 0.0013*t, 500<t<1000: ur = 0.65, 1000<t<1500: ur = 0.0013*t, 500<t<2000: ur=0.65+0.65.

I tried something like:
pr = (ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)
+ (ur+0.0013*t)*(1000<t<1500)+ (ur+0.65+0.65)*(1500<t<2000);

but it works only for the two first blocks (up to t=1000).

Please, could you help me to find an easy way to get this BC applied ?

Thank you very much,

Kamel
ps: below is the main program


function pdex1

m = 0;
x = linspace(0,1.41,50);
t = linspace(0,2000,50);

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
u = sol(:,:,1);
save -ascii relax_10pts.dat u,x,t

% A solution profile can also be illuminating.
figure
%surf(u,x,t)
plot(t,u(:,:))
%plot(0.4*diff(u))
title('Displacement vs time for each z')
xlabel('t (s)')
ylabel('U (mm)')
% --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1/(0.4*2.7e-3);
f = DuDx;
s = 0;
% --------------------------------------------------------------
function u0 = pdex1ic(x)
u0 = 0;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)+(ur+0.0013*t)*(1000<t<1500)+(ur+0.65+0.65)*(1500<t<2000);
qr = 0;
Torsten Hennig
2010-04-22 06:48:01 UTC
Permalink
Post by Kamel
Hello,
I try to impose a boundary condition to resolve a
parabolic equation using the pdepe solver. Do you
know how to impose a discrete time boundary conition
for the right BC like this (series of ramp
for t<500s, ur = 0.0013*t, 500<t<1000: ur = 0.65,
: ur=0.65+0.65.
pr = (ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)
+ (ur+0.0013*t)*(1000<t<1500)+
(ur+0.65+0.65)*(1500<t<2000);
The boundary condition for u at the right boundary
point does not change continuously at time t=1000.
ur(1000-) = 0.65 ; ur(1000+) = 1.3.
Maybe this is a mistake in your setting.
If not, you will have to restart pdepe at time t=1000
with the solution obtained so far and artificially set
u=1.3 in the routine where the initial values are
supplied.

Another point is that you have to set
pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<t<1000)
+ (ur-0.0013*t)*(1000<t<1500)+
(ur-(0.65+0.65))*(1500<t<2000);

Best wishes
Torsten.
Post by Kamel
but it works only for the two first blocks (up to
t=1000).
Please, could you help me to find an easy way to get
this BC applied ?
Thank you very much,
Kamel
ps: below is the main program
function pdex1
m = 0;
x = linspace(0,1.41,50);
t = linspace(0,2000,50);
% Extract the first solution component as u.
u = sol(:,:,1);
save -ascii relax_10pts.dat u,x,t
% A solution profile can also be illuminating.
figure
%surf(u,x,t)
plot(t,u(:,:))
%plot(0.4*diff(u))
title('Displacement vs time for each z')
xlabel('t (s)')
ylabel('U (mm)')
%
------------------------------------------------------
--------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1/(0.4*2.7e-3);
f = DuDx;
s = 0;
%
------------------------------------------------------
--------
function u0 = pdex1ic(x)
u0 = 0;
%
------------------------------------------------------
--------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr =
ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)+(ur+0.0013
*t)*(1000<t<1500)+(ur+0.65+0.65)*(1500<t<2000);
qr = 0;
Kamel
2010-04-23 11:39:08 UTC
Permalink
Thank you Torsten,

dui you know how to use the restart the calculation with the solution obtained so far ?

Thank you very much,

Kamel
Post by Torsten Hennig
Post by Kamel
Hello,
I try to impose a boundary condition to resolve a
parabolic equation using the pdepe solver. Do you
know how to impose a discrete time boundary conition
for the right BC like this (series of ramp
for t<500s, ur = 0.0013*t, 500<t<1000: ur = 0.65,
: ur=0.65+0.65.
pr = (ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)
+ (ur+0.0013*t)*(1000<t<1500)+
(ur+0.65+0.65)*(1500<t<2000);
The boundary condition for u at the right boundary
point does not change continuously at time t=1000.
ur(1000-) = 0.65 ; ur(1000+) = 1.3.
Maybe this is a mistake in your setting.
If not, you will have to restart pdepe at time t=1000
with the solution obtained so far and artificially set
u=1.3 in the routine where the initial values are
supplied.
Another point is that you have to set
pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<t<1000)
+ (ur-0.0013*t)*(1000<t<1500)+
(ur-(0.65+0.65))*(1500<t<2000);
Best wishes
Torsten.
Post by Kamel
but it works only for the two first blocks (up to
t=1000).
Please, could you help me to find an easy way to get
this BC applied ?
Thank you very much,
Kamel
ps: below is the main program
function pdex1
m = 0;
x = linspace(0,1.41,50);
t = linspace(0,2000,50);
% Extract the first solution component as u.
u = sol(:,:,1);
save -ascii relax_10pts.dat u,x,t
% A solution profile can also be illuminating.
figure
%surf(u,x,t)
plot(t,u(:,:))
%plot(0.4*diff(u))
title('Displacement vs time for each z')
xlabel('t (s)')
ylabel('U (mm)')
%
------------------------------------------------------
--------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1/(0.4*2.7e-3);
f = DuDx;
s = 0;
%
------------------------------------------------------
--------
function u0 = pdex1ic(x)
u0 = 0;
%
------------------------------------------------------
--------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr =
ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)+(ur+0.0013
*t)*(1000<t<1500)+(ur+0.65+0.65)*(1500<t<2000);
qr = 0;
Torsten Hennig
2010-04-23 13:22:14 UTC
Permalink
Post by Kamel
Thank you Torsten,
dui you know how to use the restart the calculation
with the solution obtained so far ?
Thank you very much,
Kamel
An ad-hoc solution:

Write the solution vector at t=1000 into a
global variable vector and change the value at
the right boundary point appropriately.
Then call pdepe again with different function handles
for the inital and boundary routines.
In the new routine for the initial values, supply the
values from the globally accessible solution vector
at t = 1000.

Best wishes
Torsten.
Post by Kamel
in message
athforum.org>...
Post by Torsten Hennig
Post by Kamel
Hello,
I try to impose a boundary condition to resolve a
parabolic equation using the pdepe solver. Do you
know how to impose a discrete time boundary
conition
Post by Torsten Hennig
Post by Kamel
for the right BC like this (series of ramp
for t<500s, ur = 0.0013*t, 500<t<1000: ur = 0.65,
: ur=0.65+0.65.
pr = (ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)
+ (ur+0.0013*t)*(1000<t<1500)+
(ur+0.65+0.65)*(1500<t<2000);
The boundary condition for u at the right boundary
point does not change continuously at time t=1000.
ur(1000-) = 0.65 ; ur(1000+) = 1.3.
Maybe this is a mistake in your setting.
If not, you will have to restart pdepe at time
t=1000
Post by Torsten Hennig
with the solution obtained so far and artificially
set
Post by Torsten Hennig
u=1.3 in the routine where the initial values are
supplied.
Another point is that you have to set
pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<t<1000)
+ (ur-0.0013*t)*(1000<t<1500)+
(ur-(0.65+0.65))*(1500<t<2000);
Best wishes
Torsten.
Post by Kamel
but it works only for the two first blocks (up to
t=1000).
Please, could you help me to find an easy way to
get
Post by Torsten Hennig
Post by Kamel
this BC applied ?
Thank you very much,
Kamel
ps: below is the main program
function pdex1
m = 0;
x = linspace(0,1.41,50);
t = linspace(0,2000,50);
% Extract the first solution component as u.
u = sol(:,:,1);
save -ascii relax_10pts.dat u,x,t
% A solution profile can also be illuminating.
figure
%surf(u,x,t)
plot(t,u(:,:))
%plot(0.4*diff(u))
title('Displacement vs time for each z')
xlabel('t (s)')
ylabel('U (mm)')
%
------------------------------------------------------
Post by Torsten Hennig
Post by Kamel
--------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1/(0.4*2.7e-3);
f = DuDx;
s = 0;
%
------------------------------------------------------
Post by Torsten Hennig
Post by Kamel
--------
function u0 = pdex1ic(x)
u0 = 0;
%
------------------------------------------------------
Post by Torsten Hennig
Post by Kamel
--------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr =
ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)+(ur+0.0013
Post by Torsten Hennig
Post by Kamel
*t)*(1000<t<1500)+(ur+0.65+0.65)*(1500<t<2000);
qr = 0;
Torsten Hennig
2010-04-24 10:45:41 UTC
Permalink
Post by Kamel
Post by Kamel
Thank you Torsten,
dui you know how to use the restart the
calculation
Post by Kamel
with the solution obtained so far ?
Thank you very much,
Kamel
Write the solution vector at t=1000 into a
global variable vector and change the value at
the right boundary point appropriately.
Then call pdepe again with different function
handles
for the inital and boundary routines.
In the new routine for the initial values, supply
the
values from the globally accessible solution vector
at t = 1000.
Best wishes
Torsten.
If this makes technical difficulties, you could
modify your original boundary condition for ur to
pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<=t<=1000-eps)
+ (ur-(0.65+0.65/eps*(t-(1000-eps))))*
(1000-eps<t<1000)+
+ (ur-0.0013*t)*(1000<=t<1500)
for a small quantity eps > 0.

(Same at t=1500).

Best wishes
Torsten.
Post by Kamel
Post by Kamel
in message
athforum.org>...
Post by Torsten Hennig
Post by Kamel
Hello,
I try to impose a boundary condition to resolve
a
Post by Kamel
Post by Torsten Hennig
Post by Kamel
parabolic equation using the pdepe solver. Do
you
Post by Kamel
Post by Torsten Hennig
Post by Kamel
know how to impose a discrete time boundary
conition
Post by Torsten Hennig
Post by Kamel
for the right BC like this (series of ramp
for t<500s, ur = 0.0013*t, 500<t<1000: ur =
0.65,
Post by Kamel
Post by Torsten Hennig
Post by Kamel
: ur=0.65+0.65.
pr =
(ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)
Post by Kamel
Post by Torsten Hennig
Post by Kamel
+ (ur+0.0013*t)*(1000<t<1500)+
(ur+0.65+0.65)*(1500<t<2000);
The boundary condition for u at the right
boundary
Post by Kamel
Post by Torsten Hennig
point does not change continuously at time
t=1000.
Post by Kamel
Post by Torsten Hennig
ur(1000-) = 0.65 ; ur(1000+) = 1.3.
Maybe this is a mistake in your setting.
If not, you will have to restart pdepe at time
t=1000
Post by Torsten Hennig
with the solution obtained so far and
artificially
Post by Kamel
set
Post by Torsten Hennig
u=1.3 in the routine where the initial values
are
Post by Kamel
Post by Torsten Hennig
supplied.
Another point is that you have to set
pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<t<1000)
+ (ur-0.0013*t)*(1000<t<1500)+
(ur-(0.65+0.65))*(1500<t<2000);
Best wishes
Torsten.
Post by Kamel
but it works only for the two first blocks (up
to
Post by Kamel
Post by Torsten Hennig
Post by Kamel
t=1000).
Please, could you help me to find an easy way
to
Post by Kamel
get
Post by Torsten Hennig
Post by Kamel
this BC applied ?
Thank you very much,
Kamel
ps: below is the main program
function pdex1
m = 0;
x = linspace(0,1.41,50);
t = linspace(0,2000,50);
sol =
% Extract the first solution component as u.
u = sol(:,:,1);
save -ascii relax_10pts.dat u,x,t
% A solution profile can also be illuminating.
figure
%surf(u,x,t)
plot(t,u(:,:))
%plot(0.4*diff(u))
title('Displacement vs time for each z')
xlabel('t (s)')
ylabel('U (mm)')
%
------------------------------------------------------
Post by Kamel
Post by Torsten Hennig
Post by Kamel
--------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1/(0.4*2.7e-3);
f = DuDx;
s = 0;
%
------------------------------------------------------
Post by Kamel
Post by Torsten Hennig
Post by Kamel
--------
function u0 = pdex1ic(x)
u0 = 0;
%
------------------------------------------------------
Post by Kamel
Post by Torsten Hennig
Post by Kamel
--------
function [pl,ql,pr,qr] =
pdex1bc(xl,ul,xr,ur,t)
Post by Kamel
Post by Torsten Hennig
Post by Kamel
pl = ul;
ql = 0;
pr =
ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)+(ur+0.0013
Post by Kamel
Post by Torsten Hennig
Post by Kamel
*t)*(1000<t<1500)+(ur+0.65+0.65)*(1500<t<2000);
qr = 0;
Kamel
2010-04-27 11:49:04 UTC
Permalink
Thanks! I will try it.
By the way, do you know if I can resolve simultaneously these equations using the pdepe solver:
DuDt = k1*k2*DuDx2
DpDx = (1/k1)*DuDt

where u, p are dipslacment and pressure fields, k1 and k2 constants.
The second equation looks like hyperbolic form and I don't know if we can resolve it with pdepe.

Thank you,

Kamel
Post by Torsten Hennig
Post by Kamel
Post by Kamel
Thank you Torsten,
dui you know how to use the restart the
calculation
Post by Kamel
with the solution obtained so far ?
Thank you very much,
Kamel
Write the solution vector at t=1000 into a
global variable vector and change the value at
the right boundary point appropriately.
Then call pdepe again with different function
handles
for the inital and boundary routines.
In the new routine for the initial values, supply
the
values from the globally accessible solution vector
at t = 1000.
Best wishes
Torsten.
If this makes technical difficulties, you could
modify your original boundary condition for ur to
pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<=t<=1000-eps)
+ (ur-(0.65+0.65/eps*(t-(1000-eps))))*
(1000-eps<t<1000)+
+ (ur-0.0013*t)*(1000<=t<1500)
for a small quantity eps > 0.
(Same at t=1500).
Best wishes
Torsten.
Post by Kamel
Post by Kamel
in message
athforum.org>...
Post by Torsten Hennig
Post by Kamel
Hello,
I try to impose a boundary condition to resolve
a
Post by Kamel
Post by Torsten Hennig
Post by Kamel
parabolic equation using the pdepe solver. Do
you
Post by Kamel
Post by Torsten Hennig
Post by Kamel
know how to impose a discrete time boundary
conition
Post by Torsten Hennig
Post by Kamel
for the right BC like this (series of ramp
for t<500s, ur = 0.0013*t, 500<t<1000: ur =
0.65,
Post by Kamel
Post by Torsten Hennig
Post by Kamel
: ur=0.65+0.65.
pr =
(ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)
Post by Kamel
Post by Torsten Hennig
Post by Kamel
+ (ur+0.0013*t)*(1000<t<1500)+
(ur+0.65+0.65)*(1500<t<2000);
The boundary condition for u at the right
boundary
Post by Kamel
Post by Torsten Hennig
point does not change continuously at time
t=1000.
Post by Kamel
Post by Torsten Hennig
ur(1000-) = 0.65 ; ur(1000+) = 1.3.
Maybe this is a mistake in your setting.
If not, you will have to restart pdepe at time
t=1000
Post by Torsten Hennig
with the solution obtained so far and
artificially
Post by Kamel
set
Post by Torsten Hennig
u=1.3 in the routine where the initial values
are
Post by Kamel
Post by Torsten Hennig
supplied.
Another point is that you have to set
pr = (ur-0.0013*t)*(t<500)+(ur-0.65)*(500<t<1000)
+ (ur-0.0013*t)*(1000<t<1500)+
(ur-(0.65+0.65))*(1500<t<2000);
Best wishes
Torsten.
Post by Kamel
but it works only for the two first blocks (up
to
Post by Kamel
Post by Torsten Hennig
Post by Kamel
t=1000).
Please, could you help me to find an easy way
to
Post by Kamel
get
Post by Torsten Hennig
Post by Kamel
this BC applied ?
Thank you very much,
Kamel
ps: below is the main program
function pdex1
m = 0;
x = linspace(0,1.41,50);
t = linspace(0,2000,50);
sol =
% Extract the first solution component as u.
u = sol(:,:,1);
save -ascii relax_10pts.dat u,x,t
% A solution profile can also be illuminating.
figure
%surf(u,x,t)
plot(t,u(:,:))
%plot(0.4*diff(u))
title('Displacement vs time for each z')
xlabel('t (s)')
ylabel('U (mm)')
%
------------------------------------------------------
Post by Kamel
Post by Torsten Hennig
Post by Kamel
--------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
c = 1/(0.4*2.7e-3);
f = DuDx;
s = 0;
%
------------------------------------------------------
Post by Kamel
Post by Torsten Hennig
Post by Kamel
--------
function u0 = pdex1ic(x)
u0 = 0;
%
------------------------------------------------------
Post by Kamel
Post by Torsten Hennig
Post by Kamel
--------
function [pl,ql,pr,qr] =
pdex1bc(xl,ul,xr,ur,t)
Post by Kamel
Post by Torsten Hennig
Post by Kamel
pl = ul;
ql = 0;
pr =
ur+0.0013*t)*(t<500)+(ur+0.65)*(500<t<1000)+(ur+0.0013
Post by Kamel
Post by Torsten Hennig
Post by Kamel
*t)*(1000<t<1500)+(ur+0.65+0.65)*(1500<t<2000);
qr = 0;
Torsten Hennig
2010-04-28 05:35:46 UTC
Permalink
Post by Kamel
Thanks! I will try it.
.. or you can define an event function to restart the
simulation at t = 1000.
See the pdepe documentation for more details.
Post by Kamel
By the way, do you know if I can resolve
simultaneously these equations using the pdepe
DuDt = k1*k2*DuDx2
DpDx = (1/k1)*DuDt
From
dp/dx = (du/dt)/k1 = k2*d^2u/dx^2
you get
p(x) - p(x_0) = k2 * ((du/dx)(x) - (du/dx)(x_0))
and so
p(x) = p(x_0) + k2*((du/dx)(x) - (du/dx)(x_0)) (1)

Use pdeval to compute du/dx at the end of the simulation
for the output times and 'reconstruct' p by (1).
Post by Kamel
where u, p are dipslacment and pressure fields, k1
and k2 constants.
The second equation looks like hyperbolic form and I
don't know if we can resolve it with pdepe.
Thank you,
Kamel
Best wishes
Torsten.
Kamel
2010-05-23 22:16:04 UTC
Permalink
Hello Torsten,

I would like to prescribe this expression for the right boundary condition (ur vs t) :

step_up_times= [100 400 600 800];
n = 0;
ur = zeros(1,1007);
t = 1:1007;
for k = 1:length(ur)
if t(k) <= step_up_times(end) && t(k) >= step_up_times(n+1)
n = n+1;
end
ur(k) = .1*n;
end

This is still a ramp boundary condition.
Do you know how to introduce it easily in the function pdex1bc ?

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
pl = ul;
ql = 0;
pr = "ur vs t" ;
qr = 0;

Cheers,

Kamel
Post by Torsten Hennig
Post by Kamel
Thanks! I will try it.
.. or you can define an event function to restart the
simulation at t = 1000.
See the pdepe documentation for more details.
Post by Kamel
By the way, do you know if I can resolve
simultaneously these equations using the pdepe
DuDt = k1*k2*DuDx2
DpDx = (1/k1)*DuDt
From
dp/dx = (du/dt)/k1 = k2*d^2u/dx^2
you get
p(x) - p(x_0) = k2 * ((du/dx)(x) - (du/dx)(x_0))
and so
p(x) = p(x_0) + k2*((du/dx)(x) - (du/dx)(x_0)) (1)
Use pdeval to compute du/dx at the end of the simulation
for the output times and 'reconstruct' p by (1).
Post by Kamel
where u, p are dipslacment and pressure fields, k1
and k2 constants.
The second equation looks like hyperbolic form and I
don't know if we can resolve it with pdepe.
Thank you,
Kamel
Best wishes
Torsten.
Loading...