Evolution problem error

sumantkr · April 10, 2020, 10:23am

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. : )

simon.garnotel · April 10, 2020, 6:52pm

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.

sumantkr · April 11, 2020, 3:38am

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;
    }

simon.garnotel · April 11, 2020, 7:21am

I reorganize your code:

// 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.
Is there a problem with your model instead ?

sumantkr · April 11, 2020, 9:36am

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.

simon.garnotel · April 11, 2020, 9:02pm

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
macro grad(u)[dx(u), dy(u)]//

// 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
    + kt * grad(T)' * grad(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.

sumantkr · April 12, 2020, 1:16am

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

sumantkr · April 12, 2020, 1:19am

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?

prj · April 12, 2020, 9:44am

ofstream output("output.dat");
output.precision(10);
for(int i = 0; i < u[].n; ++i)
  output << u[][i] << "\n";

sumantkr · April 12, 2020, 5:49pm

Thanks all for your help.
time variation is changing per year, so, after converting it into second would give me the desired solution.

Topic		Replies	Views
Numerical results error General Discussion	8	169	September 7, 2023
Non-Dimensional Equation error General Discussion	0	431	July 4, 2020
Operation error in problem definition General Discussion	20	2486	September 15, 2021
Error in FreeFEM Solving Equations with Time Terms General Discussion	1	72	June 24, 2024
Newton Method Convergence Error General Discussion	0	527	August 29, 2020

Evolution problem error

Related Topics