Consider a (1-dimensional) advection-diffusion process, with (constant) diffusivity D > 0, 
for some concentration u(x,t) whose initial profile is u(x,0) = uinit(x), driven by a given velocity V, 
in an (appropriately chosen) interval a ≤ x ≤ b whose ends are impermeable, during some time  0 ≤ t ≤ tmax.

1. State precisely and fully the mathematical problem modeling this process (PDE, where, IC, BCs).

2. Briefly describe two different physical problems that could be modeled by this mathematical problem.

Take as initial profile the "square bump" :   uinit(x) = 5.0  for 1 ≤ x ≤ 2,   uinit(x) = 0  otherwise.
Assume constant velocity V > 0 and diffusivity D > 0;  take tmax=4.0, dtout=2, 
and Dx = 1/MM with MM = 16 (see Note below). 

Starting at time=0, and every dtout, your code should output the time, number of time-steps, 
and the entire profile of u (including at time 0 and tmax) for plotting.  
Also, as a check of how well it is conserving, your code should compute and output the total area 
under the concentration profile at each output time (you can write your own Trapezoidal Rule routine, 
see below, or may use Matlab's if you are coding in Matlab).

First, consider the case of pure advection ( D = 0 , V > 0 ).

3. State precisely and fully the mathematical problem modeling pure advection.
   Describe what you expect to happen qualitatively and sketch (by hand) the initial profile uinit(x) 
   and the expected profile at time t=4 if V=1 and if V=5.  How far will the pulse travel in each case ?

Next, consider the case of pure diffusion ( V = 0 , D > 0 ).

4. State precisely and fully the mathematical problem modeling pure diffusion.
   Describe what you expect to happen qualitatively for low and high diffusivity and sketch some profiles by hand.

Now we want to examine how the presence of diffusion affects the advection profiles and vice versa.

For each of the cases listed bellow, do the following:
 Plot the profiles on one plot. Set the y-axis range to [0:5.2] to see the top clearly 
    (in gnuplot: set yrange [0 : 5.2] ). 
 On the plot, mark the parameter values that generated it, mark which curve is at what time  
    [by hand, if more convenient].
 Below the plot, write the computed area, and comment on what you observe (including the "also try"
    and other explorations), any surprises/unexpected behavior, what you think is happening and why, 
    how this case compares with other cases.

Here are the cases to examine.

5. Pure advection:
          (5a)  V=1., D=0., a=0., b= 7., factor=1.0
          (5b)  V=5., D=0., a=0., b=25., factor=1.0
          (5c)  V=5., D=0., a=0., b=25., factor=1.1
          (5d)  V=5., D=0., a=0., b=25., factor=0.97 (also try .99)

6. Pure diffusion:
          (6a)  V=0., D=0.1, a=0., b= 4., factor=0.99 (also try 1.0)
          (6b)  V=0., D=0.5, a=0., b= 4., factor=0.99 (also try 1.0)
          (6c) same as (6a) but with  Dx = 1/32

7. Advection-Diffusion: Low velocity:
          (7a)  V=1., D=0.1, a= 0., b=10., factor=1.0
          (7b)  V=1., D=0.1, a= 0., b=10., factor=0.98
          (7c)  V=1., D=0.5, a=-2., b=13., factor=1.0

8. Advection-Diffusion: High velocity:
          (8a)  V=5., D=0.1, a=0., b=30., factor=1.0
          (8b)  V=5., D=0.5, a=0., b=30., factor=1.0
 Note: The best way to set Dx is to read in an integer: MM = number of control volumes per unit length,
and then set Dx and M = total number of control volumes: Dx = 1.0/MM , M = (b-a)/Dx = (b-a)*MM.
 Trapezoidal Rule Quadrature: 
Here is an implementation (for our mesh!), in fortran,
which assumes the U array is indexed as U(0),U(1), ..., U(M+1).

	double precision FUNCTION TrapzRule( M, Dx, U )
		dimension U(0:M+1)
		sum = ( U(0)+3*(U(1)+U(M))+U(M+1) )/4
	 DO i= 2,M-1
		sum = sum + U(i)
	 ENDDO
	 TrapzRule = Dx * sum
	end