Lab 7

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:

(then dx=1/MM, M=(b-a)*MM), and the time-step dt should be set as:

The code should output the profile: x(i) U(i) (i=0,...,M+1) initially (at time=0), at desired times

and also at time

# 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:

u(x,0) = 0, u(0,t) = 1

whose exact solution is:

with

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:

Thus, the exact solution and the numerical solution solve the

and we can compare the two to calculate the maximum error at any time t

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:

With only 8 nodes (MM=8), the agreement won't be very good.

Nevertheless, try

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

4. To see how easy it is to go unstable, make runs with

Look at the plot and the error.

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)=Remarks - HintsF) [ 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)=_{i-1/2}F) 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_{i-1/2}declare variables, dimension arraysread in the data from(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 IF( time >= tend ) THEN print*, ' DONE, at time= ', time, ' nsteps=', nsteps print*, ' max error =', ERR STOP ENDIF ENDDO !.......... end of time-stepping ............. print*,'... out of timesteps: need bigger Nend' END 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 adatfiledatfile. e.g. here is a convenient arrangement (in a filedat):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-filedat.m:M = 8 ; tend=3. ; dtout=1. ; Nend=1000; factor=0.98 ; D = 0.1 ; a=0.0 ; b=1.0which the code reads with: dat ] These are good values to start with. On unix, you execute with:a.out < dat(output to screen), ora.out < dat > out(output to fileout) d. The OUTPUT routine, called at eachtout(everydtout) after COMPARE(...), should record thetimeandnstepsof 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) 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 2To get the curves only initially and at tend=3.0, run with dtout=3. e.in double precision, you must use DERFC for the (double precision) complementary error function.For Fortran programmers: