teaser.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python/FORTRAN interfacing teaser\n",
    "\n",
    "This document is a collection of general ideas and starting points on how to simplify the further development of ocean, atmosphere and climate models. It might serve as discussion basis for the development and application of best practices in ocean/climate model development.\n",
    "\n",
    "If you're experienced with both, Python and Fortran, you should be able to go on from here.\n",
    "\n",
    "This notebook is meant to provide examples on how to access compiled FORTRAN routines using Python on-board packages/features. It covers\n",
    "- Shared library compiling,\n",
    "- shared library function calling\n",
    "- and multi-dimensional array manipulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A subjective ocean modeler's perspective\n",
    "\n",
    "Ocean model code usually consists of a main program unit\n",
    "\n",
    "```fortran\n",
    "!--------------------------- main.f90 -----------------------------------\n",
    " PROGRAM main\n",
    "\n",
    "  ! declaration\n",
    "\n",
    "  ! initialization\n",
    "\n",
    "  ! main loop\n",
    "  do it = 1,itend\n",
    "    CALL update_external_forcing\n",
    "    CALL numerical_integration\n",
    "    CALL write_disk_output\n",
    "  end do\n",
    "\n",
    "  ! restart preparation\n",
    "    \n",
    " END PROGRAM main\n",
    "!------------------------------------------------------------------------\n",
    "```\n",
    "\n",
    "and many subroutines organized in \n",
    "\n",
    "```fortran\n",
    "!------------------------------ subroutines.f90 -------------------------\n",
    "\n",
    " ...\n",
    "\n",
    " SUBROUTINE square_root(a)\n",
    "  IMPLICIT NONE\n",
    "  real :: a\n",
    "  a = SQRT(a)\n",
    " END SUBROUTINE square_root\n",
    "\n",
    " ...\n",
    " \n",
    "!------------------------------------------------------------------------\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing during the development of ocean model code usually involves\n",
    "\n",
    "* changing the FORTRAN code\n",
    "* compiling of _complete_ model code\n",
    "* setup of a test experiment directory structure\n",
    "* job submission into an HPC queue\n",
    "* ... waiting ...\n",
    "* ... more waiting ...\n",
    "* analysing model output\n",
    "* drawing conclusions based on output\n",
    "* probably starting over ...\n",
    "\n",
    "It would be very useful if it could be shortened to\n",
    "\n",
    "* changing FORTRAN code\n",
    "* compiling only subroutine\n",
    "* using only the subroutine with data manipulation and plotting attached\n",
    "* drawing conclusions\n",
    "* starting over ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shared library compiling\n",
    "\n",
    "ask your compiler's documentation about how to do it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile subroutines.f90\n",
    "\n",
    "SUBROUTINE square(a)\n",
    "    IMPLICIT NONE\n",
    "    real :: a\n",
    "    a = a**2\n",
    "END SUBROUTINE square\n",
    "\n",
    "SUBROUTINE square_root(b)\n",
    "    IMPLICIT NONE\n",
    "    real :: b\n",
    "    b = SQRT(b)\n",
    "END SUBROUTINE square_root\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%script bash\n",
    "echo $(which gfortran) && gfortran --version\n",
    "gfortran -shared -fpic subroutines.f90 -o subroutines.so && ls -lrth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "unix tools are useful..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!readelf --symbols subroutines.so"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python/FORTRAN interfacing\n",
    "\n",
    "Don't worry, Numpy will handle most of this for you! However: https://docs.python.org/3/library/ctypes.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ctypes as ctypes\n",
    "\n",
    "fortran_routines = ctypes.CDLL('./subroutines.so')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "python_float = 4.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ctypes_float = ctypes.c_float(python_float)\n",
    "ctypes_pointer = ctypes.pointer(ctypes_float)\n",
    "\n",
    "_ = fortran_routines.square_(ctypes_pointer)\n",
    "print(ctypes_pointer.contents, ':', python_float ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens here?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = fortran_routines.square_root_(ctypes_pointer)\n",
    "print(ctypes_pointer.contents, ':', python_float ** 0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-dimensional array example\n",
    "\n",
    "* https://docs.scipy.org/doc/numpy/reference/routines.ctypeslib.html\n",
    "* https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.ctypes.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting subroutines.f90\n"
     ]
    }
   ],
   "source": [
    "%%writefile subroutines.f90\n",
    "\n",
    "SUBROUTINE square(a,ij,ik)  !! global vars are bad!\n",
    "    IMPLICIT NONE\n",
    "    integer :: ij,ik\n",
    "    real, dimension(ij,ik) :: a\n",
    "    a = a**2\n",
    "END SUBROUTINE square\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%script bash\n",
    "gfortran -shared -fpic subroutines.f90 -o subroutines.so && ls -lrth\n",
    "#readelf --symbols subroutines.so"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "fortran_routine_square = np.ctypeslib.load_library('subroutines.so', './').square_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_array = np.random.random((2, 2))\n",
    "print(target_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_original = target_array.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ensure Fortran array and ensure single precision\n",
    "target_array = np.asfortranarray(target_array, dtype=np.float32)\n",
    "print(target_array.flags)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ik = np.array(target_array.shape[0], dtype=np.int32)\n",
    "ij = np.array(target_array.shape[1], dtype=np.int32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call by reference! So no meaningful output\n",
    "fortran_routine_square(target_array.ctypes, ik.ctypes, ij.ctypes);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(target_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(target_original**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Caveats: Beware of datatypes!\n",
    "* https://docs.scipy.org/doc/numpy/user/basics.types.html\n",
    "* e.g. https://northstar-www.dartmouth.edu/doc/solaris-forte/manuals/fortran/prog_guide/11_cfort.html\n",
    "\n",
    "and use FORTRAN wrappers with ISO_C_BINDING declarations...\n",
    "\n",
    "* get your datatypes clear!\n",
    "* name mangling\n",
    "* pass-by-value functionality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting wrappers.f90\n"
     ]
    }
   ],
   "source": [
    "%%writefile wrappers.f90\n",
    "\n",
    "SUBROUTINE wrapper_square(a,ij,ik) bind(c, name='wrapper_square')\n",
    "    use iso_c_binding, only: c_int, c_float\n",
    "    IMPLICIT NONE\n",
    "    integer(c_int) :: ij,ik\n",
    "    real(c_float), dimension(ij,ik) :: a\n",
    "    call square(a,ij,ik)\n",
    "END SUBROUTINE wrapper_square\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 36K\n",
      "drwxr-xr-x 7 khoeflich khoeflich 4,0K Dez  5 11:30 RESOURCES\n",
      "-rw-r--r-- 1 khoeflich khoeflich 1,8K Dez  5 11:53 wrapper.txt\n",
      "-rw-r--r-- 1 khoeflich khoeflich  11K Dez  5 11:55 teaser.ipynb\n",
      "-rw-r--r-- 1 khoeflich khoeflich  159 Dez  5 11:55 subroutines.f90\n",
      "-rw-r--r-- 1 khoeflich khoeflich  253 Dez  5 11:55 wrappers.f90\n",
      "-rwxr-xr-x 1 khoeflich khoeflich 7,5K Dez  5 11:55 subroutines.so\n",
      "\n",
      "Symbol table '.dynsym' contains 12 entries:\n",
      "   Num:    Value          Size Type    Bind   Vis      Ndx Name\n",
      "     0: 0000000000000000     0 NOTYPE  LOCAL  DEFAULT  UND \n",
      "     1: 0000000000000000     0 NOTYPE  WEAK   DEFAULT  UND __cxa_finalize\n",
      "     2: 0000000000000000     0 NOTYPE  WEAK   DEFAULT  UND _ITM_registerTMCloneTable\n",
      "     3: 0000000000000000     0 NOTYPE  WEAK   DEFAULT  UND _ITM_deregisterTMCloneTab\n",
      "     4: 0000000000000000     0 NOTYPE  WEAK   DEFAULT  UND __gmon_start__\n",
      "     5: 0000000000201028     0 NOTYPE  GLOBAL DEFAULT   19 _edata\n",
      "     6: 0000000000201030     0 NOTYPE  GLOBAL DEFAULT   20 _end\n",
      "     7: 00000000000004a0     0 FUNC    GLOBAL DEFAULT    7 _init\n",
      "     8: 00000000000005ca   123 FUNC    GLOBAL DEFAULT   10 wrapper_square\n",
      "     9: 0000000000201028     0 NOTYPE  GLOBAL DEFAULT   20 __bss_start\n",
      "    10: 000000000000070c     0 FUNC    GLOBAL DEFAULT   11 _fini\n",
      "    11: 0000000000000645   196 FUNC    GLOBAL DEFAULT   10 square_\n",
      "\n",
      "Symbol table '.symtab' contains 51 entries:\n",
      "   Num:    Value          Size Type    Bind   Vis      Ndx Name\n",
      "     0: 0000000000000000     0 NOTYPE  LOCAL  DEFAULT  UND \n",
      "     1: 00000000000001c8     0 SECTION LOCAL  DEFAULT    1 \n",
      "     2: 00000000000001f0     0 SECTION LOCAL  DEFAULT    2 \n",
      "     3: 0000000000000230     0 SECTION LOCAL  DEFAULT    3 \n",
      "     4: 0000000000000350     0 SECTION LOCAL  DEFAULT    4 \n",
      "     5: 00000000000003e0     0 SECTION LOCAL  DEFAULT    5 \n",
      "     6: 0000000000000488     0 SECTION LOCAL  DEFAULT    6 \n",
      "     7: 00000000000004a0     0 SECTION LOCAL  DEFAULT    7 \n",
      "     8: 00000000000004c0     0 SECTION LOCAL  DEFAULT    8 \n",
      "     9: 00000000000004e0     0 SECTION LOCAL  DEFAULT    9 \n",
      "    10: 00000000000004f0     0 SECTION LOCAL  DEFAULT   10 \n",
      "    11: 000000000000070c     0 SECTION LOCAL  DEFAULT   11 \n",
      "    12: 0000000000000718     0 SECTION LOCAL  DEFAULT   12 \n",
      "    13: 0000000000000748     0 SECTION LOCAL  DEFAULT   13 \n",
      "    14: 0000000000200e50     0 SECTION LOCAL  DEFAULT   14 \n",
      "    15: 0000000000200e58     0 SECTION LOCAL  DEFAULT   15 \n",
      "    16: 0000000000200e60     0 SECTION LOCAL  DEFAULT   16 \n",
      "    17: 0000000000200fe0     0 SECTION LOCAL  DEFAULT   17 \n",
      "    18: 0000000000201000     0 SECTION LOCAL  DEFAULT   18 \n",
      "    19: 0000000000201020     0 SECTION LOCAL  DEFAULT   19 \n",
      "    20: 0000000000201028     0 SECTION LOCAL  DEFAULT   20 \n",
      "    21: 0000000000000000     0 SECTION LOCAL  DEFAULT   21 \n",
      "    22: 0000000000000000     0 FILE    LOCAL  DEFAULT  ABS crtstuff.c\n",
      "    23: 00000000000004f0     0 FUNC    LOCAL  DEFAULT   10 deregister_tm_clones\n",
      "    24: 0000000000000530     0 FUNC    LOCAL  DEFAULT   10 register_tm_clones\n",
      "    25: 0000000000000580     0 FUNC    LOCAL  DEFAULT   10 __do_global_dtors_aux\n",
      "    26: 0000000000201028     1 OBJECT  LOCAL  DEFAULT   20 completed.7697\n",
      "    27: 0000000000200e58     0 OBJECT  LOCAL  DEFAULT   15 __do_global_dtors_aux_fin\n",
      "    28: 00000000000005c0     0 FUNC    LOCAL  DEFAULT   10 frame_dummy\n",
      "    29: 0000000000200e50     0 OBJECT  LOCAL  DEFAULT   14 __frame_dummy_init_array_\n",
      "    30: 0000000000000000     0 FILE    LOCAL  DEFAULT  ABS wrappers.f90\n",
      "    31: 0000000000000000     0 FILE    LOCAL  DEFAULT  ABS subroutines.f90\n",
      "    32: 0000000000000000     0 FILE    LOCAL  DEFAULT  ABS crtstuff.c\n",
      "    33: 00000000000007ec     0 OBJECT  LOCAL  DEFAULT   13 __FRAME_END__\n",
      "    34: 0000000000000000     0 FILE    LOCAL  DEFAULT  ABS \n",
      "    35: 0000000000200e60     0 OBJECT  LOCAL  DEFAULT   16 _DYNAMIC\n",
      "    36: 0000000000201028     0 OBJECT  LOCAL  DEFAULT   19 __TMC_END__\n",
      "    37: 0000000000201020     0 OBJECT  LOCAL  DEFAULT   19 __dso_handle\n",
      "    38: 0000000000000718     0 NOTYPE  LOCAL  DEFAULT   12 __GNU_EH_FRAME_HDR\n",
      "    39: 0000000000201000     0 OBJECT  LOCAL  DEFAULT   18 _GLOBAL_OFFSET_TABLE_\n",
      "    40: 0000000000000000     0 NOTYPE  WEAK   DEFAULT  UND __cxa_finalize\n",
      "    41: 00000000000004a0     0 FUNC    GLOBAL DEFAULT    7 _init\n",
      "    42: 00000000000005ca   123 FUNC    GLOBAL DEFAULT   10 wrapper_square\n",
      "    43: 0000000000000000     0 NOTYPE  WEAK   DEFAULT  UND _ITM_registerTMCloneTable\n",
      "    44: 0000000000000000     0 NOTYPE  WEAK   DEFAULT  UND _ITM_deregisterTMCloneTab\n",
      "    45: 0000000000201028     0 NOTYPE  GLOBAL DEFAULT   20 __bss_start\n",
      "    46: 000000000000070c     0 FUNC    GLOBAL DEFAULT   11 _fini\n",
      "    47: 0000000000000645   196 FUNC    GLOBAL DEFAULT   10 square_\n",
      "    48: 0000000000201028     0 NOTYPE  GLOBAL DEFAULT   19 _edata\n",
      "    49: 0000000000201030     0 NOTYPE  GLOBAL DEFAULT   20 _end\n",
      "    50: 0000000000000000     0 NOTYPE  WEAK   DEFAULT  UND __gmon_start__\n"
     ]
    }
   ],
   "source": [
    "%%script bash\n",
    "gfortran -shared -fpic wrappers.f90 subroutines.f90 -o subroutines.so && ls -lrth\n",
    "readelf --symbols subroutines.so"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "fortran_routine_square = np.ctypeslib.load_library('subroutines.so','./').wrapper_square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_array = np.random.random((2,2))\n",
    "ik = np.array(target_array.shape[0],dtype=np.int32)\n",
    "ij = np.array(target_array.shape[1],dtype=np.int32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_original = target_array.copy()\n",
    "target_array = np.asfortranarray(target_array,dtype=np.float32)\n",
    "fortran_routine_square(target_array.ctypes,ik.ctypes,ij.ctypes);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.1843822  0.97012556]\n",
      " [0.08261189 0.784025  ]]\n",
      "[[0.18438221 0.97012553]\n",
      " [0.08261189 0.78402494]]\n"
     ]
    }
   ],
   "source": [
    "print(target_array)\n",
    "print(target_original**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Questions\n",
    "\n",
    "* possible to create thin interface to existing modules with ctypes only-use?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:py3-scipystack-20191202]",
   "language": "python",
   "name": "conda-env-py3-scipystack-20191202-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}