Industrial Math - Alexiades
                  Lab 7
      Explicit Finite Volume scheme for Diffusion

Implement the explicit numerical scheme for 1-D diffusion of a concentration u(x,t)
in an interval [a, b], during time [0, tend], with boundary conditions specified below.
Your code should accept as inputs: MM, tend, dtout, factor, a, b, D ,
(then dx=1/MM, M=(b-a)*MM), and the time-step dt should be set as: dt = factor * dtEXPL
The code should output the profile: x(i) U(i) (i=0,...,M+1) initially (at time=0), at desired times tout (every dtout),
and also at time tend, in a file for plotting, i.e. in 2 columns, like (see below):
 # Profile at time: ???? , nsteps= ????
          x(0)             U(0)
          ....		  ....
          x(M+1)	   U(M+1)  
An explicitly solvable problem on which you can debug your code is the following diffusion problem in a semi-infinite interval:
    ut = D uxx , 0 < x < Inf , t > 0
    u(x,0) = 0, u(0,t) = 1

whose exact solution is:   u(x,t) = ERFC( 0.5 x / SQRT(D t) ),
with ERFC(·) = 1 − ERF(·) the Complementary ERROR function.
Verify (by hand) that this u(x,t) satisfies the problem.

To simulate this problem on a finite interval [0,b], we simply impose the exact solution at x=b, i.e. use the boundary conditions:
    U(0) = 1.0 ,   U(M+1) = ERFC( 0.5*b / SQRT(D*time) )
Thus, the exact solution and the numerical solution solve the same IBVP problem
and we can compare the two to calculate the maximum error at any time tn:
    ERRn = MAX{ |U(i) − u(xi, tn)| : i=0,...,M+1 }.
Something like the following routine can be called at tout (or at each time-step) to find the exact solution and the error:
      subroutine COMPARE( M, D, x, time, U, uEXACT,  ERR )
!.... Compare with exact solution: u(x,t) = ERFC(0.5*x/SQRT(D*t)) 
!.... of:    u_t = D u_xx , 0<x<Inf , t>0 ; u(x,0)=0 ; u(0,t)=1 
      dimension x(0:M+1), U(0:M+1), uEXACT(0:M+1)
                ERR = 0.0
          DO i = 0, M+1
                        arg = 0.5 * x(i) / SQRT( D*time )
             uEXACT(i) = ERFC( arg )
                     ERRi = ABS( U(i) − uEXACT(i) )
                     ERR  = MAX( ERRi, ERR )
1. Debug using: a=0., b=1., D=0.1, MM=8, tend=3., dtout=1., factor=0.98
  With only 8 nodes (MM=8), the agreement won't be very good.
  Nevertheless, try factor=0.90 to see how it improves.
2. Try MM=16 to see how it improves. The plotted curves should be indistinguishable,
  and the error at tend=3 should be less than 1%.
3. Try MM=32 and MM=64 to see further improvement in the error.
  Make runs with factor = 0.90 , 0.99 to see the effect.
Now you have a debugged diffusion code !
4. To see how easy it is to go unstable, make runs with factor=1.0, and then with factor=1.01 and 1.1.
  Look at the plot and the error.

Submit ONLY the ERR at tend from each of 1., 2., 3., 4. , and your code.
		   Remarks - Hints 

a. Only x, F, U, uEXACT  need to be indexed as (spatial) arrays. 
   The natural indexing for M control volumes is:
	   a=x(0),  x(1), ..., x(M),  x(M+1)=b
	[ i=1:M are the nodes, and i=0, i=M+1 are the boundary faces]
	  U(0),  U(1), ..., U(M),  U(M+1)
	[ i=1:M are the nodal values, and i=0, i=M+1 are the boundary values]
	  F(1),  F(2), ..., F(M),  F(M+1)		( F(i)=Fi-1/2 )
	[ i=2:M are the internal fluxes, and i=1, i=M+1 are the boundary fluxes]
   [For Matlab coding, indices must be shifted to start with 1, unfortunately]

b. Should use modular programming, with different tasks in separate subprograms, like: 
	MESH	construct the mesh: x(i), i=0,...,M+1
	INIT	set initial values for U(i), i=0,...,M+1
	FLUX	compute fluxes F(i), i=1,...,M+1  ( F(i)=Fi-1/2 )
	PDE	update U(i), i=0,...,M+1
	COMPARE	compute exact solution and error
	OUTPUT	print profile at particular time (for plotting)

   The structure of the main program will be something like this

	declare variables, dimension arrays 
	read in the data from dat file (see below)
	    call MESH( ... )
	!............ set the timestep ...........
	  dtEXPL = dx*dx/(2*D)		! max timestep for stability
	  dt     = factor*dtEXPL	! factor<1 increases stability
	  Nend   = INT( (tend-t0)/dt ) + 1	!#of timesteps 
	!............... initialize ...........
	    call INIT( ... )
	  nsteps = 0
	    time = 0.0
	    tout = MAX( dtout, dt )  
	  call OUTPUT( ... )	! to print initial profile
	!............... begin time-stepping ...............
        DO nsteps = 1, Nend
	      time = nsteps*dt	
	      call FLUX( ... )
	      call PDE ( ... )
	    IF( time >= tout ) THEN
	        call COMPARE( ... )  
		call OUTPUT( ... )
		tout = tout + dtout 
	    IF( time >= tend ) THEN
	       print*, ' DONE, at time= ', time, '  nsteps=', nsteps
	       print*, '       max error =', ERR
        ENDDO !.......... end of time-stepping .............
     	print*,'... out of timesteps: need bigger Nend'

c. The input data should be read in at run time so they can be changed easily for each run. 
    They are too many to type each time, so the best way is to read them from a dat file.
    e.g. here is a convenient arrangement (in a file dat):
	8    3.   1.                  : MM, tend, dtout
	0.98    0.1  0.0   1.0    : factor,  D,  a, b 
  [ In matlab, can use a script m-file  dat.m:
      M = 8 ;  tend=3. ;  dtout=1. ;  Nend=1000; factor=0.98 ;  D = 0.1 ;  a=0.0 ; b=1.0  
     which the code reads with:  dat  ]

   These are good values to start with. On unix, you execute with:
	a.out < dat        (output to screen), or
	a.out < dat > out  (output to file out)

d. The OUTPUT routine, called at each tout (every dtout)
   after COMPARE(...), should record the time and nsteps
   of the profile, so it should be something like:

      print"( '# Profile at time:', f9.4,'  nsteps=', I6 )",time,nsteps
      print"( '# Error up to this time:', e15.6 )", ERR
      DO i = 0, M+1
           print"( f12.4, g16.8, g16.8 )", x(i), U(i), uEXACT(i)
      print"()"       !! blank line for gnuplot to break plot

   Profiles should be written to a file for plotting (e.g. "") 
   and can use gnuplot to see the numerical and exact solutions:
      gnuplot> plot "" u 1:2 w lines lt 1 lw 2, "" u 1:3 w lines lt 3 lw 2

   To get the curves only initially and at tend=3.0, run with dtout=3.

e. For Fortran programmers: in double precision, you must use DERFC 
    for the (double precision) complementary error function.