An not entirely trivial experiment with the CFFI package

Arjen Markus (10 april 2022) Having done some very simple experiments with Ashok's CFFI package (which comes with an amazingly comprehensive documentation) I decided to raise the stakes: what about using it for solving partial differential equations? Do the heavy numerics on the C or Fortran side and led Tcl do the control. The PDE of choice is a simple two-dimensional version of the advection-diffusion equation:

               2          2
   dC         d C        d C        dC
   --   = D ( ---    +   --- )  - v --
   dt           2          2        dx
              dx         dy

on a square grid (n x n cells) and simple initial and boundary conditions.

Of course, my favourite numerical language is Fortran, so let me first introduce the Fortran part of this experiment:

! advdiff.f90 --
!     Module implementing a straightforeward solver for an advection-diffusion
!     equation in two dimensions
!
!     Note:
!     I first used:  ccentre = ccentre + veldt * (...) + diffdt * (...)
!     but that gave odd results, with both gfortran and Intel Fortran oneAPI.
!     This needs to be investigated
!
module advdiff
    use iso_c_binding

    implicit none

    !
    ! Type holding the various parameters
    ! (alternatives: put them in separate arguments
    ! store them in module variables)
    !
    type, bind(C) :: param_def
        real(kind=c_double) :: diffusion      ! Diffusion coefficient (m2/s)
        real(kind=c_double) :: flow_velocity  ! Flow velocity (left to right; m/s)
        real(kind=c_double) :: deltax         ! Grid cell size (m)
        real(kind=c_double) :: deltat         ! Time step (s)
        integer             :: nreport        ! Number of steps to set before returning
    end type param_def

contains

subroutine update( n, n2, conc_in, conc_out, params ) bind(c, name = 'update')
!dec attributes dllexport :: update
    integer, value                   :: n
    integer, value                   :: n2            ! Dummy required for the Tcl side
    real(kind=c_double), intent(in)  :: conc_in(n,n)
    real(kind=c_double), intent(out) :: conc_out(n,n)
    type(param_def), intent(in)      :: params

    integer                          :: i
    real(kind=c_double)              :: diffdt, veldt

    real(kind=c_double)              :: dconc(n-2,n-2)

    !
    ! We want a compact calculation
    !
    conc_out = conc_in

    diffdt   = params%diffusion / params%deltax**2  * params%deltat
    veldt    = params%flow_velocity / params%deltax * params%deltat

    !
    ! Using associate we can formulate the difference equation in a way that avoids
    ! explicit loops and boundary conditions.
    ! Note:
    ! It is not necessarily very efficient ;)
    !
    ! By first copying the initial values and then using the output array, we
    ! avoid intermediate copying.
    !
    associate( ccentre => conc_out(2:n-1,2:n-1), &
               cleft   => conc_out(1:n-2,2:n-1), &
               cright  => conc_out(3:n  ,2:n-1), &
               cup     => conc_out(2:n-1,3:n  ), &
               cdown   => conc_out(2:n-1,1:n-2)  )

        do i = 1,params%nreport
           dconc    =           veldt  * ( cleft - ccentre ) + &
                                diffdt * ( cleft + cright + cup + cdown - 4.0_c_double * ccentre )

           ccentre  = ccentre + dconc
        enddo
    end associate

end subroutine update

end module advdiff

Note: I use the current Fortran 2008 standard, so some of the constructs may not be familiar to everyone. You can create a DLL or shared object using gfortran with this command:

gfortran -o advdiff.dll advdiff.f90 -shared

The routine update will update the incoming matrix conc_in and put the result in the outgoing matrix conc_out. To avoid unnecessary conversion steps, it does a number of steps before returning. This number is contained in the derived type params as nreport.

So far the Fortran side of things. Let us move on to the Tcl side:

  • Load the DLL or shared object via the ::cffi::Wrapper command
  • Define a convenient cffi structure for the derived type (C struct) to hold the numerical parameters
  • Define the interface of the function

This is done in the first few lines of the Tcl program. The rest is merely to:

  • Set up the matrix
  • Define a print procedure - cffi only knows about linear (non-nested) lists, so use that as the data structure
  • Do the time loop

And so we have the following Tcl program:

# model_advdiff.tcl --
#     Demonstration of cffi used on an advection-diffusion problem
#
set auto_path [list [file dirname [file dirname [file normalize [
        file join [info script] ...]]]] {*}$auto_path]
package require cffi

::cffi::Wrapper create advdiff advdiff.dll
::cffi::Struct create PARAMS {diffusion double flow_velocity double deltax double deltat double nreport int}
advdiff function update void {n int nsq int conc_in double[nsq] conc_out {double[nsq] out} params {struct.PARAMS byref} }

set params [dict create diffusion 0.1 flow_velocity 0.1 deltax 0.1 deltat 0.02 nreport 10]

#
# Set up the model area: 12x12, the outer edges are the boundaries
# Boundary: concentration = 0
# Interior: concentration = 1
#
set n 12
set nsq [expr {$n**2}]

set row [lrepeat $n 0.0]
set conc_in [lrepeat $n $row]

for {set row 1} {$row < $n-1} {incr row} {
    for {set col 1} {$col < $n-1} {incr col} {
        lset conc_in $row $col 1.0
    }
}

#
# cffi only support one-dimension arrays ...
#
set conc_in [concat {*}$conc_in]


# print_matrix --
#     Print procedure
#
# Arguments:
#     outfile      File to print to
#     n            The list represents a matrix of n by n
#     values       List of n x n values, as a single long list
#
proc print_matrix {outfile n values} {

    set i 0
    for {set row 0} {$row < $n} {incr row} {
        set line ""
        for {set col 0} {$col < $n} {incr col} {
            append line [format "%8.4f" [lindex $values $i]]
            incr i
        }
        puts $outfile $line
    }
}


# Run the calculation
set outfile [open "model_advdiff.out" w]

puts $outfile "Initial:"
print_matrix $outfile $n $conc_in

for {set time 0} {$time < 10} {incr time} {
    puts $outfile "\nTime: [expr {$time+1}]"

    update $n $nsq $conc_in conc_out $params

    print_matrix $outfile $n $conc_out

    set conc_in $conc_out
}

There are a few things to note:

  • As cffi knows nothing about n-dimensional arrays, I had to include a dummy argument nsq to define a matrix.
  • Care must be taken about the precise style of the arguments - pass by value or pass by reference etc. This is not always easy, especially with C pointers having very different roles without the syntax per se revealing that role. With Fortran similar issues appear.
  • I have been careful to use the C-Fortran interfacing features as per the Fortran 2003 standard. This allowed me to avoid all interfacing via C.
  • There is a quirk in my Fortran routine that I had to work around - see the comments. I do not entirely understand why this is required and that bugs me ;)

Just for fun, here is the output:

Initial:
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  1.0000  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 1
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.1226  0.2296  0.3045  0.3462  0.3639  0.3662  0.3547  0.3217  0.2546  0.1454  0.0000
  0.0000  0.2210  0.4137  0.5479  0.6223  0.6535  0.6576  0.6372  0.5784  0.4583  0.2619  0.0000
  0.0000  0.2843  0.5316  0.7030  0.7973  0.8365  0.8416  0.8159  0.7414  0.5881  0.3364  0.0000
  0.0000  0.3166  0.5915  0.7811  0.8848  0.9275  0.9329  0.9049  0.8230  0.6537  0.3742  0.0000
  0.0000  0.3288  0.6139  0.8101  0.9169  0.9607  0.9663  0.9375  0.8531  0.6780  0.3883  0.0000
  0.0000  0.3288  0.6139  0.8101  0.9169  0.9607  0.9663  0.9375  0.8531  0.6780  0.3883  0.0000
  0.0000  0.3166  0.5915  0.7811  0.8848  0.9275  0.9329  0.9049  0.8230  0.6537  0.3742  0.0000
  0.0000  0.2843  0.5316  0.7030  0.7973  0.8365  0.8416  0.8159  0.7414  0.5881  0.3364  0.0000
  0.0000  0.2210  0.4137  0.5479  0.6223  0.6535  0.6576  0.6372  0.5784  0.4583  0.2619  0.0000
  0.0000  0.1226  0.2296  0.3045  0.3462  0.3639  0.3662  0.3547  0.3217  0.2546  0.1454  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 2
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0624  0.1226  0.1736  0.2109  0.2324  0.2377  0.2261  0.1965  0.1480  0.0813  0.0000
  0.0000  0.1177  0.2313  0.3275  0.3977  0.4383  0.4482  0.4263  0.3705  0.2791  0.1533  0.0000
  0.0000  0.1609  0.3160  0.4474  0.5433  0.5986  0.6121  0.5822  0.5061  0.3813  0.2095  0.0000
  0.0000  0.1896  0.3724  0.5271  0.6400  0.7050  0.7208  0.6856  0.5960  0.4491  0.2468  0.0000
  0.0000  0.2037  0.4001  0.5662  0.6874  0.7572  0.7740  0.7362  0.6401  0.4823  0.2651  0.0000
  0.0000  0.2037  0.4001  0.5662  0.6874  0.7572  0.7740  0.7362  0.6401  0.4823  0.2651  0.0000
  0.0000  0.1896  0.3724  0.5271  0.6400  0.7050  0.7208  0.6856  0.5960  0.4491  0.2468  0.0000
  0.0000  0.1609  0.3160  0.4474  0.5433  0.5986  0.6121  0.5822  0.5061  0.3813  0.2095  0.0000
  0.0000  0.1177  0.2313  0.3275  0.3977  0.4383  0.4482  0.4263  0.3705  0.2791  0.1533  0.0000
  0.0000  0.0624  0.1226  0.1736  0.2109  0.2324  0.2377  0.2261  0.1965  0.1480  0.0813  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 3
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0398  0.0793  0.1146  0.1423  0.1597  0.1650  0.1571  0.1358  0.1014  0.0553  0.0000
  0.0000  0.0760  0.1515  0.2189  0.2718  0.3050  0.3151  0.3001  0.2594  0.1937  0.1056  0.0000
  0.0000  0.1056  0.2105  0.3042  0.3776  0.4237  0.4378  0.4170  0.3604  0.2691  0.1467  0.0000
  0.0000  0.1264  0.2519  0.3641  0.4520  0.5071  0.5239  0.4991  0.4313  0.3220  0.1756  0.0000
  0.0000  0.1371  0.2732  0.3949  0.4901  0.5500  0.5682  0.5412  0.4678  0.3492  0.1904  0.0000
  0.0000  0.1371  0.2732  0.3949  0.4901  0.5500  0.5682  0.5412  0.4678  0.3492  0.1904  0.0000
  0.0000  0.1264  0.2519  0.3641  0.4520  0.5071  0.5239  0.4991  0.4313  0.3220  0.1756  0.0000
  0.0000  0.1056  0.2105  0.3042  0.3776  0.4237  0.4378  0.4170  0.3604  0.2691  0.1467  0.0000
  0.0000  0.0760  0.1515  0.2189  0.2718  0.3050  0.3151  0.3001  0.2594  0.1937  0.1056  0.0000
  0.0000  0.0398  0.0793  0.1146  0.1423  0.1597  0.1650  0.1571  0.1358  0.1014  0.0553  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 4
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0272  0.0544  0.0792  0.0992  0.1121  0.1166  0.1115  0.0965  0.0721  0.0393  0.0000
  0.0000  0.0520  0.1043  0.1518  0.1901  0.2149  0.2234  0.2137  0.1850  0.1382  0.0753  0.0000
  0.0000  0.0726  0.1455  0.2119  0.2653  0.3000  0.3119  0.2982  0.2582  0.1928  0.1051  0.0000
  0.0000  0.0873  0.1749  0.2547  0.3188  0.3605  0.3748  0.3584  0.3103  0.2317  0.1263  0.0000
  0.0000  0.0949  0.1902  0.2769  0.3466  0.3920  0.4075  0.3897  0.3374  0.2519  0.1373  0.0000
  0.0000  0.0949  0.1902  0.2769  0.3466  0.3920  0.4075  0.3897  0.3374  0.2519  0.1373  0.0000
  0.0000  0.0873  0.1749  0.2547  0.3188  0.3605  0.3748  0.3584  0.3103  0.2317  0.1263  0.0000
  0.0000  0.0726  0.1455  0.2119  0.2653  0.3000  0.3119  0.2982  0.2582  0.1928  0.1051  0.0000
  0.0000  0.0520  0.1043  0.1518  0.1901  0.2149  0.2234  0.2137  0.1850  0.1382  0.0753  0.0000
  0.0000  0.0272  0.0544  0.0792  0.0992  0.1121  0.1166  0.1115  0.0965  0.0721  0.0393  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 5
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0190  0.0381  0.0556  0.0699  0.0793  0.0828  0.0794  0.0689  0.0515  0.0281  0.0000
  0.0000  0.0364  0.0730  0.1066  0.1340  0.1522  0.1587  0.1523  0.1321  0.0988  0.0539  0.0000
  0.0000  0.0508  0.1020  0.1490  0.1872  0.2126  0.2218  0.2128  0.1846  0.1381  0.0753  0.0000
  0.0000  0.0611  0.1227  0.1793  0.2253  0.2558  0.2669  0.2560  0.2221  0.1661  0.0906  0.0000
  0.0000  0.0665  0.1335  0.1951  0.2451  0.2783  0.2903  0.2785  0.2417  0.1807  0.0986  0.0000
  0.0000  0.0665  0.1335  0.1951  0.2451  0.2783  0.2903  0.2785  0.2417  0.1807  0.0986  0.0000
  0.0000  0.0611  0.1227  0.1793  0.2253  0.2558  0.2669  0.2560  0.2221  0.1661  0.0906  0.0000
  0.0000  0.0508  0.1020  0.1490  0.1872  0.2126  0.2218  0.2128  0.1846  0.1381  0.0753  0.0000
  0.0000  0.0364  0.0730  0.1066  0.1340  0.1522  0.1587  0.1523  0.1321  0.0988  0.0539  0.0000
  0.0000  0.0190  0.0381  0.0556  0.0699  0.0793  0.0828  0.0794  0.0689  0.0515  0.0281  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 6
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0133  0.0268  0.0392  0.0494  0.0562  0.0588  0.0565  0.0491  0.0368  0.0201  0.0000
  0.0000  0.0256  0.0515  0.0753  0.0948  0.1079  0.1128  0.1084  0.0942  0.0705  0.0385  0.0000
  0.0000  0.0358  0.0719  0.1052  0.1325  0.1508  0.1576  0.1515  0.1317  0.0986  0.0538  0.0000
  0.0000  0.0431  0.0866  0.1267  0.1595  0.1815  0.1897  0.1823  0.1585  0.1187  0.0648  0.0000
  0.0000  0.0469  0.0942  0.1378  0.1735  0.1974  0.2064  0.1984  0.1724  0.1291  0.0705  0.0000
  0.0000  0.0469  0.0942  0.1378  0.1735  0.1974  0.2064  0.1984  0.1724  0.1291  0.0705  0.0000
  0.0000  0.0431  0.0866  0.1267  0.1595  0.1815  0.1897  0.1823  0.1585  0.1187  0.0648  0.0000
  0.0000  0.0358  0.0719  0.1052  0.1325  0.1508  0.1576  0.1515  0.1317  0.0986  0.0538  0.0000
  0.0000  0.0256  0.0515  0.0753  0.0948  0.1079  0.1128  0.1084  0.0942  0.0705  0.0385  0.0000
  0.0000  0.0133  0.0268  0.0392  0.0494  0.0562  0.0588  0.0565  0.0491  0.0368  0.0201  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 7
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0094  0.0190  0.0278  0.0350  0.0399  0.0418  0.0402  0.0350  0.0262  0.0143  0.0000
  0.0000  0.0181  0.0364  0.0533  0.0672  0.0765  0.0801  0.0771  0.0671  0.0503  0.0275  0.0000
  0.0000  0.0253  0.0509  0.0745  0.0939  0.1070  0.1120  0.1078  0.0938  0.0703  0.0384  0.0000
  0.0000  0.0304  0.0612  0.0897  0.1130  0.1288  0.1348  0.1297  0.1128  0.0846  0.0462  0.0000
  0.0000  0.0331  0.0666  0.0976  0.1230  0.1401  0.1467  0.1411  0.1228  0.0920  0.0503  0.0000
  0.0000  0.0331  0.0666  0.0976  0.1230  0.1401  0.1467  0.1411  0.1228  0.0920  0.0503  0.0000
  0.0000  0.0304  0.0612  0.0897  0.1130  0.1288  0.1348  0.1297  0.1128  0.0846  0.0462  0.0000
  0.0000  0.0253  0.0509  0.0745  0.0939  0.1070  0.1120  0.1078  0.0938  0.0703  0.0384  0.0000
  0.0000  0.0181  0.0364  0.0533  0.0672  0.0765  0.0801  0.0771  0.0671  0.0503  0.0275  0.0000
  0.0000  0.0094  0.0190  0.0278  0.0350  0.0399  0.0418  0.0402  0.0350  0.0262  0.0143  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 8
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0067  0.0134  0.0197  0.0248  0.0283  0.0297  0.0286  0.0249  0.0186  0.0102  0.0000
  0.0000  0.0128  0.0258  0.0378  0.0476  0.0543  0.0569  0.0548  0.0477  0.0358  0.0195  0.0000
  0.0000  0.0179  0.0360  0.0528  0.0666  0.0759  0.0795  0.0766  0.0667  0.0500  0.0273  0.0000
  0.0000  0.0216  0.0434  0.0636  0.0802  0.0914  0.0957  0.0922  0.0803  0.0602  0.0329  0.0000
  0.0000  0.0235  0.0472  0.0692  0.0872  0.0994  0.1042  0.1003  0.0873  0.0655  0.0358  0.0000
  0.0000  0.0235  0.0472  0.0692  0.0872  0.0994  0.1042  0.1003  0.0873  0.0655  0.0358  0.0000
  0.0000  0.0216  0.0434  0.0636  0.0802  0.0914  0.0957  0.0922  0.0803  0.0602  0.0329  0.0000
  0.0000  0.0179  0.0360  0.0528  0.0666  0.0759  0.0795  0.0766  0.0667  0.0500  0.0273  0.0000
  0.0000  0.0128  0.0258  0.0378  0.0476  0.0543  0.0569  0.0548  0.0477  0.0358  0.0195  0.0000
  0.0000  0.0067  0.0134  0.0197  0.0248  0.0283  0.0297  0.0286  0.0249  0.0186  0.0102  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 9
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0047  0.0095  0.0140  0.0176  0.0201  0.0211  0.0203  0.0177  0.0132  0.0072  0.0000
  0.0000  0.0091  0.0183  0.0268  0.0338  0.0386  0.0404  0.0389  0.0339  0.0254  0.0139  0.0000
  0.0000  0.0127  0.0256  0.0375  0.0473  0.0539  0.0565  0.0544  0.0474  0.0355  0.0194  0.0000
  0.0000  0.0153  0.0308  0.0451  0.0569  0.0649  0.0680  0.0655  0.0570  0.0428  0.0234  0.0000
  0.0000  0.0166  0.0335  0.0490  0.0619  0.0706  0.0740  0.0713  0.0621  0.0465  0.0254  0.0000
  0.0000  0.0166  0.0335  0.0490  0.0619  0.0706  0.0740  0.0713  0.0621  0.0465  0.0254  0.0000
  0.0000  0.0153  0.0308  0.0451  0.0569  0.0649  0.0680  0.0655  0.0570  0.0428  0.0234  0.0000
  0.0000  0.0127  0.0256  0.0375  0.0473  0.0539  0.0565  0.0544  0.0474  0.0355  0.0194  0.0000
  0.0000  0.0091  0.0183  0.0268  0.0338  0.0386  0.0404  0.0389  0.0339  0.0254  0.0139  0.0000
  0.0000  0.0047  0.0095  0.0140  0.0176  0.0201  0.0211  0.0203  0.0177  0.0132  0.0072  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000

Time: 10
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.0000  0.0034  0.0068  0.0099  0.0125  0.0143  0.0150  0.0144  0.0125  0.0094  0.0051  0.0000
  0.0000  0.0064  0.0130  0.0190  0.0240  0.0274  0.0287  0.0276  0.0241  0.0181  0.0099  0.0000
  0.0000  0.0090  0.0181  0.0266  0.0335  0.0383  0.0401  0.0386  0.0337  0.0253  0.0138  0.0000
  0.0000  0.0108  0.0218  0.0320  0.0404  0.0460  0.0483  0.0465  0.0405  0.0304  0.0166  0.0000
  0.0000  0.0118  0.0237  0.0348  0.0439  0.0501  0.0525  0.0506  0.0441  0.0331  0.0181  0.0000
  0.0000  0.0118  0.0237  0.0348  0.0439  0.0501  0.0525  0.0506  0.0441  0.0331  0.0181  0.0000
  0.0000  0.0108  0.0218  0.0320  0.0404  0.0460  0.0483  0.0465  0.0405  0.0304  0.0166  0.0000
  0.0000  0.0090  0.0181  0.0266  0.0335  0.0383  0.0401  0.0386  0.0337  0.0253  0.0138  0.0000
  0.0000  0.0064  0.0130  0.0190  0.0240  0.0274  0.0287  0.0276  0.0241  0.0181  0.0099  0.0000
  0.0000  0.0034  0.0068  0.0099  0.0125  0.0143  0.0150  0.0144  0.0125  0.0094  0.0051  0.0000
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000