-- FreeFem++ v4.14 (mer. 06 mars 2024 16:59:04 CET - git v4.14-1-g2b2052ae)
   file : poiseuille.edp
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : //  poiseuille.edp
    2 : //
    3 : //  Discussion:
    4 : //
    5 : //    Solve the Navier Stokes equations defining Poiseuille flow in a channel.
    6 : //
    7 : //    This example was adapted from code at
    8 : //      https://modules.freefem.org/modules/stokes/
    9 : //
   10 : //  Licensing:
   11 : //
   12 : //    This code is distributed under the MIT license.
   13 : //
   14 : //  Modified:
   15 : //
   16 : //    12 February 2021
   17 : //
   18 : cout << "\n";
   19 : cout << "poiseuille():\n";
   20 : cout << "  FreeFem++ version.\n";
   21 : cout << "  Poiseuille flow through a channel.\n";
   22 : //
   23 : //  Parameters
   24 : //
   25 : real uMax = 10.;
   26 : 
   27 : real Mu = 1.;
   28 : 
   29 : func fx = 0.;
   30 : func fy = 0.;
   31 : //
   32 : //  Geometry
   33 : //
   34 : int nn = 4;	//Mesh quality
   35 : real L = 5.0;	//Pipe length
   36 : real D = 1.0;	//Pipe height
   37 : int Wall = 1;	//Pipe wall label
   38 : int Inlet = 2;	//Pipe inlet label
   39 : int Outlet = 3;	//Pipe outlet label
   40 : //
   41 : //                 Wall (y=D)
   42 : //           ------------------------------
   43 : //           :--->                        :
   44 : // Inlet     :----->                      :  Outlet
   45 : // (x=0)     :--->                        :  (x=L)
   46 : //           ------------------------------
   47 : //                 Wall (y=0)
   48 : //
   49 : border b1 ( t = 0.0, 1.0 ) { x=L*t;   y=0.0;   label=Wall; };
   50 : border b2 ( t = 0.0, 1.0 ) { x=L;     y=D*t;   label=Outlet; };
   51 : border b3 ( t = 0.0, 1.0 ) { x=L-L*t; y=D;     label=Wall; };
   52 : border b4 ( t = 0.0, 1.0 ) { x=0.0;   y=D-D*t; label=Inlet; };
   53 : 
   54 : int nnL = max ( 2.0, L * nn );
   55 : int nnD = max ( 2.0, D * nn );
   56 : //
   57 : //  Build the mesh.
   58 : //
   59 : mesh Th = buildmesh ( b1 ( nnL ) 
   60 :                     + b2 ( nnD ) 
   61 :                     + b3 ( nnL ) 
   62 :                     + b4 ( nnD ) );
   63 : //
   64 : //  Fespace
   65 : //  Set the finite element velocity space to piecewise quadratics.
   66 : //
   67 : fespace Uh ( Th, [ P2, P2 ] );
   68 : Uh [ux, uy];
   69 : Uh [vx, vy];
   70 : //
   71 : //  Set the finite element pressure space to piecewise linears.
   72 : //
   73 : fespace Ph(Th, P1);
   74 : Ph p;
   75 : Ph q;
   76 : //
   77 : //  Macros
   78 : //
   79 : macro Gradient(u) [dx(u), dy(u)]  )  //
   80 : macro Divergence(ux, uy) (dx(ux) + dy(uy))  )  //
   81 : //
   82 : //  Define a function for the parabolic inflow velocity.
   83 : //
   84 : func uIn = uMax * (1.0-(y-D/2.)^2 / (D/2.)^2);
   85 : //
   86 : //  Problem:
   87 : //  Solve for velocity u=(ux,uy) and pressure p,
   88 : //  using test functions v=(vx,vy) and q,
   89 : //  which zero out the following function.
   90 : //
   91 : problem S ( [ux, uy, p],[vx, vy, q])
   92 : 	= int2d(Th)(
   93 : 		  Mu * (
   94 : 			  Gradient(ux)    [dx(ux), dy(ux)] ' * Gradient(vx)    [dx(vx), dy(vx)] 
   95 : 			+ Gradient(uy)    [dx(uy), dy(uy)] ' * Gradient(vy)    [dx(vy), dy(vy)] 
   96 : 		)
   97 : 		- p * Divergence(vx, vy)     (dx(vx) + dy( vy)) 
   98 : 		- Divergence(ux, uy)     (dx(ux) + dy( uy))  * q
   99 : 	)
  100 : 	- int2d(Th)(
  101 : 		  fx*vx + fy*vy
  102 : 	)
  103 : 	+ on ( Inlet, ux = uIn, uy = 0.0 )
  104 : 	+ on ( Wall, ux = 0.0, uy = 0.0 )
  105 : 	;
  106 : //
  107 : //  Request the solution ( ux, uy, p )
  108 : //
  109 : S;
  110 : //
  111 : //  Plot the velocity vector field.
  112 : //
  113 : plot ( [ ux, uy ], nbarrow = 1, wait = true, 
  114 :   ps = "poiseuille_velocity.ps", 
  115 :   cmm = "Poiseuille Velocity" );
  116 : //
  117 : //  Plot the pressure as color contours.
  118 : //
  119 : plot ( p, nbiso = 10, fill = true, 
  120 :   ps = "poiseuille_pressure.ps", 
  121 :   cmm = "Poiseuille Pressure" );
  122 : //
  123 : //  Terminate.
  124 : //
  125 : cout << "\n";
  126 : cout << "poiseuille():\n";
  127 : cout << "  Normal end of execution.\n";
  128 : 
  129 : exit ( 0 );
  130 :  sizestack + 1024 =2064  ( 1040 )


poiseuille():
  FreeFem++ version.
  Poiseuille flow through a channel.
  --  mesh:  Nb of Triangles =    168, Nb of Vertices 109
  -- Solve : 
          min -9.04787e-15  max 10
          min -6.80489e-15  max 400

poiseuille():
  Normal end of execution.
  current line = 129
exit(0)
 err code 0 ,  mpirank 0
 CodeAlloc : nb ptr  4094,  size :531304 mpirank: 0
Ok: Normal End