/[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.4 by heimbach, Fri Oct 5 22:22:20 2001 UTC revision 1.22 by jmc, Fri Aug 27 13:09:40 2010 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 partial derivatives from an
9  partial derivatives from an existing fully non-linear prognostic code.  existing fully non-linear prognostic code.  (see \cite{gri:00}).  A
10  (see \cite{gri:00}).  software tool is used that parses and transforms source files
11  A software tool is used that parses and transforms source files  according to a set of linguistic and mathematical rules.  AD tools are
12  according to a set of linguistic and mathematical rules.  like source-to-source translators in that they parse a program code as
13  AD tools are like source-to-source translators in that  input and produce a new program code as output
14  they parse a program code as input and produce a new program code  (we restrict our discussion to source-to-source tools, ignoring
15  as output.  operator-overloading tools).  However, unlike a
16  However, unlike a pure source-to-source translation, the output program  pure source-to-source translation, the output program represents a new
17  represents a new algorithm, such as the evaluation of the  algorithm, such as the evaluation of the Jacobian, the Hessian, or
18  Jacobian, the Hessian, or higher derivative operators.  higher derivative operators.  In principle, a variety of derived
19  In principle, a variety of derived algorithms  algorithms can be generated automatically in this way.
20  can be generated automatically in this way.  
21    MITgcm has been adapted for use with the Tangent linear and Adjoint
22  The MITGCM has been adapted for use with the  Model Compiler (TAMC) and its successor TAF (Transformation of
23  Tangent linear and Adjoint Model Compiler (TAMC) and its successor TAF  Algorithms in Fortran), developed by Ralf Giering (\cite{gie-kam:98},
24  (Transformation of Algorithms in Fortran), developed  \cite{gie:99,gie:00}).  The first application of the adjoint of MITgcm
25  by Ralf Giering (\cite{gie-kam:98}, \cite{gie:99,gie:00}).  for sensitivity studies has been published by \cite{maro-eta:99}.
26  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
27  studies has been published by \cite{maro-eta:99}.  state estimation studies.  In the following we shall refer to TAMC and
28  \cite{sta-eta:97,sta-eta:01} use the MITGCM and its adjoint  TAF synonymously, except were explicitly stated otherwise.
29  for ocean state estimation studies.  
30  In the following we shall refer to TAMC and TAF synonymously,  As of mid-2007 we are also able to generate fairly efficient
31  except were explicitly stated otherwise.  adjoint code of the MITgcm using a new, open-source AD tool,
32    called OpenAD (see \cite{naum-etal:06,utke-etal:08}.
33  TAMC exploits the chain rule for computing the first  This enables us for the first time to compare adjoint models
34  derivative of a function with  generated from different AD tools, providing an additional
35  respect to a set of input variables.  accuracy check, complementary to finite-difference gradient checks.
36  Treating a given forward code as a composition of operations --  OpenAD and its application to  MITgcm is described in detail
37  each line representing a compositional element, the chain rule is  in section \ref{sec_ad_openad}.
38  rigorously applied to the code, line by line. The resulting  
39  tangent linear or adjoint code,  The AD tool exploits the chain rule for computing the first derivative of a
40  then, may be thought of as the composition in  function with respect to a set of input variables.  Treating a given
41  forward or reverse order, respectively, of the  forward code as a composition of operations -- each line representing
42  Jacobian matrices of the forward code's compositional elements.  a compositional element, the chain rule is rigorously applied to the
43    code, line by line. The resulting tangent linear or adjoint code,
44    then, may be thought of as the composition in forward or reverse
45    order, respectively, of the Jacobian matrices of the forward code's
46    compositional elements.
47    
48  %**********************************************************************  %**********************************************************************
49  \section{Some basic algebra}  \section{Some basic algebra}
50  \label{sec_ad_algebra}  \label{sec_ad_algebra}
51    \begin{rawhtml}
52    <!-- CMIREDIR:sec_ad_algebra: -->
53    \end{rawhtml}
54  %**********************************************************************  %**********************************************************************
55    
56  Let $ \cal{M} $ be a general nonlinear, model, i.e. a  Let $ \cal{M} $ be a general nonlinear, model, i.e. a
# Line 52  $\vec{u}=(u_1,\ldots,u_m)$ Line 61  $\vec{u}=(u_1,\ldots,u_m)$
61  such as forcing functions) to the $n$-dimensional space  such as forcing functions) to the $n$-dimensional space
62  $V \subset I\!\!R^n$ of  $V \subset I\!\!R^n$ of
63  model output variable $\vec{v}=(v_1,\ldots,v_n)$  model output variable $\vec{v}=(v_1,\ldots,v_n)$
64  (model state, model diagnostcs, objective function, ...)  (model state, model diagnostics, objective function, ...)
65  under consideration,  under consideration,
66  %  %
67  \begin{equation}  \begin{equation}
# Line 220  model integration, Line 229  model integration,
229  starting at step 0 and moving up to step $\Lambda$, with intermediate  starting at step 0 and moving up to step $\Lambda$, with intermediate
230  ${\cal M}_{\lambda} (\vec{u}) = \vec{v}^{(\lambda+1)}$ and final  ${\cal M}_{\lambda} (\vec{u}) = \vec{v}^{(\lambda+1)}$ and final
231  ${\cal M}_{\Lambda} (\vec{u}) = \vec{v}^{(\Lambda+1)} = \vec{v}$.  ${\cal M}_{\Lambda} (\vec{u}) = \vec{v}^{(\Lambda+1)} = \vec{v}$.
232  Let ${\cal J}$ be a cost funciton which explicitly depends on the  Let ${\cal J}$ be a cost function which explicitly depends on the
233  final state $\vec{v}$ only  final state $\vec{v}$ only
234  (this restriction is for clarity reasons only).  (this restriction is for clarity reasons only).
235  %  %
# Line 301  We note in passing that that the $\delta Line 310  We note in passing that that the $\delta
310  are the Lagrange multipliers of the model equations which determine  are the Lagrange multipliers of the model equations which determine
311  $ \vec{v}^{(\lambda)}$.  $ \vec{v}^{(\lambda)}$.
312    
313  In coponents, eq. (\ref{adjoint}) reads as follows.  In components, eq. (\ref{adjoint}) reads as follows.
314  Let  Let
315  \[  \[
316  \begin{array}{rclcrcl}  \begin{array}{rclcrcl}
# Line 322  Let Line 331  Let
331  \end{array}  \end{array}
332  \]  \]
333  denote the perturbations in $\vec{u}$ and $\vec{v}$, respectively,  denote the perturbations in $\vec{u}$ and $\vec{v}$, respectively,
334  and their adjoint varaiables;  and their adjoint variables;
335  further  further
336  \[  \[
337  M \, = \, \left(  M \, = \, \left(
# Line 468  variables $u$ Line 477  variables $u$
477  {\it all} intermediate states $ \vec{v}^{(\lambda)} $) are sought.  {\it all} intermediate states $ \vec{v}^{(\lambda)} $) are sought.
478  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
479  $ \partial {\cal J} / \partial u_{i} $ in (\ref{forward})  $ \partial {\cal J} / \partial u_{i} $ in (\ref{forward})
480  a forward calulation has to be performed for each component seperately,  a forward calculation has to be performed for each component separately,
481  i.e. $ \delta \vec{u} = \delta u_{i} {\vec{e}_{i}} $  i.e. $ \delta \vec{u} = \delta u_{i} {\vec{e}_{i}} $
482  for  the $i$-th forward calculation.  for  the $i$-th forward calculation.
483  Then, (\ref{forward}) represents the  Then, (\ref{forward}) represents the
# Line 487  M^T \left( \nabla_v {\cal J}^T \left(\de Line 496  M^T \left( \nabla_v {\cal J}^T \left(\de
496  \nabla_u {\cal J}^T \cdot \delta \vec{J}  \nabla_u {\cal J}^T \cdot \delta \vec{J}
497  \]  \]
498  where now $ \delta \vec{J} \in I\!\!R^l $ is a vector of  where now $ \delta \vec{J} \in I\!\!R^l $ is a vector of
499  dimenison $ l $.  dimension $ l $.
500  In this case $ l $ reverse simulations have to be performed  In this case $ l $ reverse simulations have to be performed
501  for each $ \delta J_{k}, \,\, k = 1, \ldots, l $.  for each $ \delta J_{k}, \,\, k = 1, \ldots, l $.
502  Then, the reverse mode is more efficient as long as  Then, the reverse mode is more efficient as long as
503  $ l < n $, otherwise the forward mode is preferable.  $ l < n $, otherwise the forward mode is preferable.
504  Stricly, the reverse mode is called adjoint mode only for  Strictly, the reverse mode is called adjoint mode only for
505  $ l = 1 $.  $ l = 1 $.
506    
507  A detailed analysis of the underlying numerical operations  A detailed analysis of the underlying numerical operations
# Line 557  Because of the local character of the de Line 566  Because of the local character of the de
566  (a derivative is defined w.r.t. a point along the trajectory),  (a derivative is defined w.r.t. a point along the trajectory),
567  the intermediate results of the model trajectory  the intermediate results of the model trajectory
568  $\vec{v}^{(\lambda+1)}={\cal M}_{\lambda}(v^{(\lambda)})$  $\vec{v}^{(\lambda+1)}={\cal M}_{\lambda}(v^{(\lambda)})$
569  are needed to evaluate the intermediate Jacobian  may be required to evaluate the intermediate Jacobian
570  $M_{\lambda}|_{\vec{v}^{(\lambda)}} \, \delta \vec{v}^{(\lambda)} $.  $M_{\lambda}|_{\vec{v}^{(\lambda)}} \, \delta \vec{v}^{(\lambda)} $.
571    This is the case e.g. for nonlinear expressions
572    (momentum advection, nonlinear equation of state), state-dependent
573    conditional statements (parameterization schemes).
574  In the forward mode, the intermediate results are required  In the forward mode, the intermediate results are required
575  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}$,
576  but in the reverse mode they are required in the reverse order.  but in the reverse mode they are required in the reverse order.
# Line 569  point of evaluation has to be recomputed Line 581  point of evaluation has to be recomputed
581    
582  A method to balance the amount of recomputations vs.  A method to balance the amount of recomputations vs.
583  storage requirements is called {\sf checkpointing}  storage requirements is called {\sf checkpointing}
584  (e.g. \cite{res-eta:98}).  (e.g. \cite{gri:92}, \cite{res-eta:98}).
585  It is depicted in \reffig{3levelcheck} for a 3-level checkpointing  It is depicted in \ref{fig:3levelcheck} for a 3-level checkpointing
586  [as an example, we give explicit numbers for a 3-day  [as an example, we give explicit numbers for a 3-day
587  integration with a 1-hourly timestep in square brackets].  integration with a 1-hourly timestep in square brackets].
588  \begin{itemize}  \begin{itemize}
# Line 580  In a first step, the model trajectory is Line 592  In a first step, the model trajectory is
592  $ {n}^{lev3} $ subsections [$ {n}^{lev3} $=3 1-day intervals],  $ {n}^{lev3} $ subsections [$ {n}^{lev3} $=3 1-day intervals],
593  with the label $lev3$ for this outermost loop.  with the label $lev3$ for this outermost loop.
594  The model is then integrated along the full trajectory,  The model is then integrated along the full trajectory,
595  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
596  [i.e. 3 times, at  [i.e. 3 times, at
597  $ i = 0,1,2 $ corresponding to $ k_{i}^{lev3} = 0, 24, 48 $].  $ i = 0,1,2 $ corresponding to $ k_{i}^{lev3} = 0, 24, 48 $].
598    In addition, the cost function is computed, if needed.
599  %  %
600  \item [$lev2$]  \item [$lev2$]
601  In a second step each subsection itself is divided into  In a second step each subsection itself is divided into
602  $ {n}^{lev2} $ sub-subsections  $ {n}^{lev2} $ subsections
603  [$ {n}^{lev2} $=4 6-hour intervals per subsection].  [$ {n}^{lev2} $=4 6-hour intervals per subsection].
604  The model picks up at the last outermost dumped state  The model picks up at the last outermost dumped state
605  $ v_{k_{n}^{lev3}} $ and is integrated forward in time along  $ v_{k_{n}^{lev3}} $ and is integrated forward in time along
606  the last subsection, with the label $lev2$ for this    the last subsection, with the label $lev2$ for this  
607  intermediate loop.  intermediate loop.
608  The model state is now stored at every $ k_{i}^{lev2} $-th  The model state is now stored to disk at every $ k_{i}^{lev2} $-th
609  timestep  timestep
610  [i.e. 4 times, at  [i.e. 4 times, at
611  $ 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 $].
# Line 600  $ i = 0,1,2,3 $ corresponding to $ k_{i} Line 613  $ i = 0,1,2,3 $ corresponding to $ k_{i}
613  \item [$lev1$]  \item [$lev1$]
614  Finally, the model picks up at the last intermediate dump state  Finally, the model picks up at the last intermediate dump state
615  $ v_{k_{n}^{lev2}} $ and is integrated forward in time along  $ v_{k_{n}^{lev2}} $ and is integrated forward in time along
616  the last sub-subsection, with the label $lev1$ for this    the last subsection, with the label $lev1$ for this  
617  intermediate loop.  intermediate loop.
618  Within this sub-subsection only, the model state is stored  Within this sub-subsection only, parts of the model state is stored
619  at every timestep  to memory at every timestep
620  [i.e. every hour $ i=0,...,5$ corresponding to  [i.e. every hour $ i=0,...,5$ corresponding to
621  $ k_{i}^{lev1} = 66, 67, \ldots, 71 $].  $ k_{i}^{lev1} = 66, 67, \ldots, 71 $].
622  Thus, the  final state $ v_n = v_{k_{n}^{lev1}} $ is reached  The  final state $ v_n = v_{k_{n}^{lev1}} $ is reached
623  and the model state of all peceeding timesteps along the last  and the model state of all preceding timesteps along the last
624  sub-subsections are available, enabling integration backwards  innermost subsection are available, enabling integration backwards
625  in time along the last sub-subsection.  in time along the last subsection.
626  Thus, the adjoint can be computed along this last  The adjoint can thus be computed along this last
627  sub-subsection $k_{n}^{lev2}$.  subsection $k_{n}^{lev2}$.
628  %  %
629  \end{itemize}  \end{itemize}
630  %  %
631  This procedure is repeated consecutively for each previous  This procedure is repeated consecutively for each previous
632  sub-subsection $k_{n-1}^{lev2}, \ldots, k_{1}^{lev2} $  subsection $k_{n-1}^{lev2}, \ldots, k_{1}^{lev2} $
633  carrying the adjoint computation to the initial time  carrying the adjoint computation to the initial time
634  of the subsection $k_{n}^{lev3}$.  of the subsection $k_{n}^{lev3}$.
635  Then, the procedure is repeated for the previous subsection  Then, the procedure is repeated for the previous subsection
# Line 627  $k_{1}^{lev3}$. Line 640  $k_{1}^{lev3}$.
640  For the full model trajectory of  For the full model trajectory of
641  $ n^{lev3} \cdot n^{lev2} \cdot n^{lev1} $ timesteps  $ n^{lev3} \cdot n^{lev2} \cdot n^{lev1} $ timesteps
642  the required storing of the model state was significantly reduced to  the required storing of the model state was significantly reduced to
643  $ n^{lev1} + n^{lev2} + n^{lev3} $  $ n^{lev2} + n^{lev3} $ to disk and roughly $ n^{lev1} $ to memory
644  [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
645  the model state was stored 13 times].  the model state was stored 7 times to disk and roughly 6 times
646    to memory].
647  This saving in memory comes at a cost of a required  This saving in memory comes at a cost of a required
648  3 full forward integrations of the model (one for each  3 full forward integrations of the model (one for each
649  checkpointing level).  checkpointing level).
650  The balance of storage vs. recomputation certainly depends  The optimal balance of storage vs. recomputation certainly depends
651  on the computing resources available.  on the computing resources available and may be adjusted by
652    adjusting the partitioning among the
653    $ n^{lev3}, \,\, n^{lev2}, \,\, n^{lev1} $.
654    
655  \begin{figure}[t!]  \begin{figure}[t!]
656  \centering  \begin{center}
657  %\psdraft  %\psdraft
658  \psfrag{v_k1^lev3}{\mathinfigure{v_{k_{1}^{lev3}}}}  %\psfrag{v_k1^lev3}{\mathinfigure{v_{k_{1}^{lev3}}}}
659  \psfrag{v_kn-1^lev3}{\mathinfigure{v_{k_{n-1}^{lev3}}}}  %\psfrag{v_kn-1^lev3}{\mathinfigure{v_{k_{n-1}^{lev3}}}}
660  \psfrag{v_kn^lev3}{\mathinfigure{v_{k_{n}^{lev3}}}}  %\psfrag{v_kn^lev3}{\mathinfigure{v_{k_{n}^{lev3}}}}
661  \psfrag{v_k1^lev2}{\mathinfigure{v_{k_{1}^{lev2}}}}  %\psfrag{v_k1^lev2}{\mathinfigure{v_{k_{1}^{lev2}}}}
662  \psfrag{v_kn-1^lev2}{\mathinfigure{v_{k_{n-1}^{lev2}}}}  %\psfrag{v_kn-1^lev2}{\mathinfigure{v_{k_{n-1}^{lev2}}}}
663  \psfrag{v_kn^lev2}{\mathinfigure{v_{k_{n}^{lev2}}}}  %\psfrag{v_kn^lev2}{\mathinfigure{v_{k_{n}^{lev2}}}}
664  \psfrag{v_k1^lev1}{\mathinfigure{v_{k_{1}^{lev1}}}}  %\psfrag{v_k1^lev1}{\mathinfigure{v_{k_{1}^{lev1}}}}
665  \psfrag{v_kn^lev1}{\mathinfigure{v_{k_{n}^{lev1}}}}  %\psfrag{v_kn^lev1}{\mathinfigure{v_{k_{n}^{lev1}}}}
666  \mbox{\epsfig{file=part5/checkpointing.eps, width=0.8\textwidth}}  %\mbox{\epsfig{file=s_autodiff/figs/checkpointing.eps, width=0.8\textwidth}}
667    \resizebox{5.5in}{!}{\includegraphics{s_autodiff/figs/checkpointing.eps}}
668  %\psfull  %\psfull
669  \caption  \end{center}
670  {Schematic view of intermediate dump and restart for  \caption{
671    Schematic view of intermediate dump and restart for
672  3-level checkpointing.}  3-level checkpointing.}
673  \label{fig:3levelcheck}  \label{fig:3levelcheck}
674  \end{figure}  \end{figure}
# Line 662  on the computing resources available. Line 680  on the computing resources available.
680  % \subsection{Error covariance estimate and Hessian matrix}  % \subsection{Error covariance estimate and Hessian matrix}
681  % \label{sec_hessian}  % \label{sec_hessian}
682    
683  \newpage  \newpage
684    
685  %**********************************************************************  %**********************************************************************
686  \section{AD-specific setup by example: sensitivity of carbon sequestration}  \section{TLM and ADM generation in general}
687  \label{sec_ad_setup_ex}  \label{sec_ad_setup_gen}
688    \begin{rawhtml}
689    <!-- CMIREDIR:sec_ad_setup_gen: -->
690    \end{rawhtml}
691  %**********************************************************************  %**********************************************************************
692    
693  The MITGCM has been adapted to enable AD using TAMC or TAF.  In this section we describe in a general fashion
694  The present description, therefore, is specific to the  the parts of the code that are relevant for automatic
695  use of TAMC or TAF as AD tool.  differentiation using the software tool TAF.
696  The following sections describe the steps which are necessary to  Modifications to use OpenAD are described in \ref{sec_ad_openad}.
697  generate a tangent linear or adjoint model of the MITGCM.  
698  We take as an example the sensitivity of carbon sequestration  \input{s_autodiff/text/doc_ad_the_model}
699  in the ocean.  
700  The AD-relevant hooks in the code are sketched in  The basic flow is depicted in \ref{fig:adthemodel}.
701  \reffig{adthemodel}, \reffig{adthemain}.  If CPP option \texttt{ALLOW\_AUTODIFF\_TAMC} is defined,
702    the driver routine
703  \subsection{Overview of the experiment}  {\it the\_model\_main}, instead of calling {\it the\_main\_loop},
704    invokes the adjoint of this routine, {\it adthe\_main\_loop}
705  We describe an adjoint sensitivity analysis of outgassing from  (case \texttt{\#define ALLOW\_ADJOINT\_RUN}), or
706  the ocean into the atmosphere of a carbon-like tracer injected  the tangent linear of this routine {\it g\_the\_main\_loop}
707  into the ocean interior (see \cite{hil-eta:01}).  (case \texttt{\#define ALLOW\_TANGENTLINEAR\_RUN}),
708    which are the toplevel routines in terms of automatic differentiation.
709  \subsubsection{Passive tracer equation}  The routines {\it adthe\_main\_loop} or {\it g\_the\_main\_loop}
710    are generated by TAF.
711  For this work the MITGCM was augmented with a thermodynamically  It contains both the forward integration of the full model, the
712  inactive tracer, $C$. Tracer residing in the ocean  cost function calculation,
713  model surface layer is outgassed according to a relaxation time scale,  any additional storing that is required for efficient checkpointing,
714  $\mu$. Within the ocean interior, the tracer is passively advected  and the reverse integration of the adjoint model.
715  by the ocean model currents. The full equation for the time evolution  
716  %  [DESCRIBE IN A SEPARATE SECTION THE WORKING OF THE TLM]
717  \begin{equation}  
718  \label{carbon_ddt}  In Fig. \ref{fig:adthemodel}
719  \frac{\partial C}{\partial t} \, = \,  the structure of {\it adthe\_main\_loop} has been strongly
720  -U\cdot \nabla C \, - \, \mu C \, + \, \Gamma(C) \,+ \, S  simplified to focus on the essentials; in particular, no checkpointing
721  \end{equation}  procedures are shown here.
722  %  Prior to the call of {\it adthe\_main\_loop}, the routine
723  also includes a source term $S$. This term  {\it ctrl\_unpack} is invoked to unpack the control vector
724  represents interior sources of $C$ such as would arise due to  or initialise the control variables.
725  direct injection.  Following the call of {\it adthe\_main\_loop},
726  The velocity term, $U$, is the sum of the  the routine {\it ctrl\_pack}
727  model Eulerian circulation and an eddy-induced velocity, the latter  is invoked to pack the control vector
728  parameterized according to Gent/McWilliams  (cf. Section \ref{section_ctrl}).
729  (\cite{gen-mcw:90, gen-eta:95}).  If gradient checks are to be performed, the option
730  The convection function, $\Gamma$, mixes $C$ vertically wherever the  {\tt ALLOW\_GRADIENT\_CHECK} is defined. In this case
731  fluid is locally statically unstable.  the driver routine {\it grdchk\_main} is called after
732    the gradient has been computed via the adjoint
733  The outgassing time scale, $\mu$, in eqn. (\ref{carbon_ddt})  (cf. Section \ref{section_grdchk}).
734  is set so that \( 1/\mu \sim 1 \ \mathrm{year} \) for the surface  
735  ocean and $\mu=0$ elsewhere. With this value, eqn. (\ref{carbon_ddt})  %------------------------------------------------------------------
736  is valid as a prognostic equation for small perturbations in oceanic  
737  carbon concentrations. This configuration provides a  \subsection{General setup
738  powerful tool for examining the impact of large-scale ocean circulation  \label{section_ad_setup}}
739  on $ CO_2 $ outgassing due to interior injections.  
740  As source we choose a constant in time injection of  In order to configure AD-related setups the following packages need
741  $ S = 1 \,\, {\rm mol / s}$.  to be enabled:
742    {\it
743  \subsubsection{Model configuration}  \begin{table}[h!]
744    \begin{tabular}{l}
745  The model configuration employed has a constant  autodiff \\
746  $4^\circ \times 4^\circ$ resolution horizontal grid and realistic  ctrl \\
747  geography and bathymetry. Twenty vertical layers are used with  cost \\
748  vertical spacing ranging  grdchk \\
749  from 50 m near the surface to 815 m at depth.  \end{tabular}
750  Driven to steady-state by climatalogical wind-stress, heat and  \end{table}
751  fresh-water forcing the model reproduces well known large-scale  }
752  features of the ocean general circulation.  The packages are enabled by adding them to your experiment-specific
753    configuration file
754  \subsubsection{Outgassing cost function}  {\it packages.conf} (see Section ???).
755    
756  To quantify and understand outgassing due to injections of $C$  The following AD-specific CPP option files need to be customized:
 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/}:  
757  %  %
758  \begin{itemize}  \begin{itemize}
759  %  %
760  \item {\it .genmakerc}  \item {\it ECCO\_CPPOPTIONS.h} \\
761  %  This header file collects CPP options for the packages
762  \item {\it COST\_CPPOPTIONS.h}  {\it autodiff, cost, ctrl} as well as AD-unrelated options for
763  %  the external forcing package {\it exf}.
764  \item {\it CPP\_EEOPTIONS.h}  \footnote{NOTE: These options are not set in their package-specific
765  %  headers such as {\it COST\_CPPOPTIONS.h}, but are instead collected
766  \item {\it CPP\_OPTIONS.h}  in the single header file {\it ECCO\_CPPOPTIONS.h}.
767  %  The package-specific header files serve as simple
768  \item {\it CTRL\_OPTIONS.h}  placeholders at this point.}
769  %  %
770  \item {\it ECCO\_OPTIONS.h}  \item {\it tamc.h} \\
771  %  This header configures the splitting of the time stepping loop
772  \item {\it SIZE.h}  w.r.t. the 3-level checkpointing (see section ???).
773  %  
 \item {\it adcommon.h}  
 %  
 \item {\it tamc.h}  
774  %  %
775  \end{itemize}  \end{itemize}
776    
777    %------------------------------------------------------------------
778    
779    \subsection{Building the AD code using TAF
780    \label{section_ad_build}}
781    
782    The build process of an AD code is very similar to building
783    the forward model. However, depending on which AD code one wishes
784    to generate, and on which AD tool is available (TAF or TAMC),
785    the following {\tt make} targets are available:
786    
787    \begin{table}[h!]
788    {\footnotesize
789    \begin{tabular}{|ccll|}
790    \hline
791    ~ & {\it AD-target} & {\it output} & {\it description} \\
792    \hline
793    \hline
794    (1) & {\tt <MODE><TOOL>only} & {\tt <MODE>\_<TOOL>\_output.f}  &
795    generates code for $<$MODE$>$ using $<$TOOL$>$ \\
796    ~ & ~ & ~ & no {\tt make} dependencies on {\tt .F .h} \\
797    ~ & ~ & ~ & useful for compiling on remote platforms \\
798    \hline
799    (2) & {\tt <MODE><TOOL>} & {\tt <MODE>\_<TOOL>\_output.f}  &
800    generates code for $<$MODE$>$ using $<$TOOL$>$ \\
801    ~ & ~ & ~ & includes {\tt make} dependencies on {\tt .F .h} \\
802    ~ & ~ & ~ & i.e. input for $<$TOOL$>$ may be re-generated \\
803    \hline
804    (3) & {\tt <MODE>all} & {\tt mitgcmuv\_<MODE>}  &
805    generates code for $<$MODE$>$ using $<$TOOL$>$ \\
806    ~ & ~ & ~ & and compiles all code \\
807    ~ & ~ & ~ & (use of TAF is set as default) \\
808    \hline
809    \end{tabular}
810    }
811    \end{table}
812  %  %
813  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:  
814  %  %
815  \begin{itemize}  \begin{itemize}
816  %  %
817  \item {\it data}  \item $<$TOOL$>$
 %  
 \item {\it data.cost}  
 %  
 \item {\it data.ctrl}  
 %  
 \item {\it data.gmredi}  
818  %  %
819  \item {\it data.grdchk}  \begin{itemize}
 %  
 \item {\it data.optim}  
 %  
 \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}  
820  %  %
821  \item {\it pickup*}  \item {\tt TAF}
822    \item {\tt TAMC}
823  %  %
824  \end{itemize}  \end{itemize}
825  %  %
826  Finally, the file to generate the adjoint code resides in  \item $<$MODE$>$
 $ adjoint/ $:  
827  %  %
828  \begin{itemize}  \begin{itemize}
829  %  %
830  \item {\it makefile}  \item {\tt ad} generates the adjoint model (ADM)
831    \item {\tt ftl} generates the tangent linear model (TLM)
832    \item {\tt svd} generates both ADM and TLM for \\
833    singular value decomposition (SVD) type calculations
834  %  %
835  \end{itemize}  \end{itemize}
836  %  %
837    \end{itemize}
838    
839  Below we describe the customisations of this files which are  For example, to generate the adjoint model using TAF after routines ({\tt .F})
840  specific to this experiment.  or headers ({\tt .h}) have been modified, but without compilation,
841    type {\tt make adtaf};
842  \subsubsection{File {\it .genmakerc}}  or, to generate the tangent linear model using TAMC without
843  This file overwrites 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 gmredi grdchk kpp )}  \\  
 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 a 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 three options specific to this experiment are \\  
 \hspace*{4ex} {\tt \#define ALLOW\_PASSIVE\_TRACER} \\  
 This flag enables the code to carry through the  
 advection/diffusion of a passive tracer along the  
 model integration. \\  
 \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}. \\  
 \hspace*{4ex} {\tt \#define ALLOW\_GRADIENT\_CHECK} \\  
 This flag enables the gradient check package.  
 After computing the unperturbed cost function and its gradient,  
 a series of computations are performed for which \\  
 $\bullet$ an element of the control vector is perturbed \\  
 $\bullet$ the cost function w.r.t. the perturbed element is  
 computed \\  
 $\bullet$ the difference between the perturbed and unperturbed  
 cost function is computed to compute the finite difference gradient \\  
 $\bullet$ the finite difference gradient is compared with the  
 adjoint-generated gradient.  
 The gradient check package is further described in Section ???.  
   
 \subsubsection{File {\it ECCO\_OPTIONS.h}}  
844    
 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 call to 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 \#define ALLOW\_DIFFKR\_CONTROL} &  
 diapycnal diffusivity \\  
 \hspace*{2ex} {\tt \#undef ALLOW\_KAPPAGM\_CONTROL} &  
 isopycnal diffusivity \\  
 \end{tabular}  
 %  
 \end{itemize}  
845    
846  \subsubsection{File {\it SIZE.h}}  A typical full build process to generate the ADM via TAF would
847    look like follows:
848    \begin{verbatim}
849    % mkdir build
850    % cd build
851    % ../../../tools/genmake2 -mods=../code_ad
852    % make depend
853    % make adall
854    \end{verbatim}
855    
856  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 /addynvars\_diffkr/} &  
 \hspace*{4ex} is related to {\it DYNVARS.h} \\  
 \hspace*{4ex} {\tt common /addynvars\_kapgm/} &  
 \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}  
857    
858  Note that if the structure of the common block changes in the  \subsection{The AD build process in detail
859  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}.  
860    
861  \subsubsection{File {\it tamc.h}}  The {\tt make <MODE>all} target consists of the following procedures:
862    
863  This routine contains the dimensions for TAMC checkpointing.  \begin{enumerate}
864  %  %
865    \item
866    A header file {\tt AD\_CONFIG.h} is generated which contains a CPP option
867    on which code ought to be generated. Depending on the {\tt make} target,
868    the contents is one of the following:
869  \begin{itemize}  \begin{itemize}
870  %  \item
871  \item {\tt \#ifdef ALLOW\_TAMC\_CHECKPOINTING} \\  {\tt \#define ALLOW\_ADJOINT\_RUN}
872  3-level checkpointing is enabled, i.e. the timestepping  \item
873  is divided into three different levels (see Section \ref{???}).  {\tt \#define ALLOW\_TANGENTLINEAR\_RUN}
874  The model state of the outermost ({\tt nchklev\_3}) and the  \item
875  intermediate ({\tt nchklev\_2}) timestepping loop are stored to file  {\tt \#define ALLOW\_ECCO\_OPTIMIZATION}
 (handled in {\it the\_main\_loop}).  
 The innermost loop ({\tt nchklev\_1})  
 avoids I/O by storing all required variables  
 to common blocks. This storing may also be necessary if  
 no checkpointing is chosen  
 (nonlinear functions, if-statements, iterative loops, ...).  
 In the present example the dimensions are chosen as follows: \\  
 \hspace*{4ex} {\tt nchklev\_1      =  36 } \\  
 \hspace*{4ex} {\tt nchklev\_2      =  30 } \\  
 \hspace*{4ex} {\tt nchklev\_3      =  60 } \\  
 To guarantee that the checkpointing intervals span the entire  
 integration period the following relation must be satisfied: \\  
 \hspace*{4ex} {\tt nchklev\_1*nchklev\_2*nchklev\_3 $ \ge $ nTimeSteps} \\  
 where {\tt nTimeSteps} is either specified in {\it data}  
 or computed via \\  
 \hspace*{4ex} {\tt nTimeSteps = (endTime-startTime)/deltaTClock }.  
 %  
 \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}.  
 %  
876  \end{itemize}  \end{itemize}
   
 The following parameters may be worth describing: \\  
877  %  %
878  \hspace*{4ex} {\tt isbyte} \\  \item
879  \hspace*{4ex} {\tt maxpass} \\  A single file {\tt <MODE>\_input\_code.f} is concatenated
880  ~  consisting of all {\tt .f} files that are part of the list {\bf AD\_FILES}
881    and all {\tt .flow} files that are part of the list {\bf AD\_FLOW\_FILES}.
882  \subsubsection{File {\it makefile}}  %
883    \item
884  This file contains all relevant paramter flags and  The AD tool is invoked with the {\tt <MODE>\_<TOOL>\_FLAGS}.
885  lists to run TAMC or TAF.  The default AD tool flags in {\tt genmake2} can be overrwritten by
886  It is assumed that TAMC is available to you, either locally,  an {\tt adjoint\_options} file (similar to the platform-specific
887  being installed on your network, or remotely through the 'TAMC Utility'.  {\tt build\_options}, see Section ???.
888  TAMC is called with the command {\tt tamc} followed by a  The AD tool writes the resulting AD code into the file
889  number of options. They are described in detail in the  {\tt <MODE>\_input\_code\_ad.f}
890  TAMC manual \cite{gie:99}.  %
891  Here we briefly discuss the main flags used in the {\it makefile}  \item
892    A short sed script {\tt adjoint\_sed} is applied to
893    {\tt <MODE>\_input\_code\_ad.f}
894    to reinstate {\bf myThid} into the CALL argument list of active file I/O.
895    The result is written to file {\tt <MODE>\_<TOOL>\_output.f}.
896    %
897    \item
898    All routines are compiled and an executable is generated
899    (see Table ???).
900  %  %
901  \begin{itemize}  \end{enumerate}
902  \item [{\tt tamc}] {\tt  
903  -input <variable names>  \subsubsection{The list AD\_FILES and {\tt .list} files}
904  -output <variable name> -r4 ... \\  
905  -toplevel <S/R name> -reverse <file names>  Not all routines are presented to the AD tool.
906  }  Routines typically hidden are diagnostics routines which
907  \end{itemize}  do not influence the cost function, but may create
908    artificial flow dependencies such as I/O of active variables.
909    
910    {\tt genmake2} generates a list (or variable) {\bf AD\_FILES}
911    which contains all routines that are shown to the AD tool.
912    This list is put together from all files with suffix {\tt .list}
913    that {\tt genmake2} finds in its search directories.
914    The list file for the core MITgcm routines is in {\tt model/src/}
915    is called {\tt model\_ad\_diff.list}.
916    Note that no wrapper routine is shown to TAF. These are either
917    not visible at all to the AD code, or hand-written AD code
918    is available (see next section).
919    
920    Each package directory contains its package-specific
921    list file {\tt <PKG>\_ad\_diff.list}. For example,
922    {\tt pkg/ptracers/} contains the file {\tt ptracers\_ad\_diff.list}.
923    Thus, enabling a package will automatically extend the
924    {\bf AD\_FILES} list of {\tt genmake2} to incorporate the
925    package-specific routines.
926    Note that you will need to regenerate the {\tt Makefile} if
927    you enable a package (e.g. by adding it to {\tt packages.conf})
928    and a {\tt Makefile} already exists.
929    
930    \subsubsection{The list AD\_FLOW\_FILES and {\tt .flow} files}
931    
932    TAMC and TAF can evaluate user-specified directives
933    that start with a specific syntax ({\tt CADJ}, {\tt C\$TAF}, {\tt !\$TAF}).
934    The main categories of directives are STORE directives and
935    FLOW directives. Here, we are concerned with flow directives,
936    store directives are treated elsewhere.
937    
938    Flow directives enable the AD tool to evaluate how it should treat
939    routines that are 'hidden' by the user, i.e. routines which are
940    not contained in the {\bf AD\_FILES} list (see previous section),
941    but which are called in part of the code that the AD tool does see.
942    The flow directive tell the AD tool
943  %  %
944  \begin{itemize}  \begin{itemize}
945  %  %
946  \item {\tt -toplevel <S/R name>} \\  \item which subroutine arguments are input/output
947  Name of the toplevel routine, with respect to which the  \item which subroutine arguments are active
948  control flow analysis is performed.  \item which subroutine arguments are required to compute the cost
949  %  \item which subroutine arguments are dependent
 \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.  
 In the discussion of Section ???  
 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  
 (or for each package which contains such routines),  
 a corresponding file {\it .flow} exists containing flow directives  
 for TAMC.  
 %  
 \item {\tt -r4} \\  
 ~  
950  %  %
951  \end{itemize}  \end{itemize}
952    %
953    The syntax for the flow directives can be found in the
954    AD tool manuals.
955    
956    {\tt genmake2} generates a list (or variable) {\bf AD\_FLOW\_FILES}
957    which contains all files with suffix{\tt .flow} that it finds
958    in its search directories.
959    The flow directives for the core MITgcm routines of
960    {\tt eesupp/src/} and {\tt model/src/}
961    reside in {\tt pkg/autodiff/}.
962    This directory also contains hand-written adjoint code
963    for the MITgcm WRAPPER (section \ref{chap:sarch}).
964    
965    Flow directives for package-specific routines are contained in
966    the corresponding package directories in the file
967    {\tt <PKG>\_ad.flow}, e.g. ptracers-specific directives are in
968    {\tt ptracers\_ad.flow}.
969    
970    \subsubsection{Store directives for 3-level checkpointing}
971    
972    The storing that is required at each period of the
973    3-level checkpointing is controled by three
974    top-level headers.
975    
976  \subsubsection{File {\it data}}  \begin{verbatim}
977    do ilev_3 = 1, nchklev_3
978  \subsubsection{File {\it data.cost}}  #  include ``checkpoint_lev3.h''
979       do ilev_2 = 1, nchklev_2
980  \subsubsection{File {\it data.ctrl}}  #     include ``checkpoint_lev2.h''
981          do ilev_1 = 1, nchklev_1
982  \subsubsection{File {\it data.gmredi}}  #        include ``checkpoint_lev1.h''
983    
984  \subsubsection{File {\it data.grdchk}}  ...
985    
986  \subsubsection{File {\it data.optim}}        end do
987       end do
988  \subsubsection{File {\it data.pkg}}  end do
989    \end{verbatim}
 \subsubsection{File {\it eedata}}  
   
 \subsubsection{File {\it topog.bin}}  
990    
991  \subsubsection{File {\it windx.bin, windy.bin}}  All files {\tt checkpoint\_lev?.h} are contained in directory
992    {\tt pkg/autodiff/}.
993    
 \subsubsection{File {\it salt.bin, theta.bin}}  
994    
995  \subsubsection{File {\it SSS.bin, SST.bin}}  \subsubsection{Changing the default AD tool flags: ad\_options files}
996    
 \subsubsection{File {\it pickup*}}  
997    
998  \subsection{Compiling the model and its adjoint}  \subsubsection{Hand-written adjoint code}
999    
1000  \newpage  %------------------------------------------------------------------
   
 %**********************************************************************  
 \section{TLM and ADM generation in general}  
 \label{sec_ad_setup_gen}  
 %**********************************************************************  
   
 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.  
   
 \begin{figure}[b!]  
 \input{part5/doc_ad_the_model}  
 \caption{~}  
 \label{fig:adthemodel}  
 \end{figure}  
   
 The basic flow is depicted in \reffig{adthemodel}.  
 If the option {\tt ALLOW\_AUTODIFF\_TAMC} is defined, the driver routine  
 {\it the\_model\_main}, instead of calling {\it the\_main\_loop},  
 invokes the adjoint of this routine, {\it adthe\_main\_loop},  
 which is the toplevel routine in terms of reverse mode computation.  
 The routine {\it adthe\_main\_loop} has been generated using TAMC.  
 It contains both the forward integration of the full model,  
 any additional storing that is required for efficient checkpointing,  
 and the reverse integration of the adjoint model.  
 The structure of {\it adthe\_main\_loop} has been strongly  
 simplified for clarification; in particular, no checkpointing  
 procedures are shown here.  
 Prior to the call of {\it adthe\_main\_loop}, the routine  
 {\it ctrl\_unpack} is invoked to unpack the control vector,  
 and following that call, the routine {\it ctrl\_pack}  
 is invoked to pack the control vector  
 (cf. Section \ref{section_ctrl}).  
 If gradient checks are to be performed, the option  
 {\tt ALLOW\_GRADIENT\_CHECK} is defined. In this case  
 the driver routine {\it grdchk\_main} is called after  
 the gradient has been computed via the adjoint  
 (cf. Section \ref{section_grdchk}).  
1001    
1002  \subsection{The cost function (dependent variable)  \subsection{The cost function (dependent variable)
1003  \label{section_cost}}  \label{section_cost}}
# Line 1165  the gradient has been computed via the a Line 1005  the gradient has been computed via the a
1005  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}.
1006  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
1007  $ {\cal J}(\vec{u}) \, = \, {\cal J}(M(\vec{u})) $.  $ {\cal J}(\vec{u}) \, = \, {\cal J}(M(\vec{u})) $.
1008  The input is referred to as the  The input are referred to as the
1009  {\sf independent variables} or {\sf control variables}.  {\sf independent variables} or {\sf control variables}.
1010  All aspects relevant to the treatment of the cost function $ {\cal J} $  All aspects relevant to the treatment of the cost function $ {\cal J} $
1011  (parameter setting, initialisation, accumulation,  (parameter setting, initialization, accumulation,
1012  final evaluation), are controlled by the package {\it pkg/cost}.  final evaluation), are controlled by the package {\it pkg/cost}.
1013    The aspects relevant to the treatment of the independent variables
1014    are controlled by the package {\it pkg/ctrl} and will be treated
1015    in the next section.
1016    
1017  \begin{figure}[h!]  \input{s_autodiff/text/doc_cost_flow}
1018  \input{part5/doc_cost_flow}  
1019  \caption{~}  \subsubsection{Enabling the package}
 \label{fig:costflow}  
 \end{figure}  
1020    
 \subsubsection{genmake and CPP options}  
 %  
 \begin{itemize}  
 %  
 \item  
1021  \fbox{  \fbox{
1022  \begin{minipage}{12cm}  \begin{minipage}{12cm}
1023  {\it genmake}, {\it CPP\_OPTIONS.h}, {\it ECCO\_CPPOPTIONS.h}  {\it packages.conf}, {\it ECCO\_CPPOPTIONS.h}
1024  \end{minipage}  \end{minipage}
1025  }  }
1026  \end{itemize}  \begin{itemize}
 %  
 The directory {\it pkg/cost} can be included to the  
 compile list in 3 different ways (cf. Section \ref{???}):  
1027  %  %
1028  \begin{enumerate}  \item
1029    The package is enabled by adding {\it cost} to your file {\it packages.conf}
1030    (see Section ???)
1031  %  %
1032  \item {\it genmake}: \\  \item
1033  Change the default settings in the file {\it genmake} by adding  
1034  {\bf cost} to the {\bf enable} list (not recommended).  
1035  %  \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}.  
1036  %  %
1037  \end{enumerate}  
1038    N.B.: In general the following packages ought to be enabled
1039    simultaneously: {\it autodiff, cost, ctrl}.
1040  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}.
1041  Each specific cost function contribution has its own option.  Each specific cost function contribution has its own option.
1042  For the present example the option is {\bf ALLOW\_COST\_TRACER}.  For the present example the option is {\bf ALLOW\_COST\_TRACER}.
1043  All cost-specific options are set in {\it ECCO\_CPPOPTIONS.h}  All cost-specific options are set in {\it ECCO\_CPPOPTIONS.h}
1044  Since the cost function is usually used in conjunction with  Since the cost function is usually used in conjunction with
1045  automatic differentiation, the CPP option  automatic differentiation, the CPP option
1046  {\bf ALLOW\_ADJOINT\_RUN} should be defined  {\bf ALLOW\_ADJOINT\_RUN} (file {\it CPP\_OPTIONS.h}) and
1047  (file {\it CPP\_OPTIONS.h}).  {\bf ALLOW\_AUTODIFF\_TAMC} (file {\it ECCO\_CPPOPTIONS.h})
1048    should be defined.
1049    
1050  \subsubsection{Initialisation}  \subsubsection{Initialization}
1051  %  %
1052  The initialisation of the {\it cost} package is readily enabled  The initialization of the {\it cost} package is readily enabled
1053  as soon as the CPP option {\bf ALLOW\_ADJOINT\_RUN} is defined.  as soon as the CPP option {\bf ALLOW\_COST} is defined.
1054  %  %
1055  \begin{itemize}  \begin{itemize}
1056  %  %
# Line 1250  Variables: {\it cost\_init} Line 1080  Variables: {\it cost\_init}
1080  }  }
1081  \\  \\
1082  This S/R  This S/R
1083  initialises the different cost function contributions.  initializes the different cost function contributions.
1084  The contribtion for the present example is {\bf objf\_tracer}  The contribution for the present example is {\bf objf\_tracer}
1085  which is defined on each tile (bi,bj).  which is defined on each tile (bi,bj).
1086  %  %
1087  \end{itemize}  \end{itemize}
# Line 1294  from each contribution and sums over all Line 1124  from each contribution and sums over all
1124  \begin{equation}  \begin{equation}
1125  {\cal J} \, = \,  {\cal J} \, = \,
1126  {\rm fc} \, = \,  {\rm fc} \, = \,
1127  {\rm mult\_tracer} \sum_{bi,\,bj}^{nSx,\,nSy}  {\rm mult\_tracer} \sum_{\text{global sum}} \sum_{bi,\,bj}^{nSx,\,nSy}
1128  {\rm objf\_tracer}(bi,bj) \, + \, ...  {\rm objf\_tracer}(bi,bj) \, + \, ...
1129  \end{equation}  \end{equation}
1130  %  %
1131  The total cost function {\bf fc} will be the  The total cost function {\bf fc} will be the
1132  'dependent' variable in the argument list for TAMC, i.e.  'dependent' variable in the argument list for TAF, i.e.
1133  \begin{verbatim}  \begin{verbatim}
1134  tamc -output 'fc' ...  taf -output 'fc' ...
1135  \end{verbatim}  \end{verbatim}
1136    
1137  %%%% \end{document}  %%%% \end{document}
1138    
1139  \begin{figure}  \input{s_autodiff/text/doc_ad_the_main}
 \input{part5/doc_ad_the_main}  
 \caption{~}  
 \label{fig:adthemain}  
 \end{figure}  
1140    
1141  \subsection{The control variables (independent variables)  \subsection{The control variables (independent variables)
1142  \label{section_ctrl}}  \label{section_ctrl}}
# Line 1327  as variable assignments. Therefore, file Line 1153  as variable assignments. Therefore, file
1153  active variables are written and from which active variables  active variables are written and from which active variables
1154  are read are called {\sf active files}.  are read are called {\sf active files}.
1155  All aspects relevant to the treatment of the control variables  All aspects relevant to the treatment of the control variables
1156  (parameter setting, initialisation, perturbation)  (parameter setting, initialization, perturbation)
1157  are controled by the package {\it pkg/ctrl}.  are controlled by the package {\it pkg/ctrl}.
1158    
1159  \begin{figure}[h!]  \input{s_autodiff/text/doc_ctrl_flow}
 \input{part5/doc_ctrl_flow}  
 \caption{~}  
 \label{fig:ctrlflow}  
 \end{figure}  
1160    
1161  \subsubsection{genmake and CPP options}  \subsubsection{genmake and CPP options}
1162  %  %
# Line 1350  are controled by the package {\it pkg/ct Line 1172  are controled by the package {\it pkg/ct
1172  %  %
1173  To enable the directory to be included to the compile list,  To enable the directory to be included to the compile list,
1174  {\bf ctrl} has to be added to the {\bf enable} list in  {\bf ctrl} has to be added to the {\bf enable} list in
1175  {\it .genmakerc} (or {\it genmake} itself).  {\it .genmakerc} or in {\it genmake} itself (analogous to {\it cost}
1176    package, cf. previous section).
1177  Each control variable is enabled via its own CPP option  Each control variable is enabled via its own CPP option
1178  in {\it ECCO\_CPPOPTIONS.h}.  in {\it ECCO\_CPPOPTIONS.h}.
1179    
1180  \subsubsection{Initialisation}  \subsubsection{Initialization}
1181  %  %
1182  \begin{itemize}  \begin{itemize}
1183  %  %
# Line 1391  and their gradients: {\it ctrl\_unpack} Line 1214  and their gradients: {\it ctrl\_unpack}
1214  \\  \\
1215  %  %
1216  Two important issues related to the handling of the control  Two important issues related to the handling of the control
1217  variables in the MITGCM need to be addressed.  variables in MITgcm need to be addressed.
1218  First, in order to save memory, the control variable arrays  First, in order to save memory, the control variable arrays
1219  are not kept in memory, but rather read from file and added  are not kept in memory, but rather read from file and added
1220  to the initial fields during the model initialisation phase.  to the initial fields during the model initialization phase.
1221  Similarly, the corresponding adjoint fields which represent  Similarly, the corresponding adjoint fields which represent
1222  the gradient of the cost function w.r.t. the control variables  the gradient of the cost function w.r.t. the control variables
1223  are written to file at the end of the adjoint integration.  are written to file at the end of the adjoint integration.
# Line 1423  and gradient are generated and initialis Line 1246  and gradient are generated and initialis
1246  %  %
1247  The dependency flow for differentiation w.r.t. the controls  The dependency flow for differentiation w.r.t. the controls
1248  starts with adding a perturbation onto the input variable,  starts with adding a perturbation onto the input variable,
1249  thus defining the independent or control variables for TAMC.  thus defining the independent or control variables for TAF.
1250  Three types of controls may be considered:  Three types of controls may be considered:
1251  %  %
1252  \begin{itemize}  \begin{itemize}
# Line 1456  u         & = \, u_{[0]} \, + \, \Delta Line 1279  u         & = \, u_{[0]} \, + \, \Delta
1279  holding the perturbation. In the case of a simple  holding the perturbation. In the case of a simple
1280  sensitivity study this array is identical to zero.  sensitivity study this array is identical to zero.
1281  However, it's specification is essential in the context  However, it's specification is essential in the context
1282  of automatic differentiation since TAMC  of automatic differentiation since TAF
1283  treats the corresponding line in the code symbolically  treats the corresponding line in the code symbolically
1284  when determining the differentiation chain and its origin.  when determining the differentiation chain and its origin.
1285  Thus, the variable names are part of the argument list  Thus, the variable names are part of the argument list
1286  when calling TAMC:  when calling TAF:
1287  %  %
1288  \begin{verbatim}  \begin{verbatim}
1289  tamc -input 'xx_tr1 ...' ...  taf -input 'xx_tr1 ...' ...
1290  \end{verbatim}  \end{verbatim}
1291  %  %
1292  Now, as mentioned above, the MITGCM avoids maintaining  Now, as mentioned above, MITgcm avoids maintaining
1293  an array for each control variable by reading the  an array for each control variable by reading the
1294  perturbation to a temporary array from file.  perturbation to a temporary array from file.
1295  To ensure the symbolic link to be recognized by TAMC, a scalar  To ensure the symbolic link to be recognized by TAF, a scalar
1296  dummy variable {\bf xx\_tr1\_dummy} is introduced  dummy variable {\bf xx\_tr1\_dummy} is introduced
1297  and an 'active read' routine of the adjoint support  and an 'active read' routine of the adjoint support
1298  package {\it pkg/autodiff} is invoked.  package {\it pkg/autodiff} is invoked.
1299  The read-procedure is tagged with the variable  The read-procedure is tagged with the variable
1300  {\bf xx\_tr1\_dummy} enabbling TAMC to recognize the  {\bf xx\_tr1\_dummy} enabling TAF to recognize the
1301  initialisation of the perturbation.  initialization of the perturbation.
1302  The modified call of TAMC thus reads  The modified call of TAF thus reads
1303  %  %
1304  \begin{verbatim}  \begin{verbatim}
1305  tamc -input 'xx_tr1_dummy ...' ...  taf -input 'xx_tr1_dummy ...' ...
1306  \end{verbatim}  \end{verbatim}
1307  %  %
1308  and the modified operation to (\ref{perturb})  and the modified operation to (\ref{perturb})
# Line 1494  in the code takes on the form Line 1317  in the code takes on the form
1317  %  %
1318  Note, that reading an active variable corresponds  Note, that reading an active variable corresponds
1319  to a variable assignment. Its derivative corresponds  to a variable assignment. Its derivative corresponds
1320  to a write statement of the adjoint variable.  to a write statement of the adjoint variable, followed by
1321    a reset.
1322  The 'active file' routines have been designed  The 'active file' routines have been designed
1323  to support active read and corresponding adjoint active write  to support active read and corresponding adjoint active write
1324  operations (and vice versa).  operations (and vice versa).
# Line 1586  variables are written to {\bf adxx\_ ... Line 1410  variables are written to {\bf adxx\_ ...
1410  \begin{itemize}  \begin{itemize}
1411  %  %
1412  \item {\bf vector\_ctrl}: the control vector \\  \item {\bf vector\_ctrl}: the control vector \\
1413  At the very beginning of the model initialisation,  At the very beginning of the model initialization,
1414  the updated compressed control vector is read (or initialised)  the updated compressed control vector is read (or initialised)
1415  and distributed to 2-dim. and 3-dim. control variable fields.  and distributed to 2-dim. and 3-dim. control variable fields.
1416  %  %
# Line 1611  at intermediate times can be written usi Line 1435  at intermediate times can be written usi
1435  {\it addummy\_in\_stepping}.  {\it addummy\_in\_stepping}.
1436  This routine is part of the adjoint support package  This routine is part of the adjoint support package
1437  {\it pkg/autodiff} (cf.f. below).  {\it pkg/autodiff} (cf.f. below).
1438    The procedure is enabled using via the CPP-option
1439    {\bf ALLOW\_AUTODIFF\_MONITOR} (file {\it ECCO\_CPPOPTIONS.h}).
1440  To be part of the adjoint code, the corresponding S/R  To be part of the adjoint code, the corresponding S/R
1441  {\it dummy\_in\_stepping} has to be called in the forward  {\it dummy\_in\_stepping} has to be called in the forward
1442  model (S/R {\it the\_main\_loop}) at the appropriate place.  model (S/R {\it the\_main\_loop}) at the appropriate place.
1443    The adjoint common blocks are extracted from the adjoint code
1444    via the header file {\it adcommon.h}.
1445    
1446  {\it dummy\_in\_stepping} is essentially empty,  {\it dummy\_in\_stepping} is essentially empty,
1447  the corresponding adjoint routine is hand-written rather  the corresponding adjoint routine is hand-written rather
# Line 1640  the common blocks Line 1468  the common blocks
1468  {\bf /adtr1\_r/}, {\bf /adffields/},  {\bf /adtr1\_r/}, {\bf /adffields/},
1469  which have been extracted from the adjoint code to enable  which have been extracted from the adjoint code to enable
1470  access to the adjoint variables.  access to the adjoint variables.
1471    
1472    {\bf WARNING:} If the structure of the common blocks
1473    {\bf /dynvars\_r/}, {\bf /dynvars\_cd/}, etc., changes
1474    similar changes will occur in the adjoint common blocks.
1475    Therefore, consistency between the TAMC-generated common blocks
1476    and those in {\it adcommon.h} have to be checked.
1477  %  %
1478  \end{itemize}  \end{itemize}
1479    
# Line 1654  The gradient $ \nabla _{u}{\cal J} |_{u_ Line 1488  The gradient $ \nabla _{u}{\cal J} |_{u_
1488  with the value of the cost function itself $ {\cal J}(u_{[k]}) $  with the value of the cost function itself $ {\cal J}(u_{[k]}) $
1489  at iteration step $ k $ serve  at iteration step $ k $ serve
1490  as input to a minimization routine (e.g. quasi-Newton method,  as input to a minimization routine (e.g. quasi-Newton method,
1491  conjugate gradient, ... \cite{gil_lem:89})  conjugate gradient, ... \cite{gil-lem:89})
1492  to compute an update in the  to compute an update in the
1493  control variable for iteration step $k+1$  control variable for iteration step $k+1$
1494  \[  \[
# Line 1785  to {\it adxx\_...$<$k$>$}, again via the Line 1619  to {\it adxx\_...$<$k$>$}, again via the
1619  Finally, {\it ctrl\_pack} collects all adjoint files  Finally, {\it ctrl\_pack} collects all adjoint files
1620  and writes them to the compressed vector file  and writes them to the compressed vector file
1621  {\bf vector\_grad\_$<$k$>$}.  {\bf vector\_grad\_$<$k$>$}.
   
 \subsection{TLM and ADM generation via TAMC}  
   
   
   
 \subsection{Flow directives and adjoint support routines \label{section_flowdir}}  
   
 \subsection{Store directives and checkpointing \label{section_checkpointing}}  
   
 \subsection{Gradient checks \label{section_grdchk}}  
   
 \subsection{Second derivative generation via TAMC}  
   
 \section{Example of adjoint code}  

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.22

  ViewVC Help
Powered by ViewVC 1.1.22