## package for solving PDEs

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:

• It is not ready for any practical use yet
• When/if it is finished it will still only deal with a particular class of PDEs
• It will use finite-volume methods on rectangular grids, so the PDE must be fit for such a treatment
• It will be limited to one temporal and two spatial dimensions for practical reasons

But on the other hand:

• You will be able to deal with systems of PDEs
• Both Dirichlet and Neumann boundary conditions are supported
• It will easy to define the particular PDEs

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              dy```

where 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"
}
}
}

#     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)"
}
}

#     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