Exercise 1

Last update: 12 June 2003

Problem description
-------------------

The assignment is concerned with the solution of the advection-diffusion equation:


dC/dt + U dC/dx = K d2C/dx2 + S

For this assignment, the source term S is 0.

The problem to be solved is the concentration C(x,t) for the following
region:

x = 0 to   x = L, with L a length given below.
t = 0 till t = T, with T a time to be given below.


To fully specify the problem, we need to give an initial condition
C(x,0) at time t=0, and boundary conditions.

The initial condition is C(x,0) = 0.

The boundary conditions are C(0,t) = Cleft, with Cleft a concentration value given below.
                            C(L,t) = Cright = 0.


The velocity U and the diffusion coefficient are given, see below.

Part of the program (in FORTRAN) is given to you. Before starting any
programming, read the Fortran language description given on the web
pages for numerical fluid dynamics, and make the FORTRAN exercise.

In the FORTRAN program, you need to:
-- set the numerical grid (routine grid)
-- implement the exact solution (routine output)
-- implement two difference schemes for the numerical solution (routine
   explicit), see below
-- calculate the error of the numerical solution (routine output)
-- use .d0 after numerical constants to retain double precision accuracy


Remember to re-compile if you change anything in the program.
The program reads the number of time steps, nstep, time step dt, and
icase from a file called param so that you can quickly change these aprameters
without re-compiling. The number icase is used to choose between
two types of discretisation (see below).
The program has output in the form of a file with three columns, namely
x, the numerical solution and the exact solution

Discretisation
--------------

The problem is to be solved numerically on an equidistant grid with N
grid cells:

   |   |   |   |   |   |   |
   i=0                     i=N

For the advection part, use two methods to discretise the problem. The methods are
upwind and central differences in space. For the
diffusion part, use central differences in space. Use forward time
differencing (NOT central time differencing!!).

Write out the discretisation in full detail before programming it. You
will have to hand in  the discretisation, written out in full detail,
anyhow.

NOTE: the schemes mentioned have a time-step limitation. Find out from
literature or by numerical experiment what the timestep limitation is.
Use a timestep which is still stable.

Check on the implementation
---------------------------
You test correct programming by solving the problem for U=0. and
a given K, which is not 0, so that you have only diffusion.
In this case, whatever the value of K, the solution is a linear function
which is Cleft at x=0 and 0 at x=L Moreover, for the case that U is not 0,
the numerical solution must closely resemble the exact solution on very
fine grids.

Values for the constants
------------------------

Set the value for L equal to 1. followed by the first number of your student number.
If the first number of your student number happens to be 0 take the next non-zero number.

Set the value for K equal to 0 followed by the second number of your student number.
If the second number of your student number happens to be 0 take the next non-zero number.

Set the value for Cleft equal to the third number of your student number.
If third number of your student number happens to be 0 take the next non-zero number.

Set the value for U equal to 1. followed by the fourth number of your student number.
If the fourth number of your student number happens to be 0 take the next non-zero number.

Example.
-------

Student number: 9035678

L: first number is 9              -->       L = 1.9
K: second numer = 0 --> next: 3   -->       K = 0.3
Cleft: third number = 3           -->   Cleft = 3.
U: fourth number = 5              -->       U = 1.5

Cases to be calculated, results expected.
-----------------------------------------

Use dt = small enough: Verify, that the time step you use is small enough for central
differences for the advection part. Check this in subroutine checkst.

Determine a number of grid points such that the the cell Peclet number is just smaller
than 1.5  If this number is (say) Ntest
Determine the solution for grids with N = Ntest, 2*Ntest, 4*Ntest

Write down the resulting L, K, Cleft, U, Ntest, Pe number and cell Pe number

Plot the solution at times t=T1 and T2 (total of 6 plots)
Here, T1 is 0.25*L/U and T2 is a time, so long that the solution has become stationary.
You have to determine yourself whether the solution has become stationary.
How can you check that? Check it. Plot all plots with the same scale, so that you can really
compare. Number the plots so you can refer to them. Also hand in
the final solution in the form of a table for N=10

Especially determine the error between the exact, stationary solution
and the solution at time=10. Determine the order in space for the two cases
(central diferences for the advection part and upwind for the advection part) numerically.
You may do it graphically, using a double log graph of error vs number of points or
error vs grid spacing.  Is the result for the order as you expect it to be?
You may need to go to g\more grid points to see what you expect (NO plots of higher
resolutions necessary!).
Why would you use the finer resolutions, not the coarser resolutions
to estimate the order of the method?

I expect you to be able to discuss the results with me (differences between
results and exact solution, differences between the numerical solutions
with the two schemes). Include the cell Pe number in your discussion. Sketch
(qualitatively) how the solution would look like for another cell Pe number, say Pe2
which you calculate as follows:
Pe2 = 2-Pe (the number 2 minus the Pe number)
HOWEVER: if Pe2 becomes < 0 then take Pe2 = 1

Hints and Directions
--------------------
-- mind the stability criterion for the scheme you use

-- do you know the stability criterion for the scheme with  upwind for advection combined with central for diffusion?

    (you do not need to derive it)

   it is not enough to have stability for advection and diffusion separately

    use: number = U/dx + 2*Kdif/dx**2

         dt  < number

    
-- be aware: if you add text in the program , the length
   of the line may go over 72 characters: this is against Fortran rules.
   Remove spaces to make the line shorter, or split the line in two using
   a & in column position 6
   especially be aware of this if you implement long formulas (diffference schemes
   and exact solution)
--  do not simply present graphs/results: discuss them (for instance, compare with results for other resolutions,

    other methods, exact solution, observe special features (spread, height of peak, wiggles etc))

--  if, for some reason, you find things which appear wrong, discuss them with the mentor

     instead of   having to make conclusions for results with errors in them

c
c the following in Fortran:
c
      x = sin(x)

c
c is equivalent to the following, splitting the line in 2 parts
c with an ampersand, &, in character position 6
c
      x =
     &  sin(x)

-- NOTE: x=0 corresponds with i=0, x=L corresponds with i=N. In the fortran
   program, be careful with choosing the index i at which you set the
   boundary condition

   The point i=N+1 is NOT necessary in this assignment

-- The exact solution to the stationary problem is:

         C = Cleft*(1.-exp(-U*(L-x)/K) )/(1.-exp(-L*U/K))

   Validate your results by comparing to this exact solution (you need it anyhow
   for calculating the error)