# Need help in plotting graph of velocity (ux) vs time

**URL:** https://community.freefem.org/t/need-help-in-plotting-graph-of-velocity-ux-vs-time/529
**Category:** General Discussion
**Created:** [August 2, 2020, 9:22am UTC](https://community.freefem.org/t/need-help-in-plotting-graph-of-velocity-ux-vs-time/529 "2020-08-02T09:22:14Z")
**Posts on this page:** 4
**Page:** 1

<div class="post-metadata">

### Author: ![Rhushikesh](https://avatars.discourse-cdn.com/v4/letter/r/e19b73/32.png) [@Rhushikesh](https://community.freefem.org/u/Rhushikesh)
#### Post date: [August 2, 2020, 9:22am UTC](https://community.freefem.org/t/need-help-in-plotting-graph-of-velocity-ux-vs-time/529/1 "2020-08-02T09:22:14Z")

</div>

```
	///////////////////////////////////////////////////
	/// UNSTEADY 2D NAVIER--STOKES ////
	///////////////////////////////////////////////////

```

///////////////////////////////////////////////////////////////////  
/////////// MESH /////////////

// Mesh Trial //  
real DD = 1.0;  
real WW = 10.0;  
// CAVITY Sipp LEBEDEV 2007 :: Trial //

///Thick line object  
border upstream1 (t =-0.6,0.0) {x = t; y = 0.0; label = 2;}; //solid thick line  
border cavityleft1 (t =0.0,-0.1) {x = 0.0; y = t; label = 2;}; //solid thick line  
border cavitybottom1 (t =0.0,1.0) {x=t;y=-0.1;label =2;}; //solid thick line  
border cavityright1 (t =-0.1,0.0) {x = 1.0; y = t; label = 2;}; // solid thick line  
border downstream1 (t =1.0, 2.50) {x = t; y = 0.0; label = 2;}; //solid thick line  
border outlet1 (t =0.0,0.1) {x = 2.5; y = t; label = 3;}; // solid thick line  
border top1 (t =2.5,-0.6) {x=t;y=0.1; label = 999;}; //solid thick line  
border inlet1 (t =0.1,0.0) {x=-0.6;y=t; label = 1;}; // solid thick line

//thin line bottom  
border cavityleft2 (t =-0.1,-0.2) {x = 0.0; y = t; label = 2;}; //solid thin line  
border cavitybottom2 (t =0.0,1.0) {x=t;y=-0.20;label =2;}; //solid thin line  
border cavityright2 (t =-0.2,-0.1) {x = 1.0; y = t; label = 2;}; //solid thin line

//dashed bottom

border cavityleft3 (t =-0.2,-0.35) {x = 0.0; y = t; label = 2;}; // dashed  
border cavitybottom3 (t =0.0,1.0) {x=t;y=-0.35;label = 2;}; //dashed  
border cavityright3 (t =-0.35,-0.2) {x = 1.0; y = t; label = 2;}; //dashed

//dotted bottom  
border cavityleft4 (t =-0.35,-1.0) {x = 0.0; y = t; label = 2;}; // dotted  
border cavitybottom4 (t =0.0,1.0) {x = t; y = -1.0; label = 2;}; //dotted  
border cavityright4 (t =-1.0,-0.35) {x = 1.0; y = t; label = 2;}; //dotted

//thin line top  
border outlet5 (t =0.1,0.15) {x = 2.5; y = t; label = 3;}; //solid thin line  
border top5 (t =2.5,-1.2) {x=t;y=0.15;label = 4;}; //solid thin line  
border inlet5 (t =0.15,0.0) {x=-1.2;y=t;label=1;}; // solid thin line  
border upstream5 (t =-1.2,-0.6) {x = t; y = 0.0; label = 2;}; //solid thin line

//dashed top  
border outlet6 (t =0.15,0.30) {x = 2.5; y=t; label = 3;}; //dashed  
border top6 (t =2.5,-1.2) {x=t;y=0.3; label = 4;}; //dashed  
border inlet6 (t =0.3,0.15) {x=-1.2; y=t; label=1;}; //dashed

//dotted top  
border outlet7 (t = 0.30,0.50) {x = 2.5; y=t; label=3;}; //dotted  
border top7 (t = 2.5, -1.2) {x = t; y = 0.5; label = 4;}; //dotted  
border inlet7 (t = 0.50,0.30) {x = -1.2; y = t; label = 1;}; //dotted

// For refinement while considering a large domain size

mesh th = buildmesh( upstream1(350)  
+ cavityleft1(350)  
+ cavitybottom1(350)  
+ cavityright1(350)  
+ downstream1(350)  
+ outlet1(350)  
+ top1(350)  
+ inlet1(350)

```
                + cavityleft2(150)
                + cavitybottom2(150)
                + cavityright2(150)

                + cavityleft3(100)
                + cavitybottom3(100)
                + cavityright3(100)

                + cavityleft4(50)
                + cavitybottom4(50)
                + cavityright4(50)

                + outlet5(200)
                + top5(200)
                + inlet5(200)
                + upstream5(200)

                + outlet6(100)
                + top6(100)
                + inlet6(100)

                + outlet7(100)
                + top7(100)
                + inlet7(100)
                
                );

```

plot(th, value = 1, wait = 1);

//////////////////////////////////////////////////////////////////  
//////////// FINITE ELEMENT SPACES ////////////////

/////////\*

fespace Uh(th, P2); //velocities  
fespace Ph(th, P1); //pressure  
fespace UhxUhxPh(th,[P2,P2,P1]); // 2 velocity components and 1 pressure field

//////////////////////////////////////////////////////////////////  
////////// NUMERICAL PARAMETERS ////////////////

real u0 = 1.0; // free-stream velocity  
real diam = 1.0; // diameter of the cylinder, length scale

real Re = 100.0; // Reynolds number  
real dt = 0.01; // time-step  
real nu = u0\*diam/Re;

int rep=0; // To re-start the simulation using a known flow field  
int stepvisu=1000; // To visualise the flow fields  
int stepsave=1000; // To save the flow fields  
int vis=0; // 0: don’t wait; 1: wait  
int iter=0; // if rep neq 0, the value of the iteration corresponding to the starting flow field

////////////////////////////////////////////////////////////////////  
/////// INLET CONDITIONS (if needed) //////////

// Thickness of the boundary layer at the inlet  
real delta=0.2;

func real uinit(real x, real y)  
{  
if(y\<0)  
return 0;  
else  
return u0\*tanh((y)/delta);  
};

//////////////////////////////////////////////////////////////////////  
/////// VELOCITY FIELDS //////////

UhxUhxPh [ux,uy,p]; //Flow field for the present time step

UhxUhxPh [vx,vy,q]; // Test functions

UhxUhxPh [upx,upy,pp]; // Flow field for the previous time step

UhxUhxPh [rhs1,rhs2,rhs3]; // Right hand side of the NS formed by the terms from the previous time step

//////////////////////////////////////////////////////////////////////  
//////// VARIATIONAL FORMULATION ////////////////

// LINEAR TERMS ///  
varf NS ([ux,uy,p],[vx,vy,q])=  
int2d(th)(  
(ux_vx+uy_vy)/dt  
+ nu\*( dx(ux)\*dx(vx)+dy(ux)\*dy(vx)  
+ dx(uy)_dx(vy)+dy(uy)dy(vy) )  
- p( dx(vx) + dy(vy) )  
+ q_( dx(ux) + dy(uy) ) )

// NO-STRESS CONDITION  
+int1d(th,3)(nu\*( (N.y_dy(ux)+N.x_dx(ux))_vx)  
- ( (N.y_vy + N.x\*vx )\*p) )

- on(1,uy=0.0,ux=u0)
- on(2,uy=0.0, ux = 0.0)
- on(3,uy=0.0)
- on(4,uy=0.0);  
// Boundary cond.

// NONLINEAR TERMS

// The nonlinear terms: upx upy stand for flow fields from the previous time step //

// you should also go through the manual where they have discussed “varf” in detail…

varf NSRITE ([ux,uy,p],[vx,vy,q])=  
int2d(th)( (upx_vx+upy_vy)/dt  
- ( upy\*dy(upy)_vy + upx_dx(upy)_vy  
+ upy_dy(upx)_vx + upx_dx(upx)\*vx )  
)

+int1d(th,3)(nu\*( (N.y_dy(ux)+N.x_dx(ux))_vx)  
- ( (N.y_vy + N.x\*vx )\*p) )

// same boundary conditions of course apply to the convective nonlinear terms too…

- on(1,uy=0.0,ux=u0)
- on(2,uy=0.0, ux = 0.0)
- on(3,uy=0.0)
- on(4,uy=0.0);  
// Boundary cond.

//////////////////////////////////////////////////////////////////////////  
///////// INITIAL FLOW FIELDS //////////////////

real tps;

// Start from a scratch solution  
tps = 0.0;

if(rep == 0)  
{  
[ux,uy,p] = [0.0,0.0,0.0]; // good starting flow field in general!  
plot(ux,value=1,wait=vis);  
plot(uy,value=1,wait=vis);  
plot(p,value=1,wait=vis);  
}  
else  
{  
// Start from a previous solution  
tps=rep\*dt;  
{  
ifstream file(“sol-20-2000\_2D\_DNS\_CYL\_dt\_0p01\_Re\_20\_mesh1.txt”);  
file \>\> ux[];  
};  
plot(ux,value=1,wait=vis);  
plot(uy,value=1,wait=vis);  
plot(p,value=1,wait=vis);  
};

////////////////////////////////////////////////////////  
//////// MATRIX FORMULATION //////////////

// varf NS ([ux,uy,p],[vx,vy,q])

matrix NSMAT = NS(UhxUhxPh,UhxUhxPh,solver=UMFPACK);

////////////////////////////////////////////////////////  
//////// TIME MARCHING //////////////////

for (int i=1; i\<=10000; i++)  
{

```
[upx,upy,pp] = [ux,uy,p];

rhs1[] = NSRITE(0,UhxUhxPh);
ux[] = NSMAT^-1*rhs1[];

tps=tps+dt;
iter = iter + 1;

cout << " i = " << i << "\n";
cout << " tps = " << tps << "\n";
cout << " UX = " << ux(0.5, 0.5) << "\n";
cout << " UY = " << uy(0.5, 0.5) << "\n";

// Vizualise the flow fields
if (!(i%stepvisu))
{
    plot(p,value=1,wait=vis);
    plot(ux,value=1,wait=vis);
    plot(uy,value=1,wait=vis);
  // plot([UX,tps],value=1,wait=vis); //
};

// Save the flow fields
if(!(i%stepsave))
{
    ofstream file("sol-"+Re+"-"+iter+"_2D_DNS_CYL_dt_0p01_Re_46p6_Mesh1.txt");
    file << ux[] << endl;
}
///// For time traces at certain points within the domain /////

ofstream uxvel1("ux_2D_DNS_CAV_TEST.dat",append);
ofstream uyvel1("uy_2D_DNS_CAV_TEST.dat",append);

uxvel1<< tps << " " << ux(0.5,0.5) <<endl ;
uyvel1<< tps << " " << uy(0.5,0.5) <<endl ;

};
// End of the time loop

plot(p,value=1,wait=vis);
plot(ux,value=1,wait=vis);
plot(uy,value=1,wait=vis);

```

///////////////////////////////////////////////  
/////// POST-PROCESSING ////////////

```
fespace PPH(th,P1);
PPH UX, UY, VORT;
UX = ux;
UY = uy;

Uh vort;
vort = dx(uy) - dy(ux);
VORT = vort;

```

//////

```
///// MATLAB ///////

include "../../ffmatlib/ffmatlib.idp"

savemesh(th,"th_Mesh1.msh");
ffSaveVh(th,PPH,"fespace_PPH_P1.txt");
ffSaveData(UX,"UX_2D_DNS_CYL_"+iter+"_dt_0p01_Re_46p6_Mesh1.txt");
ffSaveData(UY,"UY_2D_DNS_CYL_"+iter+"_dt_0p01_Re_46p6_Mesh1.txt");
ffSaveData(VORT,"VORT_2D_DNS_CYL_"+iter+"_dt_0p01_Re_46p6_Mesh1.txt");

//////////////////////////
```

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [August 2, 2020, 10:24am UTC](https://community.freefem.org/t/need-help-in-plotting-graph-of-velocity-ux-vs-time/529/2 "2020-08-02T10:24:49Z")

</div>

Use `load "PETSc"` + `savevtk(..., append = true);`, see [this](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/hpddm/laplace-adapt-dist-3d-PETSc.edp#L54) example.  
It will generate a `.pvd` file that you can open with ParaView, and you’ll have a movie.

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [August 2, 2020, 2:34pm UTC](https://community.freefem.org/t/need-help-in-plotting-graph-of-velocity-ux-vs-time/529/4 "2020-08-02T14:34:24Z")

</div>

What did you write and what is not working?

---

<div class="post-metadata">

### Author: ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)
#### Post date: [August 6, 2020, 10:17am UTC](https://community.freefem.org/t/need-help-in-plotting-graph-of-velocity-ux-vs-time/529/6 "2020-08-06T10:17:03Z")

</div>

Sorry, I don’t use MATLAB. Again, I’d suggest using ParaView.