/[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.1 by adcroft, Wed Aug 8 16:16:26 2001 UTC revision 1.19 by heimbach, Tue Aug 2 22:26:58 2005 UTC
# Line 1  Line 1 
1  % $Header$  % $Header$
2  % $Name$  % $Name$
3    
4    Author: Patrick Heimbach
5    
6  {\sf Automatic differentiation} (AD), also referred to as algorithmic  {\sf Automatic differentiation} (AD), also referred to as algorithmic
7  (or, more loosely, computational) differentiation, involves  (or, more loosely, computational) differentiation, involves
8  automatically deriving code to calculate  automatically deriving code to calculate
# Line 18  In principle, a variety of derived algor Line 20  In principle, a variety of derived algor
20  can be generated automatically in this way.  can be generated automatically in this way.
21    
22  The MITGCM has been adapted for use with the  The MITGCM has been adapted for use with the
23  Tangent linear and Adjoint Model Compiler (TAMC) and its succssor TAF  Tangent linear and Adjoint Model Compiler (TAMC) and its successor TAF
24  (Transformation of Algorithms in Fortran), developed  (Transformation of Algorithms in Fortran), developed
25  by Ralf Giering (\cite{gie-kam:98}, \cite{gie:99,gie:00}).  by Ralf Giering (\cite{gie-kam:98}, \cite{gie:99,gie:00}).
26  The first application of the adjoint of the MITGCM for senistivity  The first application of the adjoint of the MITGCM for sensitivity
27  studies has been published by \cite{maro-eta:99}.  studies has been published by \cite{maro-eta:99}.
28  \cite{sta-eta:97,sta-eta:01} use the MITGCM and its adjoint  \cite{sta-eta:97,sta-eta:01} use the MITGCM and its adjoint
29  for ocean state estimation studies.  for ocean state estimation studies.
30    In the following we shall refer to TAMC and TAF synonymously,
31    except were explicitly stated otherwise.
32    
33  TAMC exploits the chain rule for computing the first  TAMC exploits the chain rule for computing the first
34  derivative of a function with  derivative of a function with
35  respect to a set of input variables.  respect to a set of input variables.
36  Treating a given forward code as a composition of operations --  Treating a given forward code as a composition of operations --
37  each line representing a compositional element -- the chain rule is  each line representing a compositional element, the chain rule is
38  rigorously applied to the code, line by line. The resulting  rigorously applied to the code, line by line. The resulting
39  tangent linear or adjoint code,  tangent linear or adjoint code,
40  then, may be thought of as the composition in  then, may be thought of as the composition in
41  forward or reverse order, respectively, of the  forward or reverse order, respectively, of the
42  Jacobian matrices of the forward code compositional elements.  Jacobian matrices of the forward code's compositional elements.
43    
44  %**********************************************************************  %**********************************************************************
45  \section{Some basic algebra}  \section{Some basic algebra}
46  \label{sec_ad_algebra}  \label{sec_ad_algebra}
47    \begin{rawhtml}
48    <!-- CMIREDIR:sec_ad_algebra: -->
49    \end{rawhtml}
50  %**********************************************************************  %**********************************************************************
51    
52  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 57  $\vec{u}=(u_1,\ldots,u_m)$
57  such as forcing functions) to the $n$-dimensional space  such as forcing functions) to the $n$-dimensional space
58  $V \subset I\!\!R^n$ of  $V \subset I\!\!R^n$ of
59  model output variable $\vec{v}=(v_1,\ldots,v_n)$  model output variable $\vec{v}=(v_1,\ldots,v_n)$
60  (model state, model diagnostcs, objective function, ...)  (model state, model diagnostics, objective function, ...)
61  under consideration,  under consideration,
62  %  %
63  \begin{equation}  \begin{equation}
# Line 105  In contrast to the full nonlinear model Line 112  In contrast to the full nonlinear model
112  $ M $ is just a matrix  $ M $ is just a matrix
113  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
114  perturbations in  $u$,  perturbations in  $u$,
115  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
116  large-scale oceanographic application), it quickly becomes  large-scale oceanographic application), it quickly becomes
117  prohibitive to proceed directly as in (\ref{tangent_linear}),  prohibitive to proceed directly as in (\ref{tangent_linear}),
118  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 137  or a measure of some model-to-data misfi
137  \label{compo}  \label{compo}
138  \end{eqnarray}  \end{eqnarray}
139  %  %
140  The linear approximation of $ {\cal J} $,  The perturbation of $ {\cal J} $ around a fixed point $ {\cal J}_0 $,
141  \[  \[
142  {\cal J} \, \approx \, {\cal J}_0 \, + \, \delta {\cal J}  {\cal J} \, = \, {\cal J}_0 \, + \, \delta {\cal J}
143  \]  \]
144  can be expressed in both bases of $ \vec{u} $ and $ \vec{v} $  can be expressed in both bases of $ \vec{u} $ and $ \vec{v} $
145  w.r.t. their corresponding inner product  w.r.t. their corresponding inner product
# Line 152  $\left\langle \,\, , \,\, \right\rangle Line 159  $\left\langle \,\, , \,\, \right\rangle
159  \label{deljidentity}  \label{deljidentity}
160  \end{equation}  \end{equation}
161  %  %
162  (note, that the gradient $ \nabla f $ is a pseudo-vector, therefore  (note, that the gradient $ \nabla f $ is a co-vector, therefore
163  its transpose is required in the above inner product).  its transpose is required in the above inner product).
164  Then, using the representation of  Then, using the representation of
165  $ \delta {\cal J} =  $ \delta {\cal J} =
# Line 168  transpose of $ A $, Line 175  transpose of $ A $,
175  \[  \[
176  A^{\ast} \, = \, A^T  A^{\ast} \, = \, A^T
177  \]  \]
178  and from eq. (\ref{tangent_linear}), we note that  and from eq. (\ref{tangent_linear}), (\ref{deljidentity}),
179    we note that
180  (omitting $|$'s):  (omitting $|$'s):
181  %  %
182  \begin{equation}  \begin{equation}
# Line 204  the adjoint variable of the model state Line 212  the adjoint variable of the model state
212  $ \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} $.
213    
214  The {\sf reverse} nature of the adjoint calculation can be readily  The {\sf reverse} nature of the adjoint calculation can be readily
215  seen as follows. Let us decompose ${\cal J}(u)$, thus:  seen as follows.
216    Consider a model integration which consists of $ \Lambda $
217    consecutive operations
218    $ {\cal M}_{\Lambda} (  {\cal M}_{\Lambda-1} (
219    ...... ( {\cal M}_{\lambda} (
220    ......
221    ( {\cal M}_{1} ( {\cal M}_{0}(\vec{u}) )))) $,
222    where the ${\cal M}$'s could be the elementary steps, i.e. single lines
223    in the code of the model, or successive time steps of the
224    model integration,
225    starting at step 0 and moving up to step $\Lambda$, with intermediate
226    ${\cal M}_{\lambda} (\vec{u}) = \vec{v}^{(\lambda+1)}$ and final
227    ${\cal M}_{\Lambda} (\vec{u}) = \vec{v}^{(\Lambda+1)} = \vec{v}$.
228    Let ${\cal J}$ be a cost function which explicitly depends on the
229    final state $\vec{v}$ only
230    (this restriction is for clarity reasons only).
231    %
232    ${\cal J}(u)$ may be decomposed according to:
233  %  %
234  \begin{equation}  \begin{equation}
235  {\cal J}({\cal M}(\vec{u})) \, = \,  {\cal J}({\cal M}(\vec{u})) \, = \,
# Line 215  seen as follows. Let us decompose ${\cal Line 240  seen as follows. Let us decompose ${\cal
240  \label{compos}  \label{compos}
241  \end{equation}  \end{equation}
242  %  %
243  where the ${\cal M}$'s could be the elementary steps, i.e. single lines  Then, according to the chain rule, the forward calculation reads,
244  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  
245  (we've omitted the $ | $'s which, nevertheless are important  (we've omitted the $ | $'s which, nevertheless are important
246  to the aspect of {\it tangent} linearity;  to the aspect of {\it tangent} linearity;
247  note also that per definition  note also that by definition
248  $ \langle \, \nabla _{v}{\cal J}^T \, , \, \delta \vec{v} \, \rangle  $ \langle \, \nabla _{v}{\cal J}^T \, , \, \delta \vec{v} \, \rangle
249  = \nabla_v {\cal J} \cdot \delta \vec{v} $ )  = \nabla_v {\cal J} \cdot \delta \vec{v} $ )
250  %  %
# Line 259  M_{\Lambda}^T \cdot \nabla_v {\cal J}^T Line 279  M_{\Lambda}^T \cdot \nabla_v {\cal J}^T
279  %  %
280  clearly expressing the reverse nature of the calculation.  clearly expressing the reverse nature of the calculation.
281  Eq. (\ref{reverse}) is at the heart of automatic adjoint compilers.  Eq. (\ref{reverse}) is at the heart of automatic adjoint compilers.
282  The intermediate steps $\lambda$ in  If the intermediate steps $\lambda$ in
283  eqn. (\ref{compos}) -- (\ref{reverse})  eqn. (\ref{compos}) -- (\ref{reverse})
284  could represent the model state (forward or adjoint) at each  represent the model state (forward or adjoint) at each
285  intermediate time step in which case  intermediate time step as noted above, then correspondingly,
286  $ {\cal M}(\vec{v}^{(\lambda)}) = \vec{v}^{(\lambda+1)} $, and correspondingly,  $ M^T (\delta \vec{v}^{(\lambda) \, \ast}) =
287  $ M^T (\delta \vec{v}^{(\lambda) \, \ast}) = \delta \vec{v}^{(\lambda-1) \, \ast} $,  \delta \vec{v}^{(\lambda-1) \, \ast} $ for the adjoint variables.
288  but they can also be viewed more generally as  It thus becomes evident that the adjoint calculation also
289  single lines of code in the numerical algorithm.  yields the adjoint of each model state component
290  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  
291  %  %
292  \begin{equation}  \begin{equation}
293  \boxed{  \boxed{
# Line 285  M_{\Lambda}^T |_{\vec{v}^{(\lambda)}} \c Line 303  M_{\Lambda}^T |_{\vec{v}^{(\lambda)}} \c
303  %  %
304  in close analogy to eq. (\ref{adjoint})  in close analogy to eq. (\ref{adjoint})
305  We note in passing that that the $\delta \vec{v}^{(\lambda) \, \ast}$  We note in passing that that the $\delta \vec{v}^{(\lambda) \, \ast}$
306  are the Lagrange multipliers of the model state $ \vec{v}^{(\lambda)}$.  are the Lagrange multipliers of the model equations which determine
307    $ \vec{v}^{(\lambda)}$.
308    
309  In coponents, eq. (\ref{adjoint}) reads as follows.  In components, eq. (\ref{adjoint}) reads as follows.
310  Let  Let
311  \[  \[
312  \begin{array}{rclcrcl}  \begin{array}{rclcrcl}
# Line 308  Let Line 327  Let
327  \end{array}  \end{array}
328  \]  \]
329  denote the perturbations in $\vec{u}$ and $\vec{v}$, respectively,  denote the perturbations in $\vec{u}$ and $\vec{v}$, respectively,
330  and their adjoint varaiables;  and their adjoint variables;
331  further  further
332  \[  \[
333  M \, = \, \left(  M \, = \, \left(
# Line 395  and the shorthand notation for the adjoi Line 414  and the shorthand notation for the adjoi
414  $ \delta v^{(\lambda) \, \ast}_{j} = \frac{\partial}{\partial v^{(\lambda)}_{j}}  $ \delta v^{(\lambda) \, \ast}_{j} = \frac{\partial}{\partial v^{(\lambda)}_{j}}
415  {\cal J}^T $, $ j = 1, \ldots , n_{\lambda} $,  {\cal J}^T $, $ j = 1, \ldots , n_{\lambda} $,
416  for intermediate components, yielding  for intermediate components, yielding
417  \[  \begin{equation}
418  \footnotesize  \small
419    \begin{split}
420  \left(  \left(
421  \begin{array}{c}  \begin{array}{c}
422  \delta v^{(\lambda) \, \ast}_1 \\  \delta v^{(\lambda) \, \ast}_1 \\
# Line 404  for intermediate components, yielding Line 424  for intermediate components, yielding
424  \delta v^{(\lambda) \, \ast}_{n_{\lambda}} \\  \delta v^{(\lambda) \, \ast}_{n_{\lambda}} \\
425  \end{array}  \end{array}
426  \right)  \right)
427  \, = \,  \, = &
428  \left(  \left(
429  \begin{array}{ccc}  \begin{array}{ccc}
430  \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_1}  \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_1}
431  & \ldots &  & \ldots \,\, \ldots &
432  \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_1} \\  \frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_1} \\
433  \vdots & ~ & \vdots \\  \vdots & ~ & \vdots \\
434  \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_{n_{\lambda}}}  \frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_{n_{\lambda}}}
435  & \ldots  &  & \ldots \,\, \ldots  &
436  \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}}} \\
437  \end{array}  \end{array}
438  \right)  \right)
 %  
439  \cdot  \cdot
440  %  %
441    \\ ~ & ~
442    \\ ~ &
443    %
444  \left(  \left(
445  \begin{array}{ccc}  \begin{array}{ccc}
446  \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 453  for intermediate components, yielding
453  \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}}} \\
454  \end{array}  \end{array}
455  \right)  \right)
456  \cdot \ldots \ldots \cdot  \cdot \, \ldots \, \cdot
457  \left(  \left(
458  \begin{array}{c}  \begin{array}{c}
459  \delta v^{\ast}_1 \\  \delta v^{\ast}_1 \\
# Line 439  for intermediate components, yielding Line 461  for intermediate components, yielding
461  \delta v^{\ast}_{n} \\  \delta v^{\ast}_{n} \\
462  \end{array}  \end{array}
463  \right)  \right)
464  \]  \end{split}
465    \end{equation}
466    
467  Eq. (\ref{forward}) and (\ref{reverse}) are perhaps clearest in  Eq. (\ref{forward}) and (\ref{reverse}) are perhaps clearest in
468  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 473  variables $u$
473  {\it all} intermediate states $ \vec{v}^{(\lambda)} $) are sought.  {\it all} intermediate states $ \vec{v}^{(\lambda)} $) are sought.
474  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
475  $ \partial {\cal J} / \partial u_{i} $ in (\ref{forward})  $ \partial {\cal J} / \partial u_{i} $ in (\ref{forward})
476  a forward calulation has to be performed for each component seperately,  a forward calculation has to be performed for each component separately,
477  i.e. $ \delta \vec{u} = \delta u_{i} {\vec{e}_{i}} $  i.e. $ \delta \vec{u} = \delta u_{i} {\vec{e}_{i}} $
478  for  the $i$-th forward calculation.  for  the $i$-th forward calculation.
479  Then, (\ref{forward}) represents the  Then, (\ref{forward}) represents the
# Line 460  In contrast, eq. (\ref{reverse}) yields Line 483  In contrast, eq. (\ref{reverse}) yields
483  gradient $\nabla _{u}{\cal J}$ (and all intermediate gradients  gradient $\nabla _{u}{\cal J}$ (and all intermediate gradients
484  $\nabla _{v^{(\lambda)}}{\cal J}$) within a single reverse calculation.  $\nabla _{v^{(\lambda)}}{\cal J}$) within a single reverse calculation.
485    
486  Note, that in case $ {\cal J} $ is a vector-valued function  Note, that if $ {\cal J} $ is a vector-valued function
487  of dimension $ l > 1 $,  of dimension $ l > 1 $,
488  eq. (\ref{reverse}) has to be modified according to  eq. (\ref{reverse}) has to be modified according to
489  \[  \[
# Line 468  M^T \left( \nabla_v {\cal J}^T \left(\de Line 491  M^T \left( \nabla_v {\cal J}^T \left(\de
491  \, = \,  \, = \,
492  \nabla_u {\cal J}^T \cdot \delta \vec{J}  \nabla_u {\cal J}^T \cdot \delta \vec{J}
493  \]  \]
494  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
495    dimension $ l $.
496  In this case $ l $ reverse simulations have to be performed  In this case $ l $ reverse simulations have to be performed
497  for each $ \delta J_{k}, \,\, k = 1, \ldots, l $.  for each $ \delta J_{k}, \,\, k = 1, \ldots, l $.
498  Then, the reverse mode is more efficient as long as  Then, the reverse mode is more efficient as long as
499  $ l < n $, otherwise the forward mode is preferable.  $ l < n $, otherwise the forward mode is preferable.
500  Stricly, the reverse mode is called adjoint mode only for  Strictly, the reverse mode is called adjoint mode only for
501  $ l = 1 $.  $ l = 1 $.
502    
503  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 527  operator onto the $j$-th component ${\bf
527  \paragraph{Example 2:  \paragraph{Example 2:
528  $ {\cal J} = \langle \, {\cal H}(\vec{v}) - \vec{d} \, ,  $ {\cal J} = \langle \, {\cal H}(\vec{v}) - \vec{d} \, ,
529   \, {\cal H}(\vec{v}) - \vec{d} \, \rangle $} ~ \\   \, {\cal H}(\vec{v}) - \vec{d} \, \rangle $} ~ \\
530  The cost function represents the quadratic model vs.data misfit.  The cost function represents the quadratic model vs. data misfit.
531  Here, $ \vec{d} $ is the data vector and $ {\cal H} $ represents the  Here, $ \vec{d} $ is the data vector and $ {\cal H} $ represents the
532  operator which maps the model state space onto the data space.  operator which maps the model state space onto the data space.
533  Then, $ \nabla_v {\cal J} $ takes the form  Then, $ \nabla_v {\cal J} $ takes the form
# Line 534  H \cdot \left( {\cal H}(\vec{v}) - \vec{ Line 558  H \cdot \left( {\cal H}(\vec{v}) - \vec{
558    
559  We note an important aspect of the forward vs. reverse  We note an important aspect of the forward vs. reverse
560  mode calculation.  mode calculation.
561  Because of the locality of the derivative,  Because of the local character of the derivative
562    (a derivative is defined w.r.t. a point along the trajectory),
563  the intermediate results of the model trajectory  the intermediate results of the model trajectory
564  $\vec{v}^{(\lambda+1)}={\cal M}_{\lambda}(v^{(\lambda)})$  $\vec{v}^{(\lambda+1)}={\cal M}_{\lambda}(v^{(\lambda)})$
565  are needed to evaluate the intermediate Jacobian  may be required to evaluate the intermediate Jacobian
566  $M_{\lambda}|_{\vec{v}^{(\lambda)}} \, \delta \vec{v}^{(\lambda)} $.  $M_{\lambda}|_{\vec{v}^{(\lambda)}} \, \delta \vec{v}^{(\lambda)} $.
567    This is the case e.g. for nonlinear expressions
568    (momentum advection, nonlinear equation of state), state-dependent
569    conditional statements (parameterization schemes).
570  In the forward mode, the intermediate results are required  In the forward mode, the intermediate results are required
571  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}$,
572  in the reverse mode they are required in the reverse order.  but in the reverse mode they are required in the reverse order.
573  Thus, in the reverse mode the trajectory of the forward model  Thus, in the reverse mode the trajectory of the forward model
574  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
575  calculation. Alternatively, the model state would have to be  calculation. Alternatively, the complete model state up to the
576  recomputed whenever its value is required.  point of evaluation has to be recomputed whenever its value is required.
577    
578  A method to balance the amount of recomputations vs.  A method to balance the amount of recomputations vs.
579  storage requirements is called {\sf checkpointing}  storage requirements is called {\sf checkpointing}
580  (e.g. \cite{res-eta:98}).  (e.g. \cite{gri:92}, \cite{res-eta:98}).
581  It is depicted in Fig. ... for a 3-level checkpointing  It is depicted in \ref{fig:3levelcheck} for a 3-level checkpointing
582  [as concrete example, we give explicit numbers for a 3-day  [as an example, we give explicit numbers for a 3-day
583  integration with a 1-hourly timestep in square brackets].  integration with a 1-hourly timestep in square brackets].
584  \begin{itemize}  \begin{itemize}
585  %  %
# Line 559  integration with a 1-hourly timestep in Line 587  integration with a 1-hourly timestep in
587  In a first step, the model trajectory is subdivided into  In a first step, the model trajectory is subdivided into
588  $ {n}^{lev3} $ subsections [$ {n}^{lev3} $=3 1-day intervals],  $ {n}^{lev3} $ subsections [$ {n}^{lev3} $=3 1-day intervals],
589  with the label $lev3$ for this outermost loop.  with the label $lev3$ for this outermost loop.
590  The model is then integrated over the full trajectory,  The model is then integrated along the full trajectory,
591  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
592  [i.e. 3 times, at  [i.e. 3 times, at
593  $ i = 0,1,2 $ corresponding to $ k_{i}^{lev3} = 0, 24, 48 $].  $ i = 0,1,2 $ corresponding to $ k_{i}^{lev3} = 0, 24, 48 $].
594    In addition, the cost function is computed, if needed.
595  %  %
596  \item [$lev2$]  \item [$lev2$]
597  In a second step each subsection is itself divided into  In a second step each subsection itself is divided into
598  $ {n}^{lev2} $ subsubsections  $ {n}^{lev2} $ subsections
599  [$ {n}^{lev2} $=4 6-hour intervals per subsection].  [$ {n}^{lev2} $=4 6-hour intervals per subsection].
600  The model picks up at the last outermost dumped state  The model picks up at the last outermost dumped state
601  $ v_{k_{n}^{lev3}} $ and is integrated forward in time over  $ v_{k_{n}^{lev3}} $ and is integrated forward in time along
602  the last subsection, with the label $lev2$ for this    the last subsection, with the label $lev2$ for this  
603  intermediate loop.  intermediate loop.
604  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
605  timestep  timestep
606  [i.e. 4 times, at  [i.e. 4 times, at
607  $ 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 $].
608  %  %
609  \item [$lev1$]  \item [$lev1$]
610  Finally, the mode picks up at the last intermediate dump state  Finally, the model picks up at the last intermediate dump state
611  $ v_{k_{n}^{lev2}} $ and is integrated forward in time over  $ v_{k_{n}^{lev2}} $ and is integrated forward in time along
612  the last subsubsection, with the label $lev1$ for this    the last subsection, with the label $lev1$ for this  
613  intermediate loop.  intermediate loop.
614  Within this subsubsection only, the model state is stored  Within this sub-subsection only, parts of the model state is stored
615  at every timestep  to memory at every timestep
616  [i.e. every hour $ i=0,...,5$ corresponding to  [i.e. every hour $ i=0,...,5$ corresponding to
617  $ k_{i}^{lev1} = 66, 67, \ldots, 71 $].  $ k_{i}^{lev1} = 66, 67, \ldots, 71 $].
618  Thus, the  final state $ v_n = v_{k_{n}^{lev1}} $ is reached  The  final state $ v_n = v_{k_{n}^{lev1}} $ is reached
619  and the model state of all peceeding timesteps over the last  and the model state of all preceding timesteps along the last
620  subsubsections are available, enabling integration backwards  innermost subsection are available, enabling integration backwards
621  in time over the last subsubsection.  in time along the last subsection.
622  Thus, the adjoint can be computed over this last  The adjoint can thus be computed along this last
623  subsubsection $k_{n}^{lev2}$.  subsection $k_{n}^{lev2}$.
624  %  %
625  \end{itemize}  \end{itemize}
626  %  %
627  This procedure is repeated consecutively for each previous  This procedure is repeated consecutively for each previous
628  subsubsection $k_{n-1}^{lev2}, \ldots, k_{1}^{lev2} $  subsection $k_{n-1}^{lev2}, \ldots, k_{1}^{lev2} $
629  carrying the adjoint computation to the initial time  carrying the adjoint computation to the initial time
630  of the subsection $k_{n}^{lev3}$.  of the subsection $k_{n}^{lev3}$.
631  Then, the procedure is repeated for the previous subsection  Then, the procedure is repeated for the previous subsection
# Line 607  $k_{1}^{lev3}$. Line 636  $k_{1}^{lev3}$.
636  For the full model trajectory of  For the full model trajectory of
637  $ n^{lev3} \cdot n^{lev2} \cdot n^{lev1} $ timesteps  $ n^{lev3} \cdot n^{lev2} \cdot n^{lev1} $ timesteps
638  the required storing of the model state was significantly reduced to  the required storing of the model state was significantly reduced to
639  $ n^{lev1} + n^{lev2} + n^{lev3} $  $ n^{lev2} + n^{lev3} $ to disk and roughly $ n^{lev1} $ to memory
640  [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
641  the model state was stored 13 times].  the model state was stored 7 times to disk and roughly 6 times
642    to memory].
643  This saving in memory comes at a cost of a required  This saving in memory comes at a cost of a required
644  3 full forward integrations of the model (one for each  3 full forward integrations of the model (one for each
645  checkpointing level).  checkpointing level).
646  The balance of storage vs. recomputation certainly depends  The optimal balance of storage vs. recomputation certainly depends
647  on the computing resources available.  on the computing resources available and may be adjusted by
648    adjusting the partitioning among the
649    $ n^{lev3}, \,\, n^{lev2}, \,\, n^{lev1} $.
650    
651  \begin{figure}[t!]  \begin{figure}[t!]
652  \centering  \begin{center}
653  %\psdraft  %\psdraft
654  \psfrag{v_k1^lev3}{\mathinfigure{v_{k_{1}^{lev3}}}}  %\psfrag{v_k1^lev3}{\mathinfigure{v_{k_{1}^{lev3}}}}
655  \psfrag{v_kn-1^lev3}{\mathinfigure{v_{k_{n-1}^{lev3}}}}  %\psfrag{v_kn-1^lev3}{\mathinfigure{v_{k_{n-1}^{lev3}}}}
656  \psfrag{v_kn^lev3}{\mathinfigure{v_{k_{n}^{lev3}}}}  %\psfrag{v_kn^lev3}{\mathinfigure{v_{k_{n}^{lev3}}}}
657  \psfrag{v_k1^lev2}{\mathinfigure{v_{k_{1}^{lev2}}}}  %\psfrag{v_k1^lev2}{\mathinfigure{v_{k_{1}^{lev2}}}}
658  \psfrag{v_kn-1^lev2}{\mathinfigure{v_{k_{n-1}^{lev2}}}}  %\psfrag{v_kn-1^lev2}{\mathinfigure{v_{k_{n-1}^{lev2}}}}
659  \psfrag{v_kn^lev2}{\mathinfigure{v_{k_{n}^{lev2}}}}  %\psfrag{v_kn^lev2}{\mathinfigure{v_{k_{n}^{lev2}}}}
660  \psfrag{v_k1^lev1}{\mathinfigure{v_{k_{1}^{lev1}}}}  %\psfrag{v_k1^lev1}{\mathinfigure{v_{k_{1}^{lev1}}}}
661  \psfrag{v_kn^lev1}{\mathinfigure{v_{k_{n}^{lev1}}}}  %\psfrag{v_kn^lev1}{\mathinfigure{v_{k_{n}^{lev1}}}}
662  \mbox{\epsfig{file=part5/checkpointing.eps, width=0.8\textwidth}}  %\mbox{\epsfig{file=part5/checkpointing.eps, width=0.8\textwidth}}
663    \resizebox{5.5in}{!}{\includegraphics{part5/checkpointing.eps}}
664  %\psfull  %\psfull
665  \caption  \end{center}
666  {Schematic view of intermediate dump and restart for  \caption{
667    Schematic view of intermediate dump and restart for
668  3-level checkpointing.}  3-level checkpointing.}
669  \label{fig:erswns}  \label{fig:3levelcheck}
670  \end{figure}  \end{figure}
671    
672  \subsection{Optimal perturbations}  % \subsection{Optimal perturbations}
673  \label{optpert}  % \label{sec_optpert}
674    
675    
676  \subsection{Error covariance estimate and Hessian matrix}  % \subsection{Error covariance estimate and Hessian matrix}
677  \label{sec_hessian}  % \label{sec_hessian}
678    
679  \newpage  \newpage
680    
681  %**********************************************************************  %**********************************************************************
682  \section{AD-specific setup by example: sensitivity of carbon sequestration}  \section{TLM and ADM generation in general}
683  \label{sec_ad_setup_ex}  \label{sec_ad_setup_gen}
684    \begin{rawhtml}
685    <!-- CMIREDIR:sec_ad_setup_gen: -->
686    \end{rawhtml}
687  %**********************************************************************  %**********************************************************************
688    
689  The MITGCM has been adapted to enable AD using TAMC or TAF  In this section we describe in a general fashion
690  (we'll refer to TAMC and TAF interchangeably, except where  the parts of the code that are relevant for automatic
691  distinctions are explicitly mentioned).  differentiation using the software tool TAF.
692  The present description, therefore, is specific to the  
693  use of TAMC as AD tool.  \input{part5/doc_ad_the_model}
694  The following sections describe the steps which are necessary to  
695  generate a tangent linear or adjoint model of the MITGCM.  The basic flow is depicted in \ref{fig:adthemodel}.
696  We take as an example the sensitivity of carbon sequestration  If CPP option \texttt{ALLOW\_AUTODIFF\_TAMC} is defined,
697  in the ocean.  the driver routine
698  The AD-relevant hooks in the code are sketched in  {\it the\_model\_main}, instead of calling {\it the\_main\_loop},
699  \reffig{adthemodel}, \reffig{adthemain}.  invokes the adjoint of this routine, {\it adthe\_main\_loop}
700    (case \texttt{\#define ALLOW\_ADJOINT\_RUN}), or
701  \subsection{Overview of the experiment}  the tangent linear of this routine {\it g\_the\_main\_loop}
702    (case \texttt{\#define ALLOW\_TANGENTLINEAR\_RUN}),
703  We describe an adjoint sensitivity analysis of outgassing from  which are the toplevel routines in terms of automatic differentiation.
704  the ocean into the atmosphere of a carbon like tracer injected  The routines {\it adthe\_main\_loop} or {\it g\_the\_main\_loop}
705  into the ocean interior (see \cite{hil-eta:01}).  are generated by TAF.
706    It contains both the forward integration of the full model, the
707  \subsubsection{Passive tracer equation}  cost function calculation,
708    any additional storing that is required for efficient checkpointing,
709  For this work the MITGCM was augmented with a thermodynamically  and the reverse integration of the adjoint model.
710  inactive tracer, $C$. Tracer residing in the ocean  
711  model surface layer is outgassed according to a relaxation time scale,  [DESCRIBE IN A SEPARATE SECTION THE WORKING OF THE TLM]
712  $\mu$. Within the ocean interior, the tracer is passively advected  
713  by the ocean model currents. The full equation for the time evolution  In Fig. \ref{fig:adthemodel}
714  %  the structure of {\it adthe\_main\_loop} has been strongly
715  \begin{equation}  simplified to focus on the essentials; in particular, no checkpointing
716  \label{carbon_ddt}  procedures are shown here.
717  \frac{\partial C}{\partial t} \, = \,  Prior to the call of {\it adthe\_main\_loop}, the routine
718  -U\cdot \nabla C \, - \, \mu C \, + \, \Gamma(C) \,+ \, S  {\it ctrl\_unpack} is invoked to unpack the control vector
719  \end{equation}  or initialise the control variables.
720  %  Following the call of {\it adthe\_main\_loop},
721  also includes a source term $S$. This term  the routine {\it ctrl\_pack}
722  represents interior sources of $C$ such as would arise due to  is invoked to pack the control vector
723  direct injection.  (cf. Section \ref{section_ctrl}).
724  The velocity term, $U$, is the sum of the  If gradient checks are to be performed, the option
725  model Eulerian circulation and an eddy-induced velocity, the latter  {\tt ALLOW\_GRADIENT\_CHECK} is defined. In this case
726  parameterized according to Gent/McWilliams (\cite{gen:90, dan:95}).  the driver routine {\it grdchk\_main} is called after
727  The convection function, $\Gamma$, mixes $C$ vertically wherever the  the gradient has been computed via the adjoint
728  fluid is locally statically unstable.  (cf. Section \ref{section_grdchk}).
729    
730  The outgassing time scale, $\mu$, in eqn. (\ref{carbon_ddt})  %------------------------------------------------------------------
731  is set so that \( 1/\mu \sim 1 \ \mathrm{year} \) for the surface  
732  ocean and $\mu=0$ elsewhere. With this value, eqn. (\ref{carbon_ddt})  \subsection{General setup
733  is valid as a prognostic equation for small perturbations in oceanic  \label{section_ad_setup}}
734  carbon concentrations. This configuration provides a  
735  powerful tool for examining the impact of large-scale ocean circulation  In order to configure AD-related setups the following packages need
736  on $ CO_2 $ outgassing due to interior injections.  to be enabled:
737  As source we choose a constant in time injection of  {\it
738  $ S = 1 \,\, {\rm mol / s}$.  \begin{table}[h!]
739    \begin{tabular}{l}
740  \subsubsection{Model configuration}  autodiff \\
741    ctrl \\
742  The model configuration employed has a constant  cost \\
743  $4^\circ \times 4^\circ$ resolution horizontal grid and realistic  grdchk \\
744  geography and bathymetry. Twenty vertical layers are used with  \end{tabular}
745  vertical spacing ranging  \end{table}
746  from 50 m near the surface to 815 m at depth.  }
747  Driven to steady-state by climatalogical wind-stress, heat and  The packages are enabled by adding them to your experiment-specific
748  fresh-water forcing the model reproduces well known large-scale  configuration file
749  features of the ocean general circulation.  {\it packages.conf} (see Section ???).
750    
751  \subsubsection{Outgassing cost function}  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/}:  
752  %  %
753  \begin{itemize}  \begin{itemize}
754  %  %
755  \item {\it .genmakerc}  \item {\it ECCO\_CPPOPTIONS.h} \\
756  %  This header file collects CPP options for the packages
757  \item {\it COST\_CPPOPTIONS.h}  {\it autodiff, cost, ctrl} as well as AD-unrelated options for
758  %  the external forcing package {\it exf}.
759  \item {\it CPP\_EEOPTIONS.h}  \footnote{NOTE: These options are not set in their package-specific
760  %  headers such as {\it COST\_CPPOPTIONS.h}, but are instead collected
761  \item {\it CPP\_OPTIONS.h}  in the single header file {\it ECCO\_CPPOPTIONS.h}.
762  %  The package-specific header files serve as simple
763  \item {\it CTRL\_OPTIONS.h}  placeholders at this point.}
764  %  %
765  \item {\it ECCO\_OPTIONS.h}  \item {\it tamc.h} \\
766  %  This header configures the splitting of the time stepping loop
767  \item {\it SIZE.h}  w.r.t. the 3-level checkpointing (see section ???).
768  %  
 \item {\it adcommon.h}  
 %  
 \item {\it tamc.h}  
769  %  %
770  \end{itemize}  \end{itemize}
771    
772    %------------------------------------------------------------------
773    
774    \subsection{Building the AD code
775    \label{section_ad_build}}
776    
777    The build process of an AD code is very similar to building
778    the forward model. However, depending on which AD code one wishes
779    to generate, and on which AD tool is available (TAF or TAMC),
780    the following {\tt make} targets are available:
781    
782    \begin{table}[h!]
783    {\footnotesize
784    \begin{tabular}{ccll}
785    ~ & {\it AD-target} & {\it output} & {\it description} \\
786    \hline
787    \hline
788    (1) & {\tt <MODE><TOOL>only} & {\tt <MODE>\_<TOOL>\_output.f}  &
789    generates code for $<$MODE$>$ using $<$TOOL$>$ \\
790    ~ & ~ & ~ & no {\tt make} dependencies on {\tt .F .h} \\
791    ~ & ~ & ~ & useful for compiling on remote platforms \\
792    \hline
793    (2) & {\tt <MODE><TOOL>} & {\tt <MODE>\_<TOOL>\_output.f}  &
794    generates code for $<$MODE$>$ using $<$TOOL$>$ \\
795    ~ & ~ & ~ & includes {\tt make} dependencies on {\tt .F .h} \\
796    ~ & ~ & ~ & i.e. input for $<$TOOL$>$ may be re-generated \\
797    \hline
798    (3) & {\tt <MODE>all} & {\tt mitgcmuv\_<MODE>}  &
799    generates code for $<$MODE$>$ using $<$TOOL$>$ \\
800    ~ & ~ & ~ & and compiles all code \\
801    ~ & ~ & ~ & (use of TAF is set as default) \\
802    \hline
803    \hline
804    \end{tabular}
805    }
806    \end{table}
807  %  %
808  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:  
809  %  %
810  \begin{itemize}  \begin{itemize}
811  %  %
812  \item {\it data}  \item [$<$TOOL$>$]
 %  
 \item {\it data.cost}  
 %  
 \item {\it data.ctrl}  
 %  
 \item {\it data.pkg}  
 %  
 \item {\it eedata}  
813  %  %
814  \item {\it topog.bin}  \begin{itemize}
 %  
 \item {\it windx.bin, windy.bin}  
 %  
 \item {\it salt.bin, theta.bin}  
 %  
 \item {\it SSS.bin, SST.bin}  
815  %  %
816  \item {\it pickup*}  \item {\tt TAF}
817    \item {\tt TAMC}
818  %  %
819  \end{itemize}  \end{itemize}
820  %  %
821  Finally, the file to generate the adjoint code resides in  \item [$<$MODE$>$]
 $ adjoint/ $:  
822  %  %
823  \begin{itemize}  \begin{itemize}
824  %  %
825  \item {\it makefile}  \item {\tt ad} generates the adjoint model (ADM)
826    \item {\tt ftl} generates the tangent linear model (TLM)
827    \item {\tt svd} generates both ADM and TLM for \\
828    singular value decomposition (SVD) type calculations
829  %  %
830  \end{itemize}  \end{itemize}
831  %  %
832    \end{itemize}
833    
834  Below we describe the customisations of this files which are  For example, to generate the adjoint model using TAF after routines ({\tt .F})
835  specific to this experiment.  or headers ({\tt .h}) have been modified, but without compilation,
836    type {\tt make adtaf};
837  \subsubsection{File {\it .genmakerc}}  or, to generate the tangent linear model using TAMC without
838  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}}  
839    
 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}  
840    
841  \subsubsection{File {\it SIZE.h}}  A typical full build process to generate the ADM via TAF would
842    look like follows:
843    \begin{verbatim}
844    % mkdir build
845    % cd build
846    % ../../../tools/genmake2 -mods=../code_ad
847    % make depend
848    % make adall
849    \end{verbatim}
850    
851  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}  
852    
853  Note that if the structure of the common block changes in the  \subsection{The AD build process in detail
854  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}.  
855    
856  \subsubsection{File {\it tamc.h}}  The {\tt make <MODE>all} target consists of the following procedures:
857    
858  This routine contains the dimensions for TAMC checkpointing.  \begin{enumerate}
859  %  %
860    \item
861    A header file {\tt AD\_CONFIG.h} is generated which contains a CPP option
862    on which code ought to be generated. Depending on the {\tt make} target,
863    the contents is
864  \begin{itemize}  \begin{itemize}
865    \item
866    {\tt \#define ALLOW\_ADJOINT\_RUN}
867    \item
868    {\tt \#define ALLOW\_TANGENTLINEAR\_RUN}
869    \item
870    {\tt \#define ALLOW\_ECCO\_OPTIMIZATION}
871    \end{itemize}
872  %  %
873  \item {\tt \#ifdef ALLOW\_TAMC\_CHECKPOINTING} \\  \item
874  3-level checkpointing is enabled, i.e. the timestepping  A single file {\tt <MODE>\_input\_code.f} is concatenated
875  is divided into three different levels (see Section \ref{???}).  consisting of all {\tt .f} files that are part of the list {\bf AD\_FILES}
876  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}.
877  itermediate ({\tt nchklev\_2}) timestepping loop are stored to file  %
878  (handled in {\it the\_main\_loop}).  \item
879  The innermost loop ({\tt nchklev\_1})  The AD tool is invoked with the {\bf <MODE>\_<TOOL>\_FLAGS}.
880  avoids I/O by storing all required variables  The default AD tool flags in {\tt genmake2} can be overrwritten by
881  to common blocks. This storing may also be necessary if  an {\tt adjoint\_options} file (similar to the platform-specific
882  no checkpointing is chosen  {\tt build\_options}, see Section ???.
883  (nonlinear functions, if-statements, iterative loops, ...).  The AD tool writes the resulting AD code into the file
884  In the present example the dimensions are chosen as follows: \\  {\tt <MODE>\_input\_code\_ad.f}
885  \hspace*{4ex} {\tt nchklev\_1      =  36 } \\  %
886  \hspace*{4ex} {\tt nchklev\_2      =  30 } \\  \item
887  \hspace*{4ex} {\tt nchklev\_3      =  60 } \\  A short sed script {\tt adjoint\_sed} is applied to
888  To guarantee that the checkpointing intervals span the entire  {\tt <MODE>\_input\_code\_ad.f}
889  integration period the relation \\  to reinstate {\bf myThid} into the CALL argument list of active file I/O.
890  \hspace*{4ex} {\tt nchklev\_1*nchklev\_2*nchklev\_3 $ \ge $ nTimeSteps} \\  The result is written to file {\tt <MODE>\_<TOOL>\_output.f}.
891  where {\tt nTimeSteps} is either specified in {\it data}  %
892  or computed via \\  \item
893  \hspace*{4ex} {\tt nTimeSteps = (endTime-startTime)/deltaTClock }.  All routines are compiled and an executable is generated
894  %  (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}.  
895  %  %
896  \end{itemize}  \end{enumerate}
897    
898  \subsubsection{File {\it makefile}}  \subsubsection{The list AD\_FILES and {\tt .list} files}
899    
900  This file contains all relevant paramter flags and  Not all routines are presented to the AD tool.
901  lists to run TAMC.  Routines typically hidden are diagnostics routines which
902  It is assumed that TAMC is available to you, either locally,  do not influence the cost function, but may create
903  being installed on your network, or remotely through the 'TAMC Utility'.  artificial flow dependencies such as I/O of active variables.
904  TAMC is called with the command {\tt tamc} followed by a  
905  number of options. They are described in detail in the  {\tt genmake2} generates a list (or variable) {\bf AD\_FILES}
906  TAMC manual \cite{gie:99}.  which contains all routines that are shown to the AD tool.
907  Here we briefly discuss the main flags used in the {\it makefile}  This list is put together from all files with suffix {\tt .list}
908    that {\tt genmake2} finds in its search directories.
909    The list file for the core MITgcm routines is in {\tt model/src/}
910    is called {\tt model\_ad\_diff.list}.
911    Note that no wrapper routine is shown to TAF. These are either
912    not visible at all to the AD code, or hand-written AD code
913    is available (see next section).
914    
915    Each package directory contains its package-specific
916    list file {\tt <PKG>\_ad\_diff.list}. For example,
917    {\tt pkg/ptracers/} contains the file {\tt ptracers\_ad\_diff.list}.
918    Thus, enabling a package will automatically extend the
919    {\bf AD\_FILES} list of {\tt genmake2} to incorporate the
920    package-specific routines.
921    Note that you will need to regenerate the {\tt Makefile} if
922    you enable a package (e.g. by adding it to {\tt packages.conf})
923    and a {\tt Makefile} already exists.
924    
925    \subsubsection{The list AD\_FLOW\_FILES and {\tt .flow} files}
926    
927    TAMC and TAF can evaluate user-specified directives
928    that start with a specific syntax ({\tt CADJ}, {\tt C\$TAF}, {\tt !\$TAF}).
929    The main categories of directives are STORE directives and
930    FLOW directives. Here, we are concerned with flow directives,
931    store directives are treated elsewhere.
932    
933    Flow directives enable the AD tool to evaluate how it should treat
934    routines that are 'hidden' by the user, i.e. routines which are
935    not contained in the {\bf AD\_FILES} list (see previous section),
936    but which are called in part of the code that the AD tool does see.
937    The flow directive tell the AD tool
938  %  %
939  \begin{itemize}  \begin{itemize}
 \item [{\tt tamc}] {\tt  
 -input <variable names>  
 -output <variable name> ... \\  
 -toplevel <S/R name> -reverse <file names>  
 }  
 \end{itemize}  
940  %  %
941  \begin{itemize}  \item which subroutine arguments are input/output
942  %  \item which subroutine arguments are active
943  \item {\tt -toplevel <S/R name>} \\  \item which subroutine arguments are required to compute the cost
944  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.  
945  %  %
946  \end{itemize}  \end{itemize}
947    %
948    The syntax for the flow directives can be found in the
949    AD tool manuals.
950    
951    {\tt genmake2} generates a list (or variable) {\bf AD\_FLOW\_FILES}
952    which contains all files with suffix{\tt .flow} that it finds
953    in its search directories.
954    The flow directives for the core MITgcm routines of
955    {\tt eesupp/src/} and {\tt model/src/}
956    reside in {\tt pkg/autodiff/}.
957    This directory also contains hand-written adjoint code
958    for the MITgcm WRAPPER (see Section ???).
959    
960    Flow directives for package-specific routines are contained in
961    the corresponding package directories in the file
962    {\tt <PKG>\_ad.flow}, e.g. ptracers-specific directives are in
963    {\tt ptracers\_ad.flow}.
964    
965    \subsubsection{Store directives for 3-level checkpointing}
966    
967    The storing that is required at each period of the
968    3-level checkpointing is controled by three
969    top-level headers.
970    
971  \subsubsection{File {\it data}}  \begin{verbatim}
972    do ilev_3 = 1, nchklev_3
973  \subsubsection{File {\it data.cost}}  #  include ``checkpoint_lev3.h''
974       do ilev_2 = 1, nchklev_2
975  \subsubsection{File {\it data.ctrl}}  #     include ``checkpoint_lev2.h''
976          do ilev_1 = 1, nchklev_1
977  \subsubsection{File {\it data.pkg}}  #        include ``checkpoint_lev1.h''
978    
979  \subsubsection{File {\it eedata}}  ...
980    
981  \subsubsection{File {\it topog.bin}}        end do
982       end do
983  \subsubsection{File {\it windx.bin, windy.bin}}  end do
984    \end{verbatim}
 \subsubsection{File {\it salt.bin, theta.bin}}  
985    
986  \subsubsection{File {\it SSS.bin, SST.bin}}  All files {\tt checkpoint\_lev?.h} are contained in directory
987    {\tt pkg/autodiff/}.
988    
 \subsubsection{File {\it pickup*}}  
989    
990  \subsection{Compiling the model and its adjoint}  \subsubsection{Changing the default AD tool flags: ad\_options files}
991    
 \newpage  
992    
993  %**********************************************************************  \subsubsection{Hand-written adjoint code}
 \section{TLM and ADM code generation in general}  
 \label{sec_ad_setup_gen}  
 %**********************************************************************  
994    
995  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.  
996    
997  \subsection{The cost function (dependent variable)}  \subsection{The cost function (dependent variable)
998    \label{section_cost}}
999    
1000  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}.
1001  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
1002  $ {\cal J}(\vec{u}) \, = \, {\cal J}(M(\vec{u})) $.  $ {\cal J}(\vec{u}) \, = \, {\cal J}(M(\vec{u})) $.
1003  The input is referred to as the  The input are referred to as the
1004  {\sf independent variables} or {\sf control variables}.  {\sf independent variables} or {\sf control variables}.
1005  All aspects relevant to the treatment of the cost function $ {\cal J} $  All aspects relevant to the treatment of the cost function $ {\cal J} $
1006  (parameter setting, initialisation, incrementation,  (parameter setting, initialization, accumulation,
1007  final evaluation), are controled by the package {\it pkg/cost}.  final evaluation), are controlled by the package {\it pkg/cost}.
1008    The aspects relevant to the treatment of the independent variables
1009    are controlled by the package {\it pkg/ctrl} and will be treated
1010    in the next section.
1011    
1012    \input{part5/doc_cost_flow}
1013    
1014    \subsubsection{Enabling the package}
1015    
 \subsubsection{genmake and CPP options}  
 %  
 \begin{itemize}  
 %  
 \item  
1016  \fbox{  \fbox{
1017  \begin{minipage}{12cm}  \begin{minipage}{12cm}
1018  {\it genmake}, {\it CPP\_OPTIONS.h}, {\it ECCO\_CPPOPTIONS.h}  {\it packages.conf}, {\it ECCO\_CPPOPTIONS.h}
1019  \end{minipage}  \end{minipage}
1020  }  }
1021  \end{itemize}  \begin{itemize}
 %  
 The directory {\it pkg/cost} can be included to the  
 compile list in 3 different ways (cf. Section \ref{???}):  
1022  %  %
1023  \begin{enumerate}  \item
1024    The package is enabled by adding {\it cost} to your file {\it packages.conf}
1025    (see Section ???)
1026  %  %
1027  \item {\it genmake}: \\  \item
1028  Change the default settngs in the file {\it genmake} by adding  
1029  {\bf cost} to the {\bf enable} list (not recommended).  
1030  %  \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}.  
1031  %  %
1032  \end{enumerate}  
1033  Since the cost function is usually used in conjunction with  N.B.: In general the following packages ought to be enabled
1034  automatic differentiation, the CPP option  simultaneously: {\it autodiff, cost, ctrl}.
 {\bf ALLOW\_ADJOINT\_RUN} should be defined  
 (file {\it CPP\_OPTIONS.h}).  
1035  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}.
1036  Each specific cost function contribution has its own option.  Each specific cost function contribution has its own option.
1037  For the present example the option is {\bf ALLOW\_COST\_TRACER}.  For the present example the option is {\bf ALLOW\_COST\_TRACER}.
1038  All cost-specific options are set in {\it ECCO\_CPPOPTIONS.h}  All cost-specific options are set in {\it ECCO\_CPPOPTIONS.h}
1039    Since the cost function is usually used in conjunction with
1040    automatic differentiation, the CPP option
1041    {\bf ALLOW\_ADJOINT\_RUN} (file {\it CPP\_OPTIONS.h}) and
1042    {\bf ALLOW\_AUTODIFF\_TAMC} (file {\it ECCO\_CPPOPTIONS.h})
1043    should be defined.
1044    
1045  \subsubsection{Initialisation}  \subsubsection{Initialization}
1046  %  %
1047  The initialisation of the {\it cost} package is readily enabled  The initialization of the {\it cost} package is readily enabled
1048  as soon as the CPP option {\bf ALLOW\_ADJOINT\_RUN} is defined.  as soon as the CPP option {\bf ALLOW\_COST} is defined.
1049  %  %
1050  \begin{itemize}  \begin{itemize}
1051  %  %
# Line 1152  Variables: {\it cost\_init} Line 1075  Variables: {\it cost\_init}
1075  }  }
1076  \\  \\
1077  This S/R  This S/R
1078  initialises the different cost function contributions.  initializes the different cost function contributions.
1079  The contribtion for the present example is {\bf objf\_tracer}  The contribution for the present example is {\bf objf\_tracer}
1080  which is defined on each tile (bi,bj).  which is defined on each tile (bi,bj).
1081  %  %
1082  \end{itemize}  \end{itemize}
1083  %  %
1084  \subsubsection{Incrementation}  \subsubsection{Accumulation}
1085  %  %
1086  \begin{itemize}  \begin{itemize}
1087  %  %
# Line 1196  from each contribution and sums over all Line 1119  from each contribution and sums over all
1119  \begin{equation}  \begin{equation}
1120  {\cal J} \, = \,  {\cal J} \, = \,
1121  {\rm fc} \, = \,  {\rm fc} \, = \,
1122  {\rm mult\_tracer} \sum_{bi,\,bj}^{nSx,\,nSy}  {\rm mult\_tracer} \sum_{\text{global sum}} \sum_{bi,\,bj}^{nSx,\,nSy}
1123  {\rm objf\_tracer}(bi,bj) \, + \, ...  {\rm objf\_tracer}(bi,bj) \, + \, ...
1124  \end{equation}  \end{equation}
1125  %  %
# Line 1206  The total cost function {\bf fc} will be Line 1129  The total cost function {\bf fc} will be
1129  tamc -output 'fc' ...  tamc -output 'fc' ...
1130  \end{verbatim}  \end{verbatim}
1131    
1132  \begin{figure}[t!]  %%%% \end{document}
 \input{part5/doc_ad_the_model}  
 \label{fig:adthemodel}  
 \caption{~}  
 \end{figure}  
1133    
 \begin{figure}  
1134  \input{part5/doc_ad_the_main}  \input{part5/doc_ad_the_main}
 \label{fig:adthemain}  
 \caption{~}  
 \end{figure}  
1135    
1136  \subsection{The control variables (independent variables)}  \subsection{The control variables (independent variables)
1137    \label{section_ctrl}}
1138    
1139  The control variables are a subset of the model input  The control variables are a subset of the model input
1140  (initial conditions, boundary conditions, model parameters).  (initial conditions, boundary conditions, model parameters).
1141  Here we identify them with the variable $ \vec{u} $.  Here we identify them with the variable $ \vec{u} $.
1142  All intermediate variables whose derivative w.r.t. control  All intermediate variables whose derivative w.r.t. control
1143  variables don't vanish are called {\sf active variables}.  variables do not vanish are called {\sf active variables}.
1144  All subroutines whose derivative w.r.t. the control variables  All subroutines whose derivative w.r.t. the control variables
1145  don't vanish are called {\sf active routines}.  don't vanish are called {\sf active routines}.
1146  Read and write operations from and to file can be viewed  Read and write operations from and to file can be viewed
# Line 1232  as variable assignments. Therefore, file Line 1148  as variable assignments. Therefore, file
1148  active variables are written and from which active variables  active variables are written and from which active variables
1149  are read are called {\sf active files}.  are read are called {\sf active files}.
1150  All aspects relevant to the treatment of the control variables  All aspects relevant to the treatment of the control variables
1151  (parameter setting, initialisation, perturbation)  (parameter setting, initialization, perturbation)
1152  are controled by the package {\it pkg/ctrl}.  are controlled by the package {\it pkg/ctrl}.
1153    
1154    \input{part5/doc_ctrl_flow}
1155    
1156  \subsubsection{genmake and CPP options}  \subsubsection{genmake and CPP options}
1157  %  %
# Line 1249  are controled by the package {\it pkg/ct Line 1167  are controled by the package {\it pkg/ct
1167  %  %
1168  To enable the directory to be included to the compile list,  To enable the directory to be included to the compile list,
1169  {\bf ctrl} has to be added to the {\bf enable} list in  {\bf ctrl} has to be added to the {\bf enable} list in
1170  {\it .genmakerc} (or {\it genmake} itself).  {\it .genmakerc} or in {\it genmake} itself (analogous to {\it cost}
1171    package, cf. previous section).
1172  Each control variable is enabled via its own CPP option  Each control variable is enabled via its own CPP option
1173  in {\it ECCO\_CPPOPTIONS.h}.  in {\it ECCO\_CPPOPTIONS.h}.
1174    
1175  \subsubsection{Initialisation}  \subsubsection{Initialization}
1176  %  %
1177  \begin{itemize}  \begin{itemize}
1178  %  %
# Line 1293  Two important issues related to the hand Line 1212  Two important issues related to the hand
1212  variables in the MITGCM need to be addressed.  variables in the MITGCM need to be addressed.
1213  First, in order to save memory, the control variable arrays  First, in order to save memory, the control variable arrays
1214  are not kept in memory, but rather read from file and added  are not kept in memory, but rather read from file and added
1215  to the initial (or first guess) fields.  to the initial fields during the model initialization phase.
1216  Similarly, the corresponding adjoint fields which represent  Similarly, the corresponding adjoint fields which represent
1217  the gradient of the cost function w.r.t. the control variables  the gradient of the cost function w.r.t. the control variables
1218  are written to to file.  are written to file at the end of the adjoint integration.
1219  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.
1220  control variables and the gradient, a 1-dim. {\sf control vector}  control variables and the corresponding cost gradients,
1221    a 1-dim. {\sf control vector}
1222  and {\sf gradient vector} are written to file. They contain  and {\sf gradient vector} are written to file. They contain
1223  only the wet points of the control variables and the corresponding  only the wet points of the control variables and the corresponding
1224  gradient.  gradient.
1225  This leads to a significant data compression.  This leads to a significant data compression.
1226  Furthermore, the control and the gradient vector can be passed to a  Furthermore, an option is available
1227    ({\tt ALLOW\_NONDIMENSIONAL\_CONTROL\_IO}) to
1228    non-dimensionalise the control and gradient vector,
1229    which otherwise would contain different pieces of different
1230    magnitudes and units.
1231    Finally, the control and gradient vector can be passed to a
1232  minimization routine if an update of the control variables  minimization routine if an update of the control variables
1233  is sought as part of a minimization exercise.  is sought as part of a minimization exercise.
1234    
# Line 1314  and gradient are generated and initialis Line 1239  and gradient are generated and initialis
1239    
1240  \subsubsection{Perturbation of the independent variables}  \subsubsection{Perturbation of the independent variables}
1241  %  %
1242  The dependency chain for differentiation starts  The dependency flow for differentiation w.r.t. the controls
1243  with adding a perturbation onto the the input variable,  starts with adding a perturbation onto the input variable,
1244  thus defining the independent or control variables for TAMC.  thus defining the independent or control variables for TAMC.
1245  Three classes of controls may be considered:  Three types of controls may be considered:
1246  %  %
1247  \begin{itemize}  \begin{itemize}
1248  %  %
# Line 1332  Three classes of controls may be conside Line 1257  Three classes of controls may be conside
1257  Consider as an example the initial tracer distribution  Consider as an example the initial tracer distribution
1258  {\bf tr1} as control variable.  {\bf tr1} as control variable.
1259  After {\bf tr1} has been initialised in  After {\bf tr1} has been initialised in
1260  {\it ini\_tr1} (dynamical variables including  {\it ini\_tr1} (dynamical variables such as
1261  temperature and salinity are initialised in {\it ini\_fields}),  temperature and salinity are initialised in {\it ini\_fields}),
1262  a perturbation anomaly is added to the field in S/R  a perturbation anomaly is added to the field in S/R
1263  {\it ctrl\_map\_ini}  {\it ctrl\_map\_ini}
# Line 1345  u         & = \, u_{[0]} \, + \, \Delta Line 1270  u         & = \, u_{[0]} \, + \, \Delta
1270  \end{split}  \end{split}
1271  \end{equation}  \end{equation}
1272  %  %
1273  In principle {\bf xx\_tr1} is a 3-dim. global array  {\bf xx\_tr1} is a 3-dim. global array
1274  holding the perturbation. In the case of a simple  holding the perturbation. In the case of a simple
1275  sensitivity study this array is identical to zero.  sensitivity study this array is identical to zero.
1276  However, it's specification is essential since TAMC  However, it's specification is essential in the context
1277    of automatic differentiation since TAMC
1278  treats the corresponding line in the code symbolically  treats the corresponding line in the code symbolically
1279  when determining the differentiation chain and its origin.  when determining the differentiation chain and its origin.
1280  Thus, the variable names are part of the argument list  Thus, the variable names are part of the argument list
# Line 1366  dummy variable {\bf xx\_tr1\_dummy} is i Line 1292  dummy variable {\bf xx\_tr1\_dummy} is i
1292  and an 'active read' routine of the adjoint support  and an 'active read' routine of the adjoint support
1293  package {\it pkg/autodiff} is invoked.  package {\it pkg/autodiff} is invoked.
1294  The read-procedure is tagged with the variable  The read-procedure is tagged with the variable
1295  {\bf xx\_tr1\_dummy} enabbling TAMC to recognize the  {\bf xx\_tr1\_dummy} enabling TAMC to recognize the
1296  initialisation of the perturbation.  initialization of the perturbation.
1297  The modified call of TAMC thus reads  The modified call of TAMC thus reads
1298  %  %
1299  \begin{verbatim}  \begin{verbatim}
# Line 1386  in the code takes on the form Line 1312  in the code takes on the form
1312  %  %
1313  Note, that reading an active variable corresponds  Note, that reading an active variable corresponds
1314  to a variable assignment. Its derivative corresponds  to a variable assignment. Its derivative corresponds
1315  to a write statement of the adjoint variable.  to a write statement of the adjoint variable, followed by
1316    a reset.
1317  The 'active file' routines have been designed  The 'active file' routines have been designed
1318  to support active read and corresponding active write  to support active read and corresponding adjoint active write
1319  operations.  operations (and vice versa).
1320  %  %
1321  \item  \item
1322  \fbox{  \fbox{
# Line 1406  with the symbolic perturbation taking pl Line 1333  with the symbolic perturbation taking pl
1333  Note however an important difference:  Note however an important difference:
1334  Since the boundary values are time dependent with a new  Since the boundary values are time dependent with a new
1335  forcing field applied at each time steps,  forcing field applied at each time steps,
1336  the general problem may be be thought of as  the general problem may be thought of as
1337  a new control variable at each time step, i.e.  a new control variable at each time step
1338    (or, if the perturbation is averaged over a certain period,
1339    at each $ N $ timesteps), i.e.
1340  \[  \[
1341  u_{\rm forcing} \, = \,  u_{\rm forcing} \, = \,
1342  \{ \, u_{\rm forcing} ( t_n ) \, \}_{  \{ \, u_{\rm forcing} ( t_n ) \, \}_{
# Line 1432  calendar ({\it cal}~) and external forci Line 1361  calendar ({\it cal}~) and external forci
1361  %  %
1362  This routine is not yet implemented, but would proceed  This routine is not yet implemented, but would proceed
1363  proceed along the same lines as the initial value sensitivity.  proceed along the same lines as the initial value sensitivity.
1364    The mixing parameters {\bf diffkr} and {\bf kapgm}
1365    are currently added as controls in {\it ctrl\_map\_ini.F}.
1366  %  %
1367  \end{itemize}  \end{itemize}
1368  %  %
1369    
1370  \subsubsection{Output of adjoint variables and gradient}  \subsubsection{Output of adjoint variables and gradient}
1371  %  %
1372  Two ways exist to generate output of adjoint fields.  Several ways exist to generate output of adjoint fields.
1373  %  %
1374  \begin{itemize}  \begin{itemize}
1375  %  %
1376  \item  \item
1377  \fbox{  \fbox{
1378  \begin{minipage}{12cm}  \begin{minipage}{12cm}
1379  {\it ctrl\_pack}:  {\it ctrl\_map\_ini, ctrl\_map\_forcing}:
1380  \end{minipage}  \end{minipage}
1381  }  }
1382  \\  \\
 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:  
 %  
1383  \begin{itemize}  \begin{itemize}
1384  %  %
1385  \item {\bf xx\_...}: the control variable fields  \item {\bf xx\_...}: the control variable fields \\
1386    Before the forward integration, the control
1387    variables are read from file {\bf xx\_ ...} and added to
1388    the model field.
1389  %  %
1390  \item {\bf adxx\_...}: the adjoint variable fields, i.e. the gradient  \item {\bf adxx\_...}: the adjoint variable fields, i.e. the gradient
1391  $ \nabla _{u}{\cal J} $ for each control variable,  $ \nabla _{u}{\cal J} $ for each control variable \\
1392    After the adjoint integration the corresponding adjoint
1393    variables are written to {\bf adxx\_ ...}.
1394    %
1395    \end{itemize}
1396  %  %
1397  \item {\bf vector\_ctrl}: the control vector  \item
1398    \fbox{
1399    \begin{minipage}{12cm}
1400    {\it ctrl\_unpack, ctrl\_pack}:
1401    \end{minipage}
1402    }
1403    \\
1404    %
1405    \begin{itemize}
1406  %  %
1407  \item {\bf vector\_grad}: the gradient vector  \item {\bf vector\_ctrl}: the control vector \\
1408    At the very beginning of the model initialization,
1409    the updated compressed control vector is read (or initialised)
1410    and distributed to 2-dim. and 3-dim. control variable fields.
1411    %
1412    \item {\bf vector\_grad}: the gradient vector \\
1413    At the very end of the adjoint integration,
1414    the 2-dim. and 3-dim. adjoint variables are read,
1415    compressed to a single vector and written to file.
1416  %  %
1417  \end{itemize}  \end{itemize}
1418  %  %
# Line 1474  $ \nabla _{u}{\cal J} $ for each control Line 1424  $ \nabla _{u}{\cal J} $ for each control
1424  }  }
1425  \\  \\
1426  In addition to writing the gradient at the end of the  In addition to writing the gradient at the end of the
1427  forward/adjoint integration, many more adjoint variables,  forward/adjoint integration, many more adjoint variables
1428  representing the Lagrange multipliers of the model state  of the model state
1429  w.r.t. the model state  at intermediate times can be written using S/R
 at different times can be written using S/R  
1430  {\it addummy\_in\_stepping}.  {\it addummy\_in\_stepping}.
1431  This routine is part of the adjoint support package  This routine is part of the adjoint support package
1432  {\it pkg/autodiff} (cf.f. below).  {\it pkg/autodiff} (cf.f. below).
1433    The procedure is enabled using via the CPP-option
1434    {\bf ALLOW\_AUTODIFF\_MONITOR} (file {\it ECCO\_CPPOPTIONS.h}).
1435  To be part of the adjoint code, the corresponding S/R  To be part of the adjoint code, the corresponding S/R
1436  {\it dummy\_in\_stepping} has to be called in the forward  {\it dummy\_in\_stepping} has to be called in the forward
1437  model (S/R {\it the\_main\_loop}) at the appropriate place.  model (S/R {\it the\_main\_loop}) at the appropriate place.
1438    The adjoint common blocks are extracted from the adjoint code
1439    via the header file {\it adcommon.h}.
1440    
1441  {\it dummy\_in\_stepping} is essentially empty,  {\it dummy\_in\_stepping} is essentially empty,
1442  the corresponding adjoint routine is hand-written rather  the corresponding adjoint routine is hand-written rather
# Line 1491  than generated automatically. Line 1444  than generated automatically.
1444  Appropriate flow directives ({\it dummy\_in\_stepping.flow})  Appropriate flow directives ({\it dummy\_in\_stepping.flow})
1445  ensure that TAMC does not automatically  ensure that TAMC does not automatically
1446  generate {\it addummy\_in\_stepping} by trying to differentiate  generate {\it addummy\_in\_stepping} by trying to differentiate
1447  {\it dummy\_in\_stepping}, but rather takes the hand-written routine.  {\it dummy\_in\_stepping}, but instead refers to
1448    the hand-written routine.
1449    
1450  {\it dummy\_in\_stepping} is called in the forward code  {\it dummy\_in\_stepping} is called in the forward code
1451  at the beginning of each  at the beginning of each
# Line 1501  each timestep in the adjoint calculation Line 1455  each timestep in the adjoint calculation
1455  {\it addynamics}.  {\it addynamics}.
1456    
1457  {\it addummy\_in\_stepping} includes the header files  {\it addummy\_in\_stepping} includes the header files
1458  {\it adffields.h, addynamics.h, adtr1.h}.  {\it adcommon.h}.
1459  These header files are also hand-written. They contain  This header file is also hand-written. It contains
1460  the common blocks {\bf /addynvars\_r/}, {\bf /addynvars\_cd/},  the common blocks
1461    {\bf /addynvars\_r/}, {\bf /addynvars\_cd/},
1462    {\bf /addynvars\_diffkr/}, {\bf /addynvars\_kapgm/},
1463  {\bf /adtr1\_r/}, {\bf /adffields/},  {\bf /adtr1\_r/}, {\bf /adffields/},
1464  which have been extracted from the adjoint code to enable  which have been extracted from the adjoint code to enable
1465  access to the adjoint variables.  access to the adjoint variables.
1466    
1467    {\bf WARNING:} If the structure of the common blocks
1468    {\bf /dynvars\_r/}, {\bf /dynvars\_cd/}, etc., changes
1469    similar changes will occur in the adjoint common blocks.
1470    Therefore, consistency between the TAMC-generated common blocks
1471    and those in {\it adcommon.h} have to be checked.
1472  %  %
1473  \end{itemize}  \end{itemize}
1474    
# Line 1521  The gradient $ \nabla _{u}{\cal J} |_{u_ Line 1483  The gradient $ \nabla _{u}{\cal J} |_{u_
1483  with the value of the cost function itself $ {\cal J}(u_{[k]}) $  with the value of the cost function itself $ {\cal J}(u_{[k]}) $
1484  at iteration step $ k $ serve  at iteration step $ k $ serve
1485  as input to a minimization routine (e.g. quasi-Newton method,  as input to a minimization routine (e.g. quasi-Newton method,
1486  conjugate gradient, ...) to compute an update in the  conjugate gradient, ... \cite{gil-lem:89})
1487    to compute an update in the
1488  control variable for iteration step $k+1$  control variable for iteration step $k+1$
1489  \[  \[
1490  u_{[k+1]} \, = \,  u_{[0]} \, + \, \Delta u_{[k+1]}  u_{[k+1]} \, = \,  u_{[0]} \, + \, \Delta u_{[k+1]}
# Line 1535  Tab. \ref{???} sketches the flow between Line 1498  Tab. \ref{???} sketches the flow between
1498  and the minimization routine.  and the minimization routine.
1499    
1500  \begin{eqnarray*}  \begin{eqnarray*}
1501  \footnotesize  \scriptsize
1502  \begin{array}{ccccc}  \begin{array}{ccccc}
1503  u_{[0]} \,\, ,  \,\, \Delta u_{[k]}    & ~ & ~ & ~ & ~ \\  u_{[0]} \,\, ,  \,\, \Delta u_{[k]}    & ~ & ~ & ~ & ~ \\
1504  {\Big\downarrow}  {\Big\downarrow}
# Line 1552  v_{[k]} = M \left( u_{[k]} \right) & Line 1515  v_{[k]} = M \left( u_{[k]} \right) &
1515  {\cal J}_{[k]} = {\cal J} \left( M \left( u_{[k]} \right) \right)} \\  {\cal J}_{[k]} = {\cal J} \left( M \left( u_{[k]} \right) \right)} \\
1516  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1517  \hline  \hline
1518    \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~}  \\
1519    \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{{\Big\downarrow}} \\
1520    \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~}  \\
1521  \hline  \hline
1522  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1523  \multicolumn{1}{|c}{  \multicolumn{1}{|c}{
1524  \nabla_u {\cal J}_{[k]} (\delta {\cal J}) =  \nabla_u {\cal J}_{[k]} (\delta {\cal J}) =
1525  T\!\!^{\ast} \cdot \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J})} &  T^{\ast} \cdot \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J})} &
1526  \stackrel{\bf adjoint}{\mathbf \longleftarrow} &  \stackrel{\bf adjoint}{\mathbf \longleftarrow} &
1527  ad \, v_{[k]} (\delta {\cal J}) =  ad \, v_{[k]} (\delta {\cal J}) =
1528  \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J}) &  \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J}) &
# Line 1565  ad \, v_{[k]} (\delta {\cal J}) = Line 1531  ad \, v_{[k]} (\delta {\cal J}) =
1531  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
1532  \hline  \hline
1533   ~ & ~ & ~ & ~ & ~ \\   ~ & ~ & ~ & ~ & ~ \\
1534  ~ & ~ &  \hspace*{15ex}{\Bigg\downarrow}  
1535  {\cal J}_{[k]} \qquad {\Bigg\downarrow}  \qquad \nabla_u {\cal J}_{[k]}  \quad {\cal J}_{[k]}, \quad \nabla_u {\cal J}_{[k]}
1536   & ~ & ~ \\   & ~ & ~ & ~ & ~ \\
1537   ~ & ~ & ~ & ~ & ~ \\   ~ & ~ & ~ & ~ & ~ \\
1538  \hline  \hline
1539  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\  \multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\
# Line 1595  The corresponding I/O flow looks as foll Line 1561  The corresponding I/O flow looks as foll
1561    
1562  \vspace*{0.5cm}  \vspace*{0.5cm}
1563    
1564    {\scriptsize
1565  \begin{tabular}{ccccc}  \begin{tabular}{ccccc}
1566  {\bf vector\_ctrl\_$<$k$>$ } & ~ & ~ & ~ & ~ \\  {\bf vector\_ctrl\_$<$k$>$ } & ~ & ~ & ~ & ~ \\
1567  {\big\downarrow}  & ~ & ~ & ~ & ~ \\  {\big\downarrow}  & ~ & ~ & ~ & ~ \\
# Line 1605  The corresponding I/O flow looks as foll Line 1572  The corresponding I/O flow looks as foll
1572  \cline{3-3}  \cline{3-3}
1573  \multicolumn{1}{l}{\bf xx\_theta0...$<$k$>$} & ~ &  \multicolumn{1}{l}{\bf xx\_theta0...$<$k$>$} & ~ &
1574  \multicolumn{1}{|c|}{~} & ~ & ~ \\  \multicolumn{1}{|c|}{~} & ~ & ~ \\
1575  \multicolumn{1}{l}{\bf xx\_salt0...$<$k$>$} & $\longrightarrow$ &  \multicolumn{1}{l}{\bf xx\_salt0...$<$k$>$} &
1576    $\stackrel{\mbox{read}}{\longrightarrow}$ &
1577  \multicolumn{1}{|c|}{forward integration} & ~ & ~ \\  \multicolumn{1}{|c|}{forward integration} & ~ & ~ \\
1578  \multicolumn{1}{l}{\bf \vdots} & ~ & \multicolumn{1}{|c|}{~}    \multicolumn{1}{l}{\bf \vdots} & ~ & \multicolumn{1}{|c|}{~}  
1579  & ~ & ~ \\  & ~ & ~ \\
1580  \cline{3-3}  \cline{3-3}
1581  ~ & ~ & ~ & ~ & ~ \\  ~ & ~ & $\downarrow$ & ~ & ~ \\
1582  \cline{3-3}  \cline{3-3}
1583  ~ & ~ &  ~ & ~ &
1584  \multicolumn{1}{|c|}{~} & ~ &  \multicolumn{1}{|c|}{~} & ~ &
1585  \multicolumn{1}{l}{\bf adxx\_theta0...$<$k$>$}  \\  \multicolumn{1}{l}{\bf adxx\_theta0...$<$k$>$}  \\
1586  ~ & ~ & \multicolumn{1}{|c|}{adjoint integration} &  ~ & ~ & \multicolumn{1}{|c|}{adjoint integration} &
1587  $\longrightarrow$ &  $\stackrel{\mbox{write}}{\longrightarrow}$ &
1588  \multicolumn{1}{l}{\bf adxx\_salt0...$<$k$>$} \\  \multicolumn{1}{l}{\bf adxx\_salt0...$<$k$>$} \\
1589  ~ & ~ & \multicolumn{1}{|c|}{~}    ~ & ~ & \multicolumn{1}{|c|}{~}  
1590  & ~ & \multicolumn{1}{l}{\bf \vdots} \\  & ~ & \multicolumn{1}{l}{\bf \vdots} \\
# Line 1628  $\longrightarrow$ & Line 1596  $\longrightarrow$ &
1596  ~ & ~ & ~ & ~ &  {\big\downarrow} \\  ~ & ~ & ~ & ~ &  {\big\downarrow} \\
1597  ~ & ~ & ~ & ~ &  {\bf vector\_grad\_$<$k$>$ } \\  ~ & ~ & ~ & ~ &  {\bf vector\_grad\_$<$k$>$ } \\
1598  \end{tabular}  \end{tabular}
1599    }
1600    
1601  \vspace*{0.5cm}  \vspace*{0.5cm}
1602    
1603    
1604  {\it ctrl\_unpack} reads in the updated control vector  {\it ctrl\_unpack} reads the updated control vector
1605  {\bf vector\_ctrl\_$<$k$>$}.  {\bf vector\_ctrl\_$<$k$>$}.
1606  It distributes the different control variables to  It distributes the different control variables to
1607  2-dim. and 3-dim. files {\it xx\_...$<$k$>$}.  2-dim. and 3-dim. files {\it xx\_...$<$k$>$}.
1608  During the forward integration the control variables  At the start of the forward integration the control variables
1609  are read from {\it xx\_...$<$k$>$}.  are read from {\it xx\_...$<$k$>$} and added to the
1610  Correspondingly, the adjoint fields are written  field.
1611    Correspondingly, at the end of the adjoint integration
1612    the adjoint fields are written
1613  to {\it adxx\_...$<$k$>$}, again via the active file routines.  to {\it adxx\_...$<$k$>$}, again via the active file routines.
1614  Finally, {\it ctrl\_pack} collects all adjoint field files  Finally, {\it ctrl\_pack} collects all adjoint files
1615  and writes them to the compressed vector file  and writes them to the compressed vector file
1616  {\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.1  
changed lines
  Added in v.1.19

  ViewVC Help
Powered by ViewVC 1.1.22