/[MITgcm]/manual/s_autodiff/text/doc_ad_2.tex
ViewVC logotype

Diff of /manual/s_autodiff/text/doc_ad_2.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.3 by cnh, Thu Sep 27 02:00:24 2001 UTC revision 1.23 by jmc, Mon Aug 30 23:09:19 2010 UTC
# Line 1  Line 1 
1  % $Header$  % $Header$
2  % $Name$  % $Name$
3    
4    Author: Patrick Heimbach
5    \label{ask_the_author:doc_ad_2}
6    
7  {\sf Automatic differentiation} (AD), also referred to as algorithmic  {\sf Automatic differentiation} (AD), also referred to as algorithmic
8  (or, more loosely, computational) differentiation, involves  (or, more loosely, computational) differentiation, involves
9  automatically deriving code to calculate  automatically deriving code to calculate partial derivatives from an
10  partial derivatives from an existing fully non-linear prognostic code.  existing fully non-linear prognostic code.  (see \cite{gri:00}).  A
11  (see \cite{gri:00}).  software tool is used that parses and transforms source files
12  A software tool is used that parses and transforms source files  according to a set of linguistic and mathematical rules.  AD tools are
13  according to a set of linguistic and mathematical rules.  like source-to-source translators in that they parse a program code as
14  AD tools are like source-to-source translators in that  input and produce a new program code as output
15  they parse a program code as input and produce a new program code  (we restrict our discussion to source-to-source tools, ignoring
16  as output.  operator-overloading tools).  However, unlike a
17  However, unlike a pure source-to-source translation, the output program  pure source-to-source translation, the output program represents a new
18  represents a new algorithm, such as the evaluation of the  algorithm, such as the evaluation of the Jacobian, the Hessian, or
19  Jacobian, the Hessian, or higher derivative operators.  higher derivative operators.  In principle, a variety of derived
20  In principle, a variety of derived algorithms  algorithms can be generated automatically in this way.
21  can be generated automatically in this way.  
22    MITgcm has been adapted for use with the Tangent linear and Adjoint
23  The MITGCM has been adapted for use with the  Model Compiler (TAMC) and its successor TAF (Transformation of
24  Tangent linear and Adjoint Model Compiler (TAMC) and its succssor TAF  Algorithms in Fortran), developed by Ralf Giering (\cite{gie-kam:98},
25  (Transformation of Algorithms in Fortran), developed  \cite{gie:99,gie:00}).  The first application of the adjoint of MITgcm
26  by Ralf Giering (\cite{gie-kam:98}, \cite{gie:99,gie:00}).  for sensitivity studies has been published by \cite{maro-eta:99}.
27  The first application of the adjoint of the MITGCM for senistivity  \cite{stam-etal:97,stam-etal:02} use MITgcm and its adjoint for ocean
28  studies has been published by \cite{maro-eta:99}.  state estimation studies.  In the following we shall refer to TAMC and
29  \cite{sta-eta:97,sta-eta:01} use the MITGCM and its adjoint  TAF synonymously, except were explicitly stated otherwise.
30  for ocean state estimation studies.  
31    As of mid-2007 we are also able to generate fairly efficient
32  TAMC exploits the chain rule for computing the first  adjoint code of the MITgcm using a new, open-source AD tool,
33  derivative of a function with  called OpenAD (see \cite{naum-etal:06,utke-etal:08}.
34  respect to a set of input variables.  This enables us for the first time to compare adjoint models
35  Treating a given forward code as a composition of operations --  generated from different AD tools, providing an additional
36  each line representing a compositional element -- the chain rule is  accuracy check, complementary to finite-difference gradient checks.
37  rigorously applied to the code, line by line. The resulting  OpenAD and its application to  MITgcm is described in detail
38  tangent linear or adjoint code,  in section \ref{sec_ad_openad}.
39  then, may be thought of as the composition in  
40  forward or reverse order, respectively, of the  The AD tool exploits the chain rule for computing the first derivative of a
41  Jacobian matrices of the forward code compositional elements.  function with respect to a set of input variables.  Treating a given
42    forward code as a composition of operations -- each line representing
43    a compositional element, the chain rule is rigorously applied to the
44    code, line by line. The resulting tangent linear or adjoint code,
45    then, may be thought of as the composition in forward or reverse
46    order, respectively, of the Jacobian matrices of the forward code's
47    compositional elements.
48    
49  %**********************************************************************  %**********************************************************************
50  \section{Some basic algebra}  \section{Some basic algebra}
51  \label{sec_ad_algebra}  \label{sec_ad_algebra}
52    \begin{rawhtml}
53    <!-- CMIREDIR:sec_ad_algebra: -->
54    \end{rawhtml}
55  %**********************************************************************  %**********************************************************************
56    
57  Let $ \cal{M} $ be a general nonlinear, model, i.e. a  Let $ \cal{M} $ be a general nonlinear, model, i.e. a
# Line 50  $\vec{u}=(u_1,\ldots,u_m)$ Line 62  $\vec{u}=(u_1,\ldots,u_m)$
62  such as forcing functions) to the $n$-dimensional space  such as forcing functions) to the $n$-dimensional space
63  $V \subset I\!\!R^n$ of  $V \subset I\!\!R^n$ of
64  model output variable $\vec{v}=(v_1,\ldots,v_n)$  model output variable $\vec{v}=(v_1,\ldots,v_n)$
65  (model state, model diagnostcs, objective function, ...)  (model state, model diagnostics, objective function, ...)
66  under consideration,  under consideration,
67  %  %
68  \begin{equation}  \begin{equation}
69  \begin{split}  \begin{aligned}
70  {\cal M} \, : & \, U \,\, \longrightarrow \, V \\  {\cal M} \, : & \, U \,\, \longrightarrow \, V \\
71  ~      & \, \vec{u} \,\, \longmapsto \, \vec{v} \, = \,  ~      & \, \vec{u} \,\, \longmapsto \, \vec{v} \, = \,
72  {\cal M}(\vec{u})  {\cal M}(\vec{u})
73  \label{fulloperator}  \label{fulloperator}
74  \end{split}  \end{aligned}
75  \end{equation}  \end{equation}
76  %  %
77  The vectors $ \vec{u} \in U $ and $ v \in V $ may be represented w.r.t.  The vectors $ \vec{u} \in U $ and $ v \in V $ may be represented w.r.t.
# Line 105  In contrast to the full nonlinear model Line 117  In contrast to the full nonlinear model
117  $ M $ is just a matrix  $ M $ is just a matrix
118  which can readily be used to find the forward sensitivity of $\vec{v}$ to  which can readily be used to find the forward sensitivity of $\vec{v}$ to
119  perturbations in  $u$,  perturbations in  $u$,
120  but if there are very many input variables $(>>O(10^{6})$ for  but if there are very many input variables $(\gg O(10^{6})$ for
121  large-scale oceanographic application), it quickly becomes  large-scale oceanographic application), it quickly becomes
122  prohibitive to proceed directly as in (\ref{tangent_linear}),  prohibitive to proceed directly as in (\ref{tangent_linear}),
123  if the impact of each component $ {\bf e_{i}} $ is to be assessed.  if the impact of each component $ {\bf e_{i}} $ is to be assessed.
# Line 130  or a measure of some model-to-data misfi Line 142  or a measure of some model-to-data misfi
142  \label{compo}  \label{compo}
143  \end{eqnarray}  \end{eqnarray}
144  %  %
145  The linear approximation of $ {\cal J} $,  The perturbation of $ {\cal J} $ around a fixed point $ {\cal J}_0 $,
146  \[  \[
147  {\cal J} \, \approx \, {\cal J}_0 \, + \, \delta {\cal J}  {\cal J} \, = \, {\cal J}_0 \, + \, \delta {\cal J}
148  \]  \]
149  can be expressed in both bases of $ \vec{u} $ and $ \vec{v} $  can be expressed in both bases of $ \vec{u} $ and $ \vec{v} $
150  w.r.t. their corresponding inner product  w.r.t. their corresponding inner product
151  $\left\langle \,\, , \,\, \right\rangle $  $\left\langle \,\, , \,\, \right\rangle $
152  %  %
153  \begin{equation}  \begin{equation}
154  \begin{split}  \begin{aligned}
155  {\cal J} & = \,  {\cal J} & = \,
156  {\cal J} |_{\vec{u}^{(0)}} \, + \,  {\cal J} |_{\vec{u}^{(0)}} \, + \,
157  \left\langle \, \nabla _{u}{\cal J}^T |_{\vec{u}^{(0)}} \, , \, \delta \vec{u} \, \right\rangle  \left\langle \, \nabla _{u}{\cal J}^T |_{\vec{u}^{(0)}} \, , \, \delta \vec{u} \, \right\rangle
# Line 148  $\left\langle \,\, , \,\, \right\rangle Line 160  $\left\langle \,\, , \,\, \right\rangle
160  {\cal J} |_{\vec{v}^{(0)}} \, + \,  {\cal J} |_{\vec{v}^{(0)}} \, + \,
161  \left\langle \, \nabla _{v}{\cal J}^T |_{\vec{v}^{(0)}} \, , \, \delta \vec{v} \, \right\rangle  \left\langle \, \nabla _{v}{\cal J}^T |_{\vec{v}^{(0)}} \, , \, \delta \vec{v} \, \right\rangle
162  \, + \, O(\delta \vec{v}^2)  \, + \, O(\delta \vec{v}^2)
163  \end{split}  \end{aligned}
164  \label{deljidentity}  \label{deljidentity}
165  \end{equation}  \end{equation}
166  %  %
# Line 168  transpose of $ A $, Line 180  transpose of $ A $,
180  \[  \[
181  A^{\ast} \, = \, A^T  A^{\ast} \, = \, A^T
182  \]  \]
183  and from eq. (\ref{tangent_linear}), we note that  and from eq. (\ref{tangent_linear}), (\ref{deljidentity}),
184    we note that
185  (omitting $|$'s):  (omitting $|$'s):
186  %  %
187  \begin{equation}  \begin{equation}
# Line 188  the gradient $ \nabla _{u}{\cal J} $ can Line 201  the gradient $ \nabla _{u}{\cal J} $ can
201  invoking the adjoint $ M^{\ast } $ of the tangent linear model $ M $  invoking the adjoint $ M^{\ast } $ of the tangent linear model $ M $
202  %  %
203  \begin{equation}  \begin{equation}
204  \begin{split}  \begin{aligned}
205  \nabla _{u}{\cal J}^T |_{\vec{u}} &  \nabla _{u}{\cal J}^T |_{\vec{u}} &
206  = \, M^T |_{\vec{u}} \cdot \nabla _{v}{\cal J}^T |_{\vec{v}}  \\  = \, M^T |_{\vec{u}} \cdot \nabla _{v}{\cal J}^T |_{\vec{v}}  \\
207  ~ & = \, M^T |_{\vec{u}} \cdot \delta \vec{v}^{\ast} \\  ~ & = \, M^T |_{\vec{u}} \cdot \delta \vec{v}^{\ast} \\
208  ~ & = \, \delta \vec{u}^{\ast}  ~ & = \, \delta \vec{u}^{\ast}
209  \end{split}  \end{aligned}
210  \label{adjoint}  \label{adjoint}
211  \end{equation}  \end{equation}
212  %  %
# Line 204  the adjoint variable of the model state Line 217  the adjoint variable of the model state
217  $ \delta \vec{u}^{\ast} $ the adjoint variable of the control variable $ \vec{u} $.  $ \delta \vec{u}^{\ast} $ the adjoint variable of the control variable $ \vec{u} $.
218    
219  The {\sf reverse} nature of the adjoint calculation can be readily  The {\sf reverse} nature of the adjoint calculation can be readily
220  seen as follows. Let us decompose ${\cal J}(u)$, thus:  seen as follows.
221    Consider a model integration which consists of $ \Lambda $
222    consecutive operations
223    $ {\cal M}_{\Lambda} (  {\cal M}_{\Lambda-1} (
224    ...... ( {\cal M}_{\lambda} (
225    ......
226    ( {\cal M}_{1} ( {\cal M}_{0}(\vec{u}) )))) $,
227    where the ${\cal M}$'s could be the elementary steps, i.e. single lines
228    in the code of the model, or successive time steps of the
229    model integration,
230    starting at step 0 and moving up to step $\Lambda$, with intermediate
231    ${\cal M}_{\lambda} (\vec{u}) = \vec{v}^{(\lambda+1)}$ and final
232    ${\cal M}_{\Lambda} (\vec{u}) = \vec{v}^{(\Lambda+1)} = \vec{v}$.
233    Let ${\cal J}$ be a cost function which explicitly depends on the
234    final state $\vec{v}$ only
235    (this restriction is for clarity reasons only).
236    %
237    ${\cal J}(u)$ may be decomposed according to:
238  %  %
239  \begin{equation}  \begin{equation}
240  {\cal J}({\cal M}(\vec{u})) \, = \,  {\cal J}({\cal M}(\vec{u})) \, = \,
# Line 215  seen as follows. Let us decompose ${\cal Line 245  seen as follows. Let us decompose ${\cal
245  \label{compos}  \label{compos}
246  \end{equation}  \end{equation}
247  %  %
248  where the ${\cal M}$'s could be the elementary steps, i.e. single lines  Then, according to the chain rule, the forward calculation reads,
249  in the code of the model,  in terms of the Jacobi matrices
 starting at step 0 and moving up to step $\Lambda$, with intermediate  
 ${\cal M}_{\lambda} (\vec{u}) = \vec{v}^{(\lambda+1)}$ and final  
 ${\cal M}_{\Lambda} (\vec{u}) = \vec{v}^{(\Lambda+1)} = \vec{v}$  
 Then, according to the chain rule the forward calculation reads in  
 terms of the Jacobi matrices  
250  (we've omitted the $ | $'s which, nevertheless are important  (we've omitted the $ | $'s which, nevertheless are important
251  to the aspect of {\it tangent} linearity;  to the aspect of {\it tangent} linearity;
252  note also that per definition  note also that by definition
253  $ \langle \, \nabla _{v}{\cal J}^T \, , \, \delta \vec{v} \, \rangle  $ \langle \, \nabla _{v}{\cal J}^T \, , \, \delta \vec{v} \, \rangle
254  = \nabla_v {\cal J} \cdot \delta \vec{v} $ )  = \nabla_v {\cal J} \cdot \delta \vec{v} $ )
255  %  %
256  \begin{equation}  \begin{equation}
257  \begin{split}  \begin{aligned}
258  \nabla_v {\cal J} (M(\delta \vec{u})) & = \,  \nabla_v {\cal J} (M(\delta \vec{u})) & = \,
259  \nabla_v {\cal J} \cdot M_{\Lambda}  \nabla_v {\cal J} \cdot M_{\Lambda}
260  \cdot ...... \cdot M_{\lambda} \cdot ...... \cdot  \cdot ...... \cdot M_{\lambda} \cdot ...... \cdot
261  M_{1} \cdot M_{0} \cdot \delta \vec{u} \\  M_{1} \cdot M_{0} \cdot \delta \vec{u} \\
262  ~ & = \, \nabla_v {\cal J} \cdot \delta \vec{v} \\  ~ & = \, \nabla_v {\cal J} \cdot \delta \vec{v} \\
263  \end{split}  \end{aligned}
264  \label{forward}  \label{forward}
265  \end{equation}  \end{equation}
266  %  %
# Line 243  whereas in reverse mode we have Line 268  whereas in reverse mode we have
268  %  %
269  \begin{equation}  \begin{equation}
270  \boxed{  \boxed{
271  \begin{split}  \begin{aligned}
272  M^T ( \nabla_v {\cal J}^T) & = \,  M^T ( \nabla_v {\cal J}^T) & = \,
273  M_{0}^T \cdot M_{1}^T  M_{0}^T \cdot M_{1}^T
274  \cdot ...... \cdot M_{\lambda}^T \cdot ...... \cdot  \cdot ...... \cdot M_{\lambda}^T \cdot ...... \cdot
# Line 252  M_{\Lambda}^T \cdot \nabla_v {\cal J}^T Line 277  M_{\Lambda}^T \cdot \nabla_v {\cal J}^T
277  \cdot ...... \cdot  \cdot ...... \cdot
278  \nabla_{v^{(\lambda)}} {\cal J}^T \\  \nabla_{v^{(\lambda)}} {\cal J}^T \\
279  ~ & = \, \nabla_u {\cal J}^T  ~ & = \, \nabla_u {\cal J}^T
280  \end{split}  \end{aligned}
281  }  }
282  \label{reverse}  \label{reverse}
283  \end{equation}  \end{equation}
284  %  %
285  clearly expressing the reverse nature of the calculation.  clearly expressing the reverse nature of the calculation.
286  Eq. (\ref{reverse}) is at the heart of automatic adjoint compilers.  Eq. (\ref{reverse}) is at the heart of automatic adjoint compilers.
287  The intermediate steps $\lambda$ in  If the intermediate steps $\lambda$ in
288  eqn. (\ref{compos}) -- (\ref{reverse})  eqn. (\ref{compos}) -- (\ref{reverse})
289  could represent the model state (forward or adjoint) at each  represent the model state (forward or adjoint) at each
290  intermediate time step in which case  intermediate time step as noted above, then correspondingly,
291  $ {\cal M}(\vec{v}^{(\lambda)}) = \vec{v}^{(\lambda+1)} $, and correspondingly,  $ M^T (\delta \vec{v}^{(\lambda) \, \ast}) =
292  $ M^T (\delta \vec{v}^{(\lambda) \, \ast}) = \delta \vec{v}^{(\lambda-1) \, \ast} $,  \delta \vec{v}^{(\lambda-1) \, \ast} $ for the adjoint variables.
293  but they can also be viewed more generally as  It thus becomes evident that the adjoint calculation also
294  single lines of code in the numerical algorithm.  yields the adjoint of each model state component
295  In both cases it becomes evident that the adjoint calculation  $ \vec{v}^{(\lambda)} $ at each intermediate step $ \lambda $, namely
 yields at the same time the adjoint of each model state component  
 $ \vec{v}^{(\lambda)} $ at each intermediate step $ l $, namely  
296  %  %
297  \begin{equation}  \begin{equation}
298  \boxed{  \boxed{
299  \begin{split}  \begin{aligned}
300  \nabla_{v^{(\lambda)}} {\cal J}^T |_{\vec{v}^{(\lambda)}}  \nabla_{v^{(\lambda)}} {\cal J}^T |_{\vec{v}^{(\lambda)}}
301  & = \,  & = \,
302  M_{\lambda}^T |_{\vec{v}^{(\lambda)}} \cdot ...... \cdot  M_{\lambda}^T |_{\vec{v}^{(\lambda)}} \cdot ...... \cdot
303  M_{\Lambda}^T |_{\vec{v}^{(\lambda)}} \cdot \delta \vec{v}^{\ast} \\  M_{\Lambda}^T |_{\vec{v}^{(\lambda)}} \cdot \delta \vec{v}^{\ast} \\
304  ~ & = \, \delta \vec{v}^{(\lambda) \, \ast}  ~ & = \, \delta \vec{v}^{(\lambda) \, \ast}
305  \end{split}  \end{aligned}
306  }  }
307  \end{equation}  \end{equation}
308  %  %
309  in close analogy to eq. (\ref{adjoint})  in close analogy to eq. (\ref{adjoint})
310  We note in passing that that the $\delta \vec{v}^{(\lambda) \, \ast}$  We note in passing that that the $\delta \vec{v}^{(\lambda) \, \ast}$
311  are the Lagrange multipliers of the model state $ \vec{v}^{(\lambda)}$.  are the Lagrange multipliers of the model equations which determine
312    $ \vec{v}^{(\lambda)}$.
313    
314  In coponents, eq. (\ref{adjoint}) reads as follows.  In components, eq. (\ref{adjoint}) reads as follows.
315  Let  Let
316  \[  \[
317  \begin{array}{rclcrcl}  \begin{array}{rclcrcl}
# Line 308  Let Line 332  Let
332  \end{array}  \end{array}
333  \]  \]
334  denote the perturbations in $\vec{u}$ and $\vec{v}$, respectively,  denote the perturbations in $\vec{u}$ and $\vec{v}$, respectively,
335  and their adjoint varaiables;  and their adjoint variables;
336  further  further
337  \[  \[
338  M \, = \, \left(  M \, = \, \left(
# Line 395  and the shorthand notation for the adjoi Line 419  and the shorthand notation for the adjoi
419  $ \delta v^{(\lambda) \, \ast}_{j} = \frac{\partial}{\partial v^{(\lambda)}_{j}}  $ \delta v^{(\lambda) \, \ast}_{j} = \frac{\partial}{\partial v^{(\lambda)}_{j}}
420  {\cal J}^T $, $ j = 1, \ldots , n_{\lambda} $,  {\cal J}^T $, $ j = 1, \ldots , n_{\lambda} $,
421  for intermediate components, yielding  for intermediate components, yielding
422  \[  {\small
423  \footnotesize  \begin{equation}
424    \begin{aligned}
425  \left(  \left(
426  \begin{array}{c}  \begin{array}{c}
427  \delta v^{(\lambda) \, \ast}_1 \\  \delta v^{(\lambda) \, \ast}_1 \\
# Line 404  for intermediate components, yielding Line 429  for intermediate components, yielding
429  \delta v^{(\lambda) \, \ast}_{n_{\lambda}} \\  \delta v^{(\lambda) \, \ast}_{n_{\lambda}} \\
430  \end{array}  \end{array}
431  \right)  \right)
432  \, = \,  \, = &
433  \left(  \left(
434  \begin{array}{ccc}  \begin{array}{ccc}
435  \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_1}  \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_1}
436  & \ldots &  & \ldots \,\, \ldots &
437  \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_1} \\  \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_1} \\
438  \vdots & ~ & \vdots \\  \vdots & ~ & \vdots \\
439  \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_{n_{\lambda}}}  \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_{n_{\lambda}}}
440  & \ldots  &  & \ldots \,\, \ldots  &
441  \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_{n_{\lambda}}} \\  \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_{n_{\lambda}}} \\
442  \end{array}  \end{array}
443  \right)  \right)
 %  
444  \cdot  \cdot
445  %  %
446    \\ ~ & ~
447    \\ ~ &
448    %
449  \left(  \left(
450  \begin{array}{ccc}  \begin{array}{ccc}
451  \frac{\partial ({\cal M}_{\lambda+1})_1}{\partial v^{(\lambda+1)}_1}  \frac{\partial ({\cal M}_{\lambda+1})_1}{\partial v^{(\lambda+1)}_1}
# Line 431  for intermediate components, yielding Line 458  for intermediate components, yielding
458  \frac{\partial ({\cal M}_{\lambda+1})_{n_{\lambda+2}}}{\partial v^{(\lambda+1)}_{n_{\lambda+1}}} \\  \frac{\partial ({\cal M}_{\lambda+1})_{n_{\lambda+2}}}{\partial v^{(\lambda+1)}_{n_{\lambda+1}}} \\
459  \end{array}  \end{array}
460  \right)  \right)
461  \cdot \ldots \ldots \cdot  \cdot \, \ldots \, \cdot
462  \left(  \left(
463  \begin{array}{c}  \begin{array}{c}
464  \delta v^{\ast}_1 \\  \delta v^{\ast}_1 \\
# Line 439  for intermediate components, yielding Line 466  for intermediate components, yielding
466  \delta v^{\ast}_{n} \\  \delta v^{\ast}_{n} \\
467  \end{array}  \end{array}
468  \right)  \right)
469  \]  \end{aligned}
470    \end{equation}
471    }
472    
473  Eq. (\ref{forward}) and (\ref{reverse}) are perhaps clearest in  Eq. (\ref{forward}) and (\ref{reverse}) are perhaps clearest in
474  showing the advantage of the reverse over the forward mode  showing the advantage of the reverse over the forward mode
# Line 450  variables $u$ Line 479  variables $u$
479  {\it all} intermediate states $ \vec{v}^{(\lambda)} $) are sought.  {\it all} intermediate states $ \vec{v}^{(\lambda)} $) are sought.
480  In order to be able to solve for each component of the gradient  In order to be able to solve for each component of the gradient
481  $ \partial {\cal J} / \partial u_{i} $ in (\ref{forward})  $ \partial {\cal J} / \partial u_{i} $ in (\ref{forward})
482  a forward calulation has to be performed for each component seperately,  a forward calculation has to be performed for each component separately,
483  i.e. $ \delta \vec{u} = \delta u_{i} {\vec{e}_{i}} $  i.e. $ \delta \vec{u} = \delta u_{i} {\vec{e}_{i}} $
484  for  the $i$-th forward calculation.  for  the $i$-th forward calculation.
485  Then, (\ref{forward}) represents the  Then, (\ref{forward}) represents the
# Line 460  In contrast, eq. (\ref{reverse}) yields Line 489  In contrast, eq. (\ref{reverse}) yields
489  gradient $\nabla _{u}{\cal J}$ (and all intermediate gradients  gradient $\nabla _{u}{\cal J}$ (and all intermediate gradients
490  $\nabla _{v^{(\lambda)}}{\cal J}$) within a single reverse calculation.  $\nabla _{v^{(\lambda)}}{\cal J}$) within a single reverse calculation.
491    
492  Note, that in case $ {\cal J} $ is a vector-valued function  Note, that if $ {\cal J} $ is a vector-valued function
493  of dimension $ l > 1 $,  of dimension $ l > 1 $,
494  eq. (\ref{reverse}) has to be modified according to  eq. (\ref{reverse}) has to be modified according to
495  \[  \[
# Line 468  M^T \left( \nabla_v {\cal J}^T \left(\de Line 497  M^T \left( \nabla_v {\cal J}^T \left(\de
497  \, = \,  \, = \,
498  \nabla_u {\cal J}^T \cdot \delta \vec{J}  \nabla_u {\cal J}^T \cdot \delta \vec{J}
499  \]  \]
500  where now $ \delta \vec{J} \in I\!\!R $ is a vector of dimenison $ l $.  where now $ \delta \vec{J} \in I\!\!R^l $ is a vector of
501    dimension $ l $.
502  In this case $ l $ reverse simulations have to be performed  In this case $ l $ reverse simulations have to be performed
503  for each $ \delta J_{k}, \,\, k = 1, \ldots, l $.  for each $ \delta J_{k}, \,\, k = 1, \ldots, l $.
504  Then, the reverse mode is more efficient as long as  Then, the reverse mode is more efficient as long as
505  $ l < n $, otherwise the forward mode is preferable.  $ l < n $, otherwise the forward mode is preferable.
506  Stricly, the reverse mode is called adjoint mode only for  Strictly, the reverse mode is called adjoint mode only for
507  $ l = 1 $.  $ l = 1 $.
508    
509  A detailed analysis of the underlying numerical operations  A detailed analysis of the underlying numerical operations
# Line 503  operator onto the $j$-th component ${\bf Line 533  operator onto the $j$-th component ${\bf
533  \paragraph{Example 2:  \paragraph{Example 2:
534  $ {\cal J} = \langle \, {\cal H}(\vec{v}) - \vec{d} \, ,  $ {\cal J} = \langle \, {\cal H}(\vec{v}) - \vec{d} \, ,
535   \, {\cal H}(\vec{v}) - \vec{d} \, \rangle $} ~ \\   \, {\cal H}(\vec{v}) - \vec{d} \, \rangle $} ~ \\
536  The cost function represents the quadratic model vs.data misfit.  The cost function represents the quadratic model vs. data misfit.
537  Here, $ \vec{d} $ is the data vector and $ {\cal H} $ represents the  Here, $ \vec{d} $ is the data vector and $ {\cal H} $ represents the
538  operator which maps the model state space onto the data space.  operator which maps the model state space onto the data space.
539  Then, $ \nabla_v {\cal J} $ takes the form  Then, $ \nabla_v {\cal J} $ takes the form
540  %  %
541  \begin{equation*}  \begin{equation*}
542  \begin{split}  \begin{aligned}
543  \nabla_v {\cal J}^T & = \, 2 \, \, H \cdot  \nabla_v {\cal J}^T & = \, 2 \, \, H \cdot
544  \left( \, {\cal H}(\vec{v}) - \vec{d} \, \right) \\  \left( \, {\cal H}(\vec{v}) - \vec{d} \, \right) \\
545  ~          & = \, 2 \sum_{j} \left\{ \sum_k  ~          & = \, 2 \sum_{j} \left\{ \sum_k
546  \frac{\partial {\cal H}_k}{\partial v_{j}}  \frac{\partial {\cal H}_k}{\partial v_{j}}
547  \left( {\cal H}_k (\vec{v}) - d_k \right)  \left( {\cal H}_k (\vec{v}) - d_k \right)
548  \right\} \, {\vec{f}_{j}} \\  \right\} \, {\vec{f}_{j}} \\
549  \end{split}  \end{aligned}
550  \end{equation*}  \end{equation*}
551  %  %
552  where $H_{kj} = \partial {\cal H}_k / \partial v_{j} $ is the  where $H_{kj} = \partial {\cal H}_k / \partial v_{j} $ is the
# Line 534  H \cdot \left( {\cal H}(\vec{v}) - \vec{ Line 564  H \cdot \left( {\cal H}(\vec{v}) - \vec{
564    
565  We note an important aspect of the forward vs. reverse  We note an important aspect of the forward vs. reverse
566  mode calculation.  mode calculation.
567  Because of the locality of the derivative,  Because of the local character of the derivative
568    (a derivative is defined w.r.t. a point along the trajectory),
569  the intermediate results of the model trajectory  the intermediate results of the model trajectory
570  $\vec{v}^{(\lambda+1)}={\cal M}_{\lambda}(v^{(\lambda)})$  $\vec{v}^{(\lambda+1)}={\cal M}_{\lambda}(v^{(\lambda)})$
571  are needed to evaluate the intermediate Jacobian  may be required to evaluate the intermediate Jacobian
572  $M_{\lambda}|_{\vec{v}^{(\lambda)}} \, \delta \vec{v}^{(\lambda)} $.  $M_{\lambda}|_{\vec{v}^{(\lambda)}} \, \delta \vec{v}^{(\lambda)} $.
573    This is the case e.g. for nonlinear expressions
574    (momentum advection, nonlinear equation of state), state-dependent
575    conditional statements (parameterization schemes).
576  In the forward mode, the intermediate results are required  In the forward mode, the intermediate results are required
577  in the same order as computed by the full forward model ${\cal M}$,  in the same order as computed by the full forward model ${\cal M}$,
578  in the reverse mode they are required in the reverse order.  but in the reverse mode they are required in the reverse order.
579  Thus, in the reverse mode the trajectory of the forward model  Thus, in the reverse mode the trajectory of the forward model
580  integration ${\cal M}$ has to be stored to be available in the reverse  integration ${\cal M}$ has to be stored to be available in the reverse
581  calculation. Alternatively, the model state would have to be  calculation. Alternatively, the complete model state up to the
582  recomputed whenever its value is required.  point of evaluation has to be recomputed whenever its value is required.
583    
584  A method to balance the amount of recomputations vs.  A method to balance the amount of recomputations vs.
585  storage requirements is called {\sf checkpointing}  storage requirements is called {\sf checkpointing}
586  (e.g. \cite{res-eta:98}).  (e.g. \cite{gri:92}, \cite{res-eta:98}).
587  It is depicted in Fig. ... for a 3-level checkpointing  It is depicted in \ref{fig:3levelcheck} for a 3-level checkpointing
588  [as concrete example, we give explicit numbers for a 3-day  [as an example, we give explicit numbers for a 3-day
589  integration with a 1-hourly timestep in square brackets].  integration with a 1-hourly timestep in square brackets].
590  \begin{itemize}  \begin{itemize}
591  %  %
# Line 559  integration with a 1-hourly timestep in Line 593  integration with a 1-hourly timestep in
593  In a first step, the model trajectory is subdivided into  In a first step, the model trajectory is subdivided into
594  $ {n}^{lev3} $ subsections [$ {n}^{lev3} $=3 1-day intervals],  $ {n}^{lev3} $ subsections [$ {n}^{lev3} $=3 1-day intervals],
595  with the label $lev3$ for this outermost loop.  with the label $lev3$ for this outermost loop.
596  The model is then integrated over the full trajectory,  The model is then integrated along the full trajectory,
597  and the model state stored only at every $ k_{i}^{lev3} $-th timestep  and the model state stored to disk only at every $ k_{i}^{lev3} $-th timestep
598  [i.e. 3 times, at  [i.e. 3 times, at
599  $ i = 0,1,2 $ corresponding to $ k_{i}^{lev3} = 0, 24, 48 $].  $ i = 0,1,2 $ corresponding to $ k_{i}^{lev3} = 0, 24, 48 $].
600    In addition, the cost function is computed, if needed.
601  %  %
602  \item [$lev2$]  \item [$lev2$]
603  In a second step each subsection is itself divided into  In a second step each subsection itself is divided into
604  $ {n}^{lev2} $ subsubsections  $ {n}^{lev2} $ subsections
605  [$ {n}^{lev2} $=4 6-hour intervals per subsection].  [$ {n}^{lev2} $=4 6-hour intervals per subsection].
606  The model picks up at the last outermost dumped state  The model picks up at the last outermost dumped state
607  $ v_{k_{n}^{lev3}} $ and is integrated forward in time over  $ v_{k_{n}^{lev3}} $ and is integrated forward in time along
608  the last subsection, with the label $lev2$ for this    the last subsection, with the label $lev2$ for this  
609  intermediate loop.  intermediate loop.
610  The model state is now stored only at every $ k_{i}^{lev2} $-th  The model state is now stored to disk at every $ k_{i}^{lev2} $-th
611  timestep  timestep
612  [i.e. 4 times, at  [i.e. 4 times, at
613  $ i = 0,1,2,3 $ corresponding to $ k_{i}^{lev2} = 48, 54, 60, 66 $].  $ i = 0,1,2,3 $ corresponding to $ k_{i}^{lev2} = 48, 54, 60, 66 $].
614  %  %
615  \item [$lev1$]  \item [$lev1$]
616  Finally, the mode picks up at the last intermediate dump state  Finally, the model picks up at the last intermediate dump state
617  $ v_{k_{n}^{lev2}} $ and is integrated forward in time over  $ v_{k_{n}^{lev2}} $ and is integrated forward in time along
618  the last subsubsection, with the label $lev1$ for this    the last subsection, with the label $lev1$ for this  
619  intermediate loop.  intermediate loop.
620  Within this subsubsection only, the model state is stored  Within this sub-subsection only, parts of the model state is stored
621  at every timestep  to memory at every timestep
622  [i.e. every hour $ i=0,...,5$ corresponding to  [i.e. every hour $ i=0,...,5$ corresponding to
623  $ k_{i}^{lev1} = 66, 67, \ldots, 71 $].  $ k_{i}^{lev1} = 66, 67, \ldots, 71 $].
624  Thus, the  final state $ v_n = v_{k_{n}^{lev1}} $ is reached  The  final state $ v_n = v_{k_{n}^{lev1}} $ is reached
625  and the model state of all peceeding timesteps over the last  and the model state of all preceding timesteps along the last
626  subsubsections are available, enabling integration backwards  innermost subsection are available, enabling integration backwards
627  in time over the last subsubsection.  in time along the last subsection.
628  Thus, the adjoint can be computed over this last  The adjoint can thus be computed along this last
629  subsubsection $k_{n}^{lev2}$.  subsection $k_{n}^{lev2}$.
630  %  %
631  \end{itemize}  \end{itemize}
632  %  %
633  This procedure is repeated consecutively for each previous  This procedure is repeated consecutively for each previous
634  subsubsection $k_{n-1}^{lev2}, \ldots, k_{1}^{lev2} $  subsection $k_{n-1}^{lev2}, \ldots, k_{1}^{lev2} $
635  carrying the adjoint computation to the initial time  carrying the adjoint computation to the initial time
636  of the subsection $k_{n}^{lev3}$.  of the subsection $k_{n}^{lev3}$.
637  Then, the procedure is repeated for the previous subsection  Then, the procedure is repeated for the previous subsection
# Line 607  $k_{1}^{lev3}$. Line 642  $k_{1}^{lev3}$.
642  For the full model trajectory of  For the full model trajectory of
643  $ n^{lev3} \cdot n^{lev2} \cdot n^{lev1} $ timesteps  $ n^{lev3} \cdot n^{lev2} \cdot n^{lev1} $ timesteps
644  the required storing of the model state was significantly reduced to  the required storing of the model state was significantly reduced to
645  $ n^{lev1} + n^{lev2} + n^{lev3} $  $ n^{lev2} + n^{lev3} $ to disk and roughly $ n^{lev1} $ to memory
646  [i.e. for the 3-day integration with a total oof 72 timesteps  [i.e. for the 3-day integration with a total oof 72 timesteps
647  the model state was stored 13 times].  the model state was stored 7 times to disk and roughly 6 times
648    to memory].
649  This saving in memory comes at a cost of a required  This saving in memory comes at a cost of a required
650  3 full forward integrations of the model (one for each  3 full forward integrations of the model (one for each
651  checkpointing level).  checkpointing level).
652  The balance of storage vs. recomputation certainly depends  The optimal balance of storage vs. recomputation certainly depends
653  on the computing resources available.  on the computing resources available and may be adjusted by
654    adjusting the partitioning among the
655    $ n^{lev3}, \,\, n^{lev2}, \,\, n^{lev1} $.
656    
657  \begin{figure}[t!]  \begin{figure}[t!]
658  \centering  \begin{center}
659  %\psdraft  %\psdraft
660  \psfrag{v_k1^lev3}{\mathinfigure{v_{k_{1}^{lev3}}}}  %\psfrag{v_k1^lev3}{\mathinfigure{v_{k_{1}^{lev3}}}}
661  \psfrag{v_kn-1^lev3}{\mathinfigure{v_{k_{n-1}^{lev3}}}}  %\psfrag{v_kn-1^lev3}{\mathinfigure{v_{k_{n-1}^{lev3}}}}
662  \psfrag{v_kn^lev3}{\mathinfigure{v_{k_{n}^{lev3}}}}  %\psfrag{v_kn^lev3}{\mathinfigure{v_{k_{n}^{lev3}}}}
663  \psfrag{v_k1^lev2}{\mathinfigure{v_{k_{1}^{lev2}}}}  %\psfrag{v_k1^lev2}{\mathinfigure{v_{k_{1}^{lev2}}}}
664  \psfrag{v_kn-1^lev2}{\mathinfigure{v_{k_{n-1}^{lev2}}}}  %\psfrag{v_kn-1^lev2}{\mathinfigure{v_{k_{n-1}^{lev2}}}}
665  \psfrag{v_kn^lev2}{\mathinfigure{v_{k_{n}^{lev2}}}}  %\psfrag{v_kn^lev2}{\mathinfigure{v_{k_{n}^{lev2}}}}
666  \psfrag{v_k1^lev1}{\mathinfigure{v_{k_{1}^{lev1}}}}  %\psfrag{v_k1^lev1}{\mathinfigure{v_{k_{1}^{lev1}}}}
667  \psfrag{v_kn^lev1}{\mathinfigure{v_{k_{n}^{lev1}}}}  %\psfrag{v_kn^lev1}{\mathinfigure{v_{k_{n}^{lev1}}}}
668  \mbox{\epsfig{file=part5/checkpointing.eps, width=0.8\textwidth}}  %\mbox{\epsfig{file=s_autodiff/figs/checkpointing.eps, width=0.8\textwidth}}
669    \resizebox{5.5in}{!}{\includegraphics{s_autodiff/figs/checkpointing.eps}}
670  %\psfull  %\psfull
671  \caption  \end{center}
672  {Schematic view of intermediate dump and restart for  \caption{
673    Schematic view of intermediate dump and restart for
674  3-level checkpointing.}  3-level checkpointing.}
675  \label{fig:erswns}  \label{fig:3levelcheck}
676  \end{figure}  \end{figure}
677    
678  \subsection{Optimal perturbations}  % \subsection{Optimal perturbations}
679  \label{optpert}  % \label{sec_optpert}
680    
681    
682  \subsection{Error covariance estimate and Hessian matrix}  % \subsection{Error covariance estimate and Hessian matrix}
683  \label{sec_hessian}  % \label{sec_hessian}
684    
685  \newpage  \newpage
686    
687  %**********************************************************************  %**********************************************************************
688  \section{AD-specific setup by example: sensitivity of carbon sequestration}  \section{TLM and ADM generation in general}
689  \label{sec_ad_setup_ex}  \label{sec_ad_setup_gen}
690    \begin{rawhtml}
691    <!-- CMIREDIR:sec_ad_setup_gen: -->
692    \end{rawhtml}
693  %**********************************************************************  %**********************************************************************
694    
695  The MITGCM has been adapted to enable AD using TAMC or TAF  In this section we describe in a general fashion
696  (we'll refer to TAMC and TAF interchangeably, except where  the parts of the code that are relevant for automatic
697  distinctions are explicitly mentioned).  differentiation using the software tool TAF.
698  The present description, therefore, is specific to the  Modifications to use OpenAD are described in \ref{sec_ad_openad}.
699  use of TAMC as AD tool.  
700  The following sections describe the steps which are necessary to  \input{s_autodiff/text/doc_ad_the_model}
701  generate a tangent linear or adjoint model of the MITGCM.  
702  We take as an example the sensitivity of carbon sequestration  The basic flow is depicted in \ref{fig:adthemodel}.
703  in the ocean.  If CPP option \texttt{ALLOW\_AUTODIFF\_TAMC} is defined,
704  The AD-relevant hooks in the code are sketched in  the driver routine
705  \reffig{adthemodel}, \reffig{adthemain}.  {\it the\_model\_main}, instead of calling {\it the\_main\_loop},
706    invokes the adjoint of this routine, {\it adthe\_main\_loop}
707  \subsection{Overview of the experiment}  (case \texttt{\#define ALLOW\_ADJOINT\_RUN}), or
708    the tangent linear of this routine {\it g\_the\_main\_loop}
709  We describe an adjoint sensitivity analysis of outgassing from  (case \texttt{\#define ALLOW\_TANGENTLINEAR\_RUN}),
710  the ocean into the atmosphere of a carbon like tracer injected  which are the toplevel routines in terms of automatic differentiation.
711  into the ocean interior (see \cite{hil-eta:01}).  The routines {\it adthe\_main\_loop} or {\it g\_the\_main\_loop}
712    are generated by TAF.
713  \subsubsection{Passive tracer equation}  It contains both the forward integration of the full model, the
714    cost function calculation,
715  For this work the MITGCM was augmented with a thermodynamically  any additional storing that is required for efficient checkpointing,
716  inactive tracer, $C$. Tracer residing in the ocean  and the reverse integration of the adjoint model.
717  model surface layer is outgassed according to a relaxation time scale,  
718  $\mu$. Within the ocean interior, the tracer is passively advected  [DESCRIBE IN A SEPARATE SECTION THE WORKING OF THE TLM]
719  by the ocean model currents. The full equation for the time evolution  
720  %  In Fig. \ref{fig:adthemodel}
721  \begin{equation}  the structure of {\it adthe\_main\_loop} has been strongly
722  \label{carbon_ddt}  simplified to focus on the essentials; in particular, no checkpointing
723  \frac{\partial C}{\partial t} \, = \,  procedures are shown here.
724  -U\cdot \nabla C \, - \, \mu C \, + \, \Gamma(C) \,+ \, S  Prior to the call of {\it adthe\_main\_loop}, the routine
725  \end{equation}  {\it ctrl\_unpack} is invoked to unpack the control vector
726  %  or initialise the control variables.
727  also includes a source term $S$. This term  Following the call of {\it adthe\_main\_loop},
728  represents interior sources of $C$ such as would arise due to  the routine {\it ctrl\_pack}
729  direct injection.  is invoked to pack the control vector
730  The velocity term, $U$, is the sum of the  (cf. Section \ref{section_ctrl}).
731  model Eulerian circulation and an eddy-induced velocity, the latter  If gradient checks are to be performed, the option
732  parameterized according to Gent/McWilliams (\cite{gen:90, dan:95}).  {\tt ALLOW\_GRADIENT\_CHECK} is defined. In this case
733  The convection function, $\Gamma$, mixes $C$ vertically wherever the  the driver routine {\it grdchk\_main} is called after
734  fluid is locally statically unstable.  the gradient has been computed via the adjoint
735    (cf. Section \ref{sec:ad_gradient_check}).
736  The outgassing time scale, $\mu$, in eqn. (\ref{carbon_ddt})  
737  is set so that \( 1/\mu \sim 1 \ \mathrm{year} \) for the surface  %------------------------------------------------------------------
738  ocean and $\mu=0$ elsewhere. With this value, eqn. (\ref{carbon_ddt})  
739  is valid as a prognostic equation for small perturbations in oceanic  \subsection{General setup
740  carbon concentrations. This configuration provides a  \label{section_ad_setup}}
741  powerful tool for examining the impact of large-scale ocean circulation  
742  on $ CO_2 $ outgassing due to interior injections.  In order to configure AD-related setups the following packages need
743  As source we choose a constant in time injection of  to be enabled:
744  $ S = 1 \,\, {\rm mol / s}$.  {\it
745    \begin{table}[!ht]
746  \subsubsection{Model configuration}  \begin{tabular}{l}
747    autodiff \\
748  The model configuration employed has a constant  ctrl \\
749  $4^\circ \times 4^\circ$ resolution horizontal grid and realistic  cost \\
750  geography and bathymetry. Twenty vertical layers are used with  grdchk \\
751  vertical spacing ranging  \end{tabular}
752  from 50 m near the surface to 815 m at depth.  \end{table}
753  Driven to steady-state by climatalogical wind-stress, heat and  }
754  fresh-water forcing the model reproduces well known large-scale  The packages are enabled by adding them to your experiment-specific
755  features of the ocean general circulation.  configuration file
756    {\it packages.conf} (see Section ???).
757  \subsubsection{Outgassing cost function}  
758    The following AD-specific CPP option files need to be customized:
 To quantify and understand outgassing due to injections of $C$  
 in eqn. (\ref{carbon_ddt}),  
 we define a cost function $ {\cal J} $ that measures the total amount of  
 tracer outgassed at each timestep:  
 %  
 \begin{equation}  
 \label{cost_tracer}  
 {\cal J}(t=T)=\int_{t=0}^{t=T}\int_{A} \mu C \, dA \, dt  
 \end{equation}  
 %  
 Equation(\ref{cost_tracer}) integrates the outgassing term, $\mu C$,  
 from (\ref{carbon_ddt})  
 over the entire ocean surface area, $A$, and accumulates it  
 up to time $T$.  
 Physically, ${\cal J}$ can be thought of as representing the amount of  
 $CO_2$ that our model predicts would be outgassed following an  
 injection at rate $S$.  
 The sensitivity of ${\cal J}$ to the spatial location of $S$,  
 $\frac{\partial {\cal J}}{\partial S}$,  
 can be used to identify regions from which circulation  
 would cause $CO_2$ to rapidly outgas following injection  
 and regions in which $CO_2$ injections would remain effectively  
 sequesterd within the ocean.  
   
 \subsection{Code configuration}  
   
 The model configuration for this experiment resides under the  
 directory {\it verification/carbon/}.  
 The code customisation routines are in {\it verification/carbon/code/}:  
759  %  %
760  \begin{itemize}  \begin{itemize}
761  %  %
762  \item {\it .genmakerc}  \item {\it ECCO\_CPPOPTIONS.h} \\
763  %  This header file collects CPP options for the packages
764  \item {\it COST\_CPPOPTIONS.h}  {\it autodiff, cost, ctrl} as well as AD-unrelated options for
765  %  the external forcing package {\it exf}.
766  \item {\it CPP\_EEOPTIONS.h}  \footnote{NOTE: These options are not set in their package-specific
767  %  headers such as {\it COST\_CPPOPTIONS.h}, but are instead collected
768  \item {\it CPP\_OPTIONS.h}  in the single header file {\it ECCO\_CPPOPTIONS.h}.
769  %  The package-specific header files serve as simple
770  \item {\it CTRL\_OPTIONS.h}  placeholders at this point.}
771  %  %
772  \item {\it ECCO\_OPTIONS.h}  \item {\it tamc.h} \\
773  %  This header configures the splitting of the time stepping loop
774  \item {\it SIZE.h}  w.r.t. the 3-level checkpointing (see section ???).
775  %  
 \item {\it adcommon.h}  
 %  
 \item {\it tamc.h}  
776  %  %
777  \end{itemize}  \end{itemize}
778    
779    %------------------------------------------------------------------
780    
781    \subsection{Building the AD code using TAF
782    \label{section_ad_build}}
783    
784    The build process of an AD code is very similar to building
785    the forward model. However, depending on which AD code one wishes
786    to generate, and on which AD tool is available (TAF or TAMC),
787    the following {\tt make} targets are available:
788    
789    \begin{table}[!ht]
790    {\footnotesize
791    \begin{tabular}{|ccll|}
792    \hline
793    ~ & {\it AD-target} & {\it output} & {\it description} \\
794    \hline
795    \hline
796    (1) & {\tt <MODE><TOOL>only} & {\tt <MODE>\_<TOOL>\_output.f}  &
797    generates code for $<$MODE$>$ using $<$TOOL$>$ \\
798    ~ & ~ & ~ & no {\tt make} dependencies on {\tt .F .h} \\
799    ~ & ~ & ~ & useful for compiling on remote platforms \\
800    \hline
801    (2) & {\tt <MODE><TOOL>} & {\tt <MODE>\_<TOOL>\_output.f}  &
802    generates code for $<$MODE$>$ using $<$TOOL$>$ \\
803    ~ & ~ & ~ & includes {\tt make} dependencies on {\tt .F .h} \\
804    ~ & ~ & ~ & i.e. input for $<$TOOL$>$ may be re-generated \\
805    \hline
806    (3) & {\tt <MODE>all} & {\tt mitgcmuv\_<MODE>}  &
807    generates code for $<$MODE$>$ using $<$TOOL$>$ \\
808    ~ & ~ & ~ & and compiles all code \\
809    ~ & ~ & ~ & (use of TAF is set as default) \\
810    \hline
811    \end{tabular}
812    }
813    \end{table}
814  %  %
815  The runtime flag and parameters settings are contained in  Here, the following placeholders are used
 {\it verification/carbon/input/},  
 together with the forcing fields and and restart files:  
816  %  %
817  \begin{itemize}  \begin{itemize}
818  %  %
819  \item {\it data}  \item $<$TOOL$>$
820  %  %
821  \item {\it data.cost}  \begin{itemize}
 %  
 \item {\it data.ctrl}  
 %  
 \item {\it data.pkg}  
 %  
 \item {\it eedata}  
 %  
 \item {\it topog.bin}  
 %  
 \item {\it windx.bin, windy.bin}  
 %  
 \item {\it salt.bin, theta.bin}  
 %  
 \item {\it SSS.bin, SST.bin}  
822  %  %
823  \item {\it pickup*}  \item {\tt TAF}
824    \item {\tt TAMC}
825  %  %
826  \end{itemize}  \end{itemize}
827  %  %
828  Finally, the file to generate the adjoint code resides in  \item $<$MODE$>$
 $ adjoint/ $:  
829  %  %
830  \begin{itemize}  \begin{itemize}
831  %  %
832  \item {\it makefile}  \item {\tt ad} generates the adjoint model (ADM)
833    \item {\tt ftl} generates the tangent linear model (TLM)
834    \item {\tt svd} generates both ADM and TLM for \\
835    singular value decomposition (SVD) type calculations
836  %  %
837  \end{itemize}  \end{itemize}
838  %  %
839    \end{itemize}
840    
841  Below we describe the customisations of this files which are  For example, to generate the adjoint model using TAF after routines ({\tt .F})
842  specific to this experiment.  or headers ({\tt .h}) have been modified, but without compilation,
843    type {\tt make adtaf};
844  \subsubsection{File {\it .genmakerc}}  or, to generate the tangent linear model using TAMC without
845  This file overwites default settings of {\it genmake}.  re-generating the input code, type {\tt make ftltamconly}.
 In the present example it is used to switch on the following  
 packages which are related to automatic differentiation  
 and are disabled by default: \\  
 \hspace*{4ex} {\tt set ENABLE=( autodiff cost ctrl ecco )}  \\  
 Other packages which are not needed are switched off: \\  
 \hspace*{4ex} {\tt set DISABLE=( aim obcs zonal\_filt shap\_filt cal exf )}  
   
 \subsubsection{File {\it COST\_CPPOPTIONS.h,  CTRL\_OPTIONS.h}}  
   
 These files used to contain package-specific CPP-options  
 (see Section \ref{???}).  
 For technical reasons those options have been grouped together  
 in the file {\it ECCO\_OPTIONS.h}.  
 To retain the modularity, the files have been kept and contain  
 the standard include of the {\it CPP\_OPTIONS.h} file.  
   
 \subsubsection{File {\it CPP\_EEOPTIONS.h}}  
   
 This file contains 'wrapper'-specific CPP options.  
 It only needs to be changed if the code is to be run  
 in  parallel environment (see Section \ref{???}).  
   
 \subsubsection{File {\it CPP\_OPTIONS.h}}  
   
 This file contains model-specific CPP options  
 (see Section \ref{???}).  
 Most options are related to the forward model setup.  
 They are identical to the global steady circulation setup of  
 {\it verification/exp2/}.  
 The option specific to this experiment is \\  
 \hspace*{4ex} {\tt \#define ALLOW\_MIT\_ADJOINT\_RUN} \\  
 This flag enables the inclusion of some AD-related fields  
 concerning initialisation, link between control variables  
 and forward model variables, and the call to the top-level  
 forward/adjoint subroutine {\it adthe\_main\_loop}  
 instead of {\it the\_main\_loop}.  
   
 \subsubsection{File {\it ECCO\_OPTIONS.h}}  
846    
 The CPP options of several AD-related packages are grouped  
 in this file:  
 %  
 \begin{itemize}  
 %  
 \item  
 Adjoint support package: {\it pkg/autodiff/} \\  
 This package contains hand-written adjoint code such as  
 active file handling, flow directives for files which must not  
 be differentiated, and TAMC-specific header files. \\  
 \hspace*{4ex} {\tt \#define ALLOW\_AUTODIFF\_TAMC} \\  
 defines TAMC-related features in the code. \\  
 \hspace*{4ex} {\tt \#define ALLOW\_TAMC\_CHECKPOINTING} \\  
 enables the checkpointing feature of TAMC  
 (see Section \ref{???}).  
 In the present example a 3-level checkpointing is implemented.  
 The code contains the relevant store directives, common block  
 and tape initialisations, storing key computation,  
 and loop index handling.  
 The checkpointing length at each level is defined in  
 file {\it tamc.h}, cf. below.  
 %  
 \item Cost function package: {\it pkg/cost/} \\  
 This package contains all relevant routines for  
 initialising, accumulating and finalizing the cost function  
 (see Section \ref{???}). \\  
 \hspace*{4ex} {\tt \#define ALLOW\_COST} \\  
 enables all general aspects of the cost function handling,  
 in particular the hooks in the foorward code for  
 initialising, accumulating and finalizing the cost function. \\  
 \hspace*{4ex} {\tt \#define ALLOW\_COST\_TRACER} \\  
 includes the subroutine with the cost function for this  
 particular experiment, eqn. (\ref{cost_tracer}).  
 %  
 \item Control variable package: {\it pkg/ctrl/} \\  
 This package contains all relevant routines for  
 the handling of the control vector.  
 Each control variable can be enabled/disabled with its own flag: \\  
 \begin{tabular}{ll}  
 \hspace*{2ex} {\tt \#define ALLOW\_THETA0\_CONTROL} &  
 initial temperature \\  
 \hspace*{2ex} {\tt \#define ALLOW\_SALT0\_CONTROL} &  
 initial salinity \\  
 \hspace*{2ex} {\tt \#define ALLOW\_TR0\_CONTROL} &  
 initial passive tracer concentration \\  
 \hspace*{2ex} {\tt \#define ALLOW\_TAUU0\_CONTROL} &  
 zonal wind stress \\  
 \hspace*{2ex} {\tt \#define ALLOW\_TAUV0\_CONTROL} &  
 meridional wind stress \\  
 \hspace*{2ex} {\tt \#define ALLOW\_SFLUX0\_CONTROL} &  
 freshwater flux \\  
 \hspace*{2ex} {\tt \#define ALLOW\_HFLUX0\_CONTROL} &  
 heat flux \\  
 \hspace*{2ex} {\tt \#undef ALLOW\_DIFFKR\_CONTROL} &  
 diapycnal diffusivity \\  
 \hspace*{2ex} {\tt \#undef ALLOW\_KAPPAGM\_CONTROL} &  
 isopycnal diffusivity \\  
 \end{tabular}  
 %  
 \end{itemize}  
847    
848  \subsubsection{File {\it SIZE.h}}  A typical full build process to generate the ADM via TAF would
849    look like follows:
850    \begin{verbatim}
851    % mkdir build
852    % cd build
853    % ../../../tools/genmake2 -mods=../code_ad
854    % make depend
855    % make adall
856    \end{verbatim}
857    
858  The file contains the grid point dimensions of the forward  %------------------------------------------------------------------
 model. It is identical to the {\it verification/exp2/}: \\  
 \hspace*{4ex} {\tt sNx = 90} \\  
 \hspace*{4ex} {\tt sNy = 40} \\  
 \hspace*{4ex} {\tt Nr = 20} \\  
 It correpsponds to a single-tile/single-processor setup:  
 {\tt nSx = nSy = 1, nPx = nPy = 1},  
 with standard overlap dimensioning  
 {\tt OLx = OLy = 3}.  
   
 \subsubsection{File {\it adcommon.h}}  
   
 This file contains common blocks of some adjoint variables  
 that are generated by TAMC.  
 The common blocks are used by the adjoint support routine  
 {\it addummy\_in\_stepping} which needs to access those variables:  
   
 \begin{tabular}{ll}  
 \hspace*{4ex} {\tt common /addynvars\_r/} &  
 \hspace*{4ex} is related to {\it DYNVARS.h} \\  
 \hspace*{4ex} {\tt common /addynvars\_cd/} &  
 \hspace*{4ex} is related to {\it DYNVARS.h} \\  
 \hspace*{4ex} {\tt common /adtr1\_r/} &  
 \hspace*{4ex} is related to {\it TR1.h} \\  
 \hspace*{4ex} {\tt common /adffields/} &  
 \hspace*{4ex} is related to {\it FFIELDS.h}\\  
 \end{tabular}  
859    
860  Note that if the structure of the common block changes in the  \subsection{The AD build process in detail
861  above header files of the forward code, the structure  \label{section_ad_build_detail}}
 of the adjoint common blocks will change accordingly.  
 Thus, it has to be made sure that the structure of the  
 adjoint common block in the hand-written file {\it adcommon.h}  
 complies with the automatically generated adjoint common blocks  
 in {\it adjoint\_model.F}.  
862    
863  \subsubsection{File {\it tamc.h}}  The {\tt make <MODE>all} target consists of the following procedures:
864    
865  This routine contains the dimensions for TAMC checkpointing.  \begin{enumerate}
866  %  %
867    \item
868    A header file {\tt AD\_CONFIG.h} is generated which contains a CPP option
869    on which code ought to be generated. Depending on the {\tt make} target,
870    the contents is one of the following:
871  \begin{itemize}  \begin{itemize}
872    \item
873    {\tt \#define ALLOW\_ADJOINT\_RUN}
874    \item
875    {\tt \#define ALLOW\_TANGENTLINEAR\_RUN}
876    \item
877    {\tt \#define ALLOW\_ECCO\_OPTIMIZATION}
878    \end{itemize}
879  %  %
880  \item {\tt \#ifdef ALLOW\_TAMC\_CHECKPOINTING} \\  \item
881  3-level checkpointing is enabled, i.e. the timestepping  A single file {\tt <MODE>\_input\_code.f} is concatenated
882  is divided into three different levels (see Section \ref{???}).  consisting of all {\tt .f} files that are part of the list {\bf AD\_FILES}
883  The model state of the outermost ({\tt nchklev\_3}) and the  and all {\tt .flow} files that are part of the list {\bf AD\_FLOW\_FILES}.
884  itermediate ({\tt nchklev\_2}) timestepping loop are stored to file  %
885  (handled in {\it the\_main\_loop}).  \item
886  The innermost loop ({\tt nchklev\_1})  The AD tool is invoked with the {\tt <MODE>\_<TOOL>\_FLAGS}.
887  avoids I/O by storing all required variables  The default AD tool flags in {\tt genmake2} can be overrwritten by
888  to common blocks. This storing may also be necessary if  an {\tt adjoint\_options} file (similar to the platform-specific
889  no checkpointing is chosen  {\tt build\_options}, see Section ???.
890  (nonlinear functions, if-statements, iterative loops, ...).  The AD tool writes the resulting AD code into the file
891  In the present example the dimensions are chosen as follows: \\  {\tt <MODE>\_input\_code\_ad.f}
892  \hspace*{4ex} {\tt nchklev\_1      =  36 } \\  %
893  \hspace*{4ex} {\tt nchklev\_2      =  30 } \\  \item
894  \hspace*{4ex} {\tt nchklev\_3      =  60 } \\  A short sed script {\tt adjoint\_sed} is applied to
895  To guarantee that the checkpointing intervals span the entire  {\tt <MODE>\_input\_code\_ad.f}
896  integration period the relation \\  to reinstate {\bf myThid} into the CALL argument list of active file I/O.
897  \hspace*{4ex} {\tt nchklev\_1*nchklev\_2*nchklev\_3 $ \ge $ nTimeSteps} \\  The result is written to file {\tt <MODE>\_<TOOL>\_output.f}.
898  where {\tt nTimeSteps} is either specified in {\it data}  %
899  or computed via \\  \item
900  \hspace*{4ex} {\tt nTimeSteps = (endTime-startTime)/deltaTClock }.  All routines are compiled and an executable is generated
901  %  (see Table ???).
 \item {\tt \#undef ALLOW\_TAMC\_CHECKPOINTING} \\  
 No checkpointing is enabled.  
 In this case the relevant counter is {\tt nchklev\_0}.  
 Similar to above, the following relation has to be satisfied \\  
 \hspace*{4ex} {\tt nchklev\_0 $ \ge $ nTimeSteps}.  
902  %  %
903  \end{itemize}  \end{enumerate}
904    
905  \subsubsection{File {\it makefile}}  \subsubsection{The list AD\_FILES and {\tt .list} files}
906    
907  This file contains all relevant paramter flags and  Not all routines are presented to the AD tool.
908  lists to run TAMC.  Routines typically hidden are diagnostics routines which
909  It is assumed that TAMC is available to you, either locally,  do not influence the cost function, but may create
910  being installed on your network, or remotely through the 'TAMC Utility'.  artificial flow dependencies such as I/O of active variables.
911  TAMC is called with the command {\tt tamc} followed by a  
912  number of options. They are described in detail in the  {\tt genmake2} generates a list (or variable) {\bf AD\_FILES}
913  TAMC manual \cite{gie:99}.  which contains all routines that are shown to the AD tool.
914  Here we briefly discuss the main flags used in the {\it makefile}  This list is put together from all files with suffix {\tt .list}
915    that {\tt genmake2} finds in its search directories.
916    The list file for the core MITgcm routines is in {\tt model/src/}
917    is called {\tt model\_ad\_diff.list}.
918    Note that no wrapper routine is shown to TAF. These are either
919    not visible at all to the AD code, or hand-written AD code
920    is available (see next section).
921    
922    Each package directory contains its package-specific
923    list file {\tt <PKG>\_ad\_diff.list}. For example,
924    {\tt pkg/ptracers/} contains the file {\tt ptracers\_ad\_diff.list}.
925    Thus, enabling a package will automatically extend the
926    {\bf AD\_FILES} list of {\tt genmake2} to incorporate the
927    package-specific routines.
928    Note that you will need to regenerate the {\tt Makefile} if
929    you enable a package (e.g. by adding it to {\tt packages.conf})
930    and a {\tt Makefile} already exists.
931    
932    \subsubsection{The list AD\_FLOW\_FILES and {\tt .flow} files}
933    
934    TAMC and TAF can evaluate user-specified directives
935    that start with a specific syntax ({\tt CADJ}, {\tt C\$TAF}, {\tt !\$TAF}).
936    The main categories of directives are STORE directives and
937    FLOW directives. Here, we are concerned with flow directives,
938    store directives are treated elsewhere.
939    
940    Flow directives enable the AD tool to evaluate how it should treat
941    routines that are 'hidden' by the user, i.e. routines which are
942    not contained in the {\bf AD\_FILES} list (see previous section),
943    but which are called in part of the code that the AD tool does see.
944    The flow directive tell the AD tool
945  %  %
946  \begin{itemize}  \begin{itemize}
 \item [{\tt tamc}] {\tt  
 -input <variable names>  
 -output <variable name> ... \\  
 -toplevel <S/R name> -reverse <file names>  
 }  
 \end{itemize}  
947  %  %
948  \begin{itemize}  \item which subroutine arguments are input/output
949  %  \item which subroutine arguments are active
950  \item {\tt -toplevel <S/R name>} \\  \item which subroutine arguments are required to compute the cost
951  Name of the toplevel routine, with respect to which the  \item which subroutine arguments are dependent
 control flow analysis is performed.  
 %  
 \item {\tt -input <variable names>} \\  
 List of independent variables $ u $ with respect to which the  
 dependent variable $ J $ is differentiated.  
 %  
 \item {\tt -output <variable name>} \\  
 Dependent variable $ J $  which is to be differentiated.  
 %  
 \item {\tt -reverse <file names>} \\  
 Adjoint code is generated to compute the sensitivity of an  
 independent variable w.r.t.  many dependent variables.  
 The generated adjoint top-level routine computes the product  
 of the transposed Jacobian matrix $ M^T $ times  
 the gradient vector $ \nabla_v J $.  
 \\  
 {\tt <file names>} refers to the list of files {\it .f} which are to be  
 analyzed by TAMC. This list is generally smaller than the full list  
 of code to be compiled. The files not contained are either  
 above the top-level routine (some initialisations), or are  
 deliberately hidden from TAMC, either because hand-written  
 adjoint routines exist, or the routines must not (or don't have to)  
 be differentiated. For each routine which is part of the flow tree  
 of the top-level routine, but deliberately hidden from TAMC,  
 a corresponding file {\it .flow} exists containing flow directives  
 for TAMC.  
952  %  %
953  \end{itemize}  \end{itemize}
954    %
955    The syntax for the flow directives can be found in the
956    AD tool manuals.
957    
958    {\tt genmake2} generates a list (or variable) {\bf AD\_FLOW\_FILES}
959    which contains all files with suffix{\tt .flow} that it finds
960    in its search directories.
961    The flow directives for the core MITgcm routines of
962    {\tt eesupp/src/} and {\tt model/src/}
963    reside in {\tt pkg/autodiff/}.
964    This directory also contains hand-written adjoint code
965    for the MITgcm WRAPPER (section \ref{chap:sarch}).
966    
967    Flow directives for package-specific routines are contained in
968    the corresponding package directories in the file
969    {\tt <PKG>\_ad.flow}, e.g. ptracers-specific directives are in
970    {\tt ptracers\_ad.flow}.
971    
972    \subsubsection{Store directives for 3-level checkpointing}
973    
974    The storing that is required at each period of the
975    3-level checkpointing is controled by three
976    top-level headers.
977    
978  \subsubsection{File {\it data}}  \begin{verbatim}
979    do ilev_3 = 1, nchklev_3
980  \subsubsection{File {\it data.cost}}  #  include ``checkpoint_lev3.h''
981       do ilev_2 = 1, nchklev_2
982  \subsubsection{File {\it data.ctrl}}  #     include ``checkpoint_lev2.h''
983          do ilev_1 = 1, nchklev_1
984  \subsubsection{File {\it data.pkg}}  #        include ``checkpoint_lev1.h''
985    
986  \subsubsection{File {\it eedata}}  ...
987    
988  \subsubsection{File {\it topog.bin}}        end do
989       end do
990  \subsubsection{File {\it windx.bin, windy.bin}}  end do
991    \end{verbatim}
 \subsubsection{File {\it salt.bin, theta.bin}}  
992    
993  \subsubsection{File {\it SSS.bin, SST.bin}}  All files {\tt checkpoint\_lev?.h} are contained in directory
994    {\tt pkg/autodiff/}.
995    
 \subsubsection{File {\it pickup*}}  
996    
997  \subsection{Compiling the model and its adjoint}  \subsubsection{Changing the default AD tool flags: ad\_options files}
998    
 \newpage  
999    
1000  %**********************************************************************  \subsubsection{Hand-written adjoint code}
 \section{TLM and ADM code generation in general}  
 \label{sec_ad_setup_gen}  
 %**********************************************************************  
1001    
1002  In this section we describe in a general fashion  %------------------------------------------------------------------
 the parts of the code that are relevant for automatic  
 differentiation using the software tool TAMC.  
1003    
1004  \subsection{The cost function (dependent variable)}  \subsection{The cost function (dependent variable)
1005    \label{section_cost}}
1006    
1007  The cost function $ {\cal J} $ is referred to as the {\sf dependent variable}.  The cost function $ {\cal J} $ is referred to as the {\sf dependent variable}.
1008  It is a function of the input variables $ \vec{u} $ via the composition  It is a function of the input variables $ \vec{u} $ via the composition
1009  $ {\cal J}(\vec{u}) \, = \, {\cal J}(M(\vec{u})) $.  $ {\cal J}(\vec{u}) \, = \, {\cal J}(M(\vec{u})) $.
1010  The input is referred to as the  The input are referred to as the
1011  {\sf independent variables} or {\sf control variables}.  {\sf independent variables} or {\sf control variables}.
1012  All aspects relevant to the treatment of the cost function $ {\cal J} $  All aspects relevant to the treatment of the cost function $ {\cal J} $
1013  (parameter setting, initialisation, incrementation,  (parameter setting, initialization, accumulation,
1014  final evaluation), are controled by the package {\it pkg/cost}.  final evaluation), are controlled by the package {\it pkg/cost}.
1015    The aspects relevant to the treatment of the independent variables
1016    are controlled by the package {\it pkg/ctrl} and will be treated
1017    in the next section.
1018    
1019    \input{s_autodiff/text/doc_cost_flow}
1020    
1021    \subsubsection{Enabling the package}
1022    
 \subsubsection{genmake and CPP options}  
 %  
 \begin{itemize}  
 %  
 \item  
1023  \fbox{  \fbox{
1024  \begin{minipage}{12cm}  \begin{minipage}{12cm}
1025  {\it genmake}, {\it CPP\_OPTIONS.h}, {\it ECCO\_CPPOPTIONS.h}  {\it packages.conf}, {\it ECCO\_CPPOPTIONS.h}
1026  \end{minipage}  \end{minipage}
1027  }  }
1028  \end{itemize}  \begin{itemize}
 %  
 The directory {\it pkg/cost} can be included to the  
 compile list in 3 different ways (cf. Section \ref{???}):  
1029  %  %
1030  \begin{enumerate}  \item
1031    The package is enabled by adding {\it cost} to your file {\it packages.conf}
1032    (see Section ???)
1033  %  %
1034  \item {\it genmake}: \\  \item
1035  Change the default settngs in the file {\it genmake} by adding  
1036  {\bf cost} to the {\bf enable} list (not recommended).  
1037  %  \end{itemize}
 \item {\it .genmakerc}: \\  
 Customize the settings of {\bf enable}, {\bf disable} which are  
 appropriate for your experiment in the file {\it .genmakerc}  
 and add the file to your compile directory.  
 %  
 \item genmake-options: \\  
 Call {\it genmake} with the option  
 {\tt genmake -enable=cost}.  
1038  %  %
1039  \end{enumerate}  
1040  Since the cost function is usually used in conjunction with  N.B.: In general the following packages ought to be enabled
1041  automatic differentiation, the CPP option  simultaneously: {\it autodiff, cost, ctrl}.
 {\bf ALLOW\_ADJOINT\_RUN} should be defined  
 (file {\it CPP\_OPTIONS.h}).  
1042  The basic CPP option to enable the cost function is {\bf ALLOW\_COST}.  The basic CPP option to enable the cost function is {\bf ALLOW\_COST}.
1043  Each specific cost function contribution has its own option.  Each specific cost function contribution has its own option.
1044  For the present example the option is {\bf ALLOW\_COST\_TRACER}.  For the present example the option is {\bf ALLOW\_COST\_TRACER}.
1045  All cost-specific options are set in {\it ECCO\_CPPOPTIONS.h}  All cost-specific options are set in {\it ECCO\_CPPOPTIONS.h}
1046    Since the cost function is usually used in conjunction with
1047    automatic differentiation, the CPP option
1048    {\bf ALLOW\_ADJOINT\_RUN} (file {\it CPP\_OPTIONS.h}) and
1049    {\bf ALLOW\_AUTODIFF\_TAMC} (file {\it ECCO\_CPPOPTIONS.h})
1050    should be defined.
1051    
1052  \subsubsection{Initialisation}  \subsubsection{Initialization}
1053  %  %
1054  The initialisation of the {\it cost} package is readily enabled  The initialization of the {\it cost} package is readily enabled
1055  as soon as the CPP option {\bf ALLOW\_ADJOINT\_RUN} is defined.  as soon as the CPP option {\bf ALLOW\_COST} is defined.
1056  %  %
1057  \begin{itemize}  \begin{itemize}
1058  %  %
# Line 1152  Variables: {\it cost\_init} Line 1082  Variables: {\it cost\_init}
1082  }  }
1083  \\  \\
1084  This S/R  This S/R
1085  initialises the different cost function contributions.  initializes the different cost function contributions.
1086  The contribtion for the present example is {\bf objf\_tracer}  The contribution for the present example is {\bf objf\_tracer}
1087  which is defined on each tile (bi,bj).  which is defined on each tile (bi,bj).
1088  %  %
1089  \end{itemize}  \end{itemize}
1090  %  %
1091  \subsubsection{Incrementation}  \subsubsection{Accumulation}
1092  %  %
1093  \begin{itemize}  \begin{itemize}
1094  %  %
# Line 1176  Within this 'driver' routine, S/R are ca Line 1106  Within this 'driver' routine, S/R are ca
1106  the chosen cost function contributions.  the chosen cost function contributions.
1107  In the present example ({\bf ALLOW\_COST\_TRACER}),  In the present example ({\bf ALLOW\_COST\_TRACER}),
1108  S/R {\it cost\_tracer} is called.  S/R {\it cost\_tracer} is called.
1109  It accumulates {\bf objf\_tracer} according to eqn. (\ref{???}).  It accumulates {\bf objf\_tracer} according to eqn. (\ref{ask_the_author:doc_ad_2}).
1110  %  %
1111  \subsubsection{Finalize all contributions}  \subsubsection{Finalize all contributions}
1112  %  %
# Line 1196  from each contribution and sums over all Line 1126  from each contribution and sums over all
1126  \begin{equation}  \begin{equation}
1127  {\cal J} \, = \,  {\cal J} \, = \,
1128  {\rm fc} \, = \,  {\rm fc} \, = \,
1129  {\rm mult\_tracer} \sum_{bi,\,bj}^{nSx,\,nSy}  {\rm mult\_tracer} \sum_{\text{global sum}} \sum_{bi,\,bj}^{nSx,\,nSy}
1130  {\rm objf\_tracer}(bi,bj) \, + \, ...  {\rm objf\_tracer}(bi,bj) \, + \, ...
1131  \end{equation}  \end{equation}
1132  %  %
1133  The total cost function {\bf fc} will be the  The total cost function {\bf fc} will be the
1134  'dependent' variable in the argument list for TAMC, i.e.  'dependent' variable in the argument list for TAF, i.e.
1135  \begin{verbatim}  \begin{verbatim}
1136  tamc -output 'fc' ...  taf -output 'fc' ...
1137  \end{verbatim}  \end{verbatim}
1138    
 \begin{figure}[t!]  
 \input{part5/doc_ad_the_model}  
 \label{fig:adthemodel}  
 \caption{~}  
 \end{figure}  
   
1139  %%%% \end{document}  %%%% \end{document}
1140    
1141  \begin{figure}  \input{s_autodiff/text/doc_ad_the_main}
 \input{part5/doc_ad_the_main}  
 \label{fig:adthemain}  
 \caption{~}  
 \end{figure}  
1142    
1143  \subsection{The control variables (independent variables)}  \subsection{The control variables (independent variables)
1144    \label{section_ctrl}}
1145    
1146  The control variables are a subset of the model input  The control variables are a subset of the model input
1147  (initial conditions, boundary conditions, model parameters).  (initial conditions, boundary conditions, model parameters).
1148  Here we identify them with the variable $ \vec{u} $.  Here we identify them with the variable $ \vec{u} $.
1149  All intermediate variables whose derivative w.r.t. control  All intermediate variables whose derivative w.r.t. control
1150  variables don't vanish are called {\sf active variables}.  variables do not vanish are called {\sf active variables}.
1151  All subroutines whose derivative w.r.t. the control variables  All subroutines whose derivative w.r.t. the control variables
1152  don't vanish are called {\sf active routines}.  don't vanish are called {\sf active routines}.
1153  Read and write operations from and to file can be viewed  Read and write operations from and to file can be viewed
# Line 1234  as variable assignments. Therefore, file Line 1155  as variable assignments. Therefore, file
1155  active variables are written and from which active variables  active variables are written and from which active variables
1156  are read are called {\sf active files}.  are read are called {\sf active files}.
1157  All aspects relevant to the treatment of the control variables  All aspects relevant to the treatment of the control variables
1158  (parameter setting, initialisation, perturbation)  (parameter setting, initialization, perturbation)
1159  are controled by the package {\it pkg/ctrl}.  are controlled by the package {\it pkg/ctrl}.
1160    
1161    \input{s_autodiff/text/doc_ctrl_flow}
1162    
1163  \subsubsection{genmake and CPP options}  \subsubsection{genmake and CPP options}
1164  %  %
# Line 1251  are controled by the package {\it pkg/ct Line 1174  are controled by the package {\it pkg/ct
1174  %  %
1175  To enable the directory to be included to the compile list,  To enable the directory to be included to the compile list,
1176  {\bf ctrl} has to be added to the {\bf enable} list in  {\bf ctrl} has to be added to the {\bf enable} list in
1177  {\it .genmakerc} (or {\it genmake} itself).  {\it .genmakerc} or in {\it genmake} itself (analogous to {\it cost}
1178    package, cf. previous section).
1179  Each control variable is enabled via its own CPP option  Each control variable is enabled via its own CPP option
1180  in {\it ECCO\_CPPOPTIONS.h}.  in {\it ECCO\_CPPOPTIONS.h}.
1181    
1182  \subsubsection{Initialisation}  \subsubsection{Initialization}
1183  %  %
1184  \begin{itemize}  \begin{itemize}
1185  %  %
# Line 1292  and their gradients: {\it ctrl\_unpack} Line 1216  and their gradients: {\it ctrl\_unpack}
1216  \\  \\
1217  %  %
1218  Two important issues related to the handling of the control  Two important issues related to the handling of the control
1219  variables in the MITGCM need to be addressed.  variables in MITgcm need to be addressed.
1220  First, in order to save memory, the control variable arrays  First, in order to save memory, the control variable arrays
1221  are not kept in memory, but rather read from file and added  are not kept in memory, but rather read from file and added
1222  to the initial (or first guess) fields.  to the initial fields during the model initialization phase.
1223  Similarly, the corresponding adjoint fields which represent  Similarly, the corresponding adjoint fields which represent
1224  the gradient of the cost function w.r.t. the control variables  the gradient of the cost function w.r.t. the control variables
1225  are written to to file.  are written to file at the end of the adjoint integration.
1226  Second, in addition to the files holding the 2-dim. and 3-dim.  Second, in addition to the files holding the 2-dim. and 3-dim.
1227  control variables and the gradient, a 1-dim. {\sf control vector}  control variables and the corresponding cost gradients,
1228    a 1-dim. {\sf control vector}
1229  and {\sf gradient vector} are written to file. They contain  and {\sf gradient vector} are written to file. They contain
1230  only the wet points of the control variables and the corresponding  only the wet points of the control variables and the corresponding
1231  gradient.  gradient.
1232  This leads to a significant data compression.  This leads to a significant data compression.
1233  Furthermore, the control and the gradient vector can be passed to a  Furthermore, an option is available
1234    ({\tt ALLOW\_NONDIMENSIONAL\_CONTROL\_IO}) to
1235    non-dimensionalise the control and gradient vector,
1236    which otherwise would contain different pieces of different
1237    magnitudes and units.
1238    Finally, the control and gradient vector can be passed to a
1239  minimization routine if an update of the control variables  minimization routine if an update of the control variables
1240  is sought as part of a minimization exercise.  is sought as part of a minimization exercise.
1241    
# Line 1316  and gradient are generated and initialis Line 1246  and gradient are generated and initialis
1246    
1247  \subsubsection{Perturbation of the independent variables}  \subsubsection{Perturbation of the independent variables}
1248  %  %
1249  The dependency chain for differentiation starts  The dependency flow for differentiation w.r.t. the controls
1250  with adding a perturbation onto the the input variable,  starts with adding a perturbation onto the input variable,
1251  thus defining the independent or control variables for TAMC.  thus defining the independent or control variables for TAF.
1252  Three classes of controls may be considered:  Three types of controls may be considered:
1253  %  %
1254  \begin{itemize}  \begin{itemize}
1255  %  %
# Line 1334  Three classes of controls may be conside Line 1264  Three classes of controls may be conside
1264  Consider as an example the initial tracer distribution  Consider as an example the initial tracer distribution
1265  {\bf tr1} as control variable.  {\bf tr1} as control variable.
1266  After {\bf tr1} has been initialised in  After {\bf tr1} has been initialised in
1267  {\it ini\_tr1} (dynamical variables including  {\it ini\_tr1} (dynamical variables such as
1268  temperature and salinity are initialised in {\it ini\_fields}),  temperature and salinity are initialised in {\it ini\_fields}),
1269  a perturbation anomaly is added to the field in S/R  a perturbation anomaly is added to the field in S/R
1270  {\it ctrl\_map\_ini}  {\it ctrl\_map\_ini}
1271  %  %
1272    %\begin{eqnarray}
1273  \begin{equation}  \begin{equation}
1274  \begin{split}  \begin{aligned}
1275  u         & = \, u_{[0]} \, + \, \Delta u \\  u         & = \, u_{[0]} \, + \, \Delta u \\
1276  {\bf tr1}(...) & = \, {\bf tr1_{ini}}(...) \, + \, {\bf xx\_tr1}(...)  {\bf tr1}(...) & = \, {\bf tr1_{ini}}(...) \, + \, {\bf xx\_tr1}(...)
1277  \label{perturb}  \label{perturb}
1278  \end{split}  \end{aligned}
1279  \end{equation}  \end{equation}
1280    %\end{eqnarray}
1281  %  %
1282  In principle {\bf xx\_tr1} is a 3-dim. global array  {\bf xx\_tr1} is a 3-dim. global array
1283  holding the perturbation. In the case of a simple  holding the perturbation. In the case of a simple
1284  sensitivity study this array is identical to zero.  sensitivity study this array is identical to zero.
1285  However, it's specification is essential since TAMC  However, it's specification is essential in the context
1286    of automatic differentiation since TAF
1287  treats the corresponding line in the code symbolically  treats the corresponding line in the code symbolically
1288  when determining the differentiation chain and its origin.  when determining the differentiation chain and its origin.
1289  Thus, the variable names are part of the argument list  Thus, the variable names are part of the argument list
1290  when calling TAMC:  when calling TAF:
1291  %  %
1292  \begin{verbatim}  \begin{verbatim}
1293  tamc -input 'xx_tr1 ...' ...  taf -input 'xx_tr1 ...' ...
1294  \end{verbatim}  \end{verbatim}
1295  %  %
1296  Now, as mentioned above, the MITGCM avoids maintaining  Now, as mentioned above, MITgcm avoids maintaining
1297  an array for each control variable by reading the  an array for each control variable by reading the
1298  perturbation to a temporary array from file.  perturbation to a temporary array from file.
1299  To ensure the symbolic link to be recognized by TAMC, a scalar  To ensure the symbolic link to be recognized by TAF, a scalar
1300  dummy variable {\bf xx\_tr1\_dummy} is introduced  dummy variable {\bf xx\_tr1\_dummy} is introduced
1301  and an 'active read' routine of the adjoint support  and an 'active read' routine of the adjoint support
1302  package {\it pkg/autodiff} is invoked.  package {\it pkg/autodiff} is invoked.
1303  The read-procedure is tagged with the variable  The read-procedure is tagged with the variable
1304  {\bf xx\_tr1\_dummy} enabbling TAMC to recognize the  {\bf xx\_tr1\_dummy} enabling TAF to recognize the
1305  initialisation of the perturbation.  initialization of the perturbation.
1306  The modified call of TAMC thus reads  The modified call of TAF thus reads
1307  %  %
1308  \begin{verbatim}  \begin{verbatim}
1309  tamc -input 'xx_tr1_dummy ...' ...  taf -input 'xx_tr1_dummy ...' ...
1310  \end{verbatim}  \end{verbatim}
1311  %  %
1312  and the modified operation to (\ref{perturb})  and the modified operation to (\ref{perturb})
# Line 1388  in the code takes on the form Line 1321  in the code takes on the form
1321  %  %
1322  Note, that reading an active variable corresponds  Note, that reading an active variable corresponds
1323  to a variable assignment. Its derivative corresponds  to a variable assignment. Its derivative corresponds
1324  to a write statement of the adjoint variable.  to a write statement of the adjoint variable, followed by
1325    a reset.
1326  The 'active file' routines have been designed  The 'active file' routines have been designed
1327  to support active read and corresponding active write  to support active read and corresponding adjoint active write
1328  operations.  operations (and vice versa).
1329  %  %
1330  \item  \item
1331  \fbox{  \fbox{
# Line 1408  with the symbolic perturbation taking pl Line 1342  with the symbolic perturbation taking pl
1342  Note however an important difference:  Note however an important difference:
1343  Since the boundary values are time dependent with a new  Since the boundary values are time dependent with a new
1344  forcing field applied at each time steps,  forcing field applied at each time steps,
1345  the general problem may be be thought of as  the general problem may be thought of as
1346  a new control variable at each time step, i.e.  a new control variable at each time step
1347    (or, if the perturbation is averaged over a certain period,
1348    at each $ N $ timesteps), i.e.
1349  \[  \[
1350  u_{\rm forcing} \, = \,  u_{\rm forcing} \, = \,
1351  \{ \, u_{\rm forcing} ( t_n ) \, \}_{  \{ \, u_{\rm forcing} ( t_n ) \, \}_{
# Line 1434  calendar ({\it cal}~) and external forci Line 1370  calendar ({\it cal}~) and external forci
1370  %  %
1371  This routine is not yet implemented, but would proceed  This routine is not yet implemented, but would proceed
1372  proceed along the same lines as the initial value sensitivity.  proceed along the same lines as the initial value sensitivity.
1373    The mixing parameters {\bf diffkr} and {\bf kapgm}
1374    are currently added as controls in {\it ctrl\_map\_ini.F}.
1375  %  %
1376  \end{itemize}  \end{itemize}
1377  %  %
1378    
1379  \subsubsection{Output of adjoint variables and gradient}  \subsubsection{Output of adjoint variables and gradient}
1380  %  %
1381  Two ways exist to generate output of adjoint fields.  Several ways exist to generate output of adjoint fields.
1382  %  %
1383  \begin{itemize}  \begin{itemize}
1384  %  %
1385  \item  \item
1386  \fbox{  \fbox{
1387  \begin{minipage}{12cm}  \begin{minipage}{12cm}
1388  {\it ctrl\_pack}:  {\it ctrl\_map\_ini, ctrl\_map\_forcing}:
1389  \end{minipage}  \end{minipage}
1390  }  }
1391  \\  \\
 At the end of the forward/adjoint integration, the S/R  
 {\it ctrl\_pack} is called which mirrors S/R {\it ctrl\_unpack}.  
 It writes the following files:  
 %  
1392  \begin{itemize}  \begin{itemize}
1393  %  %
1394  \item {\bf xx\_...}: the control variable fields  \item {\bf xx\_...}: the control variable fields \\
1395    Before the forward integration, the control
1396    variables are read from file {\bf xx\_ ...} and added to
1397    the model field.
1398  %  %
1399  \item {\bf adxx\_...}: the adjoint variable fields, i.e. the gradient  \item {\bf adxx\_...}: the adjoint variable fields, i.e. the gradient
1400  $ \nabla _{u}{\cal J} $ for each control variable,  $ \nabla _{u}{\cal J} $ for each control variable \\
1401    After the adjoint integration the corresponding adjoint
1402    variables are written to {\bf adxx\_ ...}.
1403    %
1404    \end{itemize}
1405    %
1406    \item
1407    \fbox{
1408    \begin{minipage}{12cm}
1409    {\it ctrl\_unpack, ctrl\_pack}:
1410    \end{minipage}
1411    }
1412    \\
1413  %  %
1414  \item {\bf vector\_ctrl}: the control vector  \begin{itemize}
1415  %  %
1416  \item {\bf vector\_grad}: the gradient vector  \item {\bf vector\_ctrl}: the control vector \\
1417    At the very beginning of the model initialization,
1418    the updated compressed control vector is read (or initialised)
1419    and distributed to 2-dim. and 3-dim. control variable fields.
1420    %
1421    \item {\bf vector\_grad}: the gradient vector \\
1422    At the very end of the adjoint integration,
1423    the 2-dim. and 3-dim. adjoint variables are read,
1424    compressed to a single vector and written to file.
1425  %  %
1426  \end{itemize}  \end{itemize}
1427  %  %
# Line 1476  $ \nabla _{u}{\cal J} $ for each control Line 1433  $ \nabla _{u}{\cal J} $ for each control
1433  }  }
1434  \\  \\
1435  In addition to writing the gradient at the end of the  In addition to writing the gradient at the end of the
1436  forward/adjoint integration, many more adjoint variables,  forward/adjoint integration, many more adjoint variables
1437  representing the Lagrange multipliers of the model state  of the model state
1438  w.r.t. the model state  at intermediate times can be written using S/R
 at different times can be written using S/R  
1439  {\it addummy\_in\_stepping}.  {\it addummy\_in\_stepping}.
1440  This routine is part of the adjoint support package  This routine is part of the adjoint support package
1441  {\it pkg/autodiff} (cf.f. below).  {\it pkg/autodiff} (cf.f. below).
1442    The procedure is enabled using via the CPP-option
1443    {\bf ALLOW\_AUTODIFF\_MONITOR} (file {\it ECCO\_CPPOPTIONS.h}).
1444  To be part of the adjoint code, the corresponding S/R  To be part of the adjoint code, the corresponding S/R
1445  {\it dummy\_in\_stepping} has to be called in the forward  {\it dummy\_in\_stepping} has to be called in the forward
1446  model (S/R {\it the\_main\_loop}) at the appropriate place.  model (S/R {\it the\_main\_loop}) at the appropriate place.
1447    The adjoint common blocks are extracted from the adjoint code
1448    via the header file {\it adcommon.h}.
1449    
1450  {\it dummy\_in\_stepping} is essentially empty,  {\it dummy\_in\_stepping} is essentially empty,
1451  the corresponding adjoint routine is hand-written rather  the corresponding adjoint routine is hand-written rather
# Line 1493  than generated automatically. Line 1453  than generated automatically.
1453  Appropriate flow directives ({\it dummy\_in\_stepping.flow})  Appropriate flow directives ({\it dummy\_in\_stepping.flow})
1454  ensure that TAMC does not automatically  ensure that TAMC does not automatically
1455  generate {\it addummy\_in\_stepping} by trying to differentiate  generate {\it addummy\_in\_stepping} by trying to differentiate
1456  {\it dummy\_in\_stepping}, but rather takes the hand-written routine.  {\it dummy\_in\_stepping}, but instead refers to
1457    the hand-written routine.
1458    
1459  {\it dummy\_in\_stepping} is called in the forward code  {\it dummy\_in\_stepping} is called in the forward code
1460  at the beginning of each  at the beginning of each
# Line 1503  each timestep in the adjoint calculation Line 1464  each timestep in the adjoint calculation
1464  {\it addynamics}.  {\it addynamics}.
1465    
1466  {\it addummy\_in\_stepping} includes the header files  {\it addummy\_in\_stepping} includes the header files
1467  {\it adffields.h, addynamics.h, adtr1.h}.  {\it adcommon.h}.
1468  These header files are also hand-written. They contain  This header file is also hand-written. It contains
1469  the common blocks {\bf /addynvars\_r/}, {\bf /addynvars\_cd/},  the common blocks
1470    {\bf /addynvars\_r/}, {\bf /addynvars\_cd/},
1471    {\bf /addynvars\_diffkr/}, {\bf /addynvars\_kapgm/},
1472  {\bf /adtr1\_r/}, {\bf /adffields/},  {\bf /adtr1\_r/}, {\bf /adffields/},
1473  which have been extracted from the adjoint code to enable  which have been extracted from the adjoint code to enable
1474  access to the adjoint variables.  access to the adjoint variables.
1475    
1476    {\bf WARNING:} If the structure of the common blocks
1477    {\bf /dynvars\_r/}, {\bf /dynvars\_cd/}, etc., changes
1478    similar changes will occur in the adjoint common blocks.
1479    Therefore, consistency between the TAMC-generated common blocks
1480    and those in {\it adcommon.h} have to be checked.
1481  %  %
1482  \end{itemize}  \end{itemize}
1483    
# Line 1523  The gradient $ \nabla _{u}{\cal J} |_{u_ Line 1492  The gradient $ \nabla _{u}{\cal J} |_{u_
1492  with the value of the cost function itself $ {\cal J}(u_{[k]}) $  with the value of the cost function itself $ {\cal J}(u_{[k]}) $
1493  at iteration step $ k $ serve  at iteration step $ k $ serve
1494  as input to a minimization routine (e.g. quasi-Newton method,  as input to a minimization routine (e.g. quasi-Newton method,
1495  conjugate gradient, ...) to compute an update in the  conjugate gradient, ... \cite{gil-lem:89})
1496    to compute an update in the
1497  control variable for iteration step $k+1$  control variable for iteration step $k+1$
1498  \[  \[
1499  u_{[k+1]} \, = \,  u_{[0]} \, + \, \Delta u_{[k+1]}  u_{[k+1]} \, = \,  u_{[0]} \, + \, \Delta u_{[k+1]}
# Line 1533  u_{[k+1]} \, = \,  u_{[0]} \, + \, \Delt Line 1503  u_{[k+1]} \, = \,  u_{[0]} \, + \, \Delt
1503  $ u_{[k+1]} $ then serves as input for a forward/adjoint run  $ u_{[k+1]} $ then serves as input for a forward/adjoint run
1504  to determine $ {\cal J} $ and $ \nabla _{u}{\cal J} $ at iteration step  to determine $ {\cal J} $ and $ \nabla _{u}{\cal J} $ at iteration step
1505  $ k+1 $.  $ k+1 $.
1506  Tab. \ref{???} sketches the flow between forward/adjoint model  Tab. \ref{ask_the_author:doc_ad_2} sketches the flow between forward/adjoint model
1507  and the minimization routine.  and the minimization routine.
1508    
1509    {\scriptsize
1510  \begin{eqnarray*}  \begin{eqnarray*}
 \footnotesize  
1511  \begin{array}{ccccc}  \begin{array}{ccccc}
1512  u_{[0]} \,\, ,  \,\, \Delta u_{[k]}    & ~ & ~ & ~ & ~ \\  u_{[0]} \,\, ,  \,\, \Delta u_{[k]}    & ~ & ~ & ~ & ~ \\
1513  {\Big\downarrow}  {\Big\downarrow}
# Line 1554  v_{[k]} = M \left( u_{[k]} \right) & Line 1524  v_{[k]} = M \left( u_{[k]} \right) &
1524  {\cal J}_{[k]} = {\cal J} \left( M \left( u_{[k]} \right) \right)} \\  {\cal J}_{[k]} = {\cal J} \left( M \left( u_{[k]} \right) \right)} \\
1525  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1526  \hline  \hline
1527    \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~}  \\
1528    \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{{\Big\downarrow}} \\
1529    \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~}  \\
1530  \hline  \hline
1531  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1532  \multicolumn{1}{|c}{  \multicolumn{1}{|c}{
1533  \nabla_u {\cal J}_{[k]} (\delta {\cal J}) =  \nabla_u {\cal J}_{[k]} (\delta {\cal J}) =
1534  T\!\!^{\ast} \cdot \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J})} &  T^{\ast} \cdot \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J})} &
1535  \stackrel{\bf adjoint}{\mathbf \longleftarrow} &  \stackrel{\bf adjoint}{\mathbf \longleftarrow} &
1536  ad \, v_{[k]} (\delta {\cal J}) =  ad \, v_{[k]} (\delta {\cal J}) =
1537  \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J}) &  \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J}) &
# Line 1567  ad \, v_{[k]} (\delta {\cal J}) = Line 1540  ad \, v_{[k]} (\delta {\cal J}) =
1540  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1541  \hline  \hline
1542   ~ & ~ & ~ & ~ & ~ \\   ~ & ~ & ~ & ~ & ~ \\
1543  ~ & ~ &  \hspace*{15ex}{\Bigg\downarrow}  
1544  {\cal J}_{[k]} \qquad {\Bigg\downarrow}  \qquad \nabla_u {\cal J}_{[k]}  \quad {\cal J}_{[k]}, \quad \nabla_u {\cal J}_{[k]}
1545   & ~ & ~ \\   & ~ & ~ & ~ & ~ \\
1546   ~ & ~ & ~ & ~ & ~ \\   ~ & ~ & ~ & ~ & ~ \\
1547  \hline  \hline
1548  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
# Line 1585  ad \, v_{[k]} (\delta {\cal J}) = Line 1558  ad \, v_{[k]} (\delta {\cal J}) =
1558   ~ & ~ & ~ & ~ & \Delta u_{[k+1]} \\   ~ & ~ & ~ & ~ & \Delta u_{[k+1]} \\
1559  \end{array}  \end{array}
1560  \end{eqnarray*}  \end{eqnarray*}
1561    }
1562    
1563  The routines {\it ctrl\_unpack} and {\it ctrl\_pack} provide  The routines {\it ctrl\_unpack} and {\it ctrl\_pack} provide
1564  the link between the model and the minimization routine.  the link between the model and the minimization routine.
1565  As described in Section \ref{???}  As described in Section \ref{ask_the_author:doc_ad_2}
1566  the {\it unpack} and {\it pack} routines read and write  the {\it unpack} and {\it pack} routines read and write
1567  control and gradient {\it vectors} which are compressed  control and gradient {\it vectors} which are compressed
1568  to contain only wet points, in addition to the full  to contain only wet points, in addition to the full
# Line 1597  The corresponding I/O flow looks as foll Line 1571  The corresponding I/O flow looks as foll
1571    
1572  \vspace*{0.5cm}  \vspace*{0.5cm}
1573    
1574    {\scriptsize
1575  \begin{tabular}{ccccc}  \begin{tabular}{ccccc}
1576  {\bf vector\_ctrl\_$<$k$>$ } & ~ & ~ & ~ & ~ \\  {\bf vector\_ctrl\_$<$k$>$ } & ~ & ~ & ~ & ~ \\
1577  {\big\downarrow}  & ~ & ~ & ~ & ~ \\  {\big\downarrow}  & ~ & ~ & ~ & ~ \\
# Line 1607  The corresponding I/O flow looks as foll Line 1582  The corresponding I/O flow looks as foll
1582  \cline{3-3}  \cline{3-3}
1583  \multicolumn{1}{l}{\bf xx\_theta0...$<$k$>$} & ~ &  \multicolumn{1}{l}{\bf xx\_theta0...$<$k$>$} & ~ &
1584  \multicolumn{1}{|c|}{~} & ~ & ~ \\  \multicolumn{1}{|c|}{~} & ~ & ~ \\
1585  \multicolumn{1}{l}{\bf xx\_salt0...$<$k$>$} & $\longrightarrow$ &  \multicolumn{1}{l}{\bf xx\_salt0...$<$k$>$} &
1586    $\stackrel{\mbox{read}}{\longrightarrow}$ &
1587  \multicolumn{1}{|c|}{forward integration} & ~ & ~ \\  \multicolumn{1}{|c|}{forward integration} & ~ & ~ \\
1588  \multicolumn{1}{l}{\bf \vdots} & ~ & \multicolumn{1}{|c|}{~}    \multicolumn{1}{l}{\bf \vdots} & ~ & \multicolumn{1}{|c|}{~}  
1589  & ~ & ~ \\  & ~ & ~ \\
1590  \cline{3-3}  \cline{3-3}
1591  ~ & ~ & ~ & ~ & ~ \\  ~ & ~ & $\downarrow$ & ~ & ~ \\
1592  \cline{3-3}  \cline{3-3}
1593  ~ & ~ &  ~ & ~ &
1594  \multicolumn{1}{|c|}{~} & ~ &  \multicolumn{1}{|c|}{~} & ~ &
1595  \multicolumn{1}{l}{\bf adxx\_theta0...$<$k$>$}  \\  \multicolumn{1}{l}{\bf adxx\_theta0...$<$k$>$}  \\
1596  ~ & ~ & \multicolumn{1}{|c|}{adjoint integration} &  ~ & ~ & \multicolumn{1}{|c|}{adjoint integration} &
1597  $\longrightarrow$ &  $\stackrel{\mbox{write}}{\longrightarrow}$ &
1598  \multicolumn{1}{l}{\bf adxx\_salt0...$<$k$>$} \\  \multicolumn{1}{l}{\bf adxx\_salt0...$<$k$>$} \\
1599  ~ & ~ & \multicolumn{1}{|c|}{~}    ~ & ~ & \multicolumn{1}{|c|}{~}  
1600  & ~ & \multicolumn{1}{l}{\bf \vdots} \\  & ~ & \multicolumn{1}{l}{\bf \vdots} \\
# Line 1630  $\longrightarrow$ & Line 1606  $\longrightarrow$ &
1606  ~ & ~ & ~ & ~ &  {\big\downarrow} \\  ~ & ~ & ~ & ~ &  {\big\downarrow} \\
1607  ~ & ~ & ~ & ~ &  {\bf vector\_grad\_$<$k$>$ } \\  ~ & ~ & ~ & ~ &  {\bf vector\_grad\_$<$k$>$ } \\
1608  \end{tabular}  \end{tabular}
1609    }
1610    
1611  \vspace*{0.5cm}  \vspace*{0.5cm}
1612    
1613    
1614  {\it ctrl\_unpack} reads in the updated control vector  {\it ctrl\_unpack} reads the updated control vector
1615  {\bf vector\_ctrl\_$<$k$>$}.  {\bf vector\_ctrl\_$<$k$>$}.
1616  It distributes the different control variables to  It distributes the different control variables to
1617  2-dim. and 3-dim. files {\it xx\_...$<$k$>$}.  2-dim. and 3-dim. files {\it xx\_...$<$k$>$}.
1618  During the forward integration the control variables  At the start of the forward integration the control variables
1619  are read from {\it xx\_...$<$k$>$}.  are read from {\it xx\_...$<$k$>$} and added to the
1620  Correspondingly, the adjoint fields are written  field.
1621    Correspondingly, at the end of the adjoint integration
1622    the adjoint fields are written
1623  to {\it adxx\_...$<$k$>$}, again via the active file routines.  to {\it adxx\_...$<$k$>$}, again via the active file routines.
1624  Finally, {\it ctrl\_pack} collects all adjoint field files  Finally, {\it ctrl\_pack} collects all adjoint files
1625  and writes them to the compressed vector file  and writes them to the compressed vector file
1626  {\bf vector\_grad\_$<$k$>$}.  {\bf vector\_grad\_$<$k$>$}.
   
 \subsection{TLM and ADM generation via TAMC}  
   
   
   
 \subsection{Flow directives and adjoint support routines}  
   
 \subsection{Store directives and checkpointing}  
   
 \subsection{Gradient checks}  
   
 \subsection{Second derivative generation via TAMC}  
   
 \section{Example of adjoint code}  

Legend:
Removed from v.1.3  
changed lines
  Added in v.1.23

  ViewVC Help
Powered by ViewVC 1.1.22