Blog from October, 2011

Calling fortran from python

Calling fortran from the outside world. Using iso_c_bindings this is a bit easier. Here you can find some fortran examples and see how you can call these functions from python using ctypes and numpy.

The fortran code.

module test
use iso_c_binding
implicit none
integer, parameter :: MAXSTRLEN = 512
contains

! Utility functions
! fortran character(len=*) are not compatible with c
! To be compatible with c, strings sould be copied to a c_char array
function char_array_to_string(char_array, length)
  integer(c_int) :: length
  character(c_char) :: char_array(length)
  character(len=length) :: char_array_to_string
  integer :: i
  do i = 1, length
     char_array_to_string(i:i) = char_array(i)
  enddo
end function char_array_to_string

! C ends strings with a \0 character. Add this so it is received correctly in c compatible languages
function string_to_char_array(s, length)
  integer(c_int) :: length
  character :: s(*)
  character(c_char) :: string_to_char_array(length)
  integer :: i
  do i = 1, length
     string_to_char_array(i:i) = s(i)
  enddo
  string_to_char_array(i+1:i+1) = C_NULL_CHAR
end function string_to_char_array



! 1 int
integer(c_int) function oneint(arg1) bind(C, name="oneint")
  integer(c_int), intent(inout) :: arg1
  arg1 = 111
  oneint = 123
end function oneint

! 1 double
integer(c_int) function onedouble(arg1) bind(C, name="onedouble")
  real(c_double), intent(inout) :: arg1
  arg1 = 1.11d0
  onedouble = 123
end function onedouble


! 10by10 double
integer(c_int) function twobytwodouble(x) bind(C, name="twobytwodouble")
  real(c_double),intent(inout)  :: x(2,2)
  x = 4
  x(2,1) = 21
  x(1,2) = 12
  twobytwodouble = 123
end function twobytwodouble


! 10by10 double
integer(c_int) function twobythreedouble(x) bind(C, name="twobythreedouble")
  real(c_double),intent(inout)  :: x(2,3)
  x = 6
  x(2,1) = 21
  x(1,3) = 13
  twobythreedouble = 123
end function twobythreedouble

integer(c_int) function twobytwodoublepointer(ptr) bind(C, name="twobytwodoublepointer")
  type(c_ptr), intent(inout) :: ptr
  real(c_double), target, save  :: x(2,2)

  x = 4
  x(2,1) = 21
  x(1,2) = 12
  ptr=c_loc(x)
  twobytwodoublepointer = 123
end function twobytwodoublepointer


! 10by10 double pointer
integer(c_int) function twobythreedoublepointer(ptr) bind(C, name="twobythreedoublepointer")
  type(c_ptr), intent(inout) :: ptr
  ! Save is required here for the memory to remain available after the function call
  real(c_double), target, save  :: x(2,3)
  x = 6
  x(2,1) = 21
  x(1,3) = 13
  ptr = c_loc(x)
  twobythreedoublepointer = 123
end function twobythreedoublepointer

! character
integer(c_int) function letter(arg1) bind(C, name="letter")
  character(kind=c_char), intent(inout) :: arg1
  arg1 = 'W'
  letter = 123
end function letter

! string in (string in length is not fixed but internally you need to set a fixed string length)
integer(c_int) function stringin(arg1) bind(C, name="stringin")
  character(kind=c_char), intent(in) :: arg1(*)
  character(len=MAXSTRLEN) :: string
  string = char_array_to_string(arg1, MAXSTRLEN)
  write(*,*)string
  stringin = 123
end function stringin

! string out (requires fixed number of letters)
integer(c_int) function stringout(arg1) bind(C, name="stringout")
  ! Output string has to be fixed
  character(kind=c_char), intent(out) :: arg1(MAXSTRLEN)
  character(len=MAXSTRLEN) :: string
  string = "Hello from fortran"
  arg1 = string_to_char_array(string, len(trim(string)))
  stringout = 123
end function stringout


end module test

The corresponding python code.

#!/usr/bin/env python
import numpy as np
from ctypes import (CDLL, POINTER, ARRAY, c_void_p,
                    c_int, byref,c_double, c_char,
                    c_char_p, create_string_buffer)
from numpy.ctypeslib import ndpointer
import os

dllpath = os.path.abspath("test.dylib") # or .dll or .so
libtest = CDLL(dllpath)

# Define some extra types
# pointer to a double
c_double_p = POINTER(c_double)
# pointer to a integer
c_int_p = POINTER(c_int)

shape2x2=(2,2)
# Pointer to a 2x2 double in fortran layout
c_double2x2_c = ndpointer(shape=shape2x2, dtype="double", flags="C")
c_double2x2_f = ndpointer(shape=shape2x2, dtype="double", flags="FORTRAN")
# Pointer to a pointer to a 10x10 double in fortran layout
c_double2x2_f_p = POINTER(c_double2x2_f)
c_double2x2_c_p = POINTER(c_double2x2_c)
shape3x2=(3,2)
shape2x3=(2,3)
# Pointer to a 2x3,3x2 double in fortran layout
c_double2x3_c = ndpointer(shape=shape2x3, dtype="double", flags="C")
c_double2x3_f = ndpointer(shape=shape2x3, dtype="double", flags="FORTRAN")
c_double3x2_c = ndpointer(shape=shape3x2, dtype="double", flags="C")
c_double3x2_f = ndpointer(shape=shape3x2, dtype="double", flags="FORTRAN")
# Pointer to a pointer to a 2x3,3x2 double in fortran layout
c_double2x3_f_p = POINTER(c_double2x3_f)
c_double2x3_c_p = POINTER(c_double2x3_c)
c_double3x2_f_p = POINTER(c_double3x2_f)
c_double3x2_c_p = POINTER(c_double3x2_c)

# Pointer to a character pointer
c_char_p_p = POINTER(c_char_p)
MAXSTRLEN=512
# Character array (Fortran can only return c_char arrays in c compatible mode)
c_char_array = ARRAY(c_char,MAXSTRLEN)
# Pointer to a character array
c_char_array_p = POINTER(c_char_array)

# oneint
f = libtest.oneint
f.argtypes=[c_int_p]
arg1 = c_int(1)
rc=f(byref(arg1))
print arg1.value

# onedouble
f = libtest.onedouble
f.argtypes=[c_double_p]
arg1 = c_double(1)
rc=f(byref(arg1))
print arg1.value

# 2x2
f = libtest.twobytwodouble
f.argtypes=[c_double2x2_f]
arg1 = np.zeros(shape2x2, order="F")
rc=f(arg1)
arr = np.array(arg1)
print arr
print arr.flags


# 2x2 p 
f = libtest.twobytwodoublepointer
f.argtypes=[c_double2x2_c_p]
arg1 = c_double2x2_c()
rc=f(byref(arg1))
arr = np.array(arg1)
print arr
print arr.flags

# 2x3 
f = libtest.twobythreedouble
f.argtypes=[c_double2x3_f]
arg1 = np.zeros(shape2x3,order="F")
rc=f(arg1)
arr = np.array(arg1)
print arr
print arr.flags


# 2x3 corresponds to 3x2 p in C order, reversed from F.
f = libtest.twobythreedoublepointer
f.argtypes=[c_double3x2_c_p]
arg1 = c_double3x2_c()
rc=f(byref(arg1))
arr = np.array(arg1, order="C")
print arr
print arr.flags


# Exchange one letter
f = libtest.letter
f.argtypes=[c_char_p]
arg1 = c_char('H')
rc=f(byref(arg1))
print arg1.value

# Exchange a string (in)
f = libtest.stringin
f.argtypes=[c_char_array_p]
arg1 = create_string_buffer('Hello from python',MAXSTRLEN)
rc=f(byref(arg1))

# Exchange a string (out)
f = libtest.stringout
f.argtypes=[c_char_array_p]
arg1 = create_string_buffer('',MAXSTRLEN)
rc=f(byref(arg1))
print arg1.value


del libtest

The following paper was just published in NHESS. It was a cooperation between TUD, Deltares and TNO.

Abstract: For the design of cost-effective coastal defence a precise estimate is needed of the 1/10 000 per year storm surge. A more precise estimate requires more observations. Therefore, the three greatest storm surges that hit the northern part of the Holland Coast in the 18th century are reconstructed. The reconstructions are based on paintings, drawings, written records and shell deposits that have recently appeared. The storm-surge levels of these storms have been estimated using numerical modelling of the coastal processes. Here we show how these reconstructions can be used in combination with extreme value statistics to give a more confident estimate of low probability events.

Using 18th century storm-surge data from the Dutch Coast to improve the confidence in flood-risk estimates

This is the presentation I gave at the mini symposium. I discussed the trends in sea-level trend analysis.

  File Modified
PDF File present.pdf Presentation for the Dutch Geodetic Committee 12-10-2011 by Fedor Baart
XBeach XDune coupling

I'm currently working on the XDune XBeach coupling.

Tests

The test case, which Jaap set up, is located at https://repos.deltares.nl/repos/XBeach/branches/XD_XB_coupling/tests.

Progress

The following tasks can be separated:

  • (tick) Create an XBeach library https://repos.deltares.nl/repos/XBeach/branches/fedortmp
    I plan to merge this into the trunk, but it would be convenient for maintenance if we'd merge the libxbeach.F90 and xbeach.F90 and apply the library by default approach (xbeach.F90 depends on libxbeach.F90).
  • (tick) Create an XDune library https://repos.deltares.nl/repos/XBeach/branches/dune
    We still have to merge back changes to the main dune line, for better cooperation with the dune people. I have to check the revision and what Martijn Muller checked in.
  • (tick) Create a fortran and c compatible interface for XBeach. The libxbeach wrapper exposes all functions using iso_c_binding. The most incovenient of this is that you have to pass the string length to character functions, not sure how to prevent this.
  • (tick) Create a fortran and c compatible interface for XDune (https://repos.deltares.nl/repos/XBeach/branches/dune). I first created a simple c compatible library (wrapper.cc) but found that the trick I used for exposing pointers to c++ classes (as part of a c structure) was not compatible with fortran. I switched to a version with a global model and global parameters. This is ugly but I have not figured out a better way yet.
  • (tick) Create an ESMF wrapper for XBeach (gridded component) https://repos.deltares.nl/repos/XBeach/branches/XD_XB_coupling/esmf_xd_xb_coupler I created a gridded component earlier. It seems to work fine. I moved it from the xbeach directory to the xd_xb coupling directory to keep all esmf specific code together and out of the XBeach trunk.
  • (question) Create an ESMF wrapper for XDune (gridded component) https://repos.deltares.nl/repos/XBeach/branches/XD_XB_coupling/esmf_xd_xb_coupler. I have implemented most of the code. The only part that is missing is the grid. I have to look into how the grid is defined and map that to the ESMF types.
  • (tick) Create a coupler https://repos.deltares.nl/repos/XBeach/branches/XD_XB_coupling/esmf_xd_xb_coupler The coupler and app are running. Features that are missing are configurable time, output and timestep definitions
  • (error) Test data exchange. I'm still missing the grid of dune so I can't match the grids. I'm discussing with ESMF user mailling list on how to best implement the rotating grid.
  • (tick) Test timestep. Timesteps of both Dune and XBeach align (tested for cases where dune has timesteps >= esmf timestep)
  • Test on h4 (not done)

There was no latex style file available for the Journal of Coastal Research. Below you can find a .bst file that you can use to format your bibliography in the style of JCR.

You can use the file like this:

\bibliographystyle{jcr}	% (uses file "jcr.bst")
\bibliography{myrefs}		% expects file "myrefs.bib"

The file is made using the jcr.dbj file that is also attached. If you want to make some changes, change the jcr.dbj file and call:

latex jcr.dbj

This will result in a new jcr.bst file.

For formatting the rest of the latex file you'll also need the following macros:

\usepackage{times}
\renewcommand{\refname}{Literature cited} % custom heading for JCR
\renewcommand{\keywordsname}{ADDITIONAL INDEX WORDS:} % custom named index words

After using this your paper is pretty much well formatted, except from some details in section header formatting. I didn't know how to fix these so I did those by hand.

  File Modified
File jcr.bst Journal of Coastal Research bibtex style file 03-10-2011 by Fedor Baart
File jcr.dbj Journal of Coastal Research bibtex style docstrip source file 03-10-2011 by Fedor Baart