manual/s_autodiff/doc_ad_examples.tex

% $Header: /u/u0/gcmpack/manual/part5/doc_ad_2.tex,v 1.13 2002/01/18 22:56:45 heimbach Exp $
% $Name:  $

%**********************************************************************
\section{Sensitivity of Air-Sea Exchange to Tracer Injection Site }
\label{sect:eg-simple-tracer-adjoint}
\label{sec_ad_setup_ex}
%**********************************************************************

The MITGCM has been adapted to enable AD using TAMC or TAF.
The present description, therefore, is specific to the
use of TAMC or TAF as AD tool.
The following sections describe the steps which are necessary to 
generate a tangent linear or adjoint model of the MITGCM.
We take as an example the sensitivity of carbon sequestration
in the ocean.
The AD-relevant hooks in the code are sketched in 
\ref{fig:adthemodel}, \ref{fig:adthemain}.

\subsection{Overview of the experiment}

We describe an adjoint sensitivity analysis of out-gassing from 
the ocean into the atmosphere of a carbon-like tracer injected
into the ocean interior (see \cite{hil-eta:01}).

\subsubsection{Passive tracer equation}

For this work the MITGCM was augmented with a thermodynamically 
inactive tracer, $C$. Tracer residing in the ocean 
model surface layer is out-gassed according to a relaxation time scale, 
$\mu$. Within the ocean interior, the tracer is passively advected 
by the ocean model currents. The full equation for the time evolution
%
\begin{equation}
\label{carbon_ddt}
\frac{\partial C}{\partial t} \, = \, 
-U\cdot \nabla C \, - \, \mu C \, + \, \Gamma(C) \,+ \, S
\end{equation}
%
also includes a source term $S$. This term 
represents interior sources of $C$ such as would arise due to
direct injection.
The velocity term, $U$, is the sum of the
model Eulerian circulation and an eddy-induced velocity, the latter
parameterized according to Gent/McWilliams 
(\cite{gen-mcw:90, gen-eta:95}).
The convection function, $\Gamma$, mixes $C$ vertically wherever the 
fluid is locally statically unstable. 

The out-gassing time scale, $\mu$, in eqn. (\ref{carbon_ddt})
is set so that \( 1/\mu \sim 1 \ \mathrm{year} \) for the surface
ocean and $\mu=0$ elsewhere. With this value, eqn. (\ref{carbon_ddt})
is valid as a prognostic equation for small perturbations in oceanic 
carbon concentrations. This configuration provides a 
powerful tool for examining the impact of large-scale ocean circulation
on $ CO_2 $ out-gassing due to interior injections.
As source we choose a constant in time injection of 
$ S = 1 \,\, {\rm mol / s}$.

\subsubsection{Model configuration}

The model configuration employed has a constant 
$4^\circ \times 4^\circ$ resolution horizontal grid and realistic 
geography and bathymetry. Twenty vertical layers are used with 
vertical spacing ranging
from 50 m near the surface to 815 m at depth. 
Driven to steady-state by climatological wind-stress, heat and
fresh-water forcing the model reproduces well known large-scale
features of the ocean general circulation. 

\subsubsection{Out-gassing cost function}

To quantify and understand out-gassing due to injections of $C$
in eqn. (\ref{carbon_ddt}),
we define a cost function $ {\cal J} $ that measures the total amount of 
tracer out-gassed 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 out-gassing 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 out-gassed 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 out-gas following injection
and regions in which $CO_2$ injections would remain effectively 
sequestered within the ocean.

\subsection{Code configuration}

The model configuration for this experiment resides under the
directory {\it verification/carbon/}.
The code customization routines are in {\it verification/carbon/code/}:
%
\begin{itemize}
%
\item {\it .genmakerc}
%
\item {\it COST\_CPPOPTIONS.h}
%
\item {\it CPP\_EEOPTIONS.h}
%
\item {\it CPP\_OPTIONS.h}
%
\item {\it CTRL\_OPTIONS.h}
%
\item {\it ECCO\_OPTIONS.h}
%
\item {\it SIZE.h}
%
\item {\it adcommon.h}
%
\item {\it tamc.h}
%
\end{itemize}
%
The runtime flag and parameters settings are contained in 
{\it verification/carbon/input/},
together with the forcing fields and and restart files:
%
\begin{itemize}
%
\item {\it data}
%
\item {\it data.cost}
%
\item {\it data.ctrl}
%
\item {\it data.gmredi}
%
\item {\it data.grdchk}
%
\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}
%
\item {\it pickup*}
%
\end{itemize}
%
Finally, the file to generate the adjoint code resides in
$ adjoint/ $:
%
\begin{itemize}
%
\item {\it makefile}
%
\end{itemize}
%

Below we describe the customizations of this files which are
specific to this experiment.

\subsubsection{File {\it .genmakerc}}
This file overwrites default settings of {\it genmake}.
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 initialization, 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}}

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 initializations, 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
initializing, 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 forward code for
initializing, 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}

\subsubsection{File {\it SIZE.h}}

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 corresponds 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}

Note that if the structure of the common block changes in the 
above header files of the forward code, the structure
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}.

\subsubsection{File {\it tamc.h}}

This routine contains the dimensions for TAMC checkpointing.
%
\begin{itemize}
%
\item {\tt \#ifdef ALLOW\_TAMC\_CHECKPOINTING} \\
3-level checkpointing is enabled, i.e. the timestepping
is divided into three different levels (see Section \ref{???}).
The model state of the outermost ({\tt nchklev\_3}) and the
intermediate ({\tt nchklev\_2}) timestepping loop are stored to file
(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}.
%
\end{itemize}

The following parameters may be worth describing: \\
%
\hspace*{4ex} {\tt isbyte} \\
\hspace*{4ex} {\tt maxpass} \\
~
 
\subsubsection{File {\it makefile}}

This file contains all relevant parameter flags and
lists to run TAMC or TAF.
It is assumed that TAMC is available to you, either locally,
being installed on your network, or remotely through the 'TAMC Utility'.
TAMC is called with the command {\tt tamc} followed by a
number of options. They are described in detail in the
TAMC manual \cite{gie:99}.
Here we briefly discuss the main flags used in the {\it makefile}
%
\begin{itemize}
\item [{\tt tamc}] {\tt
-input <variable names>
-output <variable name> -r4 ... \\
-toplevel <S/R name> -reverse <file names> 
}
\end{itemize}
%
\begin{itemize}
%
\item {\tt -toplevel <S/R name>} \\
Name of the toplevel routine, with respect to which the
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.
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 initializations), 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} \\
~
%
\end{itemize}


\subsubsection{The input parameter files}

\paragraph{File {\it data}}

\paragraph{File {\it data.cost}}

\paragraph{File {\it data.ctrl}}

\paragraph{File {\it data.gmredi}}

\paragraph{File {\it data.grdchk}}

\paragraph{File {\it data.optim}}

\paragraph{File {\it data.pkg}}

\paragraph{File {\it eedata}}

\paragraph{File {\it topog.bin}}

\paragraph{File {\it windx.bin, windy.bin}}

\paragraph{File {\it salt.bin, theta.bin}}

\paragraph{File {\it SSS.bin, SST.bin}}

\paragraph{File {\it pickup*}}

\subsection{Compiling the model and its adjoint}

The built process of the adjoint model is slightly more
complex than that of compiling the forward code.
The main reason is that the adjoint code generation requires
a specific list of routines that are to be differentiated
(as opposed to the automatic generation of a list of
files to be compiled by genmake).
This list excludes routines that don't have to be or must not be
differentiated. For some of the latter routines flow directives
may be necessary, a list of which has to be given as well.
For this reason, a separate {\it makefile} is currently 
maintained in the directory {\tt adjoint/}. This
makefile is responsible for the adjoint code generation.

In the following we describe the build process step by step,
assuming you are in the directory {\tt bin/}.
A summary of steps to follow is given at the end.

\paragraph{Adjoint code generation and compilation -- step by step}

\begin{enumerate}
%
\item
{\tt ln -s ../verification/???/code/.genmakerc .} \\
{\tt ln -s ../verification/???/code/*.[Fh] .} \\
Link your customized genmake options, header files,
and modified code to the compile directory.
%
\item
{\tt ../tools/genmake -makefile} \\
Generate your Makefile (cf. Section ???).
%
\item
{\tt make depend} \\
Dependency analysis for the CPP pre-compiler (cf. Section ???).
%
\item
{\tt make small\_f} \\
This is the first difference between forward code compilation
and adjoint code generation and compilation.
Instead of going through the entire compilation process
(CPP precompiling -- {\tt .f}, object code generation -- {\tt .o},
linking of object files and libraries to generate executable),
only the CPP compiler is invoked at this stage to generate
the {\tt .f} files.
%
\item
{\tt cd ../adjoint} \\
{\tt make adtaf} or {\tt make adtamc} \\
Depending on whether you have TAF or TAMC at your disposal,
you'll choose {\tt adtaf} or {\tt adtamc} as your
make target for the {\it makefile} in the directory {\tt adjoint/}.
Several things happen at this stage.
%
\begin{enumerate}
%
\item
The initial template file {\it adjoint\_model.F} which is part
of the compiling list created by {\it genmake} is restored.
%
\item
All Fortran routines {\tt *.f} in {\tt bin/} are 
concatenated into a single file (it's current name is 
{\it tamc\_code.f}).
%
\item
Adjoint code is generated by TAMC or TAF.
The adjoint code is written to the file {\it tamc\_code\_ad.f}.
It contains all adjoint routines of the forward routines
concatenated in {\it tamc\_code.f}.
For a given forward routines {\tt subroutine routinename}
the adjoint routine is named {\tt adsubroutine routinename}
by default (that default can be changed via the flag 
{\tt -admark <markname>}).
Furthermore, it may contain modified code which
incorporates the translation of adjoint store directives
into specific Fortran code.
For a given forward routines {\tt subroutine routinename}
the modified routine is named {\tt mdsubroutine routinename}.
TAMC or TAF info is written to file
{\it tamc\_code.prot} or {\it taf.log}, respectively.
%
\end{enumerate}
%
\item
{\tt make adchange} \\
The multi-threading capability of the MITGCM requires a slight
change in the parameter list of some routines that are related to
to active file handling.
This post-processing invokes the sed script {\it adjoint\_ecco\_sed.com}
to insert the threading counter {\bf myThId} into the parameter list
of those subroutines.
The resulting code is written to file {\it tamc\_code\_sed\_ad.f}
and appended to the file {\it adjoint\_model.F}.
This concludes the adjoint code generation.
%
\item
{\tt cd ../bin} \\
{\tt make} \\
The file {\it adjoint\_model.F} now contains the full adjoint code.
All routines are now compiled.
%
\end{enumerate}

\paragraph{Adjoint code generation and compilation -- summary}
~ \\

\[
\boxed{
\begin{split}
 ~ & \mbox{\tt cd bin} \\
 ~ & \mbox{\tt ln -s ../verification/my\_experiment/code/.genmakerc .} \\
 ~ & \mbox{\tt ln -s ../verification/my\_experiment/code/*.[Fh] .} \\
 ~ & \mbox{\tt ../tools/genmake -makefile} \\
 ~ & \mbox{\tt make depend} \\
 ~ & \mbox{\tt make small\_f} \\
 ~ & \mbox{\tt cd ../adjoint} \\
 ~ & \mbox{\tt make adtaf <OR: make adtamc>} \\
 ~ & \mbox{\tt make adchange} \\
 ~ & \mbox{\tt cd ../bin} \\
 ~ & \mbox{\tt make} \\
\end{split}
}
\]
