Simple example to compute derivatives of any Fortran routine with DNAD

James C. Orr$^1$, Jean-Marie Epitalon$^2$, and James Kermode$^3$

$^1$Laboratoire des Sciences du Climat et de l'Environnement/IPSL, CEA-CNRS-UVSQ, Gif-sur-Yvette, France
$^2$Geoscientific Programming Services, Fanjeax, France
$^3$ Warwick Centre for Predictive Modelling, University of Warwick, UK

21 September 2015

Do you want to compute derivatives from an existing fortran code simply and accurately? If so, read on.

There is a simple and accurate way to compute first derivatives from an existing Fortran subroutine with minimal code modification. The method, called Dual Number Automatic Differentiation (DNAD), is described by Yu and Blair (2013). The DNAD approach yields derivatives $\partial y_j / \partial x_i$, where $x_i$ are all input variables ($i = 1,n$) and $y_j$ are all output variables ($j = 1,m$).

To do compute the derivatives, one simply needs to change the TYPE definition of the variables (and import the DNAD module). Results are as accurate as the analytical solution, i.e., to machine precision.

For this little demo, we first wrapped the fortran routines in f90wrap to access them in python. To do that, just download the files in the same directory where you found this Jupyter notebook file, and then type make. Then just execute the cells below, where we show how to run the wrapped code in python.

Specify working directory



In [13]:

    
# Working directory (change as needed) 
#cylynder_dnad_dir = "/home/my-user-name/etc/"
cylynder_dnad_dir = "."
import sys
sys.path.append(cylynder_dnad_dir)

Use print() from Python3 instead of print from Python2



In [14]:

    
from __future__ import print_function

Import `numpy` and the fortran routines (including the `cylinder` demo and the `DNAD` module)



In [ ]:

    
import numpy as np
import Example # or Example_pkg, as you prefer

Specify definitions to use later



In [ ]:

    
# Some definitions for later use (thanks to James Kermode)
d1 = Example.Dual_Num_Auto_Diff.Dual_Num()
d2 = Example.Dual_Num_Auto_Diff.Dual_Num()
d3 = Example.Mcyldnad.cyldnad(d1, d2)

Some documentation



In [16]:

    
cylin=Example.Mcyldnad.cyldnad
print(cylin.__doc__)









    



        vol = cyldnad(radius, height)
        
        
        Defined at cyldnad.fpp lines 8-16
        
        Parameters
        ----------
        radius : Dual_Num
        height : Dual_Num
        
        Returns
        -------
        vol : Dual_Num

Initialize d1 and d2 variables, each of the type dual_num (defined in the `DNAD` module):



In [18]:

    
# Specify radius (r)
d1.x_ad_ = 3
# Specify that we want dv/dr, where v is cylinder volume and r is cylinder radius
d1.xp_ad_ =np.array((1.,0.))
print("d1:", d1)

# Specify height (h)
d2.x_ad_ = 5
# Specify that we want dv/dh, where h is cylinder height
d2.xp_ad_ = np.array((0,1.))
print("d2:", d2)









    



d1: <dual_num>{
    x_ad_ : 3.0,
    xp_ad_ : array([ 1.,  0.])}
d2: <dual_num>{
    x_ad_ : 5.0,
    xp_ad_ : array([ 0.,  1.])}

Run subroutine "cylinder" to compute volume $v$, $dv/dr$, and $dv/dh$



In [20]:

    
d3 = Example.Mcyldnad.cyldnad(d1, d2)
# Print computed v, dv/dr, dv/dh (thanks to dual numbers)
print("result:", d3)









    



result: <dual_num>{
    x_ad_ : 141.3716694115407,
    xp_ad_ : array([ 94.24777961,  28.27433388])}