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 )
ENDDO
END
-----------------------------------------------------------------
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 for x and U arrays 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 ), including boundary fluxes from BCs
PDE update U(i), i=0,...,M+1 including the boundary values, consistently with the BCs
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
ENDIF
ENDDO !.......... end of time-stepping .............
IF( time >= tend ) THEN
print*, ' DONE, at time= ', time, ' nsteps=', nsteps
print*, ' max error =', ERR
ELSE
print*,'... out of timesteps: need bigger Nend'
ENDIF
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. ; 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/Linux, 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 (as comment, starting with # ),
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)
ENDDO
print"()" !! blank line for gnuplot to break plot
Profiles should be written to a file for plotting (e.g. "o.prof")
and can use gnuplot to see the numerical and exact solutions:
gnuplot> plot "o.prof" u 1:2 w lines lt 1 lw 2, "o.prof" 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.
In other languages use erfc