# 7.6.2.2. Tutorial 2 : Using Tapenade for a mountain glacier model¶

In this tutorial, we will describe the steps to use Tapenade to develop the tangent linear and adjoint codes for a simple PDE-based mountain glacier model. The model itself is inspired from the book Fundamentals of Glacier Dynamics, by CJ van der Veen and has been used for glaciology summer schools before, for example here , although they used a different tool to develop the adjoint model.

As prerequisites, ensure you have Tapenade installed and gfortran compiler available.

## 7.6.2.2.1. Mountain glacier model¶

### Model description¶

The equations for a simple, 1D mountain glacier are described here.

$\frac{\partial{H}}{\partial{t}} = -\frac{\partial}{\partial x}\left(-D(x)\frac{\partial{h}}{\partial{x}}\right) + M$

where

$D(x) = CH^{n+2}\left|\frac{\partial h}{\partial x}\right|^{n-1}$

and

$C = \frac{2A}{n+2}(\rho g)^n$

Here $$H(x,t)$$ is the thickness of the ice sheet, $$h(x,t)$$ is the height of the top surface of the ice sheet (taken from some reference height), $$b(x)$$ is the height of the base of the ice sheet (taken from the same reference height). These three parameters can be related as

$H(x,t) = h(x,t) - b(x)$

In essence it’s a highly non-linear diffusion equation with a source term. Furthermore, the boundary conditions are given by

$H_{\rm left} = 0$

and

$H_{\rm right} = 0$

### Model parameters¶

We assume that the basal topography is linear in $$x$$ such that $$\frac{\partial{b}}{\partial{x}} = -0.1$$. $$M$$ is the accumulation rate and is modeled to be constant in time, but varying linearly in space.

$M(x) = M_0 - x M_1$

where $$M_0 = 4.0 \:{\rm m/yr}, \:M_1 = 0.0002 \:{\rm yr}^{-1}$$.

The acceleration due to gravity $$g = 9.2 \:{\rm m/s}^2$$. The density of ice $$\rho = 920 \:{\rm kg/m}^3$$. The exponent $$n = 3$$ comes from Glen’s flow law. We take $$A = 10^{-16} \: {\rm Pa}^{-3} {\rm a}^{-1}$$. Our domain length is $$L = 30 \:{\rm km}$$. The total time we run the simulation for is $$T = 5000 \:{\rm years}$$.

### Discretization and Finite Volumes solution¶

The equations above cannot be solved analytically. Here we discuss a Finite Volumes method to solve them. We use a staggered grid such that $$D$$ is computed at the centres of the grid, and $$H, h$$ are computed at the vertices.

$D\left(x_{j+1/2}\right) = C\left(\frac{H_j + H_{j+1}}{2}\right)^{n+2}\left(\frac{h_{j+1}-h_j}{\Delta x}\right)^{n-1}$
$\phi_{j+1/2} = D\left(x_{j+1/2}\right)\frac{\partial{h}}{\partial{x}}, \:\text{where}\: \frac{\partial h}{\partial x} = \frac{h_{j+1}-h_j}{\Delta x}$
$\frac{\partial H_j}{\partial t} = - \frac{\phi_{j+1/2}-\phi_{j-1/2}}{\Delta x} + M(x_j)$

Here, $$\phi$$ indicates flux and is just an intermediate variable.

We step forward in time using the simple Euler forward scheme. Here we choose $$\Delta x = 1.0\:{\rm km}, \Delta t = 1.0\:{\rm month} = \frac{1}{12}\:{\rm years}$$.

### What sensitivities are we interested in?¶

Our Quantity of Interest (QoI) is the total volume of the ice sheet after 5000 years (defined as V in the code). The control variables are the spatially varying accumulation rate. We are interested in the sensitivities of the QoI to the perturbations in the accumulation rate at various locations along the ice sheet. (defined as xx in the code)

## 7.6.2.2.2. Forward Model code¶

The code for the non-linear forward model is given here.

forward.f90

module forward

implicit none
real(8), parameter :: xend = 30
real(8), parameter :: dx = 1
integer, parameter :: nx = int(xend/dx)

contains

subroutine forward_problem(xx,V)

implicit none

real(8), parameter :: rho = 920.0
real(8), parameter :: g = 9.2
real(8), parameter :: n = 3
real(8), parameter :: A = 1e-16
real(8), parameter :: dt = 1/12.0
real(8), parameter :: C = 2*A/(n+2)*(rho*g)**n*(1.e3)**n
real(8), parameter :: tend = 5000
real(8), parameter :: bx = -0.0001
real(8), parameter :: M0 = .004, M1 = 0.0002
integer, parameter :: nt = int(tend/dt)
real(8), dimension(nx+1,nt+1) :: h
real(8), dimension(nx+1,nt+1) :: h_capital
integer :: t,i
real(8), dimension(nx+1) :: xarr
real(8), dimension(nx+1) :: M
real(8), dimension(nx+1) :: b
real(8), dimension(nx) :: D, phi
real(8), intent(in), dimension(nx+1) :: xx
real(8), intent(out) :: V

xarr = (/ ((i-1)*dx, i=1,nx+1) /)
M = (/ (M0-(i-1)*dx*M1, i=1,nx+1) /)
b = (/ (1.0+bx*(i-1)*dx, i=1,nx+1) /)
M = M + xx
h(1,:) = b(1)
h(:,1) = b
h(nx+1,:) = b(nx+1)
h_capital(1,:) = h(1,:) - b(1)
h_capital(nx+1,:) = h(nx+1,:) - b(nx+1)
h_capital(:,1) = h(:,1) - b

do t = 1,nt
D(:) = C * ((h_capital(1:nx,t)+h_capital(2:nx+1,t))/2)**(n+2) * ((h(2:nx+1,t) - h(1:nx,t))/dx)**(n-1)
phi(:) = -D(:)*(h(2:nx+1,t)-h(1:nx,t))/dx
h(2:nx,t+1) = h(2:nx,t) + M(2:nx)*dt - dt/dx * (phi(2:nx)-phi(1:nx-1))

where (h(:,t+1) < b)
h(:,t+1) = b
end where

h_capital(:,t+1) = h(:,t+1) - b

end do

V = 0.

open (unit = 2, file = "results_forward_run.txt", action="write",status="replace")
write(2,*) "         #                H                h                  b"
write(2,*) "_______________________________________________________________________________"

do i = 1, size(h_capital(:,nt+1))
V = V + h_capital(i,nt+1)*dx

write(2,*) i, "    ", h_capital(i,nt+1), "    ", h(i,nt+1), "    ", b(i)
end do

close(2)
end subroutine forward_problem

end module forward


## 7.6.2.2.3. Finite Differences for validation¶

We can use finite differences with the forward differencing scheme (or alternatively, a more accurate central differencing scheme mentioned in Tutorial 1 ). to get the gradient evaluated at the given value of the control vector $$\mathbf{x}$$. For a $$N$$ dimensional control vector, the forward model needs to be called $$N+1$$ times in the forward differencing scheme to get an approximate gradient. This is because we can only evaluate directional derivatives, which we need to do $$N$$ times to get the gradient. Each directional derivative evaluation, in turn, requires 1 perturbed forward model evaluation and 1 unperturbed forward model evaluation (this can be done only once and stored).

Mathematically, this can be formulated as ($$\mathcal{J}$$ is the objective function, in our case the total volume denoted by V in the code)

$\left (\nabla_{\mathbf{x}} \mathcal{J}(\mathbf{x}), \hat{\mathbf{e}} \right) \approx \frac{\mathcal{J}\left(\mathbf{x}+\epsilon \hat{\mathbf{e}}\right) - \mathcal{J}\left(\mathbf{x}\right)}{\epsilon} + \mathcal{O}(\epsilon)$

Here $$(.,.)$$ indicates the normal inner product of two discrete vectors in $$\mathbb{R}^N$$.

The left side is the directional derivative of $$\mathcal{J}$$ at $$\mathbf{x}$$ in the direction $$\hat{\mathbf{e}}$$, while the right side is the finite differences approximation of the directional derivative.

The choice of $$\epsilon$$ can be critical, however we do not discuss that here and simply choose a value $$\epsilon = 1.e-7$$. In general, one would perform a convergence analysis with a range of values of $$\epsilon = {0.1, 0.01, 0.001, ...}$$.

The code for the finite differences can be found in the Driver Routine.

## 7.6.2.2.4. Tangent Linear model¶

The tangent linear model is described in detail in Tutorial 1. Another excellent resource is the MITgcm documentation.

### Generating Tangent Linear model code using Tapenade¶

It is pretty simple to generate the tangent linear model for our forward model.

% tapenade -tangent -tgtmodulename %_tgt -head "forward_problem(V)/(xx)" forward.f90


-tangent tells Tapenade that we want a tangent linear code.

-tgtmodulename %_tgt tells Tapenade to suffix the differentiated modules with _tgt.

-head "forward_problem(V)/(xx) tells Tapenade that the head subroutine is forward_problem, the dependent variable or the objective function is V and the independent variable is xx.

This generates 2 files -

• forward_d.f90 - This file contains the tangent linear module (called forward_tgt) as well as the original forward code module. forward_tgt) contains the forward_problem_d subroutine which is our tangent linear model. All the differential variables are suffixed with the alphabet d.

!        Generated by TAPENADE     (INRIA, Ecuador team)
!  TAPENADE 3.15 (master) - 15 Apr 2020 13:12
!
MODULE FORWARD_TGT
IMPLICIT NONE
REAL*8, PARAMETER :: xend=30
REAL*8, PARAMETER :: dx=1
INTEGER, PARAMETER :: nx=INT(xend/dx)

CONTAINS
!  Differentiation of forward_problem in forward (tangent) mode:
!   variations   of useful results: v
!   with respect to varying inputs: xx
!   RW status of diff variables: v:out xx:in
SUBROUTINE FORWARD_PROBLEM_D(xx, xxd, v, vd)
IMPLICIT NONE
REAL*8, PARAMETER :: rho=920.0
REAL*8, PARAMETER :: g=9.2
REAL*8, PARAMETER :: n=3
REAL*8, PARAMETER :: a=1e-16
REAL*8, PARAMETER :: dt=1/12.0
REAL*8, PARAMETER :: c=2*a/(n+2)*(rho*g)**n*1.e3**n
REAL*8, PARAMETER :: tend=5000
REAL*8, PARAMETER :: bx=-0.0001
REAL*8, PARAMETER :: m0=.004, m1=0.0002
INTRINSIC INT
INTEGER, PARAMETER :: nt=INT(tend/dt)
REAL*8, DIMENSION(nx+1, nt+1) :: h
REAL*8, DIMENSION(nx+1, nt+1) :: hd
REAL*8, DIMENSION(nx+1, nt+1) :: h_capital
REAL*8, DIMENSION(nx+1, nt+1) :: h_capitald
INTEGER :: t, i
REAL*8, DIMENSION(nx+1) :: xarr
REAL*8, DIMENSION(nx+1) :: m
REAL*8, DIMENSION(nx+1) :: md
REAL*8, DIMENSION(nx+1) :: b
REAL*8, DIMENSION(nx) :: d, phi
REAL*8, DIMENSION(nx) :: dd, phid
REAL*8, DIMENSION(nx+1), INTENT(IN) :: xx
REAL*8, DIMENSION(nx+1), INTENT(IN) :: xxd
REAL*8, INTENT(OUT) :: v
REAL*8, INTENT(OUT) :: vd
INTRINSIC SIZE
REAL*8, DIMENSION(nx) :: pwx1
REAL*8, DIMENSION(nx) :: pwx1d
REAL*8 :: pwy1
REAL*8, DIMENSION(nx) :: pwr1
REAL*8, DIMENSION(nx) :: pwr1d
REAL*8, DIMENSION(nx) :: pwx2
REAL*8, DIMENSION(nx) :: pwx2d
REAL*8 :: pwy2
REAL*8, DIMENSION(nx) :: pwr2
REAL*8, DIMENSION(nx) :: pwr2d
REAL*8, DIMENSION(nx) :: temp
xarr = (/((i-1)*dx, i=1,nx+1)/)
m = (/(m0-(i-1)*dx*m1, i=1,nx+1)/)
b = (/(1.0+bx*(i-1)*dx, i=1,nx+1)/)
md = xxd
m = m + xx
h(1, :) = b(1)
h(:, 1) = b
h(nx+1, :) = b(nx+1)
h_capital(1, :) = h(1, :) - b(1)
h_capital(nx+1, :) = h(nx+1, :) - b(nx+1)
h_capital(:, 1) = h(:, 1) - b
h_capitald = 0.0_8
hd = 0.0_8
DO t=1,nt
pwx1d = (h_capitald(1:nx, t)+h_capitald(2:nx+1, t))/2
pwx1 = (h_capital(1:nx, t)+h_capital(2:nx+1, t))/2
pwy1 = n + 2
WHERE (pwx1 .LE. 0.0 .AND. (pwy1 .EQ. 0.0 .OR. pwy1 .NE. INT(pwy1)&
&         ))
pwr1d = 0.0_8
ELSEWHERE
pwr1d = pwy1*pwx1**(pwy1-1)*pwx1d
END WHERE
pwr1 = pwx1**pwy1
pwx2d = (hd(2:nx+1, t)-hd(1:nx, t))/dx
pwx2 = (h(2:nx+1, t)-h(1:nx, t))/dx
pwy2 = n - 1
WHERE (pwx2 .LE. 0.0 .AND. (pwy2 .EQ. 0.0 .OR. pwy2 .NE. INT(pwy2)&
&         ))
pwr2d = 0.0_8
ELSEWHERE
pwr2d = pwy2*pwx2**(pwy2-1)*pwx2d
END WHERE
pwr2 = pwx2**pwy2
dd(:) = c*(pwr2*pwr1d+pwr1*pwr2d)
d(:) = c*pwr1*pwr2
temp = h(2:nx+1, t) - h(1:nx, t)
phid(:) = -(temp*dd(:)/dx+d(:)*(hd(2:nx+1, t)-hd(1:nx, t))/dx)
phi(:) = -(d(:)/dx*temp)
hd(2:nx, t+1) = hd(2:nx, t) + dt*md(2:nx) - dt*(phid(2:nx)-phid(1:&
&       nx-1))/dx
h(2:nx, t+1) = h(2:nx, t) + m(2:nx)*dt - dt/dx*(phi(2:nx)-phi(1:nx&
&       -1))
WHERE (h(:, t+1) .LT. b)
hd(:, t+1) = 0.0_8
h(:, t+1) = b
END WHERE
h_capitald(:, t+1) = hd(:, t+1)
h_capital(:, t+1) = h(:, t+1) - b
END DO
v = 0.
OPEN(unit=2, file='results_forward_run.txt', action='write', status=&
&  'replace')
WRITE(2, *) &
&   '         #                H                h                  b'
WRITE(2, *) '_______________________________________________________&
&________________________'
vd = 0.0_8
DO i=1,SIZE(h_capital(:, nt+1))
vd = vd + dx*h_capitald(i, nt+1)
v = v + h_capital(i, nt+1)*dx
WRITE(2, *) i, '    ', h_capital(i, nt+1), '    ', h(i, nt+1), &
&     '    ', b(i)
END DO
CLOSE(2)
END SUBROUTINE FORWARD_PROBLEM_D

SUBROUTINE FORWARD_PROBLEM(xx, v)
IMPLICIT NONE
REAL*8, PARAMETER :: rho=920.0
REAL*8, PARAMETER :: g=9.2
REAL*8, PARAMETER :: n=3
REAL*8, PARAMETER :: a=1e-16
REAL*8, PARAMETER :: dt=1/12.0
REAL*8, PARAMETER :: c=2*a/(n+2)*(rho*g)**n*1.e3**n
REAL*8, PARAMETER :: tend=5000
REAL*8, PARAMETER :: bx=-0.0001
REAL*8, PARAMETER :: m0=.004, m1=0.0002
INTRINSIC INT
INTEGER, PARAMETER :: nt=INT(tend/dt)
REAL*8, DIMENSION(nx+1, nt+1) :: h
REAL*8, DIMENSION(nx+1, nt+1) :: h_capital
INTEGER :: t, i
REAL*8, DIMENSION(nx+1) :: xarr
REAL*8, DIMENSION(nx+1) :: m
REAL*8, DIMENSION(nx+1) :: b
REAL*8, DIMENSION(nx) :: d, phi
REAL*8, DIMENSION(nx+1), INTENT(IN) :: xx
REAL*8, INTENT(OUT) :: v
INTRINSIC SIZE
REAL*8, DIMENSION(nx) :: pwx1
REAL*8 :: pwy1
REAL*8, DIMENSION(nx) :: pwr1
REAL*8, DIMENSION(nx) :: pwx2
REAL*8 :: pwy2
REAL*8, DIMENSION(nx) :: pwr2
xarr = (/((i-1)*dx, i=1,nx+1)/)
m = (/(m0-(i-1)*dx*m1, i=1,nx+1)/)
b = (/(1.0+bx*(i-1)*dx, i=1,nx+1)/)
m = m + xx
h(1, :) = b(1)
h(:, 1) = b
h(nx+1, :) = b(nx+1)
h_capital(1, :) = h(1, :) - b(1)
h_capital(nx+1, :) = h(nx+1, :) - b(nx+1)
h_capital(:, 1) = h(:, 1) - b
DO t=1,nt
pwx1 = (h_capital(1:nx, t)+h_capital(2:nx+1, t))/2
pwy1 = n + 2
pwr1 = pwx1**pwy1
pwx2 = (h(2:nx+1, t)-h(1:nx, t))/dx
pwy2 = n - 1
pwr2 = pwx2**pwy2
d(:) = c*pwr1*pwr2
phi(:) = -(d(:)*(h(2:nx+1, t)-h(1:nx, t))/dx)
h(2:nx, t+1) = h(2:nx, t) + m(2:nx)*dt - dt/dx*(phi(2:nx)-phi(1:nx&
&       -1))
WHERE (h(:, t+1) .LT. b) h(:, t+1) = b
h_capital(:, t+1) = h(:, t+1) - b
END DO
v = 0.
OPEN(unit=2, file='results_forward_run.txt', action='write', status=&
&  'replace')
WRITE(2, *) &
&   '         #                H                h                  b'
WRITE(2, *) '_______________________________________________________&
&________________________'
DO i=1,SIZE(h_capital(:, nt+1))
v = v + h_capital(i, nt+1)*dx
WRITE(2, *) i, '    ', h_capital(i, nt+1), '    ', h(i, nt+1), &
&     '    ', b(i)
END DO
CLOSE(2)
END SUBROUTINE FORWARD_PROBLEM

END MODULE FORWARD_TGT

• forward_d.msg - Contains warning messages. TAPENADE can be quite verbose with the warnings. None of the warnings below are much important.

1 Command: Procedure forward_problem understood as forward.forward_problem
2 forward_problem: (AD03) Varied variable h[i,nt+1] written by I-O to file 2
3 forward_problem: (AD03) Varied variable h_capital[i,nt+1] written by I-O to file 2
~


The adjoint model is described in detail in Tutorial 1. Another excellent resource is the MITgcm documentation.

It is pretty simple to generate the tangent linear model for our forward model.

% tapenade -reverse -head "forward_problem(V)/(xx)" forward.f90


-reverse tells Tapenade that we want a reverse i.e. adjoint code.

-head "forward_problem(V)/(xx) tells Tapenade that the head subroutine is forward_problem, the dependent variable or the objective function is V and the independent variable is xx.

This generates 2 files -

• forward_b.f90 - This file contains the adjoint module (called forward_diff) as well as the original forward code module. forward_diff contains the forward_problem_b subroutine which is our adjoint model. Notice how it requires values for xx, V which have to be obtained from a forward solve (calling the forward subroutine). We also will define Vb = 1 and propagate the adjoint sensitivities backwards to evaluate xxb. If you look at the second DO loop in forward_problem_b, it indeed runs backwards in time. Note that all the adjoint variables in the code are simply suffixed with the alphabet b. There are also a lot of PUSH and POP statements, they are essentially pushing and popping the stored variables from a stack just like we did manually in Tutorial 1.

!  TAPENADE 3.15 (master) - 15 Apr 2020 13:12
!
MODULE FORWARD_DIFF
IMPLICIT NONE
REAL*8, PARAMETER :: xend=30
REAL*8, PARAMETER :: dx=1
INTEGER, PARAMETER :: nx=INT(xend/dx)

CONTAINS
!  Differentiation of forward_problem in reverse (adjoint) mode:
!   gradient     of useful results: v
!   with respect to varying inputs: v xx
!   RW status of diff variables: v:in-zero xx:out
SUBROUTINE FORWARD_PROBLEM_B(xx, xxb, v, vb)
IMPLICIT NONE
REAL*8, PARAMETER :: rho=920.0
REAL*8, PARAMETER :: g=9.2
REAL*8, PARAMETER :: n=3
REAL*8, PARAMETER :: a=1e-16
REAL*8, PARAMETER :: dt=1/12.0
REAL*8, PARAMETER :: c=2*a/(n+2)*(rho*g)**n*1.e3**n
REAL*8, PARAMETER :: tend=5000
REAL*8, PARAMETER :: bx=-0.0001
REAL*8, PARAMETER :: m0=.004, m1=0.0002
INTRINSIC INT
INTEGER, PARAMETER :: nt=INT(tend/dt)
REAL*8, DIMENSION(nx+1, nt+1) :: h
REAL*8, DIMENSION(nx+1, nt+1) :: hb
REAL*8, DIMENSION(nx+1, nt+1) :: h_capital
REAL*8, DIMENSION(nx+1, nt+1) :: h_capitalb
INTEGER :: t, i
REAL*8, DIMENSION(nx+1) :: xarr
REAL*8, DIMENSION(nx+1) :: m
REAL*8, DIMENSION(nx+1) :: mb
REAL*8, DIMENSION(nx+1) :: b
REAL*8, DIMENSION(nx) :: d, phi
REAL*8, DIMENSION(nx) :: db, phib
REAL*8, DIMENSION(nx+1), INTENT(IN) :: xx
REAL*8, DIMENSION(nx+1) :: xxb
REAL*8 :: v
REAL*8 :: vb
INTRINSIC SIZE
REAL*8, DIMENSION(nx) :: temp
REAL*8, DIMENSION(nx) :: temp0
REAL*8, DIMENSION(nx) :: tempb
REAL*8, DIMENSION(nx) :: tempb0
REAL*8, DIMENSION(nx-1) :: tempb1
m = (/(m0-(i-1)*dx*m1, i=1,nx+1)/)
b = (/(1.0+bx*(i-1)*dx, i=1,nx+1)/)
m = m + xx
h(1, :) = b(1)
h(:, 1) = b
h(nx+1, :) = b(nx+1)
h_capital(1, :) = h(1, :) - b(1)
h_capital(nx+1, :) = h(nx+1, :) - b(nx+1)
h_capital(:, 1) = h(:, 1) - b
DO t=1,nt
CALL PUSHREAL8ARRAY(d, nx)
d(:) = c*((h_capital(1:nx, t)+h_capital(2:nx+1, t))/2)**(n+2)*((h(&
&       2:nx+1, t)-h(1:nx, t))/dx)**(n-1)
phi(:) = -(d(:)*(h(2:nx+1, t)-h(1:nx, t))/dx)
CALL PUSHREAL8ARRAY(h(2:nx, t+1), nx - 1)
h(2:nx, t+1) = h(2:nx, t) + m(2:nx)*dt - dt/dx*(phi(2:nx)-phi(1:nx&
&       -1))
mask(:) = h(:, t+1) .LT. b
CALL PUSHREAL8ARRAY(h(:, t+1), nx + 1)
WHERE (mask(:)) h(:, t+1) = b
CALL PUSHREAL8ARRAY(h_capital(:, t+1), nx + 1)
h_capital(:, t+1) = h(:, t+1) - b
END DO
DO i=1,SIZE(h_capital(:, nt+1))

END DO
CALL PUSHINTEGER4(i - 1)
h_capitalb = 0.0_8
h_capitalb(i, nt+1) = h_capitalb(i, nt+1) + dx*vb
END DO
hb = 0.0_8
mb = 0.0_8
DO t=nt,1,-1
CALL POPREAL8ARRAY(h_capital(:, t+1), nx + 1)
hb(:, t+1) = hb(:, t+1) + h_capitalb(:, t+1)
h_capitalb(:, t+1) = 0.0_8
CALL POPREAL8ARRAY(h(:, t+1), nx + 1)
phib = 0.0_8
CALL POPREAL8ARRAY(h(2:nx, t+1), nx - 1)
db = 0.0_8
temp = (h(2:nx+1, t)-h(1:nx, t))/dx
temp0 = (h_capital(1:nx, t)+h_capital(2:nx+1, t))/2
WHERE (mask(:)) hb(:, t+1) = 0.0_8
hb(2:nx, t) = hb(2:nx, t) + hb(2:nx, t+1)
mb(2:nx) = mb(2:nx) + dt*hb(2:nx, t+1)
tempb1 = -(dt*hb(2:nx, t+1)/dx)
hb(2:nx, t+1) = 0.0_8
phib(2:nx) = phib(2:nx) + tempb1
phib(1:nx-1) = phib(1:nx-1) - tempb1
db(:) = -((h(2:nx+1, t)-h(1:nx, t))*phib(:)/dx)
tempb = -(d(:)*phib(:)/dx)
hb(2:nx+1, t) = hb(2:nx+1, t) + tempb
hb(1:nx, t) = hb(1:nx, t) - tempb
CALL POPREAL8ARRAY(d, nx)
WHERE (temp0 .LE. 0.0 .AND. (n + 2 .EQ. 0.0 .OR. n + 2 .NE. INT(n &
&         + 2)))
tempb = 0.0_8
ELSEWHERE
tempb = (n+2)*temp0**(n+1)*temp**(n-1)*c*db(:)/2
END WHERE
WHERE (temp .LE. 0.0 .AND. (n - 1 .EQ. 0.0 .OR. n - 1 .NE. INT(n -&
&         1)))
tempb0 = 0.0_8
ELSEWHERE
tempb0 = (n-1)*temp**(n-2)*temp0**(n+2)*c*db(:)/dx
END WHERE
hb(2:nx+1, t) = hb(2:nx+1, t) + tempb0
hb(1:nx, t) = hb(1:nx, t) - tempb0
h_capitalb(1:nx, t) = h_capitalb(1:nx, t) + tempb
h_capitalb(2:nx+1, t) = h_capitalb(2:nx+1, t) + tempb
END DO
xxb = 0.0_8
xxb = mb
vb = 0.0_8
END SUBROUTINE FORWARD_PROBLEM_B

SUBROUTINE FORWARD_PROBLEM(xx, v)
IMPLICIT NONE
REAL*8, PARAMETER :: rho=920.0
REAL*8, PARAMETER :: g=9.2
REAL*8, PARAMETER :: n=3
REAL*8, PARAMETER :: a=1e-16
REAL*8, PARAMETER :: dt=1/12.0
REAL*8, PARAMETER :: c=2*a/(n+2)*(rho*g)**n*1.e3**n
REAL*8, PARAMETER :: tend=5000
REAL*8, PARAMETER :: bx=-0.0001
REAL*8, PARAMETER :: m0=.004, m1=0.0002
INTRINSIC INT
INTEGER, PARAMETER :: nt=INT(tend/dt)
REAL*8, DIMENSION(nx+1, nt+1) :: h
REAL*8, DIMENSION(nx+1, nt+1) :: h_capital
INTEGER :: t, i
REAL*8, DIMENSION(nx+1) :: xarr
REAL*8, DIMENSION(nx+1) :: m
REAL*8, DIMENSION(nx+1) :: b
REAL*8, DIMENSION(nx) :: d, phi
REAL*8, DIMENSION(nx+1), INTENT(IN) :: xx
REAL*8, INTENT(OUT) :: v
INTRINSIC SIZE
xarr = (/((i-1)*dx, i=1,nx+1)/)
m = (/(m0-(i-1)*dx*m1, i=1,nx+1)/)
b = (/(1.0+bx*(i-1)*dx, i=1,nx+1)/)
m = m + xx
h(1, :) = b(1)
h(:, 1) = b
h(nx+1, :) = b(nx+1)
h_capital(1, :) = h(1, :) - b(1)
h_capital(nx+1, :) = h(nx+1, :) - b(nx+1)
h_capital(:, 1) = h(:, 1) - b
DO t=1,nt
d(:) = c*((h_capital(1:nx, t)+h_capital(2:nx+1, t))/2)**(n+2)*((h(&
&       2:nx+1, t)-h(1:nx, t))/dx)**(n-1)
phi(:) = -(d(:)*(h(2:nx+1, t)-h(1:nx, t))/dx)
h(2:nx, t+1) = h(2:nx, t) + m(2:nx)*dt - dt/dx*(phi(2:nx)-phi(1:nx&
&       -1))
mask(:) = h(:, t+1) .LT. b
WHERE (mask(:)) h(:, t+1) = b
h_capital(:, t+1) = h(:, t+1) - b
END DO
v = 0.
OPEN(unit=2, file='results_forward_run.txt', action='write', status=&
&  'replace')
WRITE(2, *) &
&   '         #                H                h                  b'
WRITE(2, *) '_______________________________________________________&
&________________________'
DO i=1,SIZE(h_capital(:, nt+1))
v = v + h_capital(i, nt+1)*dx
WRITE(2, *) i, '    ', h_capital(i, nt+1), '    ', h(i, nt+1), &
&     '    ', b(i)
END DO
CLOSE(2)
END SUBROUTINE FORWARD_PROBLEM

END MODULE FORWARD_DIFF

• forward_b.msg - Contains warning messages. Tapenade can be quite verbose with the warnings. None of the warnings below are much important.

1 Command: Procedure forward_problem understood as forward.forward_problem
2 forward_problem: Command: Input variable(s) v have their derivative modified in forward_problem: added to independents
3 forward_problem: (AD03) Varied variable h[i,nt+1] written by I-O to file 2
4 forward_problem: (AD03) Varied variable h_capital[i,nt+1] written by I-O to file 2


## 7.6.2.2.6. Driver Routine¶

The driver routine evaluates the sensitivities of our cost function V to the independent variable xx in three ways - Finite Differences, Tangent Linear Model (forward mode), Adjoint Model (reverse mode). Both the TLM and adjoint methods should give sensitivities that agree within some tolerance with the Finite Differences sensitivities. One can see that the TLM and finite differences have to be called multiple times (as many times as the number of control variables) while the adjoint model only has to be called once to evaluate the entire gradient, making it much more efficient.

program driver
!!Essentially both modules had variables with the same name
!!so they had to be locally aliased for the code to compile and run
!!Obviously only one set of variables is needed so the other is useless

use forward_diff, n => nx, fp => forward_problem !! Alias to a local variable
use forward_tgt, n_useless => nx, fp_useless => forward_problem

implicit none

real(8), dimension(n+1) :: xx=0., xx_tlm =0., xxb =0., xx_fd, accuracy
real(8) :: V=0., Vb = 1., V_forward, Vd = 0.
real(8), parameter :: eps = 1.d-7
integer :: ii

call fp(xx,V)

!! Forward run
V_forward = V

xx = 0.
V = 0.
call fp(xx,V)
call forward_problem_b(xx,xxb,V,Vb)

open (unit = 1, file = "results.txt", action="write",status="replace")

write(1,*) "         #                Reverse                           FD",&
"                          Tangent                     Relative accuracy"
write(1,*) "______________________________________________________________",&
"_______________________________________________________________________"

!! Finite differences and Tangent Linear Model
do ii = 1, n+1

xx = 0.
V = 0.
xx_tlm = 0.
xx_tlm(ii) = 1.

!! TLM
call fp(xx,V)
call forward_problem_d(xx,xx_tlm,V,Vd)

!! FD
xx = 0.
V = 0.
xx(ii) = eps
call fp(xx,V)
xx_fd(ii) =  (V - V_forward)/eps

if ( xx_fd(ii).NE. 0. ) then
accuracy(ii) = 1.d0 - Vd/xx_fd(ii)
else
accuracy(ii) = 0.
end if
write(1,*) ii, "    ", xxb(ii), "    ", xx_fd(ii),"    ", Vd,"    ", accuracy(ii)
end do

close(1)

call execute_command_line('gnuplot plot.script')
end program driver


## 7.6.2.2.7. Combining all compilation commands into a Makefile¶

Makefiles are extremely useful to automate the compilation process. The Makefile shown below generates both the tangent linear and adjoint codes, and bundles everything together with the driver routine and a special file provided by the Tapenade developers, called adStack.c together into a single executable named adjoint (unfortunately, a slight misnomer since it also executes FD and TLM codes). adStack.c contains the definitions of the POP, PUSH mechanisms that Tapenade uses in its generated codes.

The adStack.c file can be found in test_ad/tapenade_supported/ADFirstAidKit sub-directory in the ad branch of the SICOPOLIS repository. The Makefile assumes that you have the entire ADFirstAidKit directory present, so it is better to copy the entire directory to work as is with this Makefile.

More information on Makefiles can be found here . Makefiles are extremely sensitive to whitespaces both at the start and end of any line, so you should be very careful with that.

• Makefile:

EXEC := adjoint
SRC  := $(wildcard *.f90) OBJ :=$(patsubst %.f90, %.o, $(SRC)) # NOTE - OBJ will not have the object files of c codes in it, this needs to be improved upon. # Options F90 := gfortran CC := gcc TAP_AD_kit:= ./ADFirstAidKit # Rules$(EXEC): $(OBJ) adStack.o forward_b.o forward_d.o$(F90) -o $@$^

%.o: %.f90
$(F90) -c$<

driver.o: forward_b.f90 forward_diff.mod forward_tgt.mod
forward_diff.mod: forward_b.o
forward_tgt.mod: forward_d.o

$(CC) -c$(TAP_AD_kit)/adStack.c

forward_b.f90: forward.f90
forward_d.f90: forward.f90
# Useful phony targets

.PHONY: clean

clean:
$(RM)$(EXEC) *.o $(MOD)$(MSG) *.msg *.mod *_db.f90 *_b.f90 *_d.f90 *~


To compile everything, just type -

% make


To clean up the compilation, just type -

% make clean


## 7.6.2.2.8. Results¶

The results are shown here. The adjoint and the tangent linear model agree with each other up to 12 decimal places in most cases. The agreement of both of these with the Finite differences is also extremely good, as seen from the Relative accuracy.

         #                Reverse                           FD                          Tangent                     Relative accuracy
_____________________________________________________________________________________________________________________________________
1        0.0000000000000000             0.0000000000000000             0.0000000000000000             0.0000000000000000
2        7.7234861280779530             7.7234756545863092             7.7234861280783651            -1.3560594380734869E-006
3        13.803334061615192             13.803310956461701             13.803334061616136            -1.6738849475395057E-006
4        19.668654270913329             19.668613369105969             19.668654270914736            -2.0795471444845504E-006
5        25.460720864192346             25.460657031572964             25.460720864194293            -2.5071081728444966E-006
6        31.276957320291768             31.276862468843092             31.276957320293935            -3.0326395730195799E-006
7        37.209906557896431             37.209769434554119             37.209906557899728            -3.6851436515661362E-006
8        43.370516019133923             43.370320437219334             43.370516019138648            -4.5095797620575695E-006
9        49.919848601374966             49.919570930256896             49.919848601380266            -5.5623699923845749E-006
10        57.138817544000176             57.138414888413536             57.138817544006628            -7.0470207107486971E-006
11        65.642858080405333             65.642249307273914             65.642858080414797            -9.2741054322775796E-006
12        77.452386462370669             77.451334217215617             77.452386462382293            -1.3585888187783723E-005
13        139.25939689848738             139.25460375929788             139.25939689851842            -3.4419969546117812E-005
14        149.07659819564697             149.07253577334245             149.07659819568059            -2.7251313040821401E-005
15        153.96500702363139             153.96133390410682             153.96500702366376            -2.3857415779593438E-005
16        156.68784327845535             156.68444133254411             156.68784327848971            -2.1712085237490797E-005
17        158.00297057747002             157.99977889585648             158.00297057750697            -2.0200545043813634E-005
18        158.20831395372016             158.20529343457679             158.20831395375782            -1.9092402760101379E-005
19        157.42238310624799             157.41950864622822             157.42238310628576            -1.8259871868764321E-005
20        155.66584913271916             155.66310167969277             155.66584913275506            -1.7649995616375591E-005
21        152.88684508060106             152.88421169046273             152.88684508064247            -1.7224735965992721E-005
22        148.96151788159625             148.95898933886542             148.96151788163439            -1.6974757818921660E-005
23        143.67549543577167             143.67306516049894             143.67549543580796            -1.6915316077614762E-005
24        136.67934044641444             136.67700816455408             136.67934044645017            -1.7064186050186336E-005
25        127.39349079030032             127.39126539429435             127.39349079032215            -1.7468984399471310E-005
26        114.79185160452661             114.78974569101297             114.79185160454820            -1.8345833267208178E-005
27        96.827067402857665             96.825136655098731             96.827067402878811            -1.9940563440234982E-005
28        68.425966236870423             68.424387524856911             68.425966236889820            -2.3072358993792008E-005
29        0.0000000000000000             0.0000000000000000             0.0000000000000000             0.0000000000000000
30        0.0000000000000000             0.0000000000000000             0.0000000000000000             0.0000000000000000
31        0.0000000000000000             0.0000000000000000             0.0000000000000000             0.0000000000000000


A plot of these results is also shown here. To the naked eye, these results are identical.

## 7.6.2.2.9. A note on Binomial Checkpointing¶

### How to use Binomial Checkpointing with Tapenade?¶

Using binomial checkpointing is extremely simple with Tapenade. First we add the following line just before our time loop starts.

!\$AD BINOMIAL-CKP nt+1 20 1
do t = 1,nt
!!!! Stuff inside the time loop.
end do


Notice that for the fortran compiler the line we just added is a comment, but when Tapenade parses the code, it knows it needs to use binomial checkpointing. The line essentially tells Tapenade the following - “Use binomial checkpointing for this loop which has nt+1 steps (if you count the exit iteration as well), checkpoint every 20 steps, with the first step numbered as 1.” Tapenade uses some subroutines whose names start with ADBINOMIAL to manage the checkpointing and thus we need to link to the adBinomial.c file in the Makefile so that we tell the compiler what these subroutines are. (essentially the same thing we did with adStack.c)