Unable to display the plot value for velocity

sumantkr · July 10, 2020, 9:59am

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


}

cmd · July 10, 2020, 12:24pm

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.

sumantkr · July 10, 2020, 12:43pm

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.

cmd · July 10, 2020, 1:16pm

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

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.

sumantkr · July 10, 2020, 1:50pm

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 ?

cmd · July 10, 2020, 3:19pm

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.

sumantkr · July 10, 2020, 4:43pm

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

cmd · July 10, 2020, 5:28pm

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.
Try using the corrected code I’ve posted below which uses an l2 norm.

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

sumantkr · July 11, 2020, 2:16am

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:

cmd · July 11, 2020, 3:20am

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.

Topic		Replies	Views
Issue in plotting data with two components (Stokes problem) General Discussion	3	240	January 28, 2022
Why the plot function doesn't plot the velocity? General Discussion	2	260	March 15, 2023
Plot of diplacement field General Discussion	2	471	April 15, 2020
Need help in plotting graph of velocity (ux) vs time General Discussion	3	631	August 6, 2020
Elasticity neo-Hookean material General Discussion	0	248	February 17, 2022

Unable to display the plot value for velocity

Related topics