Demo: A 1-D linear transport equation

This code illustrates the numerical solution of a 1-D linear transport equation using a finite difference scheme.

  • For conciseness, we skip the description of the heading of the main program. This has been presented previously. We assign values of the length of the interval and the velocity. We also declare a file where the solution is stored

       const double L = 10., c = 4.;
   ofstream pf("output.dat");

  • We now declare discretization data. nx is the number of grid sub-intervals and h is the grid size. The time step is obtained as a program argument and the final time is set to 1. These two parameters are assigned to global variables in .

       size_t nx = atoi(argv[1]);
   double h = L/double(nx);
   theTimeStep = atof(argv[2]);
   theFinalTime = 1.;

  • We compute the CFL of the numerical scheme and then declare vectors x, u and v which are respectively the vector of grid point coordinates, solution at previous time step and solution at actual one. We then define grid as a uniform one. 

       double cfl = c*theTimeStep/h;
   Vect<double> x, u(nx+1), v(nx+1);
   x.setUniform(0.,L,nx+1);

  • We define the initial condition as a Gaussian function:

       u.set("exp(-20*(x-2)^2)",x);

  • We set the time loop and implement the upwind finite difference scheme to compute the solution at current time step:

       TimeLoop {
      v(1) = u(1) - cfl*(0.5*u(2)-cfl*(u(2)-2*u(1)));
      for (size_t i=2; i <=nx+1; i++)
         v(i) = u(i) - cfl*(0.5*(u(i+1)-u(i-1))-cfl*(u(i+1)-2*u(i)+u(i-1)));
      u = v;

  • We store solution at current time step into the dedicated file and close the time loop.

          for (size_t i=1; i<=nx+1; i++)
         pf << x(i) << "   " << u(i) << endl;
      pf << endl;
   }

  • We finally print some information: 

       cout << "Number of grid intervals: " << nx << endl;
   cout << "Time step:                " << theTimeStep << endl;
   cout << "CFL:                      " << cfl << endl;
   return 0;