Why the plot function doesn't plot the velocity?

cfdyuewu · March 14, 2023, 12:47pm

the plot([uite1, uite2], pite, coef=0.2, cmm=" [Ux, Uy] - p", WindowIndex=1);
only plot the result of the first step, then it only plot the contour of the region, how to fix this?, the full code is

load "iovtk"
verbosity=0;
// Parameters
int nn = 15;
int nnPlus = 5;
real l = 1.;
real L = 15.;
real hSlope = 0.1;
real H = 6.;
real h = 0.5;

real beta = 0.01;

real eps = 9.81/303;
real nu = 1.;
real numu = nu/sqrt(0.09);
real nuep = pow(nu,1.5)/4.1;
real dt = 0.;
// Mesh
border b1(t=0, l){x=t; y=0;}
border b2(t=0, L-l){x=1.+t; y=-hSlope*t;}
border b3(t=-hSlope*(L-l), H){x=L; y=t;}
border b4(t=L, 0){x=t; y=H;}
border b5(t=H, h){x=0; y=t;}
border b6(t=h, 0){x=0; y=t;}

mesh Th=buildmesh(b1(nnPlus*nn*l) + b2(nn*sqrt((L-l)^2+(hSlope*(L-l))^2)) + b3(nn*(H + hSlope*(L-l))) + b4(nn*L) + b5(nn*(H-h)) + b6(nnPlus*nn*h));
plot(Th);

// Fespaces
fespace Vh2(Th, P1b);
Vh2 Ux, Uy,uite1,uite2,Uxst,Uyst;
Vh2 Vx, Vy,v1,v2;
Vh2 Upx, Upy=0;
fespace Vh(Th,P1);
Vh p=0, q,pp=p,pite,vp,pst;
Vh Tp, T=35,Tite,vt,Tst;
Vh k=log(0.27), kp=k,kite,vk,kst;
Vh ep=log(0.002104), epp=ep,epite,vep,epst;
 
fespace V0h(Th,P0);
 V0h muT=1;
 V0h prodk, prode;
 Vh kappa=0.25e-4, stress;
// Macro
 macro grad(u) [dx(u), dy(u)] //
 macro Grad(u1,u2) [dx(u1),dy(u1), dx(u2),dy(u2)] //
 macro Gradt(u1,u2) [dx(u1),dx(u2), dy(u1),dy(u2)] //
 macro Div(u1,u2) (dx(u1) + dy(u2)) //
macro UgradV(u1,u2,v1,v2) [[u1,u2]'*[dx(v1),dy(v1)],[u1,u2]'*[dx(v2),dy(v2)]]//
macro ugradv(u1,u2,v) [u1,u2]'*grad(v)//
// Problem
real alpha = 0.;


// Initialization with stationary solution


// Initialization
dt = 0.05; 
Tp = T - 10*((x<1)*(y<0.5) + (x>=1)*(y+0.1*(x-1)<0.5));
Upx=0+(x==0)*(y>=0)*(y<=h)*3;
plot(Upx,  cmm=" Ux ", wait=1);
plot(Upy,  cmm=" Ux ", wait=1);
// Convergence loop
real T0 = clock();
cout  << " - dt = " << dt << endl;
   alpha = 1/dt;

   // Time loop
   real t = 0.;
   for (int i = 0; i <= 500; i++){
      uite1=Upx;
      uite2=Upy;
      pite=pp;
      Tite=Tp;
      kite=kp;
      epite=epp;
      t=t+dt;
      cout << "Time step " << i << " - t = " << t << endl;
      for (int n=0;n<=4;n++){
         solve xstar ([Uxst,Uyst,pst,Tst,kst,epst],[v1,v2,vp,vt,vk,vep])
         =int2d(Th,optimize=0)(alpha*([uite1,uite2]'*[v1,v2]-Tite*vt-kite*vk-epite*vep)
         +(UgradV(uite1,uite2,uite1,uite2))'*[v1,v2]
         +(UgradV(uite1,uite2,uite1,uite2))'*[v1,v2]
         +ugradv(v1,v2,pite)
         +(0.09)*exp(2*kite-epite)*(2*kite-epite)*(Grad(uite1,uite2)+(Gradt(uite1,uite2)))'*Grad(v1,v2)
         +(0.09)*exp(2*kite-epite)*(Grad(uite1,uite2)+Gradt(uite1,uite2))'*Grad(v1,v2)
         +(9.81/303)*Tite*v2-Div(uite1,uite2)*vp-vk*(ugradv(uite1,uite2,kite)+ugradv(uite1,uite2,kite))-exp(epite-kite)*(epite-kite)*vk+(0.09)*exp(2*kite-epite)*(2*kite-epite)*[dx(kite),dy(kite)]'*[dx(kite),dy(kite)]*vk
         +(0.09)*exp(2*kite-epite)*2*[dx(kite),dy(kite)]'*[dx(kite),dy(kite)]*vk
         -(0.09)*exp(2*kite-epite)*(2*kite-epite)*[dx(vk),dy(vk)]'*[dx(kite),dy(kite)]-(0.09)*exp(2*kite-epite)*[dx(vk),dy(vk)]'*[dx(kite),dy(kite)]
         +(0.09)*exp(kite-epite)*(4*dx(uite1)*dx(uite1)+4*dy(uite2)*dy(uite2)+2*(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vk
         +(0.09)*exp(kite-epite)*(kite-epite)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vk
         -vep*(ugradv(uite1,uite2,epite)+ugradv(uite1,uite2,epite))
         -(1.92)*exp(epite-kite)*(epite-kite)*vep
         +(0.09/1.3)*exp(2*kite-epite)*(2*kite-epite)*[dx(epite),dy(epite)]'*[dx(epite),dy(epite)]*vep
         +(0.09/1.3)*exp(2*kite-epite)*2*[dx(epite),dy(epite)]'*[dx(epite),dy(epite)]*vep
         -(0.09/1.3)*exp(2*kite-epite)*(2*kite-epite)*[dx(epite),dy(epite)]'*[dx(vep),dy(vep)]
         -(0.09/1.3)*exp(2*kite-epite)*[dx(epite),dy(epite)]'*[dx(vep),dy(vep)]
         +((1.44*0.09))*exp(kite-epite)*(4*dx(uite1)*dx(uite1)+4*dy(uite2)*dy(uite2)+2*(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vep
         +((1.44*0.09))*exp(kite-epite)*(kite-epite)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vep
         -vt*(ugradv(uite1,uite2,kite)+ugradv(uite1,uite2,kite))
         -0.1*exp(2*kite-epite)*(2*kite-epite)*[dx(vt),dy(vt)]'*[dx(Tite),dy(Tite)]
         -0.1*exp(2*kite-epite)*[dx(vt),dy(vt)]'*[dx(Tite),dy(Tite)]
         )
         +int1d(Th,b1,b2,b5,optimize=0)(
            (0.09/1.3)*2.4952*exp(0.5*kite)*(0.5*kite)*vep
         )+int1d(Th,b1,b2,optimize=0)(
            0.41*((0.09)^0.25)*(exp(0.5*kite)*(0.5*kite)*uite1*v1+exp(0.5*kite)*uite1*v1)
         )-int2d(Th,optimize=0)(
            alpha*([uite1-Upx,uite2-Upy]'*[v1,v2]-(Tite-Tp)*vt-(kite-kp)*vk-(epite-epp)*vep)
            +(UgradV(uite1,uite2,uite1,uite2))'*[v1,v2]
            +(grad(pite))'*[v1,v2]
            +0.09*exp(2*kite-epite)*(Grad(uite1,uite2)+(Gradt(uite1,uite2)))'*Grad(v1,v2)
            +(9.81/303)*(Tite-35)*v2-(Div(uite1,uite2))*vp
            -(ugradv(uite1,uite2,kite))*vk
            -exp(epite-kite)*vk
            +(0.09)*vk*(exp(2*kite-epite)*grad(kite)'*grad(kite))
            -(0.09)*(exp(2*kite-epite)*grad(kite)'*grad(vk))
            +(0.09)*exp(kite-epite)*vk*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))
            -(ugradv(uite1,uite2,epite))*vep
            -1.92*exp(epite-kite)*vep
            +(0.09/1.3)*vep*(exp(2*kite-epite)*grad(epite)'*grad(epite))
            -(0.09/1.3)*(exp(2*kite-epite)*grad(epite)'*grad(vep))
            +(1.44*0.09)*exp(kite-epite)*vep*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))
            -(ugradv(uite1,uite2,Tite))*vt
            -0.1*(exp(2*kite-epite)*grad(Tite)'*grad(vt))
         )-int1d(Th,b1,b2,b5, optimize=0)(
            (0.09/2.6)*2.4952*exp(0.5*kite)*vep+(0.09/2.6)*2.4952*exp(0.5*kp)*vep)
         -int1d(Th,b1,b2, optimize=0)(0.41*((0.09)^0.25)*exp(0.5*kite)*uite1*v1)
         -int2d(Th, optimize=0)(alpha*([Uxst,Uyst]'*[v1,v2]-Tst*vt-kst*vk-epst*vep)
         +(UgradV(Uxst,Uyst,uite1,uite2))'*[v1,v2]
         +(UgradV(uite1,uite2,Uxst,Uyst))'*[v1,v2]
         +ugradv(v1,v2,pst)
         +(0.09)*exp(2*kite-epite)*(2*kst-epst)*(Grad(uite1,uite2)+(Gradt(uite1,uite2)))'*Grad(v1,v2)
         +(0.09)*exp(2*kite-epite)*(Grad(Uxst,Uyst)+Gradt(Uxst,Uyst))'*Grad(v1,v2)
         +(9.81/303)*Tst*v2-Div(Uxst,Uyst)*vp-vk*(ugradv(Uxst,Uyst,kite)+ugradv(uite1,uite2,kst))-exp(epite-kite)*(epst-kst)*vk+(0.09)*exp(2*kite-epite)*(2*kst-epst)*[dx(kite),dy(kite)]'*[dx(kite),dy(kite)]*vk
         +(0.09)*exp(2*kite-epite)*2*[dx(kite),dy(kite)]'*[dx(kst),dy(kst)]*vk
         -(0.09)*exp(2*kite-epite)*(2*kst-epst)*[dx(vk),dy(vk)]'*[dx(kite),dy(kite)]-(0.09)*exp(2*kite-epite)*[dx(vk),dy(vk)]'*[dx(kst),dy(kst)]
         +(0.09)*exp(kite-epite)*(4*dx(uite1)*dx(Uxst)+4*dy(uite2)*dy(Uyst)+2*(dy(uite1)+dx(uite2))*(dy(Uxst)+dx(Uyst)))*vk
         +(0.09)*exp(kite-epite)*(kst-epst)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vk
         -vep*(ugradv(Uxst,Uyst,epite)+ugradv(uite1,uite2,epst))
         -(1.92)*exp(epite-kite)*(epst-kst)*vep
         +(0.09/1.3)*exp(2*kite-epite)*(2*kst-epst)*[dx(epite),dy(epite)]'*[dx(epite),dy(epite)]*vep
         +(0.09/1.3)*exp(2*kite-epite)*2*[dx(epite),dy(epite)]'*[dx(epst),dy(epst)]*vep
         -(0.09/1.3)*exp(2*kite-epite)*(2*kst-epst)*[dx(epite),dy(epite)]'*[dx(vep),dy(vep)]
         -(0.09/1.3)*exp(2*kite-epite)*[dx(epst),dy(epst)]'*[dx(vep),dy(vep)]
         +((1.44*0.09))*exp(kite-epite)*(4*dx(uite1)*dx(Uxst)+4*dy(uite2)*dy(Uyst)+2*(dy(uite1)+dx(uite2))*(dy(Uxst)+dx(Uyst)))*vep
         +((1.44*0.09))*exp(kite-epite)*(kst-epst)*(2*dx(uite1)*dx(uite1)+2*dy(uite2)*dy(uite2)+(dy(uite1)+dx(uite2))*(dy(uite1)+dx(uite2)))*vep
         -vt*(ugradv(Uxst,Uyst,kite)+ugradv(uite1,uite2,kst))
         -0.1*exp(2*kite-epite)*(2*kst-epst)*[dx(vt),dy(vt)]'*[dx(Tite),dy(Tite)]
         -0.1*exp(2*kite-epite)*[dx(vt),dy(vt)]'*[dx(Tst),dy(Tst)]
         )
         -int1d(Th,b1,b2,b5, optimize=0)(
            (0.09/1.3)*2.4952*exp(0.5*kite)*(0.5*kst)*vep
         )-int1d(Th,b1,b2, optimize=0)(
            0.41*((0.09)^0.25)*(exp(0.5*kite)*(0.5*kst)*uite1*v1+exp(0.5*kite)*Uxst*v1)
         )
         +on(b5,Uxst=0,Uyst=0)
         +on(b6,Uxst=3,Uyst=0,kst=log(0.27),epst=log(0.002104),Tst=25)+on(b2,Uyst=-Upx*N.x/N.y,Tst=30)+on(b1,Uyst=0,Tst=30);
         uite1=Uxst;
         uite2=Uyst;
         pite=pst;
         Tite=Tst;
         kite=kst;
         epite=epst;
      }
      Upx=uite1;
      Upy=uite2;
      pp=pite;
      Tp=Tite;
      kp=kite;
      epp=epite;
      plot([uite1, uite2], pite, coef=0.2, cmm=" [Ux, Uy] - p", WindowIndex=1);


   }
cout << "Total Time = " << clock()-T0 << endl;

ceb · March 15, 2023, 1:24pm

Does the problem occur with a simpler (to follow) model.
I guess it is best if we include a simple example to illustrate the issue, rather than an example with lots of code lines which may not add much info to reproduce the issue.
Regards

cmd · March 15, 2023, 6:11pm

This works, but it is an extremely inefficient implementation. If you coarsen up the mesh a bit, you will see the plots you want.

Topic		Replies	Views
Unable to display the plot value for velocity General Discussion	9	502	July 11, 2020
Plot solutions functions of (x,y,t) General Discussion	0	288	June 18, 2021
Division issue on plot due to mesh? General Discussion	0	94	March 16, 2024
Plot function does't work in a loop General Discussion	0	244	November 27, 2021
Need help in plotting graph of velocity (ux) vs time General Discussion	3	631	August 6, 2020

Why the plot function doesn't plot the velocity?

Related topics