s_examples/tracer_adjsens/doc_ad_examples.tex

% $Header: /u/gcmpack/manual/s_examples/tracer_adjsens/doc_ad_examples.tex,v 1.3 2010/08/30 23:09:20 jmc Exp $
% $Name:  $

%**********************************************************************
\section{Sensitivity of Air-Sea Exchange to Tracer Injection Site }
%\label{www:tutorials}
\label{sec:eg-simple-tracer-adjoint}
\label{sec_ad_setup_ex}
\label{sec:tutorialIII}
\begin{rawhtml}
<!-- CMIREDIR:sec_ad_setup_ex: -->
\end{rawhtml}
\begin{center}
(in directory: {\it verification/tutorial\_tracer\_adjsens/})
\end{center}
%**********************************************************************

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 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}
%\label{www:tutorials}

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}
%\label{www:tutorials}

For this work 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}
%\label{www:tutorials}

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}
%\label{www:tutorials}

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}
%\label{www:tutorials}

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}}
%\label{www:tutorials}
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}}
%\label{www:tutorials}

These files used to contain package-specific CPP-options
(see Section ref:ask-the-author).
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}}
%\label{www:tutorials}

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:ask-the-author).

\subsubsection{File {\it CPP\_OPTIONS.h}}
%\label{www:tutorials}

This file contains model-specific CPP options
(see Section ref:ask-the-author).
Most options are related to the forward model setup.
They are identical to the global steady circulation setup of
{\it verification/global\_ocean.90x40x15/}.
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}}
%\label{www:tutorials}

The CPP options of several AD-related packages are grouped
in this file:
%
\begin{itemize}
%
\item 
Overall ECCO-related execution modus: \\
These determine whether a pure forward run,
a sensitivity run or an iteration of optimization is
performed. These options are not needed in the present context.
%
\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:ask-the-author).
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.
The out and intermediate loop directivs are contained
in the files {\it checkpoint\_lev3\_directives.h}, 
{\it checkpoint\_lev2\_directives.h} (package {\it pkg/autodiff}). \\
%
\hspace*{4ex} {\tt \#define ALLOW\_AUTODIFF\_MONITOOR} \\
enables the monitoring of intermediate adjoint variables
(see Section ref:ask-the-author). \\
%
\hspace*{4ex} {\tt \#define ALLOW\_DIVIDED\_ADJOINT} \\
enables adjoint dump and restart
(see Section ref:ask-the-author).
%
\item Cost function package: {\it pkg/cost/} \\
This package contains all relevant routines for
initializing, accumulating and finalizing the cost function
(see Section ref:ask-the-author). \\
\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}}
%\label{www:tutorials}

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}}
%\label{www:tutorials}

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}.
The header file is enabled via the CPP-option
{\bf ALLOW\_AUTODIFF\_MONITOR}.

\subsubsection{File {\it tamc.h}}
%\label{www:tutorials}

This routine contains the dimensions for TAMC checkpointing
and some indices relevant for storing ky computations.
%
\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:ask-the-author).
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}}
%\label{www:tutorials}

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}.
The standard output for TAF is written to file
{\it taf.log}.
%
\begin{itemize}
\item [{\tt tamc}] {\tt
-input <variable names>
-output <variable name> -i4 -r4 ... \\
-toplevel <S/R name> -reverse <file names> 
}
\item [{\tt taf}] {\tt
-input <variable names>
-output <variable name> -i4 -r4 ... \\
-toplevel <S/R name> -reverse <file names> \\
-flow taf\_flow.log -nonew\_arg
}
\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 -i4 -r4} \\
~
%
\item {\tt -flow taf\_flow.log} \\
Will cause TAF to produce a flow listing file 
named {\it taf\_flow.log} in which 
the set of active and passive variables are identified
for each subroutine.
%
\item {\tt -nonew\_arg} \\
The default in the order of the parameter list of
adjoint routines has changed.
Before TAF 1.3 the default was compatible with the
TAMC-generated list. As of TAF 1.3 the order of adjoint
routine parameter lists is no longer copatible with TAMC.
To restore compatibility when using TAF 1.3 and higher,
this argument is needed.
It is currently crucial to use since all hand-written
adjoint routines refer to the TAMC default.
%
\end{itemize}


\subsubsection{The input parameter files}
%\label{www:tutorials}

\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}} ~ \\
%
Contains two-dimendional bathymetry information

\paragraph{File {\it windx.bin, windy.bin, salt.bin, theta.bin, 
SSS.bin, SST.bin}} ~ \\
%
These contain the initial values
(salinity, temperature, {\it salt.bin, theta.bin}), 
surface boundary values (surface wind stresses,
({\it windx.bin, windy.bin}), and surface restoring fields
({\it SSS.bin, SST.bin}).

\paragraph{File {\it pickup*}} ~ \\
%
Contains model state after model spinup.

\subsection{Compiling the model and its adjoint}
%\label{www:tutorials}

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 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
{\tt make adrestore, make ftlrestore} \\
The initial template files {\it adjoint\_model.F} 
and {\it tangentlinear\_model.F} in {\it pkg/autodiff}
which are part
of the compiling list created by {\it genmake} are restored.
%
\item {\tt make depend, make small\_f} \\
The {\tt bin/} directory is brought up to date,
i.e. for recent changes in header or source code
{\it .[Fh]}, corresponding {\it .f} routines are generated
or re-generated.
Note that here, only CPP precompiling is performed;
no object code {\it .o} is generated as yet.
Precompiling is necessary for TAMC to see the full code.
%
\item
{\tt make allcode} \\
All Fortran routines {\tt *.f} in {\tt bin/} are 
concatenated into a single file called
{\it tamc\_code.f}.
%
\item
{\tt make admodeltaf/admodeltamc} \\
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.
%
\item
{\tt make adchange} \\
The multi-threading capability of 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.
%
\end{enumerate}
%
\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}

N.B.: The targets {\tt make adtaf/adtamc} now comprise a
series of targets that in previous versions had to be
invoked separarely. This was probably preferable at a more
experimental stage, but has now been dropped in favour of
a more straightforward build process.


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

{\small
\[
\boxed{
\begin{aligned}
 ~ & \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 cd ../adjoint} \\
 ~ & \mbox{\tt make adtaf <OR: make adtamc>} \\
 ~ & \hspace*{6ex} \mbox{\tt                 contains the targets:} \\
 ~ & \hspace*{6ex} \mbox{\tt                 adrestore small\_f allcode admodeltaf/admodeltamc adchange} \\
 ~ & \mbox{\tt cd ../bin} \\
 ~ & \mbox{\tt make} \\
\end{aligned}
}
\]
}