Operation error in problem definition

[sumantkr (Sumant)]

Dear all, I have defined a problem based on a non-dimensionalized coupled equation.

Unfortunately, I am getting an OPERATIONAL error in spite of having a correct operation( as per my weak formulation).

I have written the code as:

real Pr = 0.71;
real Ra = 1e5;
real pi = 3.14;
int n = 20;

mesh Th = square(n,n);
plot(Th);

fespace Uh(Th, P2);
fespace Ph(Th, P2);
fespace Kh(Th, P2);

Uh u, v, uh, vh;
Ph p, ph;
Kh T, th;



problem porous( [u, v, p, T], [uh,vh, ph, th]) =
                int2d(Th)(ph * dx(u) + ph * dy(v) )
                -int2d(Th)( u*u*dx(uh))
                -int2d(Th)( u*v*dy(uh))
                -int2d(Th)(p*dx(uh))
                +int2d(Th)(Pr*dx(u)*dx(uh) + Pr*dy(u)*dy(uh) )
                -int2d(Th)( u*v*dx(vh) )
                -int2d(Th)( v*v*dy(vh) )
                -int2d(Th)( p* dy(vh))
                +int2d(Th)( dx(v)*dx(vh)*Pr + Pr*dy(v)*dy(vh))
                +int2d(Th)( Ra*Pr*T*vh)
                +int2d(Th)(dx(t)*u*T + dy(t)*v*T)
                -int2d(Th)(dx(T)*dx(th) + dy(T)* dy(th) )
                -int1d(Th, 1)( th*dx(T)*N.x + th*dy(T)*N.y )
                +on(2,3,4, u =0, v = 0, T = 0)
                +on(1, u = 0, v= 0, T = 1/2*(1-cos(2*pi*x)) );


porous;

plot([u,v],p, T, nbiso =40, value = 1);

Error Is:

Kindly, Suggest me to fix it:

@prj @simon.garnotel

[prj]

Your problem is not linear. You need to linearize it.

[sumantkr (Sumant)]

Thanks @prj

• can you suggest me to make it linear based on the mentioned coupled equation.
  Or
  Kindly give me a link which has dealt with such non linear coupled equation.

It will give me some idea to deal with such issue.

[prj]

Here is a 2D example, adapted from this 3D code. Both should be straightforward to extend to your needs.

[sumantkr (Sumant)]

Thanks @prj
I was looking for such examples with coding.

• It has quality contents and helpful for me.

I’ll inform you after getting my code fixed with my model equation.

[frederichecht (Frédéric Hecht)]

THis is Boussnesq approximation, of a incompressible flow

Here you a have a exemple to solve this kind of equation

https://www.ljll.math.upmc.fr/hecht/ftp/ff++/2018-XVIII-las-Palmas/edp-XVIII/04Rayleigh-Benard.edp

[sumantkr (Sumant)]

Dear @frederichecht Sir,
I am glad to have an example that serves the same type of problems.
I can learn a key ideas from it.

Thanks

[sumantkr (Sumant)]

dear @prj
I have learnt the coding style for a non linear coupled equation:

few thing I am still struggling like:

• how can we control over the error
• how much iteration would be ok for a desired accuracy.

I have coded over a non linear pde:

governing equation and coding has been given below.

Results are not according to the expectation.

Kindly suggest me to figure it out( specially for a fixed value of parameters over domain)

Code:

real Pr =0.71;
real Ra = 1e3;
real Da = 1e-5;

real PrDa = Pr/Da;
int nn =10;


mesh Th = square( nn , nn );


macro grad(u) [dx(u), dy(u)] //

macro Grad(u) [ grad(u#1) , grad(u#2) ]  //

macro div(u) (dx(u#1)+ dy(u#2)) //

macro ugrad(u, v) ( [u#1, u#2]' * grad(v)) //

macro UgradV(u,v,T) ( [ugrad(u,v#1), ugrad(u,v#2), ugrad(u,T)] )  //

//defining finite element fespace

fespace Vh( Th, [P2, P2, P1, P1]);

Vh [u1, u2, p, T];  // variables tuples

Vh [ uw1, uw2, pw, Tw];  // increment in Newtons Methods

Vh [ v1,v2, q, TT];   // Test function

Vh [ vc1, vc2, pc, tc];  // computed values



//initial datd

[u1, u2, p, T] = [0,0,0, 1 ];


//plot( T, fill = 1, wait = 1, value =1);


problem porous([ uw1, uw2, pw, Tw], [ v1, v2, q, TT] ) =

int2d(Th) ( UgradV( u, uw, Tw)'* [v1, v2, TT ] + UgradV( uw, u, T)'*[ v1, v2, TT ] +
          + Pr*(Grad(uw):Grad(v)) - pw* div(v) + PrDa*[ uw1, uw2 ]'*[ v1, v2 ]
          - Ra*Pr* Tw *v1 -Ra*Pr*Tw*v2
          +grad(Tw)'* grad(TT)
          + q*div(uw)

          +1e-8*pw*q
)

+
int2d(Th)
(  UgradV(u, u, T)'*[v1, v2, TT]
- p*div(v)
+ Pr* (Grad(u):Grad(v))
+ PrDa* [u1, u2 ]'*[v1, v2]
- Ra*Pr*([T , T ]'*[v1,v2])
+ q* div(u)
+ grad(T)'* grad(TT)
+ 1e-8*p*q
)

+on( 1,2,3,4, uw1 =0, uw2 = 0)

+ on( 1, Tw = 1)
+ on(2,3,4, Tw =0);



real error = 0.0;
int count = 0;
int countmax = 5;
real err = 0;
real eps = 1e-6;

while ( (count < countmax) && (err < eps ))

count = count +1;
{
  int n;
  real err =0;
  for ( n = 0; n < 6; n++)
  {

  porous;

  u1[] = u1[] - uw1[];
  u2[] = u2[] - uw2[];
  p[] = p[] - pw[];
  T[] = T[] - Tw[];


real Lu1 = u1[].linfty, Lu2 = u2[].linfty, Lp = p[].linfty, LT = T[].linfty;

err = uw1[].linfty/Lu1 +  uw2[].linfty/Lu2 + pw[].linfty/Lp + Tw[].linfty/LT;

cout<< " Iteration ==="<<" "<< n <<" error is ===="<< err << endl;

plot([u1, u2], p, T, wait = 0 , value =1);

if (err < eps ) break;  // converge

if ( n < 3 && eps > 10. ) break; // blow up
/*
if (err< eps)

{

  plot([u1, u2], p, T, wait = 0 , value =1);

}
*/
 }//end ---> for loop
} // end  while loop on viscosity
//plot(T,cmm="Temperature Last converge of Newton Method for iteration = " + n,coef=0.3,wait=0,fill=1,value=1,nbiso=10);
plot([u1, u2], wait = 1,fill =1, value =1);
plot(u1, wait =1, fill =1, value =1);
plot( u2, wait =1, fill =1, value =1);

plot(p, wait = 1, value =1, fill =1);
plot( T, wait =1, value =1, fill =1);

After certain step, results totally get disappear.

Kindly suggest me, how to fix it.

Thanks

[prj]

Always start with a simple problem, and try to extend the complexity. Or start from something that works, and adapt it to your needs. What is not working? What is working? What did you change that broke what was previously working?

[sumantkr (Sumant)]

Thanks @prj
going to solve a simple model using Navier Stoke equation.

[frederichecht (Frédéric Hecht)]

This this a is not to hard problem,
compare to (https://www.ljll.math.upmc.fr/hecht/ftp/ff++/2018-XVIII-las-Palmas/edp-XVIII/04Rayleigh-Benard.edp) because the vortex is well defined

I just do a trivial modification on the script

Rayleigh-Benard-sumat.edp (2.2 KB)

[sumantkr (Sumant)]

Dear @frederichecht Sir,
I have thoroughly analyze the code that you have sent me.

I have one doubt:

• few terms for Newtons method under integral have not been included.
  can we exclude that linear terms ( i.e. for pressure, temperature,…)

I have rebuild that equation including all those terms but unfortunately, error terms are high ( approx 1).

Kindly suggest me to fix it.

// all coefficient used in problem

real Ra = 1e6;
real Pr = 0.71;
real Da = 1e-5;
real pRaPr = Ra*Pr;
real dPrDa = Pr/Da;

real[int] aRa = [1e3];

int nn = 20;
real  ccc = 1.5;//2.;
macro fff(x) ((1+tanh(ccc*(x*2-1))/tanh(ccc))/2)//

mesh Th=square(nn,nn,[x,(y*(1-x*0.0))]);
//plot(Th,wait=1);

// finite element fespace

fespace Vh(Th,[P2,P2,P1,P2]);

Vh [u1, u2, p, T];
Vh [uw1, uw2, pw, Tw];
Vh [v1, v2, q, TT];



//macro list used in the problem
macro div(u1, u2) ( dx(u1) + dy (u2) ) //
macro grad( u ) [ dx(u), dy(u) ] //
macro Grad(u1, u2) [ grad(u1), grad(u2) ] //
macro Ugradv(u1, u2, v) ([ u1, u2 ]'* grad(v))  //
macro UgradV(u1,u2, v1, v2 ) [Ugradv(u1, u2, v1), Ugradv(u1, u2, v2)]   //



problem porous( [ uw1, uw2, pw, Tw] , [ v1, v2, q, TT] ) =

        int2d(Th)( UgradV(uw1, uw2, u1, u2)'*[v1, v2]
                   + UgradV(u1, u2, uw1, uw2)'*[v1, v2]
                   - q*div(uw1, uw2)
                   - pw*div(v1, v2)
                   + Pr* (Grad(uw1, uw2): Grad(v1, v2))
                   + dPrDa * [uw1, uw2]'*[v1, v2]
                   - pRaPr * v2*Tw
                   +Ugradv(uw1, uw2, T)* TT
                   +Ugradv(u1, u2, Tw)* TT
                   + grad(Tw)'*grad(TT)
                   + 1e-8 * pw* q)

        - int2d(Th) ( UgradV(u1,u2,u1,u2)'* [v1, v2]

                    - q*div(u1, u2)
                    - p* div(v1, v2)
                    +Pr*(Grad(u1,u2):Grad(v1,v2))
                    + dPrDa* [u1, u2]'*[v1, v2]
                    -pRaPr*(T*v2)
                    +Ugradv(u1, u2, T)*TT
                    + grad(T)'*grad(TT)
                      + 1e-8*p*q      )
                    + on(1, 2,3,4, uw1 =0, uw2 =0, pw =0 )
                    + on(2,3, Tw = 0)
                    + on(1, Tw=sin(pi*x)); //Tw = 1);



real err;
real tol = 1e-8;
bool Newton = true;
[u1, u2, p, T] =[0, 0, 0, 1];

for(int kkk=0; kkk<aRa.n; ++kkk)
{
for (int step = 0; step <100 ; ++step)
{


  uw1[] = u1[];
  uw2[] = u2[];
  pw[] = p[];
  Tw[] = T[];

// solving problem
  porous;

  u1[]  = u1[] - uw1[];
  u2[]  = u2[] - uw2[];
  p[]  = p[] - pw[];
  T[]  = T[] - Tw[];

//plot([u1, u2], p, T, wait =1);

  real Lu1 = u1[].l2; //  , Lu2 = u2[].l2, Lp = p[].l2, LT = T[].l2;


  err = uw1[].l2/Lu1;    //+ uw2[].l2/Lu2 + pw[].l2/Lp + Tw[].l2/LT;

  cout << "err value is"<< err<<endl;

  //plot([u1,u2],T,wait=1,cmm="Ra ="+Ra);
  plot(T, wait = 1);
  if (err < tol) break;

  if ( step >3 && err >10 ) break;

}

if(kkk==0) {
   Th=square(nn,nn,[fff(x),fff(y)]);
   // u1,u2,p,T  privious mesh
   [uw1,uw2,pw,Tw]=[0,0,0,0];// // up1,up2,pp,Tp
   uw1[] = u1[];
   uw2[] = u2[];
   pw[] = p[];
   Tw[] = T[];
    //  new mesh
   [u1,u2,p,T]=[0,0,0,0];// // up1,up2,pp,Tp
   u1[]=uw1[];
   u2[] = uw2[];
   p[] = pw[];
   T[] = Tw[];
 }
plot([u1,u2],T,wait=1,cmm="Ra ="+Ra);

}

[frederichecht (Frédéric Hecht)]

Due to this writing

u1[] = u1[] - uw1[];
u2[] = u2[] - uw2[];
p[] = p[] - pw[];
T[] = T[] - Tw[];

First this is wrong because you use a coupled Vh [u1, u2, p, T];
so u1[], u2[], p[] , T[] are the same vector (same pointeur)

so you must write

u1[]  = u1[] - uw1[];

and remove

u2[]  = u2[] - uw2[];
p[]  = p[] - pw[];
T[]  = T[] - Tw[];

Secondly you use the base Newton Method so solve F(u)=0, : DF(u^n )w = F(u^1) , u^{n+1} = u^n - w.
where DF(u)w is the derivative of F in u in direction w (ie: DF(u+w) = F(u) + DF(u)w + o(w) )
The boundary condition (B.C.) must be 0 otherwise we lose B.C on u^{n+1}

This imply you have to initialise u^0 with correct B.C.

OTHER WISE YOU CAN REMARK:

DF(u^n) u^{n+1} = DF(u^n) u ^n - F(u^n)
with initial B.C.

[sumantkr (Sumant)]

Thanks @frederichecht Sir,
It really helps and you have given key information for dealing with the Newton algorithm.

• I have implemented the same that you have mentioned above.
• I have not implemented the Remark Part yet.

I have implemented your suggestion and it works fine, except for one thing.

• the error value decreasing from starting up to some iteration ( goes down up to 0.001 ) and after that it rises (approx to 0.01) and then starts decreasing ( up to 0.001 ).

Maybe I need to modify the initial guess.

[frederichecht (Frédéric Hecht)]

yes your initial given must have the Boundary condition.

[sumantkr (Sumant)]

Dear @frederichecht Sir,
After rectifying all the terms and modifying the issues, Now the error value term is decreasing.

I have one doubt for no of iteration required for achieving the sufficient accuracy.

• in my case, how much iteration required to achieve sufficient accuracy (e.g. 1.e-8).

actually, I have run the code for 300 steps and error decreases to (0.001).

but, in your case only few steps( approx 15) reached a desired accuracy ( 1e-8).

[frederichecht (Frédéric Hecht)]

possible your problem is hard, or make a bug in the computation of derivative

[sumantkr (Sumant)]

Sir, based on the above discussion and suggestion given by you, I have implemented the parallel code on a few coupled equation and results have a good agreement with the final results.

Now, I am stuck in a mixed differential type coupled problem, for which I am not getting the error reducing with the iteration. I have created a post on general discussion also.

I request you to kindly, suggest me regarding this issue.

Thansk

[sumantkr (Sumant)]

Dear @frederichecht Sir

I have been trying to solve the given governing equation with Newton Method and It is showing some bi-linearity issues.

Kindly help me to fix this issue.

Code and governing equation is attached here.

Code filenewton.edp (1.9 KB)

[sumantkr (Sumant)]

the above problem has been successfully solved and implemented according to my interest.

Thanks