arXiv:0711.4444v2 [cs.MS] 29 Nov 2007

INSTITUT NATIONAL DE RECHERCHE EN INFORMATIQUE ET EN AUTOMATIQUE

Building the Tangent and Adjoint codes of the Ocean
General Circulation Model OPA with the Automatic
Differentiation tool TAPENADE
M.H. Tber — L. Hascoët — A. Vidard — B. Dauvergne

N° 6372
Novembre 2007
Thème NUM

apport
de recherche

Building the Tangent and Adjoint odes of the
O ean General Cir ulation Model OPA with the
Automati Dierentiation tool TAPENADE
∗

M.H. Tber

, L. Has oët

∗

†

, A. Vidard

, B. Dauvergne

∗

Thème NUM  Systèmes numériques
Projet Tropi s
Rapport de re her he

Abstra t:

n° 6372  Novembre 2007  28 pages

The o ean general

ir ulation model OPA is developed by the

LODYC team at Paris VI university.

OPA has re ently undergone a major

rewriting, migrating to FORTRAN95, and its adjoint

ode needs to be re-

built. For earlier versions, the adjoint of OPA was written by hand at a high
development

ost.

We use the Automati

Dierentiation tool TAPENADE

to build me hani aly the tangent and adjoint
dierentiated

odes by

odes of OPA. We validate the

omparison with divided dieren es, and also with an

identi al twin experiment. We apply state-of-the-art methods to improve the
performan e of the adjoint
and Walther's binomial

ode.

In parti ular we implement the Griewank

he kpointing algorithm whi h gives us an optimal

trade-o between time and memory

onsumption. We apply a spe i

egy to dierentiate the iterative linear solver that

strat-

omes from the impli it time

stepping s heme.

Key-words:

OPA, general

ir ulation model, TAPENADE, Automati

ferentiation, reverse mode, Adjoint Code, Che kpointing

∗
†

INRIA Sophia Antipolis, Fran e
INRIA Grenoble, Fran e

Unité de recherche INRIA Sophia Antipolis

Dif-

Constru tion des odes Tangent et Adjoint du
modèle de ir ulation générale o éanique OPA
par l'outil de Diérentiation Automatique
TAPENADE
Résumé : Le modèle de

ir ulation générale o éanique OPA est développé

par l'équipe LODYC de l'université Paris VI. La nouvelle version 9 d'OPA
onstitue une évolution majeure, ave
TRAN95.

Les

en parti ulier une migration vers FOR-

odes Linéaire Tangent et Adjoint d'OPA, qui auparavant

étaient é rits à la main, doivent don

être redéveloppés. Nous utilisons l'outil

de Diérentiation Automatique TAPENADE pour
gent et Adjoint d'OPA 9.
raison ave

onstruire les

odes Tan-

Nous validons les dérivées obtenues par

ompa-

les Diéren es Divisées et sur deux appli ations test in luant des

expérien es jumelles. Nous utilisons le s héma de

he kpointing ré ursif bino-

mial de Griewank et Walther pour améliorer les performan es du

ode adjoint.

Nous utilisons une stratégie spé ique pour diérentier le solveur linéaire itératif provenant du s héma impli ite d'avan ement en temps.
montrent un

oût raisonnable, tant en

Nos résultats

onsommation mémoire que pour le

temps d'exé ution de l'adjoint.

Mots- lés :

OPA, Cir ulation O éanique, TAPENADE, Diérentiation Au-

tomatique, mode inverse, Code Adjoint, Che kpointing

Tangent and Adjoint dierentiation of OPA

3

1 Introdu tion
The development of tangent and adjoint models is an important step in addressing sensitivity analysis and variational data assimilation problems in O eanography. Sensitivity analysis is the study of how model output varies with
hanges in model inputs.

The sensitivity information given by the adjoint

model is used dire tly to gain an understanding of the physi al pro esses. In
data assimilation, one

onsiders a

ost fun tion whi h is a measure of the

model-data mist. The adjoint sensitivities are used to build the gradient for
des ent algorithms. Similarly the tangent model is used in the
in remental algorithms [3℄ to linearize the

ontext of the

ost fun tion around a ba kground

ontrol. For the previous version 8 of the O ean General Cir ulation Model
OPA [17℄, Weaver et al [22℄ developed the numeri al tangent and adjoint
by hand using

odes

lassi al te hniques [5, 19℄. Sin e then, the OPA model has un-

dergone a major update. Parti ularly the new versions are fully rewritten in
FORTRAN95.
adjoint

In this paper, we report on the development of tangent and

odes of OPA using the Automati

Dierentiation (AD) tool TAPE-

NADE [12℄. A brief des ription of the OPA model and the

onguration used

in this work is given in the next se tion. In se tion 3 we present the prin iples of AD and how they are ree ted into the fun tionalities of the AD tool
TAPENADE. In se tion 4 we fo us on the most interesting di ulties that we
en ountered in the appli ation of AD to su h a large

ode. Se tion 5 shows

some experiments that validate our derivatives and presents two illustrative
appli ations, fo using on

omputational aspe ts rather than impli ations for

o eanography. An outlook of further work is given in the

on lusion.

2 The O ean General Cir ulation Model OPA
Developed by the LODYC team at Paris VI university, OPA is a exible o ean
ir ulation model that
onguration.

an be used either in a regional or in a global o ean

OPA is the o ean model

omponent of NEMO (Nu leus For

European Modelling of the O ean) and is widely used in the s ienti

om-

munity. Moreover it is be oming a major a tor in operational o eanography
(Mer ator, ECMWF, UK-Met o e) Its formulation is based on the so- alled
primitive equations for the temporal evolution of o ean velo ity

RR n° 6372

urrents, tem-

Tber et al.

4

perature and salinity in its three horizontal and verti al dimensions.
equations are derived from Navier-Stokes equations

These

oupled with a state equa-

tion for water density and heat equation, under Boussinesq and hydrostati
approximations.

U

Let us introdu e the following variables:

Uh + wk

(the subs ript

tial temperature,

S

h

the velo ity ve tor,

denotes the lo al horizontal ve tor),

the salinity,

p

the pressure and

ρ

T

U =

the poten-

the in-situ density. The

ve tor invariant form of the primitive equations in an orthogonal set of unit
ve tors linked to the earth are written as follows

































∂Uh
∂t
∂p
∂z
∇·U
∂T
∂t
∂S
∂t
ρ



1
= − (∇ × U) × U + ∇ (U) 2
2



h

− f k × Uh −

1
∇h p + D U
ρ0

= −ρg
= 0
= −∇ · (T U) + D T
= −∇ · (SU) + D S
= ρ (T, S, p)

∇ is the generalized derivative ve tor operator, t the time, z the verti al
oordinate, ρ0 a referen e density, f the Coriolis a eleration, and g the gravity
U
T
S
a eleration. D , D and D are the parametrization of small s ale physi s for

where

momentum, temperature and salinity, in luding surfa e for ing terms. A full
des ription of the model basi s, dis retization, physi al and numeri al details
an be found in [17℄.
Through this paper, OPA is used in its global free surfa e
ORCA-2. In this

onguration

onguration the model uses a rotated grid with poles on

North Ameri a and Asia in order to avoid the singularity problem on the North
Pole. The spa e resolution is roughly equivalent to a geographi al mesh of 2°

by 1.3° with a meridional resolution of 0.5° near the Equator (see gure 1). The

Verti al domain, spreading from the surfa e to a depth of 5000m, is meshed
using 31 levels with levels 1 to 10 in the top 100 meters. The time step is 96
minutes so that there are 15 time steps per day. The model is for ed by heat,
freshwater, and momentum uxes from the atmosphere and/or the sea-i e.

INRIA

Tangent and Adjoint dierentiation of OPA

5

Figure 1: ORCA 2 Mesh

The solar radiation penetrates the upper layers of the o ean. Zero uxes of
heat and salt are applied through the bottom. On the lateral solid boundaries
a no-slip

ondition is also applied. Initialization of the model for temperature

and salinity is based on the Levitus et al. (1998)

limatology with a null initial

velo ity eld. For more details about the spa e time-domain and the o ean
physi s of ORCA-2, we refer to the page dedi ated to this

1

onguration in the

o ial website of NEMO-OPA .
The
to

onguration ORCA-2 is routinely used by MERCATOR/Meteo-Fran e

ompute the o eani

omponent of their seasonal fore asting system. The

size of OPA-9, 200 modules dening 800 pro edures with over 100 000 lines of
FORTRAN95, makes it the largest appli ation dierentiated by TAPENADE
to date.

The

omputational kernel whi h is a tually dierentiated a

for 330 pro edures.

1 http://www.lody

RR n° 6372

.jussieu.fr/NEMO/general/des ription/ORCA_ ong.html

ounts

Tber et al.

6

3 Prin iples of AD and the tool TAPENADE
2

TAPENADE [12℄ is an AD tool developed by the Tropi s

team at INRIA.

Given the sour e of an original program that evaluates a mathemati al fun tion, and given a sele tion of input and output variables to be dierentiated,
TAPENADE produ es a new sour e program that

omputes the partial deriva-

tives of the sele ted outputs with respe t to the sele ted inputs.
Basi ally, TAPENADE does that by inserting additional statements into a
opy of the original program. Like other AD tools, TAPENADE is based on
the fundamental observation that the original program P, whatever its size and
m
n
run time, omputes a fun tion F, X∈IR 7→ Y ∈IR whi h is the omposition
of the elementary fun tions
words if

P

omputed by ea h run-time instru tion. In other

exe utes a sequen e of elementary statements

Ik , k ∈ [1..p],

then

P

a tually evaluates

F = fp ◦ fp−1 ◦ · · · ◦ f1 ,
where ea h

fk

is the fun tion implemented by

hain rule of derivative
derivatives of ea h
Calling

X0 = X

Y
Xk = fk (Xk−1 )

omponent of

and

variables, i.e. the su

Ik .

Therefore one an apply the
F ′ , i.e. the partial

al ulus to get the Ja obian matrix

essive

with respe t to ea h
the su

omponent of

X.

essive values of all intermediate

states of the memory throughout exe

ution of

P,

we get

′
F ′ (X) = fp′ (Xp−1 ) × fp−1
(Xp−2 ) × · · · × f1′ (X0 ) .
The derivatives

fk′

(1)

of ea h elementary instru tion are easily built, and must

be inserted in the dierentiated program so that ea h of them has the values

Xk−1

dire tly available for use. This pro ess yields analyti

are exa t up to numeri al a

derivatives, that

ura y.

In pra ti e, two sorts of derivatives are of parti ular importan e in s ienti

omputing: the tangent (or dire tional) derivatives, and the adjoint (or

reverse) derivatives.

In parti ular, tangent and adjoint are the two sorts of

derivative programs required for OPA, and TAPENADE provides both. The
′
tangent derivative is the produ t Ẏ = F (X) × Ẋ of the full Ja obian times a
dire tion

Ẋ

in the input spa e. From equation (1), we nd

′
Ẏ = F ′ (X) × Ẋ = fp′ (Xp−1) × fp−1
(Xp−2 ) × · · · × f1′ (X0 ) × Ẋ
2

(2)

http://www-sop.inria.fr/tropi s/

INRIA

Tangent and Adjoint dierentiation of OPA
whi h is most
u ts are mu h

7

heaply exe uted from right to left be ause matrix×ve tor prodheaper than matrix×matrix produ ts. This is also the most

onvenient exe ution order be ause it uses the intermediate values

P builds them.
X = F ′∗ (X) × Y of

Xk

in the

same order as the program

On the other hand the adjoint

derivative is the produ t

the

transposed Ja

Y in the output spa e. The resulting X
(Y · Y ). From equation (1), we nd

weight ve tor
dot produ t

obian times a

is the gradient of the

′∗
X = F ′∗ (X) × Y = f1′∗ (X0 ) × · · · × fp−1
(Xp−2 ) × fp′∗ (Xp−1 ) × Y

(3)

whi h is also most

heaply exe uted from right to left. However, this uses the

intermediate values

Xk

in the inverse of their building order in

Regarding the runtime

Ẏ

and adjoint

X

ost only a small multiple of the original program

slowdown fa tor is less than 4 in theory. In pra ti e it
the tangent, whereas it

below. Despite its higher

m runs

m

The

an be less than 2 for

ost, the adjoint

ode is still by large the

heapest

To get the gradient with the tangent mode would

of the tangent

is independent from

P.

an rea h up to 10 for the adjoint for a reason dis ussed

way to obtain gradients.
require

P.

ost for obtaining the derivatives, both tangent

ode, one per dimension of

X,

whereas this

ost

with the adjoint mode.

The di ulty of the adjoint mode lies in the fa t that it needs the intermediate values

Xk

in reverse order. To this end, TAPENADE basi ally uses a

two-sweeps strategy,
a

alled Store-All. In the rst sweep (the forward sweep),

opy of the original program

P

is run, together with Push statements that

store intermediate values on a sta k just before they get overwritten. In the
se ond sweep (the ba kward sweep), the derivative statements ompute the
′∗
elementary derivatives fk (Xk−1 ) for k = p down to 1, using Pop statements
to restore the intermediate values as they are required. This in urs a ost in
memory spa e as the maximum sta k size needed is attained at the end of
the forward sweep, and is thus proportional to the length of the program

P.

There is also a runtime penalty for these sta k manipulations. TAPENADE
implements a number of strategies [11℄ to mitigate this
data-ow analysis of the program's
values

Xk that need to be stored.

ontrol ow graph, redu ing the number of

However for very long programs su h as OPA,

involving unsteady simulations, Store-All

an not work alone.

ombines it with a storage/re omputation trade-o

RR n° 6372

ost, based on stati

TAPENADE

alled  he kpointing.

Tber et al.

8

Che kpointing redu es the maximum sta k size at the
exe utions. Consider a pie e

C

of the original program

P.

ost of dupli ated

Che kpointing

illustrated on gure 2 means that during the main forward sweep,

C

C

as

pushes

no value on the sta k. When the ba kward sweep rea hes ba k to the pla e
where intermediate values are now missing on the sta k, it runs

C

a se ond

time, this time with the Store-All strategy i.e. pushing values on the sta k.
The ba kward sweep

an then resume safely.

To run

C

twi e requires that

enough of its input values, a snapshot, are stored but the size of a snapshot
is generally mu h less than the sta k size used by
down the adjoint program. When

C

is well

C.

hosen,

Obviously, this also slows
he kpointing

the peak size of the sta k by a fa tor of two. Che kpoints
whi h

an divide

an be nested, in

ase both the sta k's peak size and the adjoint runtime slowdown

grow as little as the logarithm of the size of
applies

he kpointing to ea h pro edure

an

P. In its default mode, TAPENADE

all.

{

C

successive
sweeps

Figure 2: Che kpointing applied to the program pie e

C.

Rightwards arrows

represent forward sweeps, thi k when they store intermediate values on the
sta k, thin otherwise.

Leftwards arrows represent ba kward sweeps.

Bla k

dots are stores, white dots are retrieves. Small dots are Push and Pops, big
dots are snapshots.

TAPENADE

apa ity to generate robust and e ient tangent and adjoint

odes has been demonstrated on several real-world test appli ations [15, 7, 1,
13, 16℄. Regarding the appli ation language, it
in FORTRAN. Taking into a

an handle programs written

ount the new programming

onstru ts provided

by FORTRAN95 has required an important programming eort in the past
few years, mostly to handle modules, stru tured data types, array notation,
pointers, and dynami

memory allo ation. Sin e the new OPA 9 is now written

in FORTRAN95, dierentiation of OPA is a very realisti

test for the new

TAPENADE 2.2.

INRIA

Tangent and Adjoint dierentiation of OPA
There exist several other AD tools.

9

Restri ting to the tools whi h, like

TAPENADE, operate by sour e transformation, provide tangent and adjoint
modes, use global program analysis to optimize the dierentiated
have demonstrated their appli ability on large industrial

odes, we

ode, and
an men-

tion TAF [4℄ a pioneer of AD for meteorology, now the standard AD tool for
the popular MIT Global Cir ulation Model. Unlike TAPENADE's, the adjoint
mode of TAF regenerates the intermediate values
an given initial point. This is

Xk

with Store-All strategy is getting blurred by nested
joint
su

by re omputation from

alled a Re ompute-All strategy. Comparison

odes grow more alike as more

he kpointing, as the ad-

he kpoints are inserted. OpenAD [20℄,

essor of ADIFOR and ADIC, uses the Store-All strategy. There are ex-

periments to also apply OpenAD to the MIT GCM. The tool Adol-C [10℄,
although using operator overloading instead of sour e transformation, is very
popular and has been applied su
adjoint mode

essfully to many industrial appli ations. Its

an be seen as an extension of the Store-All strategy: not only

the intermediate values are stored on the sta k, but also the

omputation graph

to be dierentiated. This allows the AD tool to perform further optimizations
on this graph, at the

ost of a higher memory

onsumption.

4 Applying TAPENADE to OPA
We generated working tangent and adjoint

odes for the

omputational kernel

of OPA, using TAPENADE. Depending on the nal appli ation (

f

se tion

5), the a tual fun tion to dierentiate as well as the input and output variables may be dierent, but the te hni al di ulties that we en ountered are
essentially the same. This se tion des ribes these points.

4.1

FORTRAN95

onstru ts

The new OPA 9 uses extensively the modular
had to extend the

all-graph internal representation of TAPENADE to handle

the nesting of modules and pro edures.
into the dierentiated
Be ause a module

Essentially this nesting is mirrored

ode.
an dene private

ferentiated modules do not have a

RR n° 6372

onstru ts of FORTRAN95. We

omponents, subroutines in the dif-

ess to all variables of the original module.

Tber et al.

10

Therefore the dierentiated module must ontain its own opy of all the original
module's variables, types, and pro edures. This is a
dierentiation model: the dierentiated
the original

ode; it must

dierentiated

ode

ontain its own

hange in TAPENADE's

annot just

all or use parts of

opies of those. In other words, the

ode need not be linked with the original.

The interfa e me hanism of FORTRAN95 is a way to implement overloaded pro edures.

This is stati

overloading, whi h is resolved at

ompile

time. Therefore we had to extend TAPENADE type- he king phase to
pletely solve the

alls to interfa ed pro edures.

om-

Conversely, TAPENADE is

now able to generate interfa es on the dierentiated pro edures, so that the
general stru ture of the

ode is preserved.

The array notation of FORTRAN95 is used systemati ally in OPA. At the
same time, dierentiation requires that many

alls to intrinsi

fun tions be

split to propagate the derivatives. When these fun tions are used on arrays
("elemental" intrinsi s) TAPENADE must generate a

ode whi h is far from

trivial. For instan e the single statement from OPA:

zws(:,:,:) = SQRT(ABS(psal(:,:,:)))
generates in the adjoint mode

abs1 = ABS(psal(:,:,:))
mask = (psal(:,:,:) .GT. 0.0)
...
WHERE (abs1 .EQ. 0.0)
abs1b = 0.0
ELSEWHERE
abs1b = zwsb(:, :, :)/(2.0*SQRT(abs1))
END WHERE
WHERE (.NOT.mask(:, :, :))
psalb(:, :, :) = psalb(:, :, :) - abs1b
ELSEWHERE
psalb(:, :, :) = psalb(:, :, :) + abs1b
END WHERE
Without going in too mu h detail into the adjoint dierentiation model, we
observe that the test that is needed to prote t the dierentiated

ode against

INRIA

Tangent and Adjoint dierentiation of OPA
the non-dierentiability of
dierentiation of

ABS,

SQRT

11

at 0, as well as the test that

have been turned into

WHERE

ontrols the

onstru ts to keep the

runtime benets of array notation. Some temporary variables are introdu ed
automati ally to store

ontrol-ow de isions (e.g.

abs1

and

mask),

although

TAPENADE still doesn't do this in an optimal way on the example.
OPA uses pointers and dynami

DEALLOCATE).

memory allo ation ( alls to

ALLOCATE

and

This is an appli ation for the pointer analysis now available

in TAPENADE, nding whether a variable has a derivative, even when this
variable is a

essed through a pointer. Unfortunately, dynami

allo ation is

handled partly, i.e. only in the tangent mode of TAPENADE. In the adjoint
mode, we have no general strategy for memory allo ation and TAPENADE
sometimes

annot produ e a working

an allo ate should be a

ode. We understand that the adjoint of

DEALLOCATE,

and vi e-versa, but some

be made by hand on the dierentiated

4.2

hanges must

ode to make it work.

Che kpointing and hidden variables

OPA reads and writes several data les, not only during the pre- and postpro essing stages, but also during the

omputational kernel itself. Sour e terms

su h as the wind stress are being read at intermediate time steps. Also, some
modules and pro edures dene private
but

annot be a

SAVE variables, whose value is preserved

essed from outside.

are just examples of a

Although unrelated, these two points

ommon problem: they

reentrant".
If a

alled pro edure modies an internal

sible from the outside

alling

ontext to

an identi al result. Similarly if the

SAVE

an make a pro edure non
variable, it be omes impos-

all the pro edure a se ond time with

alled pro edure reads from a previously

opened le, and just moves the read pointer further in the le, then it be omes
impossible to

all the pro edure twi e and obtain the same values read.

Non reentrant pro edures are a problem for the
the adjoint mode. We saw in se tion 3 that

he kpointing strategy of

he kpointing relies on

alling the

he kpointed pie e twi e, in su h a way that the se ond

all is equivalent to the

rst. To this end, a su ient subset of the exe ution

ontext, the snapshot,

must be saved and restored. Hidden variables like an internal
the read pointer inside an opened le

RR n° 6372

SAVE variable or

annot be saved nor restored in general.

Tber et al.

12

When
then

he kpointing would require hidden variables to be put in the snapshot,
he kpointing should be forbidden.

Similarly, when a pro edure only allo ates some memory, the allo ation
must not be done twi e.

If this pro edure is

he kpointed, then one must

deallo ate the memory when restoring the snapshot before the dupli ate

all.

TAPENADE is not yet able to do this automati ally.
TAPENADE has some fun tionalities to
problem, but in all

ope with this hidden variables

ases intera tion with the user is ne essary. First, TAPE-

NADE issues a warning message when a subroutine
be ause of a private

SAVE

annot be

he kpointed

variable. The message is issued only when this vari-

able would be part of the snapshot for this pro edure. When this happened
for OPA, we just turned by hand the variables in question into publi
variables in the original

ode. In prin iple this

global

ould also be done automati-

ally. However there are only a handful su h variables, thus developing this is
not our priority.
When a subroutine is not reentrant be ause of I/O le pointers or be ause
of isolated memory allo ation or deallo ation, then TAPENADE lets the user
label the subroutine so that it must not be

he kpointed. For OPA, we took

another strategy: we modied the main I/O subroutines so that they always
rst make sure that the le is opened and then only use dire t read into the
le without using a read I/O pointer. Thus all I/O subroutines are reentrant.

4.3

Binomial Che kpointing

Automati

Dierentiation of OPA is one of the most ambitious appli ations

of TAPENADE so far. It means building the adjoint of a pie e of

ode that

performs an unsteady nonlinear simulation over a very large number of time
steps. Ea h time step

omputes a new state whose size ranges in the hundreds

of megabytes. In adjoint mode if no

he kpointing was applied, whi h means

that all intermediate values were to be stored on a sta k, we

ould exe ute

only a handful of time steps before we run out of memory even on our largest
workstation. Che kpointing is

ompulsory to

ompute the adjoint over several

thousands of time steps, whi h is our goal.
We saw in se tion 3 that TAPENADE applies
of subroutine

alls, i.e. ea h

all is

he kpointing at the level

he kpointed. This easy strategy is often

INRIA

Tangent and Adjoint dierentiation of OPA
far from optimal. On one hand several

13

alls are better not

TAPENADE now oers the option to mark sele ted
On the other hand,

he kpointed, and

alls for

not

he kpointing.

he kpointing should be applied at other lo ations. For

example at the top level of the simulation program is a loop over many time
steps.

We denitely need an e ient

he kpointing s heme applied at this

level of time iterations.
One

lassi al solution used by TAF on the MIT GCM

multi level re ursive

he kpointing. Basi ally, it splits the

ode [14℄, is

alled

omplete time inter-

val into a small number of equidistant intervals, then apply the same strategy
to ea h of the sub-intervals.

For instan e 64 time steps

an be split into 4

large intervals of 4 small intervals of 4 time steps, as sket hed on gure 3. This
onsumes a maximum of 9 simultaneous snapshots, and the average number
of dupli ate exe utions for a time step is
1000 time steps

2.25.

In a more realisti

situation,

an be split into 10 large intervals of 10 small intervals of

10 time steps, and one

an gure out that this

onsumes a maximum of 27

simultaneous snapshots, and the average number of dupli ate exe utions for a
time step is

2.7.

0

16

32

48

52

56

60 61 62
57 58

53 54
49 50
36

40

44 45 46
41 42

37 38
33 34
20

24

28 29 30
25 26

21 22
17 18
4

8

12 13 14
9 10

5
1

6

2

Figure 3: Three-levels
ward

he kpointing with 64 time steps and 9 snapshots. For-

omputations go right, adjoint

omputations go left. Bla k

resent writing/taking a snapshot, white
snapshot.

RR n° 6372

ir les rep-

ir les represent reading an available

Tber et al.

14

However, it was shown in [21℄ that this strategy is not optimal. Under the
reasonable assumptions that all time steps

ost the same run time, and that

the snapshot needed to run again from time step
from step

n

to any later step

n + x,

the optimal distribution of nested

n to n+1 is the same as to run

Griewank and Walther have

hara terized

he kpoints, whi h follows a binomial law.

With this optimal strategy, both spatial and temporal omplexity of the adjoint
ode grow logarithmi ally with respe t to the number of time steps of the
original simulation.

In other words, both the slowdown fa tor whi h grows

like the number of times ea h time step is exe uted, and the memory whi h
grows like the number of simultaneous snapshots, grow logarithmi ally with
the total number of time steps.
In real appli ations, run-time and memory spa e do not behave symmetrially. One

an always wait a little longer for the result, whereas the memory

spa e is bounded. Therefore the maximum number of snapshots

d

that

an

be stored simultaneously is xed.

Then [8℄ shows that the optimal strategy
th
gives a slowdown fa tor that grows only like the d root of the total number of

time steps, whi h is still very good. Figure 4 shows the optimal

he kpointing

strategy for the same problem as gure 3 i.e. 64 time steps with memory for
9 snapshots.
only

2.

The average number of dupli ate exe utions for a time step is

For the more realisti

situation (1000 time steps and memory for 27

2.57.

snapshots) the average number of dupli ate exe utions is only
We implemented this optimal strategy in the adjoint

ode of OPA. We

made our rst experiments by hand modi ation of the adjoint

ode produ ed

by TAPENADE. Still, TAPENADE produ ed automati ally the pro edures
that store and retrieve the snapshot, and therefore the hand modi ation was
benign: given the number of time steps, a general pro edure

3

s hedules the

optimal sequen e of a tions (store snapshot, retrieve snapshot, run time step,
run adjoint time step) to dierentiate the

omplete simulation. Further ver-

sions of TAPENADE will fully automate this pro ess.

Figure 5 shows the

performan es on OPA. They are in good agreement with the theory. Noti e in
parti ular the two small ine tion points on the
and 800 iterations.

urve around 150 iterations

Going ba k to the optimality proof in [8℄, we see that

the optimal strategy is parti ularly e ient when the number of time steps is
3A

FORTRAN95 implementation of this s heduling pro edure an be found in www.inriasop/tropi s/ftp/Hi ham_Tber/

INRIA

Tangent and Adjoint dierentiation of OPA

0

16

30

15

38

45

51

56

60

62

57 58
52 53 54
46 47 48 49

39 40 41 42 43

31 32 33 34 35 36

19

22

24 25 26 27 28

20
17
3

6

9

10 11 12 13 14

7
4
1

Figure 4: Optimal binomial

he kpointing with 64 time steps and 9 snapshots

7.5

Slowdown factor

7
6.5
6
5.5
5
4.5

200

400

600

800

1000

1200

Total number of time steps
Figure 5: Optimal binomial he kpointing with 15 snapshots: slowdown fa tor
as a fun tion of the total length of the initial simulation. The slowdown fa tor
is the run-time ratio of the adjoint

RR n° 6372

ode

ompared to the original

ode.

Tber et al.

16

exa tly

(d + t)!
d!t!
where d is the number of snapshots and t is the number of dupli ate exe utions
allowed per time step. For our target ma hine d = 15 and we nd η(15, 2) =
136 and η(15, 3) = 816, whi h orresponds to the ine tion points of gure 5.
η(d, t) =

For the previous version OPA 8, the adjoint was written by hand. Nevertheless, even a hand-written adjoint must implement strategies to retrieve intermediate states in reverse order that is, something very

lose to

Looking at this hand-written adjoint, we rst observe that the

he kpointing.
he kpointing

strategy is neither multi level nor optimal binomial. It is more like a single level
strategy, with one snapshot stored every xed number of time steps. During the
reverse sweep, states between two stored snapshots are rebuilt approximately
using linear interpolation. The advantage is that few time steps are evaluated
twi e, and therefore the slowdown fa tor remains well below 4. We

an see at

least two drawba ks. First, this hand manipulation requires deep knowledge of
the original program and of the underlying equations. This method does not
blend easily with Automati

Dierentiation. It is not yet automated in any

AD tool and therefore tedious and error-prone

ode manipulations would still

be ne essary. Se ond, this introdu es approximation errors into the

omputed

derivatives, whose mathemati al behavior is un lear. The gradient obtained in
the end is used in

omplex optimizations or data-assimilation loops, and small

errors may result in poor

onvergen e. In any

ase, for very large numbers of

time steps, we believe a trade-o between exa t binomial

he kpointing and

approximate interpolation is worth experimenting. Interpolation is probably
good enough for many variables that vary very slowly, and whi h

ould be

designated by the end-user, and only the other variables would need to be
stored.

4.4

Iterative linear solver

The OPA model solves an ellipti

equation at the end of ea h time step, using

an iterative method that generates a sequen e of approximations of the exa t
solution. The me hani al appli ation of AD on this kind of methods gives a
sequen e of derivatives of the approximate solutions with the same number

INRIA

Tangent and Adjoint dierentiation of OPA

17

of iterations as the original solver. The reason is that AD keeps the ow of
ontrol of the original program in the dierentiated program. In parti ular the
onvergen e tests are still based only on the non-dierentiated variables. Naturally, one may ask whether and how AD-produ ed derivatives are reasonable
approximations to the desired derivative of the exa t solution. The issues of
derivative

onvergen e for iterative solvers in relation to AD are dis ussed in

[6, 9, 2℄.
OPA provides two alternative algorithms to solve the ellipti
PCG for Pre onditioned Conjugate Gradient, and SOR for Su
Relaxation method. Both algorithms give

equation:

essive Over-

orre t results for the original

ode,

but PCG is generally preferred thanks to its e ien y and ve torization properties. However, the AD-dierentiated
algorithms. Figure 6

tives obtained by divided dieren es.
with the SOR algorithm remain
reases.
be ome

On the

ode gives dierent results using the two

ompares the AD-derivatives with approximate derivaWe see that the derivatives obtained

orre t when the number of time steps in-

ontrary, the derivatives obtained with the PCG algorithm

ompletely wrong after 80 time steps. Noti e that this o

urs in tan-

gent mode as well as in adjoint mode: the derivatives obtained with PCG,
although wrong, remain identi al in tangent and adjoint. Our explanation is
that ea h iteration of PCG involves the

omputation of s alar produ ts of vari-

ables that depend on the state ve tor, thus making the numeri al algorithm
nonlinear even though the ellipti

equation is linear. In [6℄, Gilbert has shown

that the appli ation of AD to a xed point iteration gives a derivative xed
point iteration that
in the

onverges R-linearly to the desired derivative in parti ular

ase of a large

this is not the

ontra tive iterate or se ant updating.

Unfortunately

ase for quasi-Newton iterative solvers su h as PCG, for whi h

there is no similar

onvergen e result to our knowledge.

To solve this problem for the tangent-dierentiated OPA we exploit the
linearity of the ellipti

system, and for the adjoint-dierentiated OPA we ex-

ploit the self adjointness property of the ellipti

operator [22℄. We

an thus

use the original PCG routine itself to solve for the dierentiated linear systems. Pra ti ally, we do this using the so- alled bla k-box feature provided
in TAPENADE. Figure 6 shows that (here for the tangent mode) the PCG
gives the same a

RR n° 6372

ura y as the SOR solver.

18

Tber et al.

Figure 6: evolution of the relative error between tangent derivative and divided
dieren es, for the three strategies: SOR and straightforward AD, PCG and
straightforward AD, PCG with the bla k box strategy

INRIA

Tangent and Adjoint dierentiation of OPA

19

In another experiment, we tried to use straightforward AD with the PCG
solver, but this time xing the number of PCG iterations to some very high
value. We observed that the derivatives be ome
dieren es. This
expensive and the

oherent again with divided

ould be another way to solve our problem, but it is

ertainly

hoi e of the high iteration number is deli ate. This prob-

lem denitely deserves further study, and

onrms the general re ommendation

not to dierentiate solvers of a nonlinear kind, and use a bla k-box strategy
instead.

5 Validation Experiments
5.1
The

Corre tness test
lassi al way to

gent and adjoint

he k for

orre tness of the automati ally generated tan-

odes is as follows:

1. Choose an arbitrary input

X

and and arbitrary dire tion

Ẋ .

Compute

the Divided Dieren e

DD =
for a good enough small

F (X + εẊ) − F (X)
ε

ε.

2. Using the tangent dierentiated program,
3. Using the adjoint dierentiated program,
and nally
test for the
derivative

he k that

Ẏ = F ′ (X) × Ẋ .

ompute
ompute

X = F ′∗ (X) × Ẏ

(DD · DD) = (Ẏ · Ẏ ) = (X · Ẋ).

We performed this

omplete global ORCA-2 simulation on 1000 time steps and its
odes.

The results are shown in table 1.

The values mat h, and

Table 1: Dot produ t test for 1000 time steps

(DD · DD) (ε = 10−7 )
(Ẏ · Ẏ )
(X · Ẋ)

RR n° 6372

4.405352760987440e+08
4.405346876439977e+08
4.405346876439867e+08

Tber et al.

20

(Ẏ · Ẏ )

and

(X · Ẋ)

mat h very well, up to the last few digits, whi h shows

that the tangent and adjoint
in a dierent
of

odes really

ompute the same derivatives, only

omputation order as shown by equations (2) and (3). The values

(DD · DD)

and

(Ẏ · Ẏ )

don't mat h so well, be ause of the weakness of the

Divided Dieren es approximation. Figure 7 shows this weakness: For a small
0

10

−2

10

−4

10

−6

10 −12
10

−10

−8

10

10

−6

ε

10

−4

−2

10

10

Figure 7: Relative error of Divided Dieren es with respe t to AD-generated
derivatives,

value of
of

ε,

ε,

omputed for various values of of the step size

the dominant error is due to ma hine a

ura y. For a large value

the dominant error is due to the se ond derivatives of

minimizes both errors, but

5.2

ε

annot eliminate them

F.

The best

ε

ompletely.

Sensitivity analysis on a long simulation

One of the main appli ation of adjoint models is the sensitivity analysis i.e. the
study of how model output varies with

hanges in model inputs. Using dire t

or statisti al methods would require many integration of the non linear model
while one adjoint model integration is enough to ompute this sensitivity. As an
example, gure 8 shows the output map of the sensitivity of the North Atlanti
meridional heat ux at 29°N to

hanges in the initial sea surfa e temperature

(SSTt0 ) over one year integration period, starting January 1, 1998.

This is

INRIA

Tangent and Adjoint dierentiation of OPA
done by

omputing the gradient respe t to

J=

Z

tN

t0

where

Ω

v

SSTt0

of

T.v dxdzdt

Ω

ross se tion at 29°N in the North Atlanti ,

is the zonal

temperature and

ZZ

21

is the meridional

T

is the

urrent velo ity.

Contours in gure 8 show where variation of initial SST would ee t the
most upon heat transport at 29°N. It shows large s ale patterns mainly lo ated

north of the 29°N parallel and in the Caribbean sea with a strong spot o
Moro

o. These results are

onsistent with those obtained by Marotzke et al.

([18℄)
This map was

omputed by the TAPENADE-generated adjoint of OPA on

the global ORCA2 grid, over 5475 time steps (1 year). This experiment was
done with the SOR algorithm as the iterative linear solver. The TAPENADEgenerated adjoint

omputed this sensitivity map in a time that is only 8.03

times that of the original simulation.

5.3

Data Assimilation

For further validation of the automati ally generated derivatives, we
out a data assimilation experiment.

arried

This was done in a so- alled twin ex-

periment framework whereby the dire t model traje tory is used to generate
syntheti

observations.

The initial sea surfa e temperature is perturbed by

a white noise and it has to be re overed using variational data assimilation
te hniques. Syntheti

observation are given by the sea surfa e height (SSH)

and the sea surfa e salinity (SSS) generated from the model's original outputs
starting from the unperturbed SST.
The

ost fun tion to be minimised is

J(SST (t0 )) =

Z

tN

k SSH(t) − SSH o(t) k2 + k SSS(t) − SSS o (t) k2 dt

(4)

t0

Where the supers ript

SSS(t)
For

o

stands for syntheti

observation and

SSH(t)

and

are model output.
omputing

ost issues, only the Antar ti

zoom of ORCA2 is

onsid-

ered, the minimisation is done an iterative gradient sear h algorithm where

RR n° 6372

Tber et al.

22

Figure 8: Sensitivity map of the North Atlanti
line), with respe t to

heat transport at 29°N (dotted

hanges in the initial surfa e temperature

INRIA

Tangent and Adjoint dierentiation of OPA
the gradient of

J

is

23

umputed using adjoint te hniques. Figure 9 illustrates

the performan e of the optimization loop for an integration period of 1 month
i.e. 450 time steps. The

ost fun tion de reases by two orders of magnitude.

Figure 10 indi ates that the true solution (top panel) is re overed with a
good approximation (bottom panel) from the randomly perturbed one(middle
panel), showing the quality of the derivatives obtained.

6

Cost function

10

5

10

4

10 0
10

1

2

10

10

Iterations

Figure 9: Twin experiment: Convergen e of the

ost fun tion

6 Con lusion and Outlook
The eort to build the tangent and adjoint

odes for the previous version 8

of the OPA o ean General Cir ulation Model has

ost several months devel-

opment from an experien ed resear her. For the new version OPA 9 written
in FORTRAN95, the use of the AD tool TAPENADE signi antly redu es
this eort. Our rst numeri al appli ations show the quality of the derivatives
obtained. This works validates the

hoi e of AD as the strategy to obtain the

tangent and adjoint for OPA 9, and for the versions to

ome.

At the same time, OPA is the largest FORTRAN95 appli ation dierentiated with TAPENADE. This work has pointed at a number of limitations of

RR n° 6372

24

Tber et al.

Figure 10: Twin experiment: True eld (top), Initial perturbed eld (middle)
and identied optimal sea surfa e temperatures (bottom)

INRIA

Tangent and Adjoint dierentiation of OPA

25

TAPENADE that have been lifted. Other limitations remain, su h as the nonreentrant pro edures, whi h need to be addressed in future work. Su
dierentiation of OPA denitely in reases our

essful

onden e in TAPENADE.

This works is also an additional illustration of the superiority of the binomial

he kpointing strategy,

ompared to multi level

he kpointing. By the

standards of other appli ation elds, e.g. CFD, a slowdown of the adjoint
of only 7 for a nonsteady simulation on 1000 time steps would be

ode

onsidered

very good. By the standards of weather simulation or o ean modeling however, s ientists expe t yet faster adjoints, at the
tion. Even if we

ost of a radi al approxima-

onsider that these approximations

hange the mathemati al

nature of the optimization pro ess, we understand they are ne essary and we
shall study how they

an be proposed as an option by the AD tool.

This work has underlined several dire tions for further resear h in AD and
AD tools. Some of them are already being studied by resear hers in our groups.
Considering the appli ation language, two

onstru ts need to be dierentiated

better:



The next experiment to be made very soon is to apply TAPENADE to
the parallelized version of OPA. This is ne essary before the generated
tangent and adjoint



odes

an be used in produ tion

ontext.

The OPA sour e makes extensive use of the prepro essor dire tives su h
as

#IFDEF. TAPENADE does not deal with these dire

tives be ause they

do not respe t the synta ti stru ture of a ode. Handling these dire tives
in the AD tool is in our opinion hopeless. What might be done though,
is to generate dierentiated

odes for ea h possible prepro essed

and devise a tool to put the dire tives ba k into the dierentiated
This is made easier if the dierentiated
of the original, as is the

ode

ode,
odes.

losely follows the stru ture

ase with TAPENADE.

Considering spe i ally adjoint dierentiation, we hope to obtain more e ient
ode through a more systemati
ellipti

operator.

exploitation of self-adjointness, e.g.

We also hope to optimize the

its present version, TAPENADE applies
Using proling information, we believe we
whi h

he kpointing is useless or

he kpointing strategy.

he kpointing to ea h pro edure
an dete t several pro edure

In
all.

alls for

ounter-produ tive. TAPENADE is already

able to use this information to produ e a better adjoint.

RR n° 6372

of the

Tber et al.

26

Referen es
[1℄ W. Castaings, D. Dartus, M. Honnorat, F.-X. Le Dimet, Y. Loukili, and
J. Monnier. Automati

dierentiation: A tool for variational data assimi-

lation and adjoint sensitivity analysis for ood modeling. In Bü ker et al,
editor,

Automati Dierentiation: Appli ations, Theory, and Implemen-

tations, LNCSE, pages 249262. Springer, 2005.
[2℄ B. Christianson. Reverse a

Opti-

umulation and attra tive xed points.

mization Methods and Software, 3:311326, 1994.

[3℄ P. Courtier, J.-N. Thépaut, and A. Hollingsworth. A strategy for operational implementation of 4d-var, using an in remental approa h.

Meteorol. So ., 120:13671387, 1994.

[4℄ R. Giering.

Q. J. R.

Tangent linear and Adjoint Model Compiler, Users manual,

1997. http://www.autodi. om/tam .
[5℄ R. Giering and T. Kaminski. Re ipes for adjoint

ode

ACM

onstru tion.

Transa tions on Mathemati al Software, 24(4):437474, 1998.

[6℄ J.C. Gilbert. Automati

dierentiation and iterative pro esses.

tion Methods and Software, 1:1321, 1992.

[7℄ M.B. Giles, D. Ghate, and M.C. Duta. Using automati
adjoint CFD

ode development. In Uthup et al, editor,

Aerospa e Design and Optimization,

Optimiza-

dierentiation for

Re ent Trends in

pages 363373. Tata-M Graw Hill,

New Delhi, 2005. Post-SAROD-2005, Bangalore, India.
[8℄ A. Griewank. A hieving logarithmi
plexity in reverse automati

Software, 1:3554, 1992.

growth of temporal and spatial

dierentiation.

om-

Optimization Methods and

[9℄ A. Griewank, C. Bis hof, G. Corliss, A. Carle, and K. Williamson. Derivative

onvergen e for iterative equation solvers.

Software, 2:321355, 1993.

Optimization Methods and

INRIA

Tangent and Adjoint dierentiation of OPA

27

[10℄ A. Griewank, D. Juedes, J. Srinivasan, and C. Tyner. ADOL-C, a pa kage
for the automati

dierentiation of algorithms written in C/C++.

Trans. Math. Software, 22(2):131167, 1996.
[11℄ L. Has oët and M. Araya-Polo.

The adjoint data-ow analyses:

malization, properties, and appli ations.

Appli ations, Theory, and Tools, Le

In

ACM
For-

Automati Dierentiation:

ture Notes in Computational S ien e

and Engineering. Springer, 2005. Sele ted papers from AD2004 Chi ago.
[12℄ L. Has oët and V. Pas ual. Tapenade 2.1 user's guide. Te hni al Report
0300, INRIA, 2004.
[13℄ L. Has oët, M. Vázquez, and A. Dervieux. Automati
optimum design, applied to soni

dierentiation for

boom redu tion. In Kumar et al., editor,

Pro eedings of ICCSA'03, Montreal, Canada,

pages 8594. LNCS 2668,

Springer, 2003.
[14℄ Heimba h, P. and Hill, C. and Giering, R.
of the parallel MIT general
dierentiation.

An e ient exa t adjoint

ir ulation model, generated via automati

Future Generation Computer Systems,

21(8):13561371,

2005.
[15℄ J. Kim, E. Hunke, and W. Lips omb. Sensitivity analysis and parameter
tuning s heme for global sea-i e modeling.

O ean Modeling Journal, 14(1

2):6180, 2006.
[16℄ C. Lauvernet, F. Baret, L. Has oët, and F.-X. LeDimet. Improved estimates of vegetation biophysi al variables from meris toa images by using
spatial and temporal

Beijing, China, 2005.

onstraints.

In

pro eedings of the 9th ISPMSRS,

[17℄ G. Made , P. Dele luse, M. Imbard, and C. Levy. Opa8.1 o ean general
ir ulation model referen e manual. Te hni al report, Pole de Modelisation, IPSL, 1998.
[18℄ J. Marotzke, R. Giering, K.Q. Zhang, D. Stammer, C. Hill, and T. Lee.
Constru tion of the adjoint MIT o ean General Cir ulation Model and

RR n° 6372

Tber et al.

28

appli ation to atlanti

heat transport sensitivity.

J. Geophys. Res.,

(104(C12)):29,52929,548, 1999.
[19℄ O. Talagrand. The use of adjoint equations in numeri al modeling of the

Automati Dierentiation of Algorithms: Theory, Implementation and Appliation, pages 169180, 1991. Philadelphia, Penn: SIAM.
atmospheri

ir ulation. In A. Griewank and G. Corliss, editors,

[20℄ J. Utke, U. Naumann, M. Fagan, N. Tallent, Strout M., P. Heimba h,
C. Hill, and C. Wuns h.
Automati

OpenAD/F: A modular, open-sour e tool for

Dierentiation of Fortran

odes. Te hni al report ANL/MCS-

P1230-0205, Argonne National Laboratory, 2006.

Submitted to ACM

TOMS.
[21℄ A. Walther and A. Griewank. Advantages of binomial

he kpointing for

Numeri al
Mathemati s and Advan ed Appli ations, pages 834843. Springer, Berlin,

memory-redu ed adjoint al ulations. In Feistauer et al, editor,
2003. Pro eeding of ENUMATH 2003.
[22℄ A. Weaver, J. Vialard, and D. Anderson.

Three- and four-dimensional

variational assimilation with an o ean general
tropi al pa i
he ks.

ir ulation model of the

o ean: I. formulation, internal diagnosti s and

Monthly Weather Review, 131(7):13601378, 2003.

onsisten y

INRIA

Unité de recherche INRIA Sophia Antipolis
2004, route des Lucioles - BP 93 - 06902 Sophia Antipolis Cedex (France)
Unité de recherche INRIA Futurs : Parc Club Orsay Université - ZAC des Vignes
4, rue Jacques Monod - 91893 ORSAY Cedex (France)
Unité de recherche INRIA Lorraine : LORIA, Technopôle de Nancy-Brabois - Campus scientifique
615, rue du Jardin Botanique - BP 101 - 54602 Villers-lès-Nancy Cedex (France)
Unité de recherche INRIA Rennes : IRISA, Campus universitaire de Beaulieu - 35042 Rennes Cedex (France)
Unité de recherche INRIA Rhône-Alpes : 655, avenue de l’Europe - 38334 Montbonnot Saint-Ismier (France)
Unité de recherche INRIA Rocquencourt : Domaine de Voluceau - Rocquencourt - BP 105 - 78153 Le Chesnay Cedex (France)

Éditeur
INRIA - Domaine de Voluceau - Rocquencourt, BP 105 - 78153 Le Chesnay Cedex (France)
http://www.inria.fr

ISSN 0249-6399

