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
Plot solutions functions of (x,y,t) General Discussion	240	June 18, 2021
Division issue on plot due to mesh? General Discussion	46	March 16, 2024
Plot function does't work in a loop General Discussion	216	November 27, 2021
Trying to understand difference between 3D and 2D plot General Discussion	119	October 26, 2023
Plot an function General Discussion	342	June 16, 2020

Why the plot function doesn't plot the velocity?

Related Topics