Arjen Markus (14 october 2008)

The *method of manufactured solutions* (cf. [1 ]) is a method for testing the solution algorithms for (partial) differential equations. (The method is also known as *forced solutions*, it seems) The idea is very simple:

Given a PDE like:

dC d2C dC -- = D --- - u -- + S dt dx2 dx

the advection-diffusion equation with constant coefficients and a source term S.

Select a function f(t,x) that should fulfill this equation, but determine the source term like this:

df d2f df S = -- - D --- + u -- dt dx2 dx

So: simply fill in the solution and see what source term and boundary conditions we need to fulfill the equation.

We can then use the program that needs to be tested to see how accurate it can approximate the solution.

The main indicator for correctness is whether the error depends on the schematisation (delta-t and delta-x in this case) in the same way as the theory claims: error ~ dt**2 for second-order methods for instance.

Alas, in the above case, we need to deal with both delta-t, the time step, and delta-x, the finite size of the grid cells used to discretise the above equation. Furthermore, the relative magnitude of the advective and diffusive term may matter - the truncation errors involved in the discretisation may have different orders!

To see how we can deal with this problem, let us just experiment...

In the program below the manufactured solution is:

- kt f(t,x) = (1 - e ) x(L-x)/L**2

and boundary conditions:

x = 0, L: C = 0 t = 0: C = 0

The source term becomes:

-kt -kt S = (k e x(L-x) + (1 - e )(u(L-2x) - 2D))/L**2

Parameters:

- k, L, D, u
- dx = L/nC, dt = 10/knT (nC= number of cells, nT= number of time steps)

# experiment_mms.tcl -- # Experiment with Manufactured Solutions # # exactSolution -- # Determine the exact solution at time t # # Arguments: # t Time in question # # Result: # Values of the exact solution at the cell midpoints # proc exactSolution {t} { global k global nC global L set result {} set tfactor [expr {1.0 - exp(-$k * $t)}] for { set i 0 } { $i < $nC } { incr i } { set x [expr {($i + 0.5) * $L / $nC}] set value [expr {$tfactor * $x * ($L-$x)/$L/$L}] lappend result $value } return $result } # accumulateError -- # Determine the accumulative error # # Arguments: # error_ Error from previous time (name of variable) # solution Numerical solution # exact Exact solution # # Result: # None # Side effect: # error increased # proc accumulateError {error_ solution exact} { upvar 1 $error_ error foreach s $solution e $exact { set error [expr {$error + abs($s-$e)}] } } # sourceTerm -- # Compute the value of the source term at the given time # # Arguments: # t New time # Result: # Source term per cell # proc sourceTerm {t} { global k global nC global L global u global D set result {} set tfactor [expr {exp(-$k * $t)}] for { set i 0 } { $i < $nC } { incr i } { set x [expr {($i + 0.5) * $L / $nC}] set value [expr {$k*$tfactor * $x * ($L-$x) + (1.0-$tfactor)*($u*($L-2.0*$x) - 2.0*$D)}] set value [expr {$value/$L/$L}] lappend result $value } return $result } # nextTime -- # Compute the numerical solution at the next time # # Arguments: # solution_ Solution so far (name of variable) # t New time # dt Time step # Result: # None # Side effect: # solution updated # Note: # Upwind differences used, u >= 0 # proc nextTime {solution_ t dt} { global L global nC global D global u upvar 1 $solution_ solution set source [sourceTerm $t] set dx [expr {$L/($nC-1.0)}] set extended_solution [concat 0.0 $solution 0.0] set solution {} for {set i 1} {$i <= $nC} {incr i} { set cc [lindex $extended_solution $i] set cl [lindex $extended_solution [expr {$i-1}]] set cr [lindex $extended_solution [expr {$i+1}]] set src [lindex $source [expr {$i-1}]] set dcdt [expr {$D*($cl+$cr-2.0*$cc)/$dx/$dx - $u*($cc-$cl)/$dx + $src}] lappend solution [expr {$cc + $dcdt * $dt}] } } # initialCondition -- # Make the initial condition # # Arguments: # None # Result: # List representing the initial condition # proc initialCondition {} { global nC set condition {} for {set i 0} {$i < $nC} {incr i} { lappend condition 0.0 } return $condition } # determineError -- # Determine the error for a given set of constants # # Arguments: # length Length of the interval (m) # velocity Flow velocity (m/s) # diffusion Diffusion/dispersion coefficient (m2/s) # decay Decay coefficient for transient part # numberCells Number of cells (>1) # numberTimes Number of times (>1) # Result: # Absolute error, accumulated over the times # proc determineError {length velocity diffusion decay numberCells numberTimes} { global L global u global D global k global nC global nT set L $length set u $velocity set D $diffusion set k $decay set nC $numberCells set nT $numberTimes set error 0.0 set solution [initialCondition] set dt [expr {10.0/($k*$nT)}] # # Note: the 0.5 in the computation of the time makes a large # difference in the results! # for {set t 0} {$t < $nT} {incr t} { set time [expr {$dt * ($t+0.5)}] # set time [expr {$dt * ($t)}] nextTime solution $time $dt accumulateError error $solution [exactSolution [expr {$dt*($t+1)}]] } return $error } # main -- # Perform the experiments # # Decrease the time step: linear decrease in the error expected, as the # method is first-order in time # # # Checking the code ... # set L 10000.0 set u 1.0 set D 0.01 set k 0.1 set nC 100 set nT 100 set solution [initialCondition] set time [expr {10.0/($k*$nT)}] set exact [exactSolution $time] nextTime solution [expr {$time/2.0}] $time foreach s $solution e $exact { puts "[expr {$s-$e}] -- $s -- $e" } # # More systematic tests # puts "Varying time step:" foreach nT {100 300 1000 3000} { set error [determineError 10000.0 1.0 0.01 0.1 100 $nT] puts "$nT\t$error\t[expr {$error/$nT}]" } puts "Varying cell size:" foreach nC {100 300 1000 3000} { set error [determineError 10000.0 1.0 0.01 0.1 $nC 100] puts "$nC\t$error\t[expr {$error/$nC}]" } # # Residence time: 10000.0/1.0 = 1e4 seconds - based on velocity # 10000.0**2/0.01 = 1e10 seconds - based on dispersion # Period: 10.0/0.1 = 100 seconds # # Maximum time step: 1.0*dt < 10000.0, so with 100 time steps, # the period should not be longer than, say, 1.0e6 seconds or # a decay rate of less than 10.0*100/10000.0 = 0.1 # # Otherwise: instability! # # The consequence of a lazy design ;)

AM (22 october 2008) Here is an experiment with PDEs and the Tensor package.