Initial condition not used in parablic equation

General Discussion
1

I’m still struggling with my system of equations but simplified it to 1 parabolic equation to find the issue and feal like I lost some basics.

In below equation value of initial condition doesn’t influence results, looks like only boundary condition changes function values. What’s more function values don’t change much in following time steps.

What I’m missing?


border C01(t=0.06, 2.08){x=t; y=0.002; }
border C02(t=0.002, 1.002){x=2.08; y=t;}
border C03(t=2.08, 0.06){x=t; y=1.002;}
border C04(t=1.002, 0.002){x=0.06; y=t;}  
int n = 100;
mesh Th = buildmesh(C01(n) + C02(n/2) + C03(n) + C04(n/2));


fespace Vh1(Th,P1);
// Parameters
real Dgamma = 1000;
real beta = 0.85;
real gamma0 = 0.85;
real w1 = 60; //boundary condition

real t=0, T=5, dt=1;
 
Vh1 gamma = gamma0, // initial condition
gammaold,
o2,
uv1 = 1,
uv2 = 1;

macro wektoru(w1,w2) [w1,w2]	//EOM
macro grad(ff) [dx(ff),dy(ff)]	//EOM


problem Probgamma(gamma,o2,solver=LU)=
int2d(Th,qft=qf1pTlump)(gamma/dt*o2)
-int2d(Th,qft=qf1pTlump)(gammaold/dt*o2)
-int2d(Th)(Dgamma*(dx(gamma)*dx(o2) + dy(gamma)*dy(o2))) 
-int2d(Th)(wektoru(uv1, uv2)'*grad(gamma)*o2) 
+int2d(Th)(beta*gamma*o2)

+on(C04,   gamma = w1) 
;

for(real t=t;t<=T;t+=dt)
{
cout<<"===================================Solution at time: "<<t<<endl;

gammaold=gamma;
Probgamma; 

plot(gamma,value=1,fill=1,wait=0, cmm="solution gamma at time: "+t);

}
2

I didn’t try to write it out but just looking at the problem the laplacian coef is
large suggesting second der is going to be “small.” It looks like the time derivative goes
to zero when gamma is just the boundary condition plus an exponential ( almost linear )
change in the x direction. I played with solver and took out quadrature rule not
much change although setting “init=1” seems to change things a little… Beta and 1/dt
seem to be similar magnitude although again you should probably just write it
out and look at the numbers. I should know the damped thing pretty well but its late here :slight_smile: There are cases when a condition can be ignored, for example setting
initial value and initial derivative without a mixed formulation but I’m not
sure any case like that exists here. What result were your expecting or what equation
were you trying to implement?

load "medit"
border C01(t=0.06, 2.08){x=t; y=0.002; }
border C02(t=0.002, 1.002){x=2.08; y=t;}
border C03(t=2.08, 0.06){x=t; y=1.002;}
border C04(t=1.002, 0.002){x=0.06; y=t;}  
int n = 100;
mesh Th = buildmesh(C01(n) + C02(n/2) + C03(n) + C04(n/2));


fespace Vh1(Th,P1);
// Parameters
real Dgamma = 1000;
real beta = 0.85;
real gamma0 = 0.85;
real w1 = 60; //boundary condition

real t=0, T=5*100, dt=1;
 
Vh1 gamma = gamma0, // initial condition
gammaold,
o2,
uv1 = 1,
uv2 = 1;

macro wektoru(w1,w2) [w1,w2]	//EOM
macro grad(ff) [dx(ff),dy(ff)]	//EOM


//problem Probgamma(gamma,o2,solver=LU,init=1)=
problem Probgamma(gamma,o2,solver=sparsesolver,init=0)=
//int2d(Th,qft=qf1pTlump)(gamma/dt*o2)
int2d(Th)(gamma/dt*o2)
//-int2d(Th,qft=qf1pTlump)(gammaold/dt*o2)
-int2d(Th)(gammaold/dt*o2)
-int2d(Th)(Dgamma*(dx(gamma)*dx(o2) + dy(gamma)*dy(o2))) 
-int2d(Th)(wektoru(uv1, uv2)'*grad(gamma)*o2) 
+int2d(Th)(beta*gamma*o2)

+on(C04,   gamma = w1) 
;

for(real t=t;t<=T;t+=dt)
{
cout<<"===================================Solution at time: "<<t<<endl;

gammaold=gamma;
Probgamma; 

plot(gamma,value=1,fill=1,wait=0, cmm="solution gamma at time: "+t);
//medit("soln ",Th,gamma);
}
3

Thanks for answering. I’ll look at it deeply again after the weekend.

Checking shortly at your code changes I see that initial condition influences results – good and expected. However, the solution doesn’t change in time like the gamma is constant – so the time derivative is not influencing.

If the initial condition is beta = 0.85 and boundary condition w1 = 60 I would expected values ~between these two in initial step and changing in the time. This equation describes oxygen flow in tissues so expect oxygen from boundary migrates to the whole tissue.

4

What eqn are you trying to model?
AFAICT, as written, for the steady state soln characteristic eqn in x direction roots I got (.5/D)(1 +/- sqrt(1-4BD)).
Your x values are somewhere around 0-2 or so and D is 1000 so I don’t
expect the soln to change much but in any case it should not be hard to put in
analytical soln and compare.

5

I’m trying to simulate tumor growth, this is oxygen density equation.
Probably you gave me important clue which I missed as I’m not so fluent in equations. Original area was rectangle 100-200 not 1-2, I rescaled it due to some other issues I didn’t realize how it influences current version.
After changing mesh size (to ~100-200) function values start to change however become more unstable in time (growing rapidly or take negative values) – need to verify this deeply.

6

If you are just trying to simulate 1 dimensional diffusion, there are a lot of examples of
that- heat diffuses and so do charge carriers. Getting a useful result will
probably require a more complicated model or better boundary conditions.
Its good to start simple and your current model is so simple there are ways
to approach it analytically which is always a good sanity check on a computer.
I’m an engineer so my description of the math may be vulgar but hopefully
useful :slight_smile: Probably also you are interested in steady state or quasi static
situation first where the time derivatives don’t matter so much. But, you
may want to see analytically what happens to an initial guess of the form
exp(ikx) or exp(kx) whatever is easier. or for complex k. Then you get
something like dn/dt=a(k)*n and you can see if there are conditions where it grows
indefinitely. If seems to go to zero analytically and the simulation explodes,
it could be a bug or the time or space scales are too large.