[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 water quality processes
   * Add waste loads
   * 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

     #
     # Add the waste loads
     #
     # 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]
]]