I started watching the CFD lectures provided by Boston University via iTunes, and thought I could share some bad code with the world. The lectures apparently consist of 12 steps to solving the Navier-Stokes equations, where the first one is to solve a 1D convection equation. The class used Python, but I will be using C++.

Our aim is to solve

\displaystyle \frac{\partial u(x,t)}{\partial t} + c\frac{\partial u(x,t)}{\partial x}=0, \ \ \ \ \ (1)

subject to some initial condition {u(x,0)=u_0}.The analytical solution to this problem is {u(x,t)=u_0(x-ct).}
We will approximate the left-hand side by a forward difference, and the right hand side by a backward difference. Read the post on finite-differences if the topic is unfamiliar to you.

The finite-differences we will use, together with their truncation errors are:

\displaystyle \frac{\partial u(x,t)}{\partial t} = \frac{1}{k}\left(u(x,t+k)-u(x,t)\right) + O(k)

and

\displaystyle \frac{\partial u(x,t)}{\partial x} = \frac{1}{h}\left(u(x,t)-u(x-h,t)\right) + O(h).

To simplify the final expression, we introduce the notation {u(x,t)=v_{i,j}}. Using a different letter for the computed solution {(v)} indicates that the computed solution is different from the real one, {(u)}. We are now able to rewrite equation (1) in its finite-difference form:

\displaystyle \frac{1}{k}\left(v_{i,j+1}-v_{i,j}\right) = \frac{c}{h}\left( v_{i,j}-v_{i-1,j}\right).

The initial condition makes sure that {v_{i,0}} is known for all {i}. We can therefore solve for {v_{i,j+1}} to get,

\displaystyle v_{i,j+1} = v_{i,j} - \frac{ck}{h}\left(v_{i,j}-v_{i-1,j}\right).

I won’t go into details about the implementation, because it is rather simple. The initial condition is chosen to be a step-function. I think the reason Prof. Barba has chosen this function is that it is not continuous, and solving using finite-differences will give pretty bad results. The output is saved to a .dat file. I import the .dat file into MATLAB to generate plots. As expected, the results are not great. It’s quite instructive to play around with the step sizes, number of meshpoints and the initial condition. I apologize for not naming my variables as they are named in the above derivation.

/*
 * Solves 1D convection equation, also known as the one-way wave equation.
 */

#include <iostream>
#include <fstream>
#include <vector>
using namespace std;

void saveArray(vector< vector<double> > &u, int nx, int nt, double dx, double dt, int c);

int main()
{

  int nx = 20;
  int nt = 50;
  double dx = double(2)/(nx-1);
  double dt = 0.01;
  int c = 5;

  vector< vector<double> > u(nx,vector<double> (nt));
  vector< vector<double> > un(nx,vector<double> (nt));

  // Initial condition:
  for (int i=0; i <= nx-1; i++)
    {
      if (i*dx >= 0.5 && i*dx <= 1)
	{
	  u[i][0] = 2;
	}
      else
	{
	  u[i][0] = 1;
	}
    }

  // Finite-difference loop:
  for (int it=1; it<=nt-1; it++)
    {
      for (int k=0; k<=nx-1; k++)
	{
	  un[k][it-1] = u[k][it-1];
	}
      for (int i=1; i<=nx-1; i++)
	{
	  u[0][it] = un[1][it-1];
	  u[i][it] = un[i][it-1] - c*dt/dx*(un[i][it-1]-un[i-1][it-1]);
	}
    }

  saveArray(u, nx, nt, dx, dt, c);

  return 0;
}

// Save array to file:
void saveArray(vector< vector<double> > &u, int nx, int nt, double dx, double dt, int c)
{
  ofstream myfile;
  myfile.open("1d_convection02.dat");
  myfile << "%\t" << "nx = "  << nx  << "\tnt = " << nt <<
    "\tdt = " << dt << "\tc = " << c << endl;

  for (int i=0; i<=nx-1; i++)
    {
      myfile << i*dx << "\t\t";
      for (int j=0; j<=nt-1; j++)
	{
	  myfile << u[i][j] << "\t\t";
	}
      myfile << endl;
    }
}

About these ads