Arjen Markus (7 december 2010) After a demonstration of "OpenFOAM" yesterday I knew I could not postpone this anymore: I have wanted to experiment with the Tensor package and NAP for a long time now, to see how easy it is to solve partial differential equations with them.

Here is my small-scale experiment with the Tensor package: solving the one-dimensional wave equation. The boundary conditions: left a prescribed excitation (h=h(t)) and right a fully reflexive wall (u=0). The function h(t) describes a compact wave train.

Not all is settled yet, but I thought I'd wikify it anyway.

# waves.tcl -- # Simulate a one-dimensional wave packet: # - Two quantities related via first-order # equations # - Initial condition: a block # - The boundary conditions are fully reflexive # - implemented via u=0 on the boundary # But the left boundary has a prescribed height # - The grid is staggered (u and h on different points) # # Note: # I am not fully satisfied with the implementation of the boundary conditions. # This requires further attention! # package require Tk package require Plotchart package require Tensor # nextTime -- # Compute the result for the next time # # Arguments: # None # # Side effects: # Tensor u and h updated # proc nextTime {} { set ucoeff [expr {$::velocity * $::deltt / $::deltx}] set hcoeff [expr {$::deltt / $::deltx}] set rfrict [expr {1.0 - $::friction * $::deltt}] set end [expr {$::number-1}] set endm1 [expr {$::number-2}] set endm2 [expr {$::number-3}] set endp1 $::number # # Bounds depend on the type of boundary condition - unfortunately # ::tensor::expr "u(i=1:$endm1) = $rfrict * u(i=1:$endm1) - $ucoeff * (h(i=1:$endm1) - h(i=0:$endm2))" ::tensor::expr "h(i=1:$end) = h(i=1:$end) - $hcoeff * (u(i=2:$endp1) - u(i=1:$end))" } # setBoundaryConditions -- # Set the boundary conditions - almost trivial in this case # # Arguments: # time Time in the computation # # Side effects: # Tensor u is updated # proc setBoundaryConditions {time} { set end $::number set endm1 [expr {$end-1}] # # A block boundary condition like this is too severe for the # simple scheme I implemented. So use a smoother one instead # if { $time < 2.0 } { set hvalue 1.0 } else { set hvalue 0.0 } if { $time < 2.0 } { set hvalue [expr {0.5*(1.0-cos(3.1415926*$time)) * sin(3.0*3.1415926*$time)}] } else { set hvalue 0.0 } h section 0 = scalar $hvalue h section $end = scalar [h section $endm1] u section 0 = scalar 0.0 u section $end = scalar 0.0 } # plotData -- # PLot the result # # Arguments: # time Time in the computation # values List of values # # Side effects: # Graph is updated # proc plotData {time values} { .c delete data .c delete title $::p title "Time = [format "%5.2f" $time]" set x [expr {0.5 * $::deltx}] foreach v $values { set x [expr {$x + $::deltx}] $::p plot data $x $v } $::p plot data {} {} ;# Break the line of data after 100 { set go 1 } vwait go } # main -- # Create the tensors, store a few # computational paramters and then # run the time loop # set velocity 1.0 set friction 0.0 set number 1000 set length 10.0 set deltx [expr {$length/$number}] set deltt 0.001 tensor::create u -size [expr {$number+1}] -type double tensor::create h -size [expr {$number+1}] -type double # # Initial condition # u = scalar 0.0 h = scalar 0.0 # # Loop over time # pack [canvas .c -bg white -width 500 -height 300] set p [::Plotchart::createXYPlot .c {0 10.0 1.0} {-2.0 2.0 0.5}] tkwait visibility .c set t 0 while { $t < 100000} { set time [expr {$deltt * $t}] if { $t % 200 == 0 } { plotData $time [h] } setBoundaryConditions $time nextTime incr t }

Here is a snapshot:

TV Yah yah, I know I've been into some interesting materials.

AM And here is another example: Solving a 2D diffusion problem with the Tensor package