Arjen Markus (20 february 2007) I would like to claim that the package presented below is a general package for solving partial differential equations, but that would be a gross exaggeration:
But on the other hand:
Well, the type of PDE it can solve can be expressed like this:
du d2u d2u du du
---- = f(x,y,u) * --- + g(x,y,u) * --- + h(x,y,u) * ---- + j(x,y,u) ---- + k(x,y,u)
dt dx2 dy2 dx dywhere u is the dependent variable or a set of such variables.
Lars H: Is there a reason for not allowing mixed second order derivatives? Even if they're not used in the basic PDEs, I think they tend to pop out if you have to do a nontrivial change of variables.
AM Hm, you are right about that. Mixed derivatives do have a tendency though to mess up the discretisation. I will need to think about it. They will certainly even be in the basic PDE if there is anisotropy.
Without further ado, here is the code I have so far (limited to 1D, no first-order derivatives yet). It relies on NAP for the heavy-duty computations.
# pde.tcl --
# Implementation of a basic solver for partial differential
# equations (PDE), using NAP
#
# This solver is meant for quasilinear parabolic PDEs, but
# it can handle systems of PDes. Via a simple transformation
# a hyperbolic equation, like the wave equation, can be
# written as a system of two or more equations that are first-order
# in the time.
#
# The current implementation is limited to equidistant grids in
# one spatial dimension. The types of boundary conditions are
# - Dirichlet (value given at the boundary)
# - Neumann (value of the flux given at the boundary)
# The same type is assumed for all equations, but the boundaries
# are independent.
#
# Because of the amount of information involved in solving PDEs
# we store everything in a separate namespace (this is partly
# required because of the properties of NAP). This is one of the
# tasks taken care of in the initialisation phase.
#
# TODO:
# nextTime
# getValues
#
package require nap
namespace eval PDE {
namespace import ::NAP::*
variable pdeCount 0
#
# Public procedures
#
namespace export initPde1D boundary timestep values initial
}
# initPde1D --
# Create a namespace to solve a one-dimensional PDE
# Arguments:
# None
# Result:
# The name of the namespace dedicated to the PDE
# Note:
# The result is required for accessing the PDE solution system.
# Pass it as the first argument.
#
proc ::PDE::initPde1D {} {
variable pdeCount
namespace eval $pdeCount [list \
variable _dimension 1 \
variable _pde ::PDE::$pdeCount \
]
set ns "::PDE::$pdeCount"
incr pdeCount
return $ns
}
# timestep --
# Set the value of the time step for integration
# Arguments:
# pde Name of the PDE solution system
# deltt Time step to be used
# Result:
# None
#
proc ::PDE::timestep {pde deltt} {
set [set pde]::_deltt $deltt
}
# printPde --
# Print the PDE information array
# Arguments:
# pde Name of the array holding all information
# Result:
# None
# Side effects:
# Prints the information contained in the array on stdout
#
proc ::PDE::printPde {pde} {
set [set pde]::_pde $pde
namespace eval $pde {
puts "PDE: $_pde"
puts " -- grid coordinates (cell boundaries) --"
puts "[$_grid value]"
puts " -- volumes of cells --"
puts "[$_volume value]"
puts " -- grid distances --"
puts "[$_dx value]"
if { [info exists _statevars] } {
puts " -- state variables --"
foreach v $_statevars {
puts "$v: [[set $v] value]"
}
}
}
}
# values --
# Return the values of a state variables as a list
# Arguments:
# pde Name of the array holding all information
# name Name of the state variable
# Result:
# None
#
proc ::PDE::values {pde name} {
::PDE::HasParameter $_pde $name
set [set pde]::_name $name
namespace eval $pde {
return [[set $_name] value]
}
}
# HasParameter --
# Check that a correct parameter name is given
# Arguments:
# pde Name of the PDE system
# name Parameter name (empty: any parameter)
# hasvalue Does it have a value? (optional)
# Result:
# 1 if there is a parameter by that name, error otherwise
#
proc ::PDE::HasParameter {pde name {hasvalue 0}} {
set [set pde]::_name $name
set [set pde]::_hasvalue $hasvalue
namespace eval $pde {
if { ! [info exists _statevars] } {
return -code error "No state variables defined yet"
}
if { $_name != {} } {
if { [lsearch $_statevars $_name] < 0 } {
return -code error "No such state variable: $_name"
}
if { $_hasvalue && ! [info exists $_name] } {
return -code error "State variable $_name has no values yet"
}
}
}
}
# HasGrid --
# Check that the PDE system has a grid
# Arguments:
# pde Name of the PDE system
# Result:
# 1 if there is a grid, error otherwise
#
proc ::PDE::HasGrid {pde} {
namespace eval $pde {
if { ! [info exists _grid] } {
return -code error "No grid defined yet"
}
}
}
# boundary --
# Set the type and value for the boundaries
# Arguments:
# pde Name of the array holding all information
# name Which parameter
# bound Which boundary (left, right, top or bottom)
# type Type of condition (dirichlet or neumann)
# value Value of the variable or the derivative
# Result:
# None
#
proc ::PDE::boundary {pde name bound type value} {
switch -re $bound {
"l|left" { set b left }
"r|right" { set b right }
"t|top" { set b top }
"b|bottom" { set b bottom }
default {
return -code error "Boundary must be left, right, top or bottom"
}
}
switch -re $type {
"d|dirichlet" { set t dirichlet }
"n|neumann" { set n neumann }
default {
return -code error "Type must be dirichlet or neumann"
}
}
set [set pde]::_btype($b,$name) $t
set [set pde]::_bvalue($b,$name) $value
}
# SetBoundaryValues --
# Set the correct values for the boundary elements
# Arguments:
# pde Name of the PDE system
# Result:
# None
#
proc ::PDE::SetBoundaryValues {pde} {
#
# TODO: Neumann boundaries!
#
namespace eval $pde {
set _maxcellp1 [expr {$_maxcell + 1}]
foreach _v $_statevars {
foreach _b {left right} _idx [list 0 $_maxcellp1] {
if { $_btype($_b,$_v) == "dirichlet" } {
[set $_v] set value $_bvalue($_b,$_v) $_idx
}
}
}
}
}
# initial --
# Set the initial value for a particular state variable
# Arguments:
# pde Name of the array holding all information
# name Which variable
# value Initial value
# Result:
# None
# Notes:
# Useful only when the grid has been defined and the state variables
#
proc ::PDE::initial {pde name value} {
set [set pde]::_name $name
set [set pde]::_value $value
HasParameter $pde $name
HasGrid $pde
namespace eval $pde {
set $_name $_value
nap "$_name = $_value + 0.0 * $_volume"
}
}
# zeroDerivs --
# Set the derivatives for all state variables to zero
# Arguments:
# pde Name of the array holding all information
# Result:
# None
#
proc ::PDE::zeroDerivs {pde} {
HasParameter $pde {}
HasGrid $pde
namespace eval $pde {
foreach name $_statevars {
nap "_d_$name = 0.0 * $_volume"
}
}
}
# addTerm --
# Add a term to the derivative of a particular state variable
# Arguments:
# pde Name of the array holding all information
# name Which variable
# term The value in each cell to be added
# Result:
# None
#
proc ::PDE::addTerm {pde name term} {
HasParameter $pde $name 1
HasGrid $pde
set [set pde]::_name $name
set [set pde]::_term $term
namespace eval $pde {
[set _d_$_name] \
set value [string map [list Name $_name] "_d_Name(1 .. _maxcell) + _term"] \
"(1 .. _maxcell)"
}
}
# addSecondDeriv --
# Compute the second order derivative of a particular state variable,
# multiplied by a factor (for convenience)
# Arguments:
# pde Name of the array holding all information
# name Which variable
# factor Expression for multiplication factor
# dname Which variable to take the derivative of
# Result:
# None
#
proc ::PDE::addSecondDeriv {pde name factor dname} {
set [set pde]::_name $name
set [set pde]::_dname $dname
set [set pde]::_factor $factor
HasParameter $pde $name 1
HasParameter $pde $dname 1
HasGrid $pde
SetBoundaryValues $pde
namespace eval $pde {
set _maxcellm1 [expr {$_maxcell-1}]
nap [string map [list Name $_dname] \
"_term1 = (Name(1 .. _maxcell) - Name(0 .. _maxcellm1)) /
_dx(1 .. _maxcell)"]
nap [string map [list Name $_dname] \
"_term2 = (Name(2 .. _nocells) - Name(1 .. _maxcell)) /
_dx(2 .. _nocells)"]
nap "_term = $_factor * (_term2 - _term1)"
[set _d_$_name] \
set value [string map [list Name $_name] "_d_Name(1 .. _maxcell) + _term"] \
"(1 .. _maxcell)"
}
}
# nextStep --
# Compute for all state variables the values in the next time step
# Arguments:
# pde Name of the array holding all information
# Result:
# None
#
proc ::PDE::nextStep {pde} {
HasParameter $pde {} 0
HasGrid $pde
namespace eval $pde {
foreach _v $_statevars {
nap "$_v = $_v + $_deltt * _d_$_v"
}
}
}
# statevars --
# Set the names of the state variables
# Arguments:
# pde Name of the array holding all information
# args Names of the variables
# Result:
# None
#
proc ::PDE::statevars {pde args} {
foreach name $args {
if { [string index $name 0] == "_" } {
return -code error "Names starting with an underscore are reserved - do not use them for state variables"
}
}
set [set pde]::_statevars $args
}
# grid1D --
# Create a one-dimensional grid
# Arguments:
# pde Name of the array holding all information
# nocells Number of grid cells
# xmin Minimum x-coordinate (left boundary)
# xmax Maximum x-coordinate (right boundary)
# Result:
# None
#
proc ::PDE::grid1D {pde nocells xmin xmax} {
#
# Create a namespace to hold the NAP variables
#
namespace eval $pde {
namespace import ::NAP::nap
}
#
# Check the consistency
#
if { $nocells < 1 } {
return -code error "Number of cells must be at least 1"
}
if { $xmin >= $xmax } {
return -code error "Minimum x-coordinate must be smaller than maximum"
}
set [set pde]::_nocells $nocells
set [set pde]::_maxcell [expr {$nocells-1}]
set [set pde]::_xmax $xmax
set [set pde]::_xmin $xmin
namespace eval $pde {
set _deltx [expr {($_xmax-$_xmin)/double($_nocells)}]
nap "_grid = _xmin + _deltx * (0 .. _nocells)"
nap "_dx = _deltx + 0.0 * (0 .. _maxcell)"
nap "_volume = _deltx + 0.0 * (0 .. _maxcell)"
}
}
# main --
# Testing:
# - Generation of the grid must work properly
#
set g [::PDE::initPde1D]
::PDE::grid1D $g 20 -10 10
::PDE::statevars $g u
::PDE::initial $g u 2.1
::PDE::boundary $g u left dirichlet 0
::PDE::boundary $g u right dirichlet 0
::PDE::timestep $g 0.1
::PDE::printPde $g
set i 0
while { $i < 10 } {
::PDE::zeroDerivs $g
::PDE::addTerm $g u 1.0
::PDE::addSecondDeriv $g u 1.0 u
::PDE::nextStep $g
incr i
}
::PDE::printPde $g