# Evolution problem error

Dear All,
I am trying to solve transient problem for porous media.
Unfortunately, the time-dependent function is not updating after time increment, which gives us an undesired output value.
Kindly suggest, where it went wrong.

code is as mentioned below/////////////

``````//real uold=473.14;
//func uo = 473.15-10*t;
border C1(t1=-1,1) { x = 10*t1; y= -10; }
border C2(t1=-1,1) { x = 10; y = 10*t1;  }
border C3(t1=-1,1) { x = -10*t1; y =10; }
border C4(t1=-1,1) { x = -10; y = -10*t1;}
mesh Th= buildmesh(C1(20)+C2(20)+C3(20)+C4(20));
//plot(Th, wait=1);

fespace Vh(Th, P1);
Vh u, v, uu, f, g=473.15;

real dt=0.01;
real phio= 0.1003, rhof=1000.0, cvf = 4200, kts=2.49, cvs=1000.0, rhos=2910.2;
real cpf=4200, ktf=0.6, qs =0;
real kt= (1-phio)*kts+phio*ktf;
real RhoCv= (1-phio)*rhos*cvs+ phio*rhof*cvf;

problem  pht(u,v) = int2d(Th) ( dt*kt*(dx(u)*dx(v)+dy(u)*dy(v))  )
+ int2d(Th) (RhoCv*v*u)
+ int2d(Th)(-RhoCv*v*uu)
+ int2d(Th) (-qs*v)
+ on(C2, u=f)
+ on(C4, u= g);

real t =0;                          // initial time
uu = 473.15;                   // initial condition
for ( int m = 0; m <=3/dt; m++ )
{

t = t+dt;
cout<<"time is ="<< t<<endl;
f = 473.15 - 10*t;
cout<< "func is f ="<< f[]<<endl;

pht;
plot(u, wait = 0, value =1, fill= 1);
uu = u;
//cout<< "result at center="<< u(0,0)<<endl;
}
``````

///////////////////////////////////////////

Kindly suggest me how to rectify it. : )

Why did you use an fespace function to define `f` ?
I think use `func f = 473.15 - 10*t;` outside the loop do the job.

• Thanks for your suggestion but it is not working.
• I have taken the reference for solving my problem is Evolution Problem given in documentation of FreeFem++. link is
=========================================
https://doc.freefem.org/models/evolution-problems.html
==========================================

while taking function f outside an error occurs " identifier t does not exist".

• I put real t =0; before, then it’s getting fixed but again the f value is not updating after time loop increasing.

kindly help me.

Here is the boundary condition imposed over geometry.

also, the rectified code is as mentioned below.

``````//real uold=473.14;
//func uo = 473.15-10*t;
border C1(t1=-1,1) { x = 10*t1; y= -10; }
border C2(t1=-1,1) { x = 10; y = 10*t1;  }
border C3(t1=-1,1) { x = -10*t1; y =10; }
border C4(t1=-1,1) { x = -10; y = -10*t1;}
mesh Th= buildmesh(C1(20)+C2(20)+C3(20)+C4(20));
//plot(Th, wait=1);

real g=473.15;
fespace Vh(Th, P1);
Vh u, v,uu;

real dt=0.01;
real phio= 0.1003, rhof=1000.0, cvf = 4200, kts=2.49, cvs=1000.0, rhos=2910.2;
real cpf=4200, ktf=0.6, qs =0;
real kt= (1-phio)*kts+phio*ktf;
real RhoCv= (1-phio)*rhos*cvs+ phio*rhof*cvf;

real t =0;                          // after suggestion
func f = 473.15 -10*t;              // after suggestion

problem  pht(u,v) = int2d(Th) ( dt*kt*(dx(u)*dx(v)+dy(u)*dy(v))  )
+ int2d(Th) (RhoCv*v*u)
+ int2d(Th)(-RhoCv*v*uu)
+ int2d(Th) (-qs*v)
+ on(C2, u=f)
+ on(C4, u= g);

//real t =0;    // initial time

uu = 473.15;  // initial condition
for ( int m = 0; m <=3/dt; m++ )
{

t = t+dt;
cout<<"time is ="<< t<<endl;
// f = 473.15 - 10*t;
// cout<< "func is f ="<< f[]<<endl;

pht;
plot(u, wait = 0, value =1, fill= 1);
uu = u;
cout<< "result at center="<< u(0,0)<<endl;
}
``````

``````// Parameters
real g = 473.15;
real dt = 0.01;
real phio = 0.1003, rhof=1000.0, cvf = 4200, kts=2.49, cvs=1000.0, rhos=2910.2;
real cpf = 4200, ktf = 0.6, qs = 0;
real kt = (1-phio)*kts+phio*ktf;
real RhoCv = (1-phio)*rhos*cvs+ phio*rhof*cvf;

real t = 0.;
func f = 473.15 - 10*t;

// Mesh
border C1(t1=-1,1) { x = 10*t1; y= -10; }
border C2(t1=-1,1) { x = 10; y = 10*t1;  }
border C3(t1=-1,1) { x = -10*t1; y =10; }
border C4(t1=-1,1) { x = -10; y = -10*t1;}
mesh Th = buildmesh(C1(20)+C2(20)+C3(20)+C4(20));

// FE space
fespace Vh(Th, P1);
Vh u, v, uu;

// Problem
problem  pht(u,v)
= int2d(Th) ( dt*kt*(dx(u)*dx(v)+dy(u)*dy(v))  )
+ int2d(Th) (RhoCv*v*u)
+ int2d(Th)(-RhoCv*v*uu)
+ int2d(Th) (-qs*v)
+ on(C2, u=f)
+ on(C4, u= g);

// Initialization
uu = 473.15;

// Time loop
for (int m = 0; m <= 3./dt; m++) {
// Update
t = t+dt;
cout << "time is =" << t << endl;

// SOlve
pht;

// Plot
plot(u, wait = 0, value =1, fill= 1);

// Display
// cout<< "result at center="<< u(0,0)<<endl;
cout << "Value on C2 = " << u(10, 0) << endl;

// Update
uu = u;
}
``````

You can see in the `cout` ‘value on C2’ that the value of `g` is updated at each time step.

Dear Simon,
Thank you for your suggestion and now the f value is updating but still that the temperature variation should be gradient type behavior like in time evolution problem (like the above result).

but my output result from this code is totally different from the original output result (pic below).

Also, the output at the center is supposed to be u(0,0) = 458.053, which I am not getting.

where I am getting wrong according to the problem, I am not getting.

kindly, help me.

Thanks again. It just depend of your parameters, for example try this code:

``````// Parameters
real fT = 3.;
real dt = 1.e-3;

real kt = 1.;
real rho = 1.;
real cv = 1.e-3;

real T0 = 473.15;

// Mesh
border C1(t1=-1,1) { x = 10*t1; y= -10; }
border C2(t1=-1,1) { x = 10; y = 10*t1;  }
border C3(t1=-1,1) { x = -10*t1; y =10; }
border C4(t1=-1,1) { x = -10; y = -10*t1;}
mesh Mesh = buildmesh(C1(20)+C2(20)+C3(20)+C4(20));

// FE space
fespace Vh(Mesh, P1);
Vh T, Th, Tp;

// Macro

// Boundary
real t = 0.;
real T1 = 473.15;
func T2 = 473.15 - 10.*t;

// Problem
problem pht (T, Th)
= int2d(Mesh)(
(rho*cv/dt) * T * Th
)
- int2d(Mesh)(
(rho*cv/dt) * Tp * Th
)
+ on(C2, T=T2)
+ on(C4, T=T1)
;

// Initialization
T = T0;

// Time loop
int nbIter = fT / dt;
for (int m = 0; m < nbIter; ++m) {
// Update
t = t + dt;
Tp = T;
cout << "time is =" << t << endl;

// Solve
pht;

// Plot
plot(T, value=1, fill=1);
}
``````

The result looks like what you want, no?

I think you have to take care of unities of your different quantities.

Thanks Simon,
It looks like the same result that I am expecting.

I would like to tell you, I have solved the steady problem (excluding time derivative term) which has the same output result with the paper result.

But, this transient problem with the same parameter is creating a bug result.

You helped me a lot, I’ll fix the parameter once again and inform you after getting solved.

Thanks.

S

• Can you suggest me, how to get the output results up to 10 digits after the decimal place and save it in a dat file format?
``````ofstream output("output.dat");
output.precision(10);
for(int i = 0; i < u[].n; ++i)
output << u[][i] << "\n";
``````
1 Like