-- FreeFem++ v4.6 (Thu Apr  2 15:47:38 CEST 2020 - git v4.6)
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : //  traffic.edp
    2 : //
    3 : //  Discussion:
    4 : //
    5 : //    Pseudo two−dimensional vehicle trafic flow with road branching.
    6 : //
    7 : //    d/dt rho + del ( rho V(rho) n ) - nu * del^2 rho = 0
    8 : //    V(rho) = q(rho) / rho 
    9 : //    q(rho) = rho ∗ Vfree ∗ (1−rho/rhoc)
   10 : //
   11 : //    Here "n" is the unit normal vector.
   12 : //
   13 : //  Location:
   14 : //
   15 : //    http://people.sc.fsu.edu/~jburkardt/freefem_src/traffic/traffic.edp
   16 : //
   17 : //  Modified:
   18 : //
   19 : //    27 June 2015
   20 : //
   21 : //  Author:
   22 : //
   23 : //    Florian De Vuyst
   24 : //
   25 : //  Reference:
   26 : //
   27 : //    Florian De Vuyst,
   28 : //    Numerical modeling of transport problems using freefem++ software -
   29 : //    with examples in biology, CFD, traffic flow and energy transfer,
   30 : //    HAL id: cel-00842234
   31 : //    https://cel.archives-ouvertes.fr/cel-00842234
   32 : //
   33 : //  Parameters:
   34 : //
   35 : //    real DT, the time step.
   36 : //
   37 : //    real NU, the viscosity.
   38 : //
   39 : //    real T, current time.
   40 : //
   41 : cout << "\n";
   42 : cout << "traffic\n";
   43 : cout << "  FreeFem++ version\n";
   44 : cout << "  Simulate traffic flow.\n";
   45 : 
   46 : real [int] A(2), B(2), C(2), D(2), E(2);
   47 : real [int] F(2), G(2), H(2), I(2), J(2);
   48 : real [int] K(2), L(2), M(2), Ng(2);
   49 : A = [0,0]; B = [9, 1.5]; C = [12, 1.5]; D = [21,0];
   50 : E = [6, 2]; F = [15, 2]; G = [0, 4]; H = [9, 2.5];
   51 : I = [12, 2.5]; J = [21, 4]; K = [0, 1]; L = [0, 3];
   52 : M = [21, 1]; Ng = [21, 3];
   53 : 
   54 : border c1(t=0,1) {x=(1-t)*A[0]+t*B[0]; y=(1-t)*A[1]+t*B[1];}
   55 : border c2(t=0,1) {x=(1-t)*B[0]+t*C[0]; y=(1-t)*B[1]+t*C[1];}
   56 : border c3(t=0,1) {x=(1-t)*C[0]+t*D[0]; y=(1-t)*C[1]+t*D[1];}
   57 : border c4(t=0,1) {x=(1-t)*G[0]+t*L[0]; y=(1-t)*G[1]+t*L[1];}
   58 : border c5(t=0,1) {x=(1-t)*L[0]+t*E[0]; y=(1-t)*L[1]+t*E[1];}
   59 : border c6(t=0,1) {x=(1-t)*E[0]+t*K[0]; y=(1-t)*E[1]+t*K[1];}
   60 : border c7(t=0,1) {x=(1-t)*K[0]+t*A[0]; y=(1-t)*K[1]+t*A[1];}
   61 : border c8(t=0,1) {x=(1-t)*D[0]+t*M[0]; y=(1-t)*D[1]+t*M[1];}
   62 : border c9(t=0,1) {x=(1-t)*M[0]+t*F[0]; y=(1-t)*M[1]+t*F[1];}
   63 : border c10(t=0,1) {x=(1-t)*F[0]+t*Ng[0]; y=(1-t)*F[1]+t*Ng[1];}
   64 : border c11(t=0,1) {x=(1-t)*Ng[0]+t*J[0]; y=(1-t)*Ng[1]+t*J[1];}
   65 : border c12(t=0,1) {x=(1-t)*J[0]+t*I[0]; y=(1-t)*J[1]+t*I[1];}
   66 : border c13(t=0,1) {x=(1-t)*I[0]+t*H[0]; y=(1-t)*I[1]+t*H[1];}
   67 : border c14(t=0,1) {x=(1-t)*H[0]+t*G[0]; y=(1-t)*H[1]+t*G[1];}
   68 : //
   69 : //  Define the mesh.
   70 : //
   71 : mesh Th = buildmesh ( 
   72 :     c1(60) + c2(30) + c3(60) + c4(6) + c5(50) + c6(50) + c7(6) + c8(6) 
   73 :   + c9(50) + c10(50) + c11(6) + c12(60) + c13(30) + c14(60) );
   74 : 
   75 : plot ( Th, ps = "traffic_mesh.ps" );
   76 : //
   77 : //  Define the finite element spaces.
   78 : //
   79 : fespace Vh(Th, P1);
   80 : Vh rho;
   81 : Vh rhoold;
   82 : Vh sigma;
   83 : Vh sigmaold, rhoh, sigmah, p, ph;
   84 : Vh n1, n2, u1, u2, u1visu, u2visu, vrho;
   85 : fespace Wh(Th, P2);
   86 : Wh phi, phih;
   87 : //
   88 : real rho0 = 25; // [nb cars / km]
   89 : real rhoc = 300; // critical density
   90 : 
   91 : real vfree = 110; // [km / h]
   92 : real dt = 0.004; 
   93 : real nu = 10.0;
   94 : real t = 0.0;
   95 : //
   96 : // Initial field
   97 : //
   98 : rho = rho0;
   99 : //
  100 : // Step 0. Define a velocity unit vector by solving an independent Laplace problem. 
  101 : //
  102 : problem Laplace ( phi, phih ) =
  103 :   int2d (Th) ( dx(phi)*dx(phih) + dy(phi)*dy(phih) )
  104 : + int1d (Th, c4) (phih) 
  105 : + int1d (Th, c7) (phih)
  106 : + on ( c8, c11, phi=0 );
  107 : //
  108 : //  Solve the problem for PHI.
  109 : //
  110 : Laplace;
  111 : //
  112 : //  Get the unit vector [N1,N2] of the gradient of the solution PHI.
  113 : //
  114 : n1 = dx(phi) / sqrt ( dx(phi)^2 + dy(phi)^2 );
  115 : n2 = dy(phi) / sqrt ( dx(phi)^2 + dy(phi)^2 );
  116 : 
  117 : rhoold = rho;
  118 : vrho = vfree * (1.0-rho / rhoc);
  119 : u1 = vrho * n1;
  120 : u2 = vrho * n2;
  121 : 
  122 : problem step1 ( rho, rhoh ) =
  123 :   int2d (Th) ( rho*rhoh/dt )
  124 : - int2d (Th) ( convect([u1,u2], -dt, rhoold)*rhoh/dt )
  125 : + int2d (Th) ( nu*dx(rho)*dx(rhoh) + nu*dy(rho)*dy(rhoh) )
  126 : + int2d (Th) ( rho*dx(u1)*rhoh + rho*dy(u2)*rhoh )
  127 : + on (c4, rho=rho0 )
  128 : + on (c7, rho=rhoc );
  129 : 
  130 : for ( int it = 0 ; it < 200 ; it++ )
  131 : {
  132 :   t = t + dt;
  133 :   step1;
  134 : 
  135 :   rhoold = rho;
  136 :   vrho = vfree * ( 1.0 - rho / rhoc);
  137 :   u1 = vrho * n1;
  138 :   u2 = vrho * n2;
  139 : 
  140 :   plot ( rho, nbiso = 50, fill = 1, value = 1, wait = 0 );
  141 : }
  142 : 
  143 : cout << "Final time = " << t << endl;
  144 : 
  145 : plot ( rho, nbiso = 40, fill=1, value=1, ps = "traffic_rho.ps", 
  146 :   cmm = "traffic density" );
  147 : u1visu = u1 / vfree;
  148 : u2visu = u2 / vfree;
  149 : plot ( [ u1visu, u2visu ], value=1, ps = "traffic_velocity.ps",
  150 :   cmm = "traffic velocity" );
  151 : //
  152 : //  Terminate.
  153 : //
  154 : cout << "\n";
  155 : cout << "traffic:\n";
  156 : cout << "  Normal end of execution.\n";
  157 :  sizestack + 1024 =8216  ( 7192 )


traffic
  FreeFem++ version
  Simulate traffic flow.
  --  mesh:  Nb of Triangles =   4132, Nb of Vertices 2329
  -- Solve : 
          min -26.7999  max -1.55491e-32
  -- Solve : 
          min 22.3024  max 300
  -- Solve : 
          min 20.4422  max 300
  -- Solve : 
          min 18.8348  max 300
  -- Solve : 
          min 17.4596  max 300
  -- Solve : 
          min 16.2402  max 300
  -- Solve : 
          min 15.0757  max 300
  -- Solve : 
          min 13.9872  max 300
  -- Solve : 
          min 13.162  max 300
  -- Solve : 
          min 12.5872  max 300
  -- Solve : 
          min 12.2255  max 300
  -- Solve : 
          min 12.0376  max 300
  -- Solve : 
          min 11.9942  max 300
  -- Solve : 
          min 12.0951  max 300
  -- Solve : 
          min 12.3519  max 300
  -- Solve : 
          min 12.7549  max 300
  -- Solve : 
          min 12.8045  max 300
  -- Solve : 
          min 12.8652  max 300
  -- Solve : 
          min 12.9192  max 300
  -- Solve : 
          min 12.9774  max 300
  -- Solve : 
          min 13.0325  max 300
  -- Solve : 
          min 13.0898  max 300
  -- Solve : 
          min 13.1448  max 300
  -- Solve : 
          min 13.1898  max 300
  -- Solve : 
          min 12.5899  max 300
  -- Solve : 
          min 12.1341  max 300
  -- Solve : 
          min 11.8662  max 300
  -- Solve : 
          min 11.7729  max 300
  -- Solve : 
          min 11.8386  max 300
  -- Solve : 
          min 12.0459  max 300
  -- Solve : 
          min 12.3761  max 300
  -- Solve : 
          min 12.8097  max 300
  -- Solve : 
          min 13.3264  max 300
  -- Solve : 
          min 13.9061  max 300
  -- Solve : 
          min 14.5299  max 300
  -- Solve : 
          min 15.1806  max 300
  -- Solve : 
          min 15.8431  max 300
  -- Solve : 
          min 16.5047  max 300
  -- Solve : 
          min 17.1551  max 300
  -- Solve : 
          min 17.7863  max 300
  -- Solve : 
          min 18.3925  max 300
  -- Solve : 
          min 18.9702  max 300
  -- Solve : 
          min 19.5176  max 300
  -- Solve : 
          min 20.0352  max 300
  -- Solve : 
          min 20.5251  max 300
  -- Solve : 
          min 20.9908  max 300
  -- Solve : 
          min 21.4374  max 300
  -- Solve : 
          min 21.8705  max 300
  -- Solve : 
          min 22.2963  max 300
  -- Solve : 
          min 22.7208  max 300
  -- Solve : 
          min 23.1495  max 300
  -- Solve : 
          min 23.5869  max 300
  -- Solve : 
          min 24.0365  max 300
  -- Solve : 
          min 24.5004  max 300
  -- Solve : 
          min 24.9798  max 300
  -- Solve : 
          min 24.9995  max 300
  -- Solve : 
          min 24.9996  max 300
  -- Solve : 
          min 24.9996  max 300
  -- Solve : 
          min 24.9996  max 300
  -- Solve : 
          min 24.9996  max 300
  -- Solve : 
          min 24.9996  max 300
  -- Solve : 
          min 24.9997  max 300
  -- Solve : 
          min 24.9997  max 300
  -- Solve : 
          min 24.9997  max 300
  -- Solve : 
          min 24.9997  max 300
  -- Solve : 
          min 24.9997  max 300
  -- Solve : 
          min 24.9998  max 300
  -- Solve : 
          min 24.9998  max 300
  -- Solve : 
          min 24.9998  max 300
  -- Solve : 
          min 24.9998  max 300
  -- Solve : 
          min 24.9999  max 300
  -- Solve : 
          min 24.9999  max 300
  -- Solve : 
          min 24.9999  max 300
  -- Solve : 
          min 24.9999  max 300
  -- Solve : 
          min 24.9999  max 300
  -- Solve : 
          min 24.9999  max 300
  -- Solve : 
          min 24.9999  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
  -- Solve : 
          min 25  max 300
Final time = 0.8

traffic:
  Normal end of execution.
times: compile 0.007901s, execution 161.733s,  mpirank:0
 CodeAlloc : nb ptr  4660,  size :516624 mpirank: 0
Ok: Normal End