Unable to display the plot value for velocity

Hello Everyone,
I am not getting the plot for [u1, u2].

kindly, suggest me to fix it.

code has been given as:

/*
coupled equation is
u_x + v_y = 0
((u,v).grad(u,v))* + grad(p) -1/Re * ( lap(u, v) ) = 0
(u, v). ( grad(T) ) - lap( T) = 0


B. C.
over square domain,
: all velocity components are zero
: temperature value is zero at left and right edge
and sin(pi*x) at the bottom edge


*/


int nn = 20;

real err;
real tol = 1.e-8;
real Re = 2000;

//verbosity =0;

mesh Th = square(nn, nn );


macro grad(u) [ dx(u), dy(u)]                         //
macro div(u1, u2) (dx(u1) + dy(u2) )                //
macro Grad(u1, u2) [ grad(u1), grad(u2) ]            //
macro Ugradv( u1, u2, v1)  ([u1, u2]'*grad(v1) )        //
macro UgradV(u1, u2, v1, v2) [Ugradv(u1, u2, v1), Ugradv(u1, u2, v2) ]  //


// finite element spaces

fespace Vh( Th, [P2, P2, P1b, P2]);
Vh [u1, u2, p, T];  // variables
Vh [du1, du2, dp, dT];   // incremenet
Vh [v1, v2, q, TT];    // test function


problem heat([du1, du2, dp, dT], [v1, v2, q, TT]) =

                int2d(Th) ( div(du1, du2)*q
                        + [v1, v2]'*UgradV( du1, du2, u1, u2)
                        + [v1, v2]'*UgradV( u1, u2, du1, du2)
                        + 1./Re*(Grad(du1, du2): Grad(v1, v2))
                        - dp* div(v1, v2)
                        + Ugradv(du1, du2, T)*TT
                        + Ugradv(u1, u2, dT)*TT
                        + grad(dT)'*grad(TT)
                        - q* div(du1, du2)
                        - 1.e-8*dp*q
                  )

                  - int2d(Th) (
                     [v1, v2]'*UgradV(u1, u2, u1, u2)
                    +[v1, v2]'*UgradV(u1, u2, u1, u2)
                    + 1./Re*( Grad(u1, u2): Grad(v1, v2) )
                    - p* div(v1, v1)
                    + TT* Ugradv(u1, u2, T)
                    + grad(T)'*grad(TT)
                    - q*div(u1, u2)
                    - 1.e-8*p*q
                    )
                    + on( 1, 2,3,4, du1=0, du2 = 0)
                    + on( 2, 4, dT = 0)
                    + on(1, dT = sin(3.14*x));


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

for (int step = 0; step < 15; step++)
{
  heat;

  u1[] = u1[] - du1[];

  real solerr = u1[].linfty;
  real stperr = du1[].linfty;

  err = stperr/solerr;

  cout<< "step is == "<< step<<endl;
  cout<< "error is == "<< err<<endl;

  if(err< tol)  break;
  if( step > 3 && err > 1.e3) break;
  plot([u1, u2],p,T, wait = 0, value =1 );
  //plot( T, wait =0, fill =1,  value =1);


}

Hi there!

I don’t think this solver is going to work as you intend it to as is, but you can fix the plotting issue by replacing your code as follows:

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

Remember that you can view previous calls to plot by pressing “p” on the plot viewer (hit “?” while viewing the plot window for a full list of keyboard shortcuts). Your old code was covering the other plots with the plot for T.

1 Like

Dear @cmd,
Thanks for your reply. One more thing I would like to ask for the solution procedure,

  • is this process is fine for the mentioned coupled eq.

  • how much iterations are generally ok for achieving a desired accuracy

( Generally I have seen merely few steps is sufficient for achieving a desired accuracy ( 1e-8))

Kindly suggest.

  1. I think your linearization is incorrect. On the linear part you have div(du1, du2)*q ... - q* div(du1, du2) which cancels out. You also wrote the nonlinear advection term twice in the residual. Lastly, your BC are not right after the first iteration. I think you meant what I have put below. Beware of your choice of elements as I don’t believe this is a stable choice of mixed elements, though it is probably ok with the stabilization term.

    problem heat([du1, du2, dp, dT], [v1, v2, q, TT]) =
             int2d(Th) (
                       [v1, v2]'*UgradV( du1, du2, u1, u2)
                     + [v1, v2]'*UgradV( u1, u2, du1, du2)
                     + 1./Re*(Grad(du1, du2): Grad(v1, v2))
                     - dp* div(v1, v2)
                     + Ugradv(du1, du2, T)*TT
                     + Ugradv(u1, u2, dT)*TT
                     + grad(dT)'*grad(TT)
                     - q* div(du1, du2)
                     - 1.e-8*dp*q
               )
    
               - int2d(Th) (
                  [v1, v2]'*UgradV(u1, u2, u1, u2)
                 + 1./Re*( Grad(u1, u2): Grad(v1, v2) )
                 - p* div(v1, v1)
                 + TT* Ugradv(u1, u2, T)
                 + grad(T)'*grad(TT)
                 - q*div(u1, u2)
                 - 1.e-8*p*q
                 )
                 + on( 1, 2,3,4, du1=u1, du2 = u2)
                 + on( 2, 4, dT = T)
                 + on(1, dT = T-  sin(3.14*x));
    
  2. Because have homogeneous BC’s for the velocity field and no coupling terms, in this case, the (homogeneous) initial guess is already a solution of the velocity/pressure part. In this case, the temperature part is linear and only one iteration is required. In the more general case, the number of Newton iterations required for convergence (and whether they converge at all!) depends on the strength of the nonlinearity in the problem.

1 Like

Thanks @cmd,
yes, I did a mistake that you have rectified.

you have suggested about the choice of mixed elements

  • How can we combine these elements
    and in my case, the choice of elements are optimal ?

Taylor-Hood (P2-P1) or MINI (P1b-P1) elements are both suitable pairs. Taylor Hood is second order, but MINI is more economical. You can find others in the literature but these are usually a good place to start.

1 Like

Thanks @cmd,
Please clarify one more doubt

-a) in my problem, can we adopt any method or solver to solve such type of non linear coupled .

-b) rate of decrease of error is very slow ( for my problem) and it takes approx ( 3000 steps to reduce an error to 0.0001.

  • c) and, what is wrong in my vector field plotting.
    Only, temperature field, is showing after inserting code.

Specially (b) part creating a trouble for solution.

Please, suggest me to fix this error.

Thanks

  1. The Newton method you have implemented should handle the nonlinear problem just fine. If convergence is an issue you may be able to just use a continuation on Re to solve your problem.

  2. Try using the corrected code I’ve posted below which uses an l2 norm.

  3. See my earlier comment. press “p” to see older plots.

    /*
    coupled equation is
    u_x + v_y = 0
    ((u,v).grad(u,v))* + grad(p) -1/Re * ( lap(u, v) ) = 0
    (u, v). ( grad(T) ) - lap( T) = 0
    
    B. C.
    over square domain,
    : all velocity components are zero
    : temperature value is zero at left and right edge
    and sin(pi*x) at the bottom edge
    */
    int nn = 20;
    real tol = 1.e-8;
    real Re = 2000;
    //verbosity =0;
    mesh Th = square(nn, nn);
    
    macro grad(u)[dx(u), dy(u)]//
    macro div(u1, u2)(dx(u1) + dy(u2))//
    macro Grad(u1, u2)[grad(u1), grad(u2)]//
    macro Ugradv( u1, u2, v1)([u1, u2]'*grad(v1))//
    macro UgradV(u1, u2, v1, v2)[ Ugradv(u1, u2, v1), Ugradv(u1, u2, v2)]//
    
    // finite element spaces
    fespace Vh( Th, [P2, P2, P1, P2]);
    Vh [u1, u2, p, T];  // variables
    Vh [du1, du2, dp, dT];   // incremenet
    Vh [v1, v2, q, TT];    // test function
    [u1, u2, p, T] = [0,0,0,0];
    
    // problem
     problem heat([du1, du2, dp, dT], [v1, v2, q, TT])
    = int2d(Th)(
      [v1, v2]'*UgradV( du1, du2, u1, u2)
    + [v1, v2]'*UgradV( u1, u2, du1, du2)
    + 1./Re*(Grad(du1, du2): Grad(v1, v2))
    - dp*div(v1, v2)
    + Ugradv(du1, du2, T)*TT
    + Ugradv(u1, u2, dT)*TT
    + grad(dT)'*grad(TT)
    - q*div(du1, du2)
    - 1e-8*dp*q
    )
    - int2d(Th) (
      [v1, v2]'*UgradV(u1, u2, u1, u2)
    + 1./Re*( Grad(u1, u2): Grad(v1, v2) )
    - p*div(v1, v1)
    + Ugradv(u1, u2, T)*TT
    + grad(T)'*grad(TT)
    - q*div(u1, u2)
    - 1e-8*p*q
    )
    + on(1, 2, 3, 4, du1 = u1, du2 = u2)
    + on(2, 4, dT = T)
    + on(1, dT = T - sin(3.14*x));
    
    // Newton loop
    for (int step = 0; step < 15; step++){
      heat;
      u1[] -= du1[];
      real err = du1[].l2/u1[].l2;
      cout<< "step "<< step << " error = " << err << endl;
      if(err < tol)  break;
      if( step > 3 && err > 1.e3) break;
      plot([u1, u2], wait = 0, value = 1, cmm = "[u1, u2]");
      plot(p, wait = 0, value = 1, cmm = "p");
      plot(T, wait = 0, value = 1, cmm = "T");
    }
1 Like

Thanks @cmd, your updated code almost solved my convergence issue.

Only, the output results plots for pressure and velocity.
Its values are almost zero everywhere on the domain.

-Is it my Model Issue?

pressure plot is:

Velocity Plot is:

Your model has homogeneous BC and no coupling terms, so homogeneous velocity/pressure fields are the correct solution for the problem as posed… If you intended to couple the flow field to the temperature through a Boussinesq type coupling you’ll need to add that.