Arjen Markus (3 march 2006) I have been wanting to bring in some physics into the Wiki for some time now. The problem is that solving physical puzzles always carries a lot of context with it. The script below shows how to solve one particular puzzle that is not too difficult.

The actual problem is artistically based on things I deal with for a living: it considers the spreading of some pollutant in a river system. The river is a slowly flowing body of water, 100 km long. It has a main channel and a shallow bank along it (which is almost stagnant water). The pollutant is spilled at time zero (so I can use the initial condition to simulate it). Then the pollutant is followed over ten days and the concentration at three points along the river is monitored.

There are plenty of things to do before this becomes something even remotely useful:

• Add a graphical facility for displaying the results
• Try and see if a different data structure (lists instead of arrays) can speed up things

Still, it was fun to implement and test :)

``` # wqmodel.tcl --
#     Initial set-up for a toy water quality model:
#     - Define a schematisation with finite volumes
#     - Define the substances to be modelled
#     - Solve the advection-diffusion equations
#
#     Note:
#     Using arrays may not be very efficient, it
#     does keep the code simple for the moment
#

namespace eval ::WQModel {
variable segments   {}
variable substances {}
variable volume
variable surface
variable exchanges
variable conc
variable mass
variable deriv

namespace export segment substances initialise timeframe exchange \
initial-values boundary-values getconc nextstep
}

# segment --
#     Define a segment by its ID and volume and surface area
#
# Arguments:
#     id          (Numerical) ID of the segment
#     vol         Volume of the segment
#     surf        Surface area (to compute the depth for instance)
#
# Result:
#     None
#
proc ::WQModel::segment {id vol surf} {
variable segments
variable volume
variable surface

lappend segments \$id
set volume(\$id)  \$vol
set surface(\$id) \$surf
}

# substances --
#     Define the names (and order) of the substances
#
# Arguments:
#     args        List of substance names
#
# Result:
#     None
#
proc ::WQModel::substances {args} {
variable substances

set substances \$args
}

# initial-values --
#     Define the initial values of the substances per segment
#
# Arguments:
#     id          Segment ID
#     values      List of initial values
#
# Result:
#     None
#
proc ::WQModel::initial-values {id values} {
variable conc

set conc(\$id) \$values
}

# boundary-values --
#     Define the boundary values of the substances per
#     boundary segment
#
# Arguments:
#     id          Segment ID (starting with "B")
#     values      List of boundary values
#
# Result:
#     None
#
proc ::WQModel::boundary-values {id values} {
variable segments
variable conc
variable deriv

set conc(\$id)  \$values
set deriv(\$id) \$values  ;# They must simply exist
}

# exchange --
#     Define the exchanges between segments
#
# Arguments:
#     from        ID of segment "from"
#     to          ID of segment "to"
#     flow        Flow rate
#     disp        Dispersion coefficient
#     area        Area of exchange surface
#     distance    Representative distance between centres
#
# Result:
#     None
#
# Note:
#     Boundary segments should have an ID starting
#     with "B"
#
proc ::WQModel::exchange {from to flow disp area distance} {
variable exchanges

lappend exchanges \$from \$to \$flow \$disp \$area \$distance
}

# timeframe --
#     Define the start and stop time and time step
#
# Arguments:
#     start       Start time
#     stop        Stop time
#     timestep    Time step
#
# Result:
#     None
#
# Note:
#     Time in days (flow rate in m3/s)
#
proc ::WQModel::timeframe {start stop timestep} {
variable start_time
variable stop_time
variable time_step

set start_time \$start
set stop_time \$stop
set time_step \$timestep
}

# initialise_derivs --
#     Initialise the derivatives
#
# Arguments:
#     None
#
# Result:
#     None
#
proc ::WQModel::initialise_derivs {} {
variable segments
variable substances
variable deriv

set dv {}
foreach s \$substances {
lappend dv 0.0
}

foreach id \$segments {
set deriv(\$id) \$dv
}
}

# initialise --
#     Initialise the computation
#
# Arguments:
#     None
#
# Result:
#     None
#
proc ::WQModel::initialise {} {
variable segments
variable substances
variable conc
variable volume
variable mass
variable time
variable start_time

set time \$start_time

foreach id \$segments {
set mass(\$id) {}
set vol       \$volume(\$id)
foreach s \$substances c \$conc(\$id) {
lappend mass(\$id) [expr {\$c*\$vol}]
}
}
}

# new_concentrations --
#     Compute the new concentrations
#
# Arguments:
#     None
#
# Result:
#     None
#
proc ::WQModel::new_concentrations {} {
variable segments
variable substances
variable conc
variable mass
variable volume
variable deriv
variable time_step

set dv {}
foreach s \$substances {
lappend dv 0.0
}

foreach id \$segments {
set newconc {}
set newmass {}
foreach m \$mass(\$id) dv \$deriv(\$id) v \$volume(\$id) {
set m [expr {\$m + \$dv * \$time_step}]
lappend newmass \$m
lappend newconc [expr {\$m / \$v}]
}
set conc(\$id) \$newconc
set mass(\$id) \$newmass
}
}

# transport_processes --
#     Compute the derivatives due to transport
#
# Arguments:
#     None
#
# Result:
#     None
#
proc ::WQModel::transport_processes {} {
variable segments
variable conc
variable deriv
variable exchanges

foreach {from to flow disp area distance} \$exchanges {
#
# Straightforward backward differences ... not accurate, but simple!
#
if { \$flow > 0.0 } {
set qfrom [expr {86400.0*\$flow}]
set qto   0.0
} else {
set qfrom 0.0
set qto   [expr {-86400.0*\$flow}]
}
set disprate [expr {86400.0*\$disp*\$area/\$distance}]

set newfrom  {}
set newto    {}
foreach cfrom \$conc(\$from) dvfrom \$deriv(\$from) \
cto   \$conc(\$to)   dvto   \$deriv(\$to)   {
lappend newfrom [expr {\$dvfrom - \$cfrom * \$qfrom + \$cto   *\$qto   - \$disprate * (\$cfrom-\$cto)}]
lappend newto   [expr {\$dvto   + \$cfrom * \$qfrom - \$cto   *\$qto   - \$disprate * (\$cto-\$cfrom)}]
}
set deriv(\$from) \$newfrom
set deriv(\$to)   \$newto
}
}

# nextstep --
#     Compute the concentrations for the next step
#
# Arguments:
#     None
#
# Result:
#     1 if there is a next step, 0 if the computation has ended
#
proc ::WQModel::nextstep {} {
variable time
variable time_step
variable stop_time

#
# Prepare the step
#
initialise_derivs

#
# Deal with water quality processes
#
# TODO

#
#
# TODO

#
# Transport processes
#
transport_processes

#
# Compute the new concentrations
#
new_concentrations

set time [expr {\$time + \$time_step}]
set ::time \$time
return [expr {\$time <= \$stop_time}]
}

# getconc --
#     Return the concentrations for a particular segment
#
# Arguments:
#     id         The ID of the segment
#
# Result:
#     List of concentrations
#
proc ::WQModel::getconc {id} {
variable conc

return \$conc(\$id)
}

# main --
#     A stretch of river with a calamitous spill
#     For more realism: shallow sides with almost stagnant water
#
namespace import ::WQModel::*

substances pesticide

#
# The segments and their volumes
#
set length 1000.0  ;# segment is 1 kilometer long
set width   200.0  ;# 200 m wide
set depth     5.0  ;# 5 m deep

for { set i 0 } { \$i < 200 } { incr i 2 } {
set j [expr {\$i+1}]
segment \$i [expr {\$length*\$width*\$depth}]     [expr {\$length*\$width}]
segment \$j [expr {0.1*\$length*\$width*\$depth}] [expr {0.3*\$length*\$width}]

initial-values \$i 0.0
initial-values \$j 0.0
}

#
# The exchanges between main channel and shallow sides
#
set area     [expr {\$length*0.1/0.3}]
set distance [expr {0.5*(\$width+0.3*\$width)}]
set flow     0.0    ;# m3/s
set disp     0.001  ;# m2/s

for { set i 0 } { \$i < 200 } { incr i 2 } {
set j [expr {\$i+1}]
exchange \$i \$j \$flow \$disp \$area \$distance
}

#
# The exchanges along main channel
#
set area     [expr {\$width*\$depth}]
set distance \$length
set flow     100.0 ;# m3/s
set disp     0.001  ;# m2/s

for { set i -2 } { \$i < 200 } { incr i 2 } {
set j [expr {\$i+2}]

set from \$i
set to   \$j
if { \$i == -2 } {
set from "B1"
}
if { \$j >= 200 } {
set to   "B2"
}

exchange \$from \$to \$flow \$disp \$area \$distance
}

#
# Define the spill
#
boundary-values B1 0.0
boundary-values B2 0.0
initial-values 10  1.0

#
# Now compute the solution
#
global time
timeframe 0.0 10.0 0.01
initialise

while { [nextstep] } {
if { abs(10.0*\$time-int(10.0*\$time+0.0001)) < 0.0001 } {
puts "\$time [getconc 40] [getconc 80] [getconc 160]"
}
}```

 Category Mathematics Category Physics Category Science