Math578 - Alexiades

                  Lab 2:   Diffusion / Conduction
                  (always read the entire Lab before you start working on it...)

Consider a (1-dimensional) diffusion process, with (constant) diffusivity D > 0, for the concentration u(x,t)
in an interval a ≤ x ≤ b, during some time 0 ≤ t ≤ tend.
Initialy, the concentration profile is (some specified) uinit(x), a ≤ x ≤ b, and the ends are impermeable.
Part 1:
1. State precisely and fully the mathematical problem (PDE, IC, BCs) modeling this process.

2. Implement the explicit Finite Volume scheme for this problem, for a mesh of M uniform control volumes in the interval [a,b], up to a time tend.
  Your code should read values for MM, tend, factor, dtout, D, a, b from a data file,
  where MM = number of nodes per unit length. Then dx = 1.0/MM and M = (b-a)*MM = (b-a)/dx .
  For stability of the explicit scheme, the time-step dt is determined as: dt = factor * dtexpl, where dtexpl = dx2/(2D).

3. Take a=0.0, b=4.0, and consider the initial profile ("square bump"): uinit(x) = 5 for 1 ≤ x ≤ 2, uinit(x) = 0 otherwise .
  3a. Using MM=8 (so M=4*8=32), D=0.1, factor=0.95, tend=6, dtout=2,
    calculate and plot the profiles at times: 0., 2., 4., 6. on the same plot (mark which curve is at what time). What is happening to the "bump" ?
  3b. Calculate the area under the profile at time=0. and time=6. (can use Trapezoidal Rule Quadrature, see below).
    Do they agree ? should they ? What does this area represent physically ?
  (In the next Lab, you will test the code against an exact solution).

4. To speed things up... submit this part NOW, as follows: Create a single text file Lab2_1.txt containing:
  Name, Date, Lab2_part1
  =============================================== (separator line)
  Your answers to the questions in parts 1, 3a,3b
  =============================================== (separator line)
  Your code, with all routines.
  =============================================== (separator line)

  Upload to Canvas, under Lab2_1 by Jan.30. Then continue with the rest.
Part 2:
5. Now let's examine the effect of various parameters. For each of the runs listed bellow, do the following:
  • plot the initial and final (at time=tend) profiles on one plot;
  • on the plot, mark the parameter values that generated it;
  • comment as to what you think is happening and why, how it compares with the other cases, observations, remarks.
  Here are the cases to examine.
  5a. Effect of the time-step: dx = 1/16 (i.e. MM=16), D = 0.1, tend = 1., dtout=1 ,
     (5a1)  factor = 0.90
     (5a2)  factor = 0.99
     (5a3)  factor = 1.0
     (5a4)  factor = 1.05
  5b. Effect of the mesh: D=0.1, tend=1., dtout=1, factor=0.99,
     (5b1)  dx = 1/ 8 (i.e. MM= 8)
     (5b2)  dx = 1/16 (i.e. MM=16)
     (5b3)  dx = 1/32 (i.e. MM=32)
  5c. Effect of diffusivity: dx=1/8, tend=1., dtout=1, factor=0.90,
     (5c1)  D = 0.1
     (5c2)  D = 0.5

6. Create a single PDF file Lab2_2.pdf containing:
  Name, Date, Lab2_part2
  =============================================== (separator line)
  The plots and comments for 5a,5b,5c, each properly labeled ( (a1),(a2),... ).
  =============================================== (separator line)
  Your code
  =============================================== (separator line)

Upload "Lab2_2.pdf" on Canvas (under Lab2_2) by Feb.6
                CODING SUGGESTIONS
  In scientific computing, CLARITY and EFFICIENCY are paramount !
          so keep it simple, understandable, and efficient !
  • Use double precision for all real variables (not 'float').
  • Use natural indexing, short and meaningful variable names, as close to those used in the numerical scheme as possible.
      e.g. running index for space: i or j or k, running index for time: n
      PLEASE USE the same names for routines and variables as in the algorith presented in class
  • Natural indexing and array dimensions for a mesh of M control volumes are:
        x(0:M+1), with x(0)=a and x(M+1)=b the boundaries;
        U(0:M+1), with U(0) and U(M+1) the boundary values at a and b;
        F(1:M+1), with F(1)=F½ and F(M+1)=FM+½ the boundary fluxes at a and b.
  • Efficiency: Minimize operations (without sacrificing clarity), e.g.:
      * since M+1 is used repeatedly, set Mp1=M+1 and use Mp1 everafter.
      * always evaluate polynomials in nested form (Horner's Rule), e.g. ax²+bx+c = c + x*( b + a*x )
      * x*x*x is MUCH cheaper than x3 (even worse would be x3.0 which forces exp and log evaluations!)
      * Minimize storage: only x, F, U need to be indexed as (spatial) arrays.
  • Trapezoidal Rule Quadrature: Here is an efficient implementation ( for our mesh! ) in Fortran,
      which assumes the U array is indexed as U(0),U(1), ..., U(M+1). (check that it is correct!) 
    		double precision FUNCTION TrapzRule( M, Dx, U )
		dimension U(0:M+1)
	
		  Area = ( U(0) - U(1) - U(M) + U(M+1) )/4
		  DO i = 1,M
			Area = Area + U(i)
		  ENDDO
		  TrapzRule = Area * Dx
		end
  • Code outline for time-stepping scheme: 
          declare variables, dimension arrays
      read data file: MM, tend, t0, dtout, factor, D, a, b
      dx = 1.0/MM ;   M = (b-a)*MM     !! =(b-a)/dx
      dtEXPL = dx*dx/(2*D) ;  dt = factor*dtEXPL
      call MESH(...)
     !-----initialize-----!
      nsteps = 0 ; time = 0.0 ; tout = dtout ; MaxSteps=INT(tend/dt)+1
      call INIT(...)
      call OUTPUT(...)	!!profile at time=0
     !-----begin timestepping-----!
      DO nsteps = 1, MaxSteps
	 call FLUX(...)
	 call  PDE(...)
	 time = nsteps*dt	!! =time + dt
	 if( time >= tout ) then
		call OUTPUT(...)
		tout = tout + dtout
	 endif
      ENDDO
     !-----end of timestepping-----!
      print: 'DONE, exiting at t=', time ,' after ', nsteps ,' steps'
    PLEASE USE the above names for routines and variables, as in the algorith presented in class.

  • Output for plotting: 
       *  Lines starting with # are comments for gnuplot (and most other plotters).
   *  Blank line starts a new curve.
   *  It is best to print run-time info to stdout (to screen or redirect to a file: a.out < dat > o.out ),
       and to print values for plotting in separate file(s)
   *  I find it convenient to name all output files as "o.name", so profiles in  o.prof, histories in o.hist, etc.
   *  Profiles are pairs:  x(i), U(i), i=0,...,M+1, at specific desired times, printed in 2 columns.
   *  You can print all profiles in one file, say o.prof, separated by blank line(s), 
       each preceded by a comment line specifying the time of that profile. But when plotted, all will be of the same color. 
   * Or each in a separate file, e.g. o.prof0, o.prof2, o.prof4, o.prof6, which you can gnuplot with:
	   plot "o.prof0"  u 1:2 t 'time=0'  w lines lt -1 lw 2, \
		"o.prof2"  u 1:2 t 'time=2'  w lines lt 2 lw 3, \
		"o.prof4"  u 1:2 t 'time=4'  w lines lt 3 lw 3, \
		"o.prof6"  u 1:2 t 'time=6'  w lines lt 1 lw 3
       each of different color, and with a tag, much nicer!
  Caution: Do not go overboard trying to code up output at all these times!  
          just make a run for each with appropriate parameters (tend, dtout) and rename the file!
   * See Gnuplot links on main page for how to print, etc.