Math578 - Alexiades
Diffusion / Conduction
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 u0(x), a ≤ x ≤ b, and the ends are impermeable.
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/MM and M = (b-a)/dx = (b-a)*MM.
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"):
u0(x) = 5 for 1 ≤ x ≤ 2, u0(x) = 0 otherwise.
a. 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" ?
b. 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 ?
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 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.
c. Effect of the time-step: dx = 1/16 (i.e. MM=16), D = 0.1, tend = 1., dtout=1 ,
(c1) factor = 0.90
(c2) factor = 0.99
(c3) factor = 1.0
(c4) factor = 1.05
d. Effect of the mesh: D=0.1, tend=1., dtout=1, factor=0.99,
(d1) dx = 1/ 8 (i.e. MM= 8)
(d2) dx = 1/16 (i.e. MM=16)
(d3) dx = 1/32 (i.e. MM=32)
e. Effect of diffusivity: dx=1/8, tend=1., dtout=1, factor=0.90,
(e1) D = 0.1
(e2) D = 0.5
What to submit:
1. By M Sep 4 (midnight): Put your code into a txt file "lab2.txt", and submit as email ATTACHMENT to
email@example.com , (and to yourself) with Subject: M578 lab2 code.
2. On R Sep 7: Submit on PAPER: your answers to the questions in parts 1, 3a,b, the plots and comments for 3c,d,e, and your code.
(In the next Lab, you will test the code against an exact solution).
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
Natural indexing and array dimensions for a mesh of M control volumes:
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 )
Area = ( U(0) - U(1) - U(M) + U(M+1) )/4
DO i = 1,M
Area = Area + U(i)
TrapzRule = Area * Dx
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
nsteps = 0 ; time = 0.0 ; tout = dtout ; MaxSteps=INT(tend/dt)+1
call OUTPUT(...) !!profile at time=0
DO nsteps = 1, MaxSteps
time = nsteps*dt !! =time + dt
if( time >= tout ) then
tout = tout + dtout
!-----end of timestepping-----!
print: 'DONE, exiting at t=', time ,' after ', nsteps ,' steps'
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.