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

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

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

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

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

# 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 ;)