Method of Manufactured Solutions

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.