18 |
can be generated automatically in this way. |
can be generated automatically in this way. |
19 |
|
|
20 |
The MITGCM has been adapted for use with the |
The MITGCM has been adapted for use with the |
21 |
Tangent linear and Adjoint Model Compiler (TAMC) and its succssor TAF |
Tangent linear and Adjoint Model Compiler (TAMC) and its successor TAF |
22 |
(Transformation of Algorithms in Fortran), developed |
(Transformation of Algorithms in Fortran), developed |
23 |
by Ralf Giering (\cite{gie-kam:98}, \cite{gie:99,gie:00}). |
by Ralf Giering (\cite{gie-kam:98}, \cite{gie:99,gie:00}). |
24 |
The first application of the adjoint of the MITGCM for senistivity |
The first application of the adjoint of the MITGCM for sensitivity |
25 |
studies has been published by \cite{maro-eta:99}. |
studies has been published by \cite{maro-eta:99}. |
26 |
\cite{sta-eta:97,sta-eta:01} use the MITGCM and its adjoint |
\cite{sta-eta:97,sta-eta:01} use the MITGCM and its adjoint |
27 |
for ocean state estimation studies. |
for ocean state estimation studies. |
28 |
|
In the following we shall refer to TAMC and TAF synonymously, |
29 |
|
except were explicitly stated otherwise. |
30 |
|
|
31 |
TAMC exploits the chain rule for computing the first |
TAMC exploits the chain rule for computing the first |
32 |
derivative of a function with |
derivative of a function with |
33 |
respect to a set of input variables. |
respect to a set of input variables. |
34 |
Treating a given forward code as a composition of operations -- |
Treating a given forward code as a composition of operations -- |
35 |
each line representing a compositional element -- the chain rule is |
each line representing a compositional element, the chain rule is |
36 |
rigorously applied to the code, line by line. The resulting |
rigorously applied to the code, line by line. The resulting |
37 |
tangent linear or adjoint code, |
tangent linear or adjoint code, |
38 |
then, may be thought of as the composition in |
then, may be thought of as the composition in |
39 |
forward or reverse order, respectively, of the |
forward or reverse order, respectively, of the |
40 |
Jacobian matrices of the forward code compositional elements. |
Jacobian matrices of the forward code's compositional elements. |
41 |
|
|
42 |
%********************************************************************** |
%********************************************************************** |
43 |
\section{Some basic algebra} |
\section{Some basic algebra} |
52 |
such as forcing functions) to the $n$-dimensional space |
such as forcing functions) to the $n$-dimensional space |
53 |
$V \subset I\!\!R^n$ of |
$V \subset I\!\!R^n$ of |
54 |
model output variable $\vec{v}=(v_1,\ldots,v_n)$ |
model output variable $\vec{v}=(v_1,\ldots,v_n)$ |
55 |
(model state, model diagnostcs, objective function, ...) |
(model state, model diagnostics, objective function, ...) |
56 |
under consideration, |
under consideration, |
57 |
% |
% |
58 |
\begin{equation} |
\begin{equation} |
107 |
$ M $ is just a matrix |
$ M $ is just a matrix |
108 |
which can readily be used to find the forward sensitivity of $\vec{v}$ to |
which can readily be used to find the forward sensitivity of $\vec{v}$ to |
109 |
perturbations in $u$, |
perturbations in $u$, |
110 |
but if there are very many input variables $(>>O(10^{6})$ for |
but if there are very many input variables $(\gg O(10^{6})$ for |
111 |
large-scale oceanographic application), it quickly becomes |
large-scale oceanographic application), it quickly becomes |
112 |
prohibitive to proceed directly as in (\ref{tangent_linear}), |
prohibitive to proceed directly as in (\ref{tangent_linear}), |
113 |
if the impact of each component $ {\bf e_{i}} $ is to be assessed. |
if the impact of each component $ {\bf e_{i}} $ is to be assessed. |
132 |
\label{compo} |
\label{compo} |
133 |
\end{eqnarray} |
\end{eqnarray} |
134 |
% |
% |
135 |
The linear approximation of $ {\cal J} $, |
The perturbation of $ {\cal J} $ around a fixed point $ {\cal J}_0 $, |
136 |
\[ |
\[ |
137 |
{\cal J} \, \approx \, {\cal J}_0 \, + \, \delta {\cal J} |
{\cal J} \, = \, {\cal J}_0 \, + \, \delta {\cal J} |
138 |
\] |
\] |
139 |
can be expressed in both bases of $ \vec{u} $ and $ \vec{v} $ |
can be expressed in both bases of $ \vec{u} $ and $ \vec{v} $ |
140 |
w.r.t. their corresponding inner product |
w.r.t. their corresponding inner product |
154 |
\label{deljidentity} |
\label{deljidentity} |
155 |
\end{equation} |
\end{equation} |
156 |
% |
% |
157 |
(note, that the gradient $ \nabla f $ is a pseudo-vector, therefore |
(note, that the gradient $ \nabla f $ is a co-vector, therefore |
158 |
its transpose is required in the above inner product). |
its transpose is required in the above inner product). |
159 |
Then, using the representation of |
Then, using the representation of |
160 |
$ \delta {\cal J} = |
$ \delta {\cal J} = |
170 |
\[ |
\[ |
171 |
A^{\ast} \, = \, A^T |
A^{\ast} \, = \, A^T |
172 |
\] |
\] |
173 |
and from eq. (\ref{tangent_linear}), we note that |
and from eq. (\ref{tangent_linear}), (\ref{deljidentity}), |
174 |
|
we note that |
175 |
(omitting $|$'s): |
(omitting $|$'s): |
176 |
% |
% |
177 |
\begin{equation} |
\begin{equation} |
207 |
$ \delta \vec{u}^{\ast} $ the adjoint variable of the control variable $ \vec{u} $. |
$ \delta \vec{u}^{\ast} $ the adjoint variable of the control variable $ \vec{u} $. |
208 |
|
|
209 |
The {\sf reverse} nature of the adjoint calculation can be readily |
The {\sf reverse} nature of the adjoint calculation can be readily |
210 |
seen as follows. Let us decompose ${\cal J}(u)$, thus: |
seen as follows. |
211 |
|
Consider a model integration which consists of $ \Lambda $ |
212 |
|
consecutive operations |
213 |
|
$ {\cal M}_{\Lambda} ( {\cal M}_{\Lambda-1} ( |
214 |
|
...... ( {\cal M}_{\lambda} ( |
215 |
|
...... |
216 |
|
( {\cal M}_{1} ( {\cal M}_{0}(\vec{u}) )))) $, |
217 |
|
where the ${\cal M}$'s could be the elementary steps, i.e. single lines |
218 |
|
in the code of the model, or successive time steps of the |
219 |
|
model integration, |
220 |
|
starting at step 0 and moving up to step $\Lambda$, with intermediate |
221 |
|
${\cal M}_{\lambda} (\vec{u}) = \vec{v}^{(\lambda+1)}$ and final |
222 |
|
${\cal M}_{\Lambda} (\vec{u}) = \vec{v}^{(\Lambda+1)} = \vec{v}$. |
223 |
|
Let ${\cal J}$ be a cost function which explicitly depends on the |
224 |
|
final state $\vec{v}$ only |
225 |
|
(this restriction is for clarity reasons only). |
226 |
|
% |
227 |
|
${\cal J}(u)$ may be decomposed according to: |
228 |
% |
% |
229 |
\begin{equation} |
\begin{equation} |
230 |
{\cal J}({\cal M}(\vec{u})) \, = \, |
{\cal J}({\cal M}(\vec{u})) \, = \, |
235 |
\label{compos} |
\label{compos} |
236 |
\end{equation} |
\end{equation} |
237 |
% |
% |
238 |
where the ${\cal M}$'s could be the elementary steps, i.e. single lines |
Then, according to the chain rule, the forward calculation reads, |
239 |
in the code of the model, |
in terms of the Jacobi matrices |
|
starting at step 0 and moving up to step $\Lambda$, with intermediate |
|
|
${\cal M}_{\lambda} (\vec{u}) = \vec{v}^{(\lambda+1)}$ and final |
|
|
${\cal M}_{\Lambda} (\vec{u}) = \vec{v}^{(\Lambda+1)} = \vec{v}$ |
|
|
Then, according to the chain rule the forward calculation reads in |
|
|
terms of the Jacobi matrices |
|
240 |
(we've omitted the $ | $'s which, nevertheless are important |
(we've omitted the $ | $'s which, nevertheless are important |
241 |
to the aspect of {\it tangent} linearity; |
to the aspect of {\it tangent} linearity; |
242 |
note also that per definition |
note also that by definition |
243 |
$ \langle \, \nabla _{v}{\cal J}^T \, , \, \delta \vec{v} \, \rangle |
$ \langle \, \nabla _{v}{\cal J}^T \, , \, \delta \vec{v} \, \rangle |
244 |
= \nabla_v {\cal J} \cdot \delta \vec{v} $ ) |
= \nabla_v {\cal J} \cdot \delta \vec{v} $ ) |
245 |
% |
% |
274 |
% |
% |
275 |
clearly expressing the reverse nature of the calculation. |
clearly expressing the reverse nature of the calculation. |
276 |
Eq. (\ref{reverse}) is at the heart of automatic adjoint compilers. |
Eq. (\ref{reverse}) is at the heart of automatic adjoint compilers. |
277 |
The intermediate steps $\lambda$ in |
If the intermediate steps $\lambda$ in |
278 |
eqn. (\ref{compos}) -- (\ref{reverse}) |
eqn. (\ref{compos}) -- (\ref{reverse}) |
279 |
could represent the model state (forward or adjoint) at each |
represent the model state (forward or adjoint) at each |
280 |
intermediate time step in which case |
intermediate time step as noted above, then correspondingly, |
281 |
$ {\cal M}(\vec{v}^{(\lambda)}) = \vec{v}^{(\lambda+1)} $, and correspondingly, |
$ M^T (\delta \vec{v}^{(\lambda) \, \ast}) = |
282 |
$ M^T (\delta \vec{v}^{(\lambda) \, \ast}) = \delta \vec{v}^{(\lambda-1) \, \ast} $, |
\delta \vec{v}^{(\lambda-1) \, \ast} $ for the adjoint variables. |
283 |
but they can also be viewed more generally as |
It thus becomes evident that the adjoint calculation also |
284 |
single lines of code in the numerical algorithm. |
yields the adjoint of each model state component |
285 |
In both cases it becomes evident that the adjoint calculation |
$ \vec{v}^{(\lambda)} $ at each intermediate step $ \lambda $, namely |
|
yields at the same time the adjoint of each model state component |
|
|
$ \vec{v}^{(\lambda)} $ at each intermediate step $ l $, namely |
|
286 |
% |
% |
287 |
\begin{equation} |
\begin{equation} |
288 |
\boxed{ |
\boxed{ |
298 |
% |
% |
299 |
in close analogy to eq. (\ref{adjoint}) |
in close analogy to eq. (\ref{adjoint}) |
300 |
We note in passing that that the $\delta \vec{v}^{(\lambda) \, \ast}$ |
We note in passing that that the $\delta \vec{v}^{(\lambda) \, \ast}$ |
301 |
are the Lagrange multipliers of the model state $ \vec{v}^{(\lambda)}$. |
are the Lagrange multipliers of the model equations which determine |
302 |
|
$ \vec{v}^{(\lambda)}$. |
303 |
|
|
304 |
In coponents, eq. (\ref{adjoint}) reads as follows. |
In components, eq. (\ref{adjoint}) reads as follows. |
305 |
Let |
Let |
306 |
\[ |
\[ |
307 |
\begin{array}{rclcrcl} |
\begin{array}{rclcrcl} |
322 |
\end{array} |
\end{array} |
323 |
\] |
\] |
324 |
denote the perturbations in $\vec{u}$ and $\vec{v}$, respectively, |
denote the perturbations in $\vec{u}$ and $\vec{v}$, respectively, |
325 |
and their adjoint varaiables; |
and their adjoint variables; |
326 |
further |
further |
327 |
\[ |
\[ |
328 |
M \, = \, \left( |
M \, = \, \left( |
409 |
$ \delta v^{(\lambda) \, \ast}_{j} = \frac{\partial}{\partial v^{(\lambda)}_{j}} |
$ \delta v^{(\lambda) \, \ast}_{j} = \frac{\partial}{\partial v^{(\lambda)}_{j}} |
410 |
{\cal J}^T $, $ j = 1, \ldots , n_{\lambda} $, |
{\cal J}^T $, $ j = 1, \ldots , n_{\lambda} $, |
411 |
for intermediate components, yielding |
for intermediate components, yielding |
412 |
\[ |
\begin{equation} |
413 |
\footnotesize |
\small |
414 |
|
\begin{split} |
415 |
\left( |
\left( |
416 |
\begin{array}{c} |
\begin{array}{c} |
417 |
\delta v^{(\lambda) \, \ast}_1 \\ |
\delta v^{(\lambda) \, \ast}_1 \\ |
419 |
\delta v^{(\lambda) \, \ast}_{n_{\lambda}} \\ |
\delta v^{(\lambda) \, \ast}_{n_{\lambda}} \\ |
420 |
\end{array} |
\end{array} |
421 |
\right) |
\right) |
422 |
\, = \, |
\, = & |
423 |
\left( |
\left( |
424 |
\begin{array}{ccc} |
\begin{array}{ccc} |
425 |
\frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_1} |
\frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_1} |
426 |
& \ldots & |
& \ldots \,\, \ldots & |
427 |
\frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_1} \\ |
\frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_1} \\ |
428 |
\vdots & ~ & \vdots \\ |
\vdots & ~ & \vdots \\ |
429 |
\frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_{n_{\lambda}}} |
\frac{\partial ({\cal M}_{\lambda})_1}{\partial v^{(\lambda)}_{n_{\lambda}}} |
430 |
& \ldots & |
& \ldots \,\, \ldots & |
431 |
\frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_{n_{\lambda}}} \\ |
\frac{\partial ({\cal M}_{\lambda})_{n_{\lambda+1}}}{\partial v^{(\lambda)}_{n_{\lambda}}} \\ |
432 |
\end{array} |
\end{array} |
433 |
\right) |
\right) |
|
% |
|
434 |
\cdot |
\cdot |
435 |
% |
% |
436 |
|
\\ ~ & ~ |
437 |
|
\\ ~ & |
438 |
|
% |
439 |
\left( |
\left( |
440 |
\begin{array}{ccc} |
\begin{array}{ccc} |
441 |
\frac{\partial ({\cal M}_{\lambda+1})_1}{\partial v^{(\lambda+1)}_1} |
\frac{\partial ({\cal M}_{\lambda+1})_1}{\partial v^{(\lambda+1)}_1} |
448 |
\frac{\partial ({\cal M}_{\lambda+1})_{n_{\lambda+2}}}{\partial v^{(\lambda+1)}_{n_{\lambda+1}}} \\ |
\frac{\partial ({\cal M}_{\lambda+1})_{n_{\lambda+2}}}{\partial v^{(\lambda+1)}_{n_{\lambda+1}}} \\ |
449 |
\end{array} |
\end{array} |
450 |
\right) |
\right) |
451 |
\cdot \ldots \ldots \cdot |
\cdot \, \ldots \, \cdot |
452 |
\left( |
\left( |
453 |
\begin{array}{c} |
\begin{array}{c} |
454 |
\delta v^{\ast}_1 \\ |
\delta v^{\ast}_1 \\ |
456 |
\delta v^{\ast}_{n} \\ |
\delta v^{\ast}_{n} \\ |
457 |
\end{array} |
\end{array} |
458 |
\right) |
\right) |
459 |
\] |
\end{split} |
460 |
|
\end{equation} |
461 |
|
|
462 |
Eq. (\ref{forward}) and (\ref{reverse}) are perhaps clearest in |
Eq. (\ref{forward}) and (\ref{reverse}) are perhaps clearest in |
463 |
showing the advantage of the reverse over the forward mode |
showing the advantage of the reverse over the forward mode |
468 |
{\it all} intermediate states $ \vec{v}^{(\lambda)} $) are sought. |
{\it all} intermediate states $ \vec{v}^{(\lambda)} $) are sought. |
469 |
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 |
470 |
$ \partial {\cal J} / \partial u_{i} $ in (\ref{forward}) |
$ \partial {\cal J} / \partial u_{i} $ in (\ref{forward}) |
471 |
a forward calulation has to be performed for each component seperately, |
a forward calculation has to be performed for each component separately, |
472 |
i.e. $ \delta \vec{u} = \delta u_{i} {\vec{e}_{i}} $ |
i.e. $ \delta \vec{u} = \delta u_{i} {\vec{e}_{i}} $ |
473 |
for the $i$-th forward calculation. |
for the $i$-th forward calculation. |
474 |
Then, (\ref{forward}) represents the |
Then, (\ref{forward}) represents the |
478 |
gradient $\nabla _{u}{\cal J}$ (and all intermediate gradients |
gradient $\nabla _{u}{\cal J}$ (and all intermediate gradients |
479 |
$\nabla _{v^{(\lambda)}}{\cal J}$) within a single reverse calculation. |
$\nabla _{v^{(\lambda)}}{\cal J}$) within a single reverse calculation. |
480 |
|
|
481 |
Note, that in case $ {\cal J} $ is a vector-valued function |
Note, that if $ {\cal J} $ is a vector-valued function |
482 |
of dimension $ l > 1 $, |
of dimension $ l > 1 $, |
483 |
eq. (\ref{reverse}) has to be modified according to |
eq. (\ref{reverse}) has to be modified according to |
484 |
\[ |
\[ |
486 |
\, = \, |
\, = \, |
487 |
\nabla_u {\cal J}^T \cdot \delta \vec{J} |
\nabla_u {\cal J}^T \cdot \delta \vec{J} |
488 |
\] |
\] |
489 |
where now $ \delta \vec{J} \in I\!\!R $ is a vector of dimenison $ l $. |
where now $ \delta \vec{J} \in I\!\!R^l $ is a vector of |
490 |
|
dimension $ l $. |
491 |
In this case $ l $ reverse simulations have to be performed |
In this case $ l $ reverse simulations have to be performed |
492 |
for each $ \delta J_{k}, \,\, k = 1, \ldots, l $. |
for each $ \delta J_{k}, \,\, k = 1, \ldots, l $. |
493 |
Then, the reverse mode is more efficient as long as |
Then, the reverse mode is more efficient as long as |
494 |
$ l < n $, otherwise the forward mode is preferable. |
$ l < n $, otherwise the forward mode is preferable. |
495 |
Stricly, the reverse mode is called adjoint mode only for |
Strictly, the reverse mode is called adjoint mode only for |
496 |
$ l = 1 $. |
$ l = 1 $. |
497 |
|
|
498 |
A detailed analysis of the underlying numerical operations |
A detailed analysis of the underlying numerical operations |
522 |
\paragraph{Example 2: |
\paragraph{Example 2: |
523 |
$ {\cal J} = \langle \, {\cal H}(\vec{v}) - \vec{d} \, , |
$ {\cal J} = \langle \, {\cal H}(\vec{v}) - \vec{d} \, , |
524 |
\, {\cal H}(\vec{v}) - \vec{d} \, \rangle $} ~ \\ |
\, {\cal H}(\vec{v}) - \vec{d} \, \rangle $} ~ \\ |
525 |
The cost function represents the quadratic model vs.data misfit. |
The cost function represents the quadratic model vs. data misfit. |
526 |
Here, $ \vec{d} $ is the data vector and $ {\cal H} $ represents the |
Here, $ \vec{d} $ is the data vector and $ {\cal H} $ represents the |
527 |
operator which maps the model state space onto the data space. |
operator which maps the model state space onto the data space. |
528 |
Then, $ \nabla_v {\cal J} $ takes the form |
Then, $ \nabla_v {\cal J} $ takes the form |
553 |
|
|
554 |
We note an important aspect of the forward vs. reverse |
We note an important aspect of the forward vs. reverse |
555 |
mode calculation. |
mode calculation. |
556 |
Because of the locality of the derivative, |
Because of the local character of the derivative |
557 |
|
(a derivative is defined w.r.t. a point along the trajectory), |
558 |
the intermediate results of the model trajectory |
the intermediate results of the model trajectory |
559 |
$\vec{v}^{(\lambda+1)}={\cal M}_{\lambda}(v^{(\lambda)})$ |
$\vec{v}^{(\lambda+1)}={\cal M}_{\lambda}(v^{(\lambda)})$ |
560 |
are needed to evaluate the intermediate Jacobian |
are needed to evaluate the intermediate Jacobian |
561 |
$M_{\lambda}|_{\vec{v}^{(\lambda)}} \, \delta \vec{v}^{(\lambda)} $. |
$M_{\lambda}|_{\vec{v}^{(\lambda)}} \, \delta \vec{v}^{(\lambda)} $. |
562 |
In the forward mode, the intermediate results are required |
In the forward mode, the intermediate results are required |
563 |
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}$, |
564 |
in the reverse mode they are required in the reverse order. |
but in the reverse mode they are required in the reverse order. |
565 |
Thus, in the reverse mode the trajectory of the forward model |
Thus, in the reverse mode the trajectory of the forward model |
566 |
integration ${\cal M}$ has to be stored to be available in the reverse |
integration ${\cal M}$ has to be stored to be available in the reverse |
567 |
calculation. Alternatively, the model state would have to be |
calculation. Alternatively, the complete model state up to the |
568 |
recomputed whenever its value is required. |
point of evaluation has to be recomputed whenever its value is required. |
569 |
|
|
570 |
A method to balance the amount of recomputations vs. |
A method to balance the amount of recomputations vs. |
571 |
storage requirements is called {\sf checkpointing} |
storage requirements is called {\sf checkpointing} |
572 |
(e.g. \cite{res-eta:98}). |
(e.g. \cite{res-eta:98}). |
573 |
It is depicted in Fig. ... for a 3-level checkpointing |
It is depicted in \ref{fig:3levelcheck} for a 3-level checkpointing |
574 |
[as concrete example, we give explicit numbers for a 3-day |
[as an example, we give explicit numbers for a 3-day |
575 |
integration with a 1-hourly timestep in square brackets]. |
integration with a 1-hourly timestep in square brackets]. |
576 |
\begin{itemize} |
\begin{itemize} |
577 |
% |
% |
579 |
In a first step, the model trajectory is subdivided into |
In a first step, the model trajectory is subdivided into |
580 |
$ {n}^{lev3} $ subsections [$ {n}^{lev3} $=3 1-day intervals], |
$ {n}^{lev3} $ subsections [$ {n}^{lev3} $=3 1-day intervals], |
581 |
with the label $lev3$ for this outermost loop. |
with the label $lev3$ for this outermost loop. |
582 |
The model is then integrated over the full trajectory, |
The model is then integrated along the full trajectory, |
583 |
and the model state stored only at every $ k_{i}^{lev3} $-th timestep |
and the model state stored only at every $ k_{i}^{lev3} $-th timestep |
584 |
[i.e. 3 times, at |
[i.e. 3 times, at |
585 |
$ i = 0,1,2 $ corresponding to $ k_{i}^{lev3} = 0, 24, 48 $]. |
$ i = 0,1,2 $ corresponding to $ k_{i}^{lev3} = 0, 24, 48 $]. |
586 |
% |
% |
587 |
\item [$lev2$] |
\item [$lev2$] |
588 |
In a second step each subsection is itself divided into |
In a second step each subsection itself is divided into |
589 |
$ {n}^{lev2} $ subsubsections |
$ {n}^{lev2} $ sub-subsections |
590 |
[$ {n}^{lev2} $=4 6-hour intervals per subsection]. |
[$ {n}^{lev2} $=4 6-hour intervals per subsection]. |
591 |
The model picks up at the last outermost dumped state |
The model picks up at the last outermost dumped state |
592 |
$ v_{k_{n}^{lev3}} $ and is integrated forward in time over |
$ v_{k_{n}^{lev3}} $ and is integrated forward in time along |
593 |
the last subsection, with the label $lev2$ for this |
the last subsection, with the label $lev2$ for this |
594 |
intermediate loop. |
intermediate loop. |
595 |
The model state is now stored only at every $ k_{i}^{lev2} $-th |
The model state is now stored at every $ k_{i}^{lev2} $-th |
596 |
timestep |
timestep |
597 |
[i.e. 4 times, at |
[i.e. 4 times, at |
598 |
$ 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 $]. |
599 |
% |
% |
600 |
\item [$lev1$] |
\item [$lev1$] |
601 |
Finally, the mode picks up at the last intermediate dump state |
Finally, the model picks up at the last intermediate dump state |
602 |
$ v_{k_{n}^{lev2}} $ and is integrated forward in time over |
$ v_{k_{n}^{lev2}} $ and is integrated forward in time along |
603 |
the last subsubsection, with the label $lev1$ for this |
the last sub-subsection, with the label $lev1$ for this |
604 |
intermediate loop. |
intermediate loop. |
605 |
Within this subsubsection only, the model state is stored |
Within this sub-subsection only, the model state is stored |
606 |
at every timestep |
at every timestep |
607 |
[i.e. every hour $ i=0,...,5$ corresponding to |
[i.e. every hour $ i=0,...,5$ corresponding to |
608 |
$ k_{i}^{lev1} = 66, 67, \ldots, 71 $]. |
$ k_{i}^{lev1} = 66, 67, \ldots, 71 $]. |
609 |
Thus, the final state $ v_n = v_{k_{n}^{lev1}} $ is reached |
Thus, the final state $ v_n = v_{k_{n}^{lev1}} $ is reached |
610 |
and the model state of all peceeding timesteps over the last |
and the model state of all proceeding timesteps along the last |
611 |
subsubsections are available, enabling integration backwards |
sub-subsections are available, enabling integration backwards |
612 |
in time over the last subsubsection. |
in time along the last sub-subsection. |
613 |
Thus, the adjoint can be computed over this last |
Thus, the adjoint can be computed along this last |
614 |
subsubsection $k_{n}^{lev2}$. |
sub-subsection $k_{n}^{lev2}$. |
615 |
% |
% |
616 |
\end{itemize} |
\end{itemize} |
617 |
% |
% |
618 |
This procedure is repeated consecutively for each previous |
This procedure is repeated consecutively for each previous |
619 |
subsubsection $k_{n-1}^{lev2}, \ldots, k_{1}^{lev2} $ |
sub-subsection $k_{n-1}^{lev2}, \ldots, k_{1}^{lev2} $ |
620 |
carrying the adjoint computation to the initial time |
carrying the adjoint computation to the initial time |
621 |
of the subsection $k_{n}^{lev3}$. |
of the subsection $k_{n}^{lev3}$. |
622 |
Then, the procedure is repeated for the previous subsection |
Then, the procedure is repeated for the previous subsection |
637 |
on the computing resources available. |
on the computing resources available. |
638 |
|
|
639 |
\begin{figure}[t!] |
\begin{figure}[t!] |
640 |
\centering |
\begin{center} |
641 |
%\psdraft |
%\psdraft |
642 |
\psfrag{v_k1^lev3}{\mathinfigure{v_{k_{1}^{lev3}}}} |
%\psfrag{v_k1^lev3}{\mathinfigure{v_{k_{1}^{lev3}}}} |
643 |
\psfrag{v_kn-1^lev3}{\mathinfigure{v_{k_{n-1}^{lev3}}}} |
%\psfrag{v_kn-1^lev3}{\mathinfigure{v_{k_{n-1}^{lev3}}}} |
644 |
\psfrag{v_kn^lev3}{\mathinfigure{v_{k_{n}^{lev3}}}} |
%\psfrag{v_kn^lev3}{\mathinfigure{v_{k_{n}^{lev3}}}} |
645 |
\psfrag{v_k1^lev2}{\mathinfigure{v_{k_{1}^{lev2}}}} |
%\psfrag{v_k1^lev2}{\mathinfigure{v_{k_{1}^{lev2}}}} |
646 |
\psfrag{v_kn-1^lev2}{\mathinfigure{v_{k_{n-1}^{lev2}}}} |
%\psfrag{v_kn-1^lev2}{\mathinfigure{v_{k_{n-1}^{lev2}}}} |
647 |
\psfrag{v_kn^lev2}{\mathinfigure{v_{k_{n}^{lev2}}}} |
%\psfrag{v_kn^lev2}{\mathinfigure{v_{k_{n}^{lev2}}}} |
648 |
\psfrag{v_k1^lev1}{\mathinfigure{v_{k_{1}^{lev1}}}} |
%\psfrag{v_k1^lev1}{\mathinfigure{v_{k_{1}^{lev1}}}} |
649 |
\psfrag{v_kn^lev1}{\mathinfigure{v_{k_{n}^{lev1}}}} |
%\psfrag{v_kn^lev1}{\mathinfigure{v_{k_{n}^{lev1}}}} |
650 |
\mbox{\epsfig{file=part5/checkpointing.eps, width=0.8\textwidth}} |
%\mbox{\epsfig{file=part5/checkpointing.eps, width=0.8\textwidth}} |
651 |
|
\resizebox{5.5in}{!}{\includegraphics{part5/checkpointing.eps}} |
652 |
%\psfull |
%\psfull |
653 |
\caption |
\end{center} |
654 |
{Schematic view of intermediate dump and restart for |
\caption{ |
655 |
|
Schematic view of intermediate dump and restart for |
656 |
3-level checkpointing.} |
3-level checkpointing.} |
657 |
\label{fig:erswns} |
\label{fig:3levelcheck} |
658 |
\end{figure} |
\end{figure} |
659 |
|
|
660 |
\subsection{Optimal perturbations} |
% \subsection{Optimal perturbations} |
661 |
\label{optpert} |
% \label{sec_optpert} |
662 |
|
|
663 |
|
|
664 |
\subsection{Error covariance estimate and Hessian matrix} |
% \subsection{Error covariance estimate and Hessian matrix} |
665 |
\label{sec_hessian} |
% \label{sec_hessian} |
666 |
|
|
667 |
\newpage |
\newpage |
668 |
|
|
671 |
\label{sec_ad_setup_ex} |
\label{sec_ad_setup_ex} |
672 |
%********************************************************************** |
%********************************************************************** |
673 |
|
|
674 |
The MITGCM has been adapted to enable AD using TAMC or TAF |
The MITGCM has been adapted to enable AD using TAMC or TAF. |
|
(we'll refer to TAMC and TAF interchangeably, except where |
|
|
distinctions are explicitly mentioned). |
|
675 |
The present description, therefore, is specific to the |
The present description, therefore, is specific to the |
676 |
use of TAMC as AD tool. |
use of TAMC or TAF as AD tool. |
677 |
The following sections describe the steps which are necessary to |
The following sections describe the steps which are necessary to |
678 |
generate a tangent linear or adjoint model of the MITGCM. |
generate a tangent linear or adjoint model of the MITGCM. |
679 |
We take as an example the sensitivity of carbon sequestration |
We take as an example the sensitivity of carbon sequestration |
680 |
in the ocean. |
in the ocean. |
681 |
The AD-relevant hooks in the code are sketched in |
The AD-relevant hooks in the code are sketched in |
682 |
\reffig{adthemodel}, \reffig{adthemain}. |
\ref{fig:adthemodel}, \ref{fig:adthemain}. |
683 |
|
|
684 |
\subsection{Overview of the experiment} |
\subsection{Overview of the experiment} |
685 |
|
|
686 |
We describe an adjoint sensitivity analysis of outgassing from |
We describe an adjoint sensitivity analysis of out-gassing from |
687 |
the ocean into the atmosphere of a carbon like tracer injected |
the ocean into the atmosphere of a carbon-like tracer injected |
688 |
into the ocean interior (see \cite{hil-eta:01}). |
into the ocean interior (see \cite{hil-eta:01}). |
689 |
|
|
690 |
\subsubsection{Passive tracer equation} |
\subsubsection{Passive tracer equation} |
691 |
|
|
692 |
For this work the MITGCM was augmented with a thermodynamically |
For this work the MITGCM was augmented with a thermodynamically |
693 |
inactive tracer, $C$. Tracer residing in the ocean |
inactive tracer, $C$. Tracer residing in the ocean |
694 |
model surface layer is outgassed according to a relaxation time scale, |
model surface layer is out-gassed according to a relaxation time scale, |
695 |
$\mu$. Within the ocean interior, the tracer is passively advected |
$\mu$. Within the ocean interior, the tracer is passively advected |
696 |
by the ocean model currents. The full equation for the time evolution |
by the ocean model currents. The full equation for the time evolution |
697 |
% |
% |
706 |
direct injection. |
direct injection. |
707 |
The velocity term, $U$, is the sum of the |
The velocity term, $U$, is the sum of the |
708 |
model Eulerian circulation and an eddy-induced velocity, the latter |
model Eulerian circulation and an eddy-induced velocity, the latter |
709 |
parameterized according to Gent/McWilliams (\cite{gen:90, dan:95}). |
parameterized according to Gent/McWilliams |
710 |
|
(\cite{gen-mcw:90, gen-eta:95}). |
711 |
The convection function, $\Gamma$, mixes $C$ vertically wherever the |
The convection function, $\Gamma$, mixes $C$ vertically wherever the |
712 |
fluid is locally statically unstable. |
fluid is locally statically unstable. |
713 |
|
|
714 |
The outgassing time scale, $\mu$, in eqn. (\ref{carbon_ddt}) |
The out-gassing time scale, $\mu$, in eqn. (\ref{carbon_ddt}) |
715 |
is set so that \( 1/\mu \sim 1 \ \mathrm{year} \) for the surface |
is set so that \( 1/\mu \sim 1 \ \mathrm{year} \) for the surface |
716 |
ocean and $\mu=0$ elsewhere. With this value, eqn. (\ref{carbon_ddt}) |
ocean and $\mu=0$ elsewhere. With this value, eqn. (\ref{carbon_ddt}) |
717 |
is valid as a prognostic equation for small perturbations in oceanic |
is valid as a prognostic equation for small perturbations in oceanic |
718 |
carbon concentrations. This configuration provides a |
carbon concentrations. This configuration provides a |
719 |
powerful tool for examining the impact of large-scale ocean circulation |
powerful tool for examining the impact of large-scale ocean circulation |
720 |
on $ CO_2 $ outgassing due to interior injections. |
on $ CO_2 $ out-gassing due to interior injections. |
721 |
As source we choose a constant in time injection of |
As source we choose a constant in time injection of |
722 |
$ S = 1 \,\, {\rm mol / s}$. |
$ S = 1 \,\, {\rm mol / s}$. |
723 |
|
|
728 |
geography and bathymetry. Twenty vertical layers are used with |
geography and bathymetry. Twenty vertical layers are used with |
729 |
vertical spacing ranging |
vertical spacing ranging |
730 |
from 50 m near the surface to 815 m at depth. |
from 50 m near the surface to 815 m at depth. |
731 |
Driven to steady-state by climatalogical wind-stress, heat and |
Driven to steady-state by climatological wind-stress, heat and |
732 |
fresh-water forcing the model reproduces well known large-scale |
fresh-water forcing the model reproduces well known large-scale |
733 |
features of the ocean general circulation. |
features of the ocean general circulation. |
734 |
|
|
735 |
\subsubsection{Outgassing cost function} |
\subsubsection{Out-gassing cost function} |
736 |
|
|
737 |
To quantify and understand outgassing due to injections of $C$ |
To quantify and understand out-gassing due to injections of $C$ |
738 |
in eqn. (\ref{carbon_ddt}), |
in eqn. (\ref{carbon_ddt}), |
739 |
we define a cost function $ {\cal J} $ that measures the total amount of |
we define a cost function $ {\cal J} $ that measures the total amount of |
740 |
tracer outgassed at each timestep: |
tracer out-gassed at each timestep: |
741 |
% |
% |
742 |
\begin{equation} |
\begin{equation} |
743 |
\label{cost_tracer} |
\label{cost_tracer} |
744 |
{\cal J}(t=T)=\int_{t=0}^{t=T}\int_{A} \mu C \, dA \, dt |
{\cal J}(t=T)=\int_{t=0}^{t=T}\int_{A} \mu C \, dA \, dt |
745 |
\end{equation} |
\end{equation} |
746 |
% |
% |
747 |
Equation(\ref{cost_tracer}) integrates the outgassing term, $\mu C$, |
Equation(\ref{cost_tracer}) integrates the out-gassing term, $\mu C$, |
748 |
from (\ref{carbon_ddt}) |
from (\ref{carbon_ddt}) |
749 |
over the entire ocean surface area, $A$, and accumulates it |
over the entire ocean surface area, $A$, and accumulates it |
750 |
up to time $T$. |
up to time $T$. |
751 |
Physically, ${\cal J}$ can be thought of as representing the amount of |
Physically, ${\cal J}$ can be thought of as representing the amount of |
752 |
$CO_2$ that our model predicts would be outgassed following an |
$CO_2$ that our model predicts would be out-gassed following an |
753 |
injection at rate $S$. |
injection at rate $S$. |
754 |
The sensitivity of ${\cal J}$ to the spatial location of $S$, |
The sensitivity of ${\cal J}$ to the spatial location of $S$, |
755 |
$\frac{\partial {\cal J}}{\partial S}$, |
$\frac{\partial {\cal J}}{\partial S}$, |
756 |
can be used to identify regions from which circulation |
can be used to identify regions from which circulation |
757 |
would cause $CO_2$ to rapidly outgas following injection |
would cause $CO_2$ to rapidly out-gas following injection |
758 |
and regions in which $CO_2$ injections would remain effectively |
and regions in which $CO_2$ injections would remain effectively |
759 |
sequesterd within the ocean. |
sequestered within the ocean. |
760 |
|
|
761 |
\subsection{Code configuration} |
\subsection{Code configuration} |
762 |
|
|
763 |
The model configuration for this experiment resides under the |
The model configuration for this experiment resides under the |
764 |
directory {\it verification/carbon/}. |
directory {\it verification/carbon/}. |
765 |
The code customisation routines are in {\it verification/carbon/code/}: |
The code customization routines are in {\it verification/carbon/code/}: |
766 |
% |
% |
767 |
\begin{itemize} |
\begin{itemize} |
768 |
% |
% |
798 |
% |
% |
799 |
\item {\it data.ctrl} |
\item {\it data.ctrl} |
800 |
% |
% |
801 |
|
\item {\it data.gmredi} |
802 |
|
% |
803 |
|
\item {\it data.grdchk} |
804 |
|
% |
805 |
|
\item {\it data.optim} |
806 |
|
% |
807 |
\item {\it data.pkg} |
\item {\it data.pkg} |
808 |
% |
% |
809 |
\item {\it eedata} |
\item {\it eedata} |
830 |
\end{itemize} |
\end{itemize} |
831 |
% |
% |
832 |
|
|
833 |
Below we describe the customisations of this files which are |
Below we describe the customizations of this files which are |
834 |
specific to this experiment. |
specific to this experiment. |
835 |
|
|
836 |
\subsubsection{File {\it .genmakerc}} |
\subsubsection{File {\it .genmakerc}} |
837 |
This file overwites default settings of {\it genmake}. |
This file overwrites default settings of {\it genmake}. |
838 |
In the present example it is used to switch on the following |
In the present example it is used to switch on the following |
839 |
packages which are related to automatic differentiation |
packages which are related to automatic differentiation |
840 |
and are disabled by default: \\ |
and are disabled by default: \\ |
841 |
\hspace*{4ex} {\tt set ENABLE=( autodiff cost ctrl ecco )} \\ |
\hspace*{4ex} {\tt set ENABLE=( autodiff cost ctrl ecco gmredi grdchk kpp )} \\ |
842 |
Other packages which are not needed are switched off: \\ |
Other packages which are not needed are switched off: \\ |
843 |
\hspace*{4ex} {\tt set DISABLE=( aim obcs zonal\_filt shap\_filt cal exf )} |
\hspace*{4ex} {\tt set DISABLE=( aim obcs zonal\_filt shap\_filt cal exf )} |
844 |
|
|
855 |
|
|
856 |
This file contains 'wrapper'-specific CPP options. |
This file contains 'wrapper'-specific CPP options. |
857 |
It only needs to be changed if the code is to be run |
It only needs to be changed if the code is to be run |
858 |
in parallel environment (see Section \ref{???}). |
in a parallel environment (see Section \ref{???}). |
859 |
|
|
860 |
\subsubsection{File {\it CPP\_OPTIONS.h}} |
\subsubsection{File {\it CPP\_OPTIONS.h}} |
861 |
|
|
864 |
Most options are related to the forward model setup. |
Most options are related to the forward model setup. |
865 |
They are identical to the global steady circulation setup of |
They are identical to the global steady circulation setup of |
866 |
{\it verification/exp2/}. |
{\it verification/exp2/}. |
867 |
The option specific to this experiment is \\ |
The three options specific to this experiment are \\ |
868 |
|
\hspace*{4ex} {\tt \#define ALLOW\_PASSIVE\_TRACER} \\ |
869 |
|
This flag enables the code to carry through the |
870 |
|
advection/diffusion of a passive tracer along the |
871 |
|
model integration. \\ |
872 |
\hspace*{4ex} {\tt \#define ALLOW\_MIT\_ADJOINT\_RUN} \\ |
\hspace*{4ex} {\tt \#define ALLOW\_MIT\_ADJOINT\_RUN} \\ |
873 |
This flag enables the inclusion of some AD-related fields |
This flag enables the inclusion of some AD-related fields |
874 |
concerning initialisation, link between control variables |
concerning initialization, link between control variables |
875 |
and forward model variables, and the call to the top-level |
and forward model variables, and the call to the top-level |
876 |
forward/adjoint subroutine {\it adthe\_main\_loop} |
forward/adjoint subroutine {\it adthe\_main\_loop} |
877 |
instead of {\it the\_main\_loop}. |
instead of {\it the\_main\_loop}. \\ |
878 |
|
\hspace*{4ex} {\tt \#define ALLOW\_GRADIENT\_CHECK} \\ |
879 |
|
This flag enables the gradient check package. |
880 |
|
After computing the unperturbed cost function and its gradient, |
881 |
|
a series of computations are performed for which \\ |
882 |
|
$\bullet$ an element of the control vector is perturbed \\ |
883 |
|
$\bullet$ the cost function w.r.t. the perturbed element is |
884 |
|
computed \\ |
885 |
|
$\bullet$ the difference between the perturbed and unperturbed |
886 |
|
cost function is computed to compute the finite difference gradient \\ |
887 |
|
$\bullet$ the finite difference gradient is compared with the |
888 |
|
adjoint-generated gradient. |
889 |
|
The gradient check package is further described in Section ???. |
890 |
|
|
891 |
\subsubsection{File {\it ECCO\_OPTIONS.h}} |
\subsubsection{File {\it ECCO\_OPTIONS.h}} |
892 |
|
|
907 |
(see Section \ref{???}). |
(see Section \ref{???}). |
908 |
In the present example a 3-level checkpointing is implemented. |
In the present example a 3-level checkpointing is implemented. |
909 |
The code contains the relevant store directives, common block |
The code contains the relevant store directives, common block |
910 |
and tape initialisations, storing key computation, |
and tape initializations, storing key computation, |
911 |
and loop index handling. |
and loop index handling. |
912 |
The checkpointing length at each level is defined in |
The checkpointing length at each level is defined in |
913 |
file {\it tamc.h}, cf. below. |
file {\it tamc.h}, cf. below. |
914 |
% |
% |
915 |
\item Cost function package: {\it pkg/cost/} \\ |
\item Cost function package: {\it pkg/cost/} \\ |
916 |
This package contains all relevant routines for |
This package contains all relevant routines for |
917 |
initialising, accumulating and finalizing the cost function |
initializing, accumulating and finalizing the cost function |
918 |
(see Section \ref{???}). \\ |
(see Section \ref{???}). \\ |
919 |
\hspace*{4ex} {\tt \#define ALLOW\_COST} \\ |
\hspace*{4ex} {\tt \#define ALLOW\_COST} \\ |
920 |
enables all general aspects of the cost function handling, |
enables all general aspects of the cost function handling, |
921 |
in particular the hooks in the foorward code for |
in particular the hooks in the forward code for |
922 |
initialising, accumulating and finalizing the cost function. \\ |
initializing, accumulating and finalizing the cost function. \\ |
923 |
\hspace*{4ex} {\tt \#define ALLOW\_COST\_TRACER} \\ |
\hspace*{4ex} {\tt \#define ALLOW\_COST\_TRACER} \\ |
924 |
includes the subroutine with the cost function for this |
includes the call to the cost function for this |
925 |
particular experiment, eqn. (\ref{cost_tracer}). |
particular experiment, eqn. (\ref{cost_tracer}). |
926 |
% |
% |
927 |
\item Control variable package: {\it pkg/ctrl/} \\ |
\item Control variable package: {\it pkg/ctrl/} \\ |
943 |
freshwater flux \\ |
freshwater flux \\ |
944 |
\hspace*{2ex} {\tt \#define ALLOW\_HFLUX0\_CONTROL} & |
\hspace*{2ex} {\tt \#define ALLOW\_HFLUX0\_CONTROL} & |
945 |
heat flux \\ |
heat flux \\ |
946 |
\hspace*{2ex} {\tt \#undef ALLOW\_DIFFKR\_CONTROL} & |
\hspace*{2ex} {\tt \#define ALLOW\_DIFFKR\_CONTROL} & |
947 |
diapycnal diffusivity \\ |
diapycnal diffusivity \\ |
948 |
\hspace*{2ex} {\tt \#undef ALLOW\_KAPPAGM\_CONTROL} & |
\hspace*{2ex} {\tt \#undef ALLOW\_KAPPAGM\_CONTROL} & |
949 |
isopycnal diffusivity \\ |
isopycnal diffusivity \\ |
958 |
\hspace*{4ex} {\tt sNx = 90} \\ |
\hspace*{4ex} {\tt sNx = 90} \\ |
959 |
\hspace*{4ex} {\tt sNy = 40} \\ |
\hspace*{4ex} {\tt sNy = 40} \\ |
960 |
\hspace*{4ex} {\tt Nr = 20} \\ |
\hspace*{4ex} {\tt Nr = 20} \\ |
961 |
It correpsponds to a single-tile/single-processor setup: |
It corresponds to a single-tile/single-processor setup: |
962 |
{\tt nSx = nSy = 1, nPx = nPy = 1}, |
{\tt nSx = nSy = 1, nPx = nPy = 1}, |
963 |
with standard overlap dimensioning |
with standard overlap dimensioning |
964 |
{\tt OLx = OLy = 3}. |
{\tt OLx = OLy = 3}. |
975 |
\hspace*{4ex} is related to {\it DYNVARS.h} \\ |
\hspace*{4ex} is related to {\it DYNVARS.h} \\ |
976 |
\hspace*{4ex} {\tt common /addynvars\_cd/} & |
\hspace*{4ex} {\tt common /addynvars\_cd/} & |
977 |
\hspace*{4ex} is related to {\it DYNVARS.h} \\ |
\hspace*{4ex} is related to {\it DYNVARS.h} \\ |
978 |
|
\hspace*{4ex} {\tt common /addynvars\_diffkr/} & |
979 |
|
\hspace*{4ex} is related to {\it DYNVARS.h} \\ |
980 |
|
\hspace*{4ex} {\tt common /addynvars\_kapgm/} & |
981 |
|
\hspace*{4ex} is related to {\it DYNVARS.h} \\ |
982 |
\hspace*{4ex} {\tt common /adtr1\_r/} & |
\hspace*{4ex} {\tt common /adtr1\_r/} & |
983 |
\hspace*{4ex} is related to {\it TR1.h} \\ |
\hspace*{4ex} is related to {\it TR1.h} \\ |
984 |
\hspace*{4ex} {\tt common /adffields/} & |
\hspace*{4ex} {\tt common /adffields/} & |
1003 |
3-level checkpointing is enabled, i.e. the timestepping |
3-level checkpointing is enabled, i.e. the timestepping |
1004 |
is divided into three different levels (see Section \ref{???}). |
is divided into three different levels (see Section \ref{???}). |
1005 |
The model state of the outermost ({\tt nchklev\_3}) and the |
The model state of the outermost ({\tt nchklev\_3}) and the |
1006 |
itermediate ({\tt nchklev\_2}) timestepping loop are stored to file |
intermediate ({\tt nchklev\_2}) timestepping loop are stored to file |
1007 |
(handled in {\it the\_main\_loop}). |
(handled in {\it the\_main\_loop}). |
1008 |
The innermost loop ({\tt nchklev\_1}) |
The innermost loop ({\tt nchklev\_1}) |
1009 |
avoids I/O by storing all required variables |
avoids I/O by storing all required variables |
1015 |
\hspace*{4ex} {\tt nchklev\_2 = 30 } \\ |
\hspace*{4ex} {\tt nchklev\_2 = 30 } \\ |
1016 |
\hspace*{4ex} {\tt nchklev\_3 = 60 } \\ |
\hspace*{4ex} {\tt nchklev\_3 = 60 } \\ |
1017 |
To guarantee that the checkpointing intervals span the entire |
To guarantee that the checkpointing intervals span the entire |
1018 |
integration period the relation \\ |
integration period the following relation must be satisfied: \\ |
1019 |
\hspace*{4ex} {\tt nchklev\_1*nchklev\_2*nchklev\_3 $ \ge $ nTimeSteps} \\ |
\hspace*{4ex} {\tt nchklev\_1*nchklev\_2*nchklev\_3 $ \ge $ nTimeSteps} \\ |
1020 |
where {\tt nTimeSteps} is either specified in {\it data} |
where {\tt nTimeSteps} is either specified in {\it data} |
1021 |
or computed via \\ |
or computed via \\ |
1029 |
% |
% |
1030 |
\end{itemize} |
\end{itemize} |
1031 |
|
|
1032 |
|
The following parameters may be worth describing: \\ |
1033 |
|
% |
1034 |
|
\hspace*{4ex} {\tt isbyte} \\ |
1035 |
|
\hspace*{4ex} {\tt maxpass} \\ |
1036 |
|
~ |
1037 |
|
|
1038 |
\subsubsection{File {\it makefile}} |
\subsubsection{File {\it makefile}} |
1039 |
|
|
1040 |
This file contains all relevant paramter flags and |
This file contains all relevant parameter flags and |
1041 |
lists to run TAMC. |
lists to run TAMC or TAF. |
1042 |
It is assumed that TAMC is available to you, either locally, |
It is assumed that TAMC is available to you, either locally, |
1043 |
being installed on your network, or remotely through the 'TAMC Utility'. |
being installed on your network, or remotely through the 'TAMC Utility'. |
1044 |
TAMC is called with the command {\tt tamc} followed by a |
TAMC is called with the command {\tt tamc} followed by a |
1049 |
\begin{itemize} |
\begin{itemize} |
1050 |
\item [{\tt tamc}] {\tt |
\item [{\tt tamc}] {\tt |
1051 |
-input <variable names> |
-input <variable names> |
1052 |
-output <variable name> ... \\ |
-output <variable name> -r4 ... \\ |
1053 |
-toplevel <S/R name> -reverse <file names> |
-toplevel <S/R name> -reverse <file names> |
1054 |
} |
} |
1055 |
\end{itemize} |
\end{itemize} |
1070 |
\item {\tt -reverse <file names>} \\ |
\item {\tt -reverse <file names>} \\ |
1071 |
Adjoint code is generated to compute the sensitivity of an |
Adjoint code is generated to compute the sensitivity of an |
1072 |
independent variable w.r.t. many dependent variables. |
independent variable w.r.t. many dependent variables. |
1073 |
The generated adjoint top-level routine computes the product |
In the discussion of Section ??? |
1074 |
|
the generated adjoint top-level routine computes the product |
1075 |
of the transposed Jacobian matrix $ M^T $ times |
of the transposed Jacobian matrix $ M^T $ times |
1076 |
the gradient vector $ \nabla_v J $. |
the gradient vector $ \nabla_v J $. |
1077 |
\\ |
\\ |
1078 |
{\tt <file names>} refers to the list of files {\it .f} which are to be |
{\tt <file names>} refers to the list of files {\it .f} which are to be |
1079 |
analyzed by TAMC. This list is generally smaller than the full list |
analyzed by TAMC. This list is generally smaller than the full list |
1080 |
of code to be compiled. The files not contained are either |
of code to be compiled. The files not contained are either |
1081 |
above the top-level routine (some initialisations), or are |
above the top-level routine (some initializations), or are |
1082 |
deliberately hidden from TAMC, either because hand-written |
deliberately hidden from TAMC, either because hand-written |
1083 |
adjoint routines exist, or the routines must not (or don't have to) |
adjoint routines exist, or the routines must not (or don't have to) |
1084 |
be differentiated. For each routine which is part of the flow tree |
be differentiated. For each routine which is part of the flow tree |
1085 |
of the top-level routine, but deliberately hidden from TAMC, |
of the top-level routine, but deliberately hidden from TAMC |
1086 |
|
(or for each package which contains such routines), |
1087 |
a corresponding file {\it .flow} exists containing flow directives |
a corresponding file {\it .flow} exists containing flow directives |
1088 |
for TAMC. |
for TAMC. |
1089 |
% |
% |
1090 |
|
\item {\tt -r4} \\ |
1091 |
|
~ |
1092 |
|
% |
1093 |
\end{itemize} |
\end{itemize} |
1094 |
|
|
1095 |
|
|
1096 |
\subsubsection{File {\it data}} |
\subsubsection{The input parameter files} |
1097 |
|
|
1098 |
|
\paragraph{File {\it data}} |
1099 |
|
|
1100 |
|
\paragraph{File {\it data.cost}} |
1101 |
|
|
1102 |
|
\paragraph{File {\it data.ctrl}} |
1103 |
|
|
1104 |
|
\paragraph{File {\it data.gmredi}} |
1105 |
|
|
1106 |
|
\paragraph{File {\it data.grdchk}} |
1107 |
|
|
1108 |
|
\paragraph{File {\it data.optim}} |
1109 |
|
|
1110 |
|
\paragraph{File {\it data.pkg}} |
1111 |
|
|
1112 |
\subsubsection{File {\it data.cost}} |
\paragraph{File {\it eedata}} |
1113 |
|
|
1114 |
\subsubsection{File {\it data.ctrl}} |
\paragraph{File {\it topog.bin}} |
1115 |
|
|
1116 |
\subsubsection{File {\it data.pkg}} |
\paragraph{File {\it windx.bin, windy.bin}} |
1117 |
|
|
1118 |
\subsubsection{File {\it eedata}} |
\paragraph{File {\it salt.bin, theta.bin}} |
1119 |
|
|
1120 |
\subsubsection{File {\it topog.bin}} |
\paragraph{File {\it SSS.bin, SST.bin}} |
1121 |
|
|
1122 |
\subsubsection{File {\it windx.bin, windy.bin}} |
\paragraph{File {\it pickup*}} |
1123 |
|
|
1124 |
\subsubsection{File {\it salt.bin, theta.bin}} |
\subsection{Compiling the model and its adjoint} |
1125 |
|
|
1126 |
\subsubsection{File {\it SSS.bin, SST.bin}} |
The built process of the adjoint model is slightly more |
1127 |
|
complex than that of compiling the forward code. |
1128 |
|
The main reason is that the adjoint code generation requires |
1129 |
|
a specific list of routines that are to be differentiated |
1130 |
|
(as opposed to the automatic generation of a list of |
1131 |
|
files to be compiled by genmake). |
1132 |
|
This list excludes routines that don't have to be or must not be |
1133 |
|
differentiated. For some of the latter routines flow directives |
1134 |
|
may be necessary, a list of which has to be given as well. |
1135 |
|
For this reason, a separate {\it makefile} is currently |
1136 |
|
maintained in the directory {\tt adjoint/}. This |
1137 |
|
makefile is responsible for the adjoint code generation. |
1138 |
|
|
1139 |
\subsubsection{File {\it pickup*}} |
In the following we describe the build process step by step, |
1140 |
|
assuming you are in the directory {\tt bin/}. |
1141 |
|
A summary of steps to follow is given at the end. |
1142 |
|
|
1143 |
\subsection{Compiling the model and its adjoint} |
\paragraph{Adjoint code generation and compilation -- step by step} |
1144 |
|
|
1145 |
|
\begin{enumerate} |
1146 |
|
% |
1147 |
|
\item |
1148 |
|
{\tt ln -s ../verification/???/code/.genmakerc .} \\ |
1149 |
|
{\tt ln -s ../verification/???/code/*.[Fh] .} \\ |
1150 |
|
Link your customized genmake options, header files, |
1151 |
|
and modified code to the compile directory. |
1152 |
|
% |
1153 |
|
\item |
1154 |
|
{\tt ../tools/genmake -makefile} \\ |
1155 |
|
Generate your Makefile (cf. Section ???). |
1156 |
|
% |
1157 |
|
\item |
1158 |
|
{\tt make depend} \\ |
1159 |
|
Dependency analysis for the CPP pre-compiler (cf. Section ???). |
1160 |
|
% |
1161 |
|
\item |
1162 |
|
{\tt make small\_f} \\ |
1163 |
|
This is the first difference between forward code compilation |
1164 |
|
and adjoint code generation and compilation. |
1165 |
|
Instead of going through the entire compilation process |
1166 |
|
(CPP precompiling -- {\tt .f}, object code generation -- {\tt .o}, |
1167 |
|
linking of object files and libraries to generate executable), |
1168 |
|
only the CPP compiler is invoked at this stage to generate |
1169 |
|
the {\tt .f} files. |
1170 |
|
% |
1171 |
|
\item |
1172 |
|
{\tt cd ../adjoint} \\ |
1173 |
|
{\tt make adtaf} or {\tt make adtamc} \\ |
1174 |
|
Depending on whether you have TAF or TAMC at your disposal, |
1175 |
|
you'll choose {\tt adtaf} or {\tt adtamc} as your |
1176 |
|
make target for the {\it makefile} in the directory {\tt adjoint/}. |
1177 |
|
Several things happen at this stage. |
1178 |
|
% |
1179 |
|
\begin{enumerate} |
1180 |
|
% |
1181 |
|
\item |
1182 |
|
The initial template file {\it adjoint\_model.F} which is part |
1183 |
|
of the compiling list created by {\it genmake} is restored. |
1184 |
|
% |
1185 |
|
\item |
1186 |
|
All Fortran routines {\tt *.f} in {\tt bin/} are |
1187 |
|
concatenated into a single file (it's current name is |
1188 |
|
{\it tamc\_code.f}). |
1189 |
|
% |
1190 |
|
\item |
1191 |
|
Adjoint code is generated by TAMC or TAF. |
1192 |
|
The adjoint code is written to the file {\it tamc\_code\_ad.f}. |
1193 |
|
It contains all adjoint routines of the forward routines |
1194 |
|
concatenated in {\it tamc\_code.f}. |
1195 |
|
For a given forward routines {\tt subroutine routinename} |
1196 |
|
the adjoint routine is named {\tt adsubroutine routinename} |
1197 |
|
by default (that default can be changed via the flag |
1198 |
|
{\tt -admark <markname>}). |
1199 |
|
Furthermore, it may contain modified code which |
1200 |
|
incorporates the translation of adjoint store directives |
1201 |
|
into specific Fortran code. |
1202 |
|
For a given forward routines {\tt subroutine routinename} |
1203 |
|
the modified routine is named {\tt mdsubroutine routinename}. |
1204 |
|
TAMC or TAF info is written to file |
1205 |
|
{\it tamc\_code.prot} or {\it taf.log}, respectively. |
1206 |
|
% |
1207 |
|
\end{enumerate} |
1208 |
|
% |
1209 |
|
\item |
1210 |
|
{\tt make adchange} \\ |
1211 |
|
The multi-threading capability of the MITGCM requires a slight |
1212 |
|
change in the parameter list of some routines that are related to |
1213 |
|
to active file handling. |
1214 |
|
This post-processing invokes the sed script {\it adjoint\_ecco\_sed.com} |
1215 |
|
to insert the threading counter {\bf myThId} into the parameter list |
1216 |
|
of those subroutines. |
1217 |
|
The resulting code is written to file {\it tamc\_code\_sed\_ad.f} |
1218 |
|
and appended to the file {\it adjoint\_model.F}. |
1219 |
|
This concludes the adjoint code generation. |
1220 |
|
% |
1221 |
|
\item |
1222 |
|
{\tt cd ../bin} \\ |
1223 |
|
{\tt make} \\ |
1224 |
|
The file {\it adjoint\_model.F} now contains the full adjoint code. |
1225 |
|
All routines are now compiled. |
1226 |
|
% |
1227 |
|
\end{enumerate} |
1228 |
|
|
1229 |
|
\paragraph{Adjoint code generation and compilation -- summary} |
1230 |
|
~ \\ |
1231 |
|
|
1232 |
|
\[ |
1233 |
|
\boxed{ |
1234 |
|
\begin{split} |
1235 |
|
~ & \mbox{\tt cd bin} \\ |
1236 |
|
~ & \mbox{\tt ln -s ../verification/my\_experiment/code/.genmakerc .} \\ |
1237 |
|
~ & \mbox{\tt ln -s ../verification/my\_experiment/code/*.[Fh] .} \\ |
1238 |
|
~ & \mbox{\tt ../tools/genmake -makefile} \\ |
1239 |
|
~ & \mbox{\tt make depend} \\ |
1240 |
|
~ & \mbox{\tt make small\_f} \\ |
1241 |
|
~ & \mbox{\tt cd ../adjoint} \\ |
1242 |
|
~ & \mbox{\tt make adtaf <OR: make adtamc>} \\ |
1243 |
|
~ & \mbox{\tt make adchange} \\ |
1244 |
|
~ & \mbox{\tt cd ../bin} \\ |
1245 |
|
~ & \mbox{\tt make} \\ |
1246 |
|
\end{split} |
1247 |
|
} |
1248 |
|
\] |
1249 |
|
|
1250 |
\newpage |
\newpage |
1251 |
|
|
1252 |
%********************************************************************** |
%********************************************************************** |
1253 |
\section{TLM and ADM code generation in general} |
\section{TLM and ADM generation in general} |
1254 |
\label{sec_ad_setup_gen} |
\label{sec_ad_setup_gen} |
1255 |
%********************************************************************** |
%********************************************************************** |
1256 |
|
|
1258 |
the parts of the code that are relevant for automatic |
the parts of the code that are relevant for automatic |
1259 |
differentiation using the software tool TAMC. |
differentiation using the software tool TAMC. |
1260 |
|
|
1261 |
\subsection{The cost function (dependent variable)} |
\input{part5/doc_ad_the_model} |
1262 |
|
|
1263 |
|
The basic flow is depicted in \ref{fig:adthemodel}. |
1264 |
|
If the option {\tt ALLOW\_AUTODIFF\_TAMC} is defined, the driver routine |
1265 |
|
{\it the\_model\_main}, instead of calling {\it the\_main\_loop}, |
1266 |
|
invokes the adjoint of this routine, {\it adthe\_main\_loop}, |
1267 |
|
which is the toplevel routine in terms of reverse mode computation. |
1268 |
|
The routine {\it adthe\_main\_loop} has been generated using TAMC. |
1269 |
|
It contains both the forward integration of the full model, |
1270 |
|
any additional storing that is required for efficient checkpointing, |
1271 |
|
and the reverse integration of the adjoint model. |
1272 |
|
The structure of {\it adthe\_main\_loop} has been strongly |
1273 |
|
simplified for clarification; in particular, no checkpointing |
1274 |
|
procedures are shown here. |
1275 |
|
Prior to the call of {\it adthe\_main\_loop}, the routine |
1276 |
|
{\it ctrl\_unpack} is invoked to unpack the control vector, |
1277 |
|
and following that call, the routine {\it ctrl\_pack} |
1278 |
|
is invoked to pack the control vector |
1279 |
|
(cf. Section \ref{section_ctrl}). |
1280 |
|
If gradient checks are to be performed, the option |
1281 |
|
{\tt ALLOW\_GRADIENT\_CHECK} is defined. In this case |
1282 |
|
the driver routine {\it grdchk\_main} is called after |
1283 |
|
the gradient has been computed via the adjoint |
1284 |
|
(cf. Section \ref{section_grdchk}). |
1285 |
|
|
1286 |
|
\subsection{The cost function (dependent variable) |
1287 |
|
\label{section_cost}} |
1288 |
|
|
1289 |
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}. |
1290 |
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 |
1292 |
The input is referred to as the |
The input is referred to as the |
1293 |
{\sf independent variables} or {\sf control variables}. |
{\sf independent variables} or {\sf control variables}. |
1294 |
All aspects relevant to the treatment of the cost function $ {\cal J} $ |
All aspects relevant to the treatment of the cost function $ {\cal J} $ |
1295 |
(parameter setting, initialisation, incrementation, |
(parameter setting, initialization, accumulation, |
1296 |
final evaluation), are controled by the package {\it pkg/cost}. |
final evaluation), are controlled by the package {\it pkg/cost}. |
1297 |
|
|
1298 |
|
\begin{figure}[h!] |
1299 |
|
\input{part5/doc_cost_flow} |
1300 |
|
\caption{~} |
1301 |
|
\label{fig:costflow} |
1302 |
|
\end{figure} |
1303 |
|
|
1304 |
\subsubsection{genmake and CPP options} |
\subsubsection{genmake and CPP options} |
1305 |
% |
% |
1319 |
\begin{enumerate} |
\begin{enumerate} |
1320 |
% |
% |
1321 |
\item {\it genmake}: \\ |
\item {\it genmake}: \\ |
1322 |
Change the default settngs in the file {\it genmake} by adding |
Change the default settings in the file {\it genmake} by adding |
1323 |
{\bf cost} to the {\bf enable} list (not recommended). |
{\bf cost} to the {\bf enable} list (not recommended). |
1324 |
% |
% |
1325 |
\item {\it .genmakerc}: \\ |
\item {\it .genmakerc}: \\ |
1332 |
{\tt genmake -enable=cost}. |
{\tt genmake -enable=cost}. |
1333 |
% |
% |
1334 |
\end{enumerate} |
\end{enumerate} |
|
Since the cost function is usually used in conjunction with |
|
|
automatic differentiation, the CPP option |
|
|
{\bf ALLOW\_ADJOINT\_RUN} should be defined |
|
|
(file {\it CPP\_OPTIONS.h}). |
|
1335 |
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}. |
1336 |
Each specific cost function contribution has its own option. |
Each specific cost function contribution has its own option. |
1337 |
For the present example the option is {\bf ALLOW\_COST\_TRACER}. |
For the present example the option is {\bf ALLOW\_COST\_TRACER}. |
1338 |
All cost-specific options are set in {\it ECCO\_CPPOPTIONS.h} |
All cost-specific options are set in {\it ECCO\_CPPOPTIONS.h} |
1339 |
|
Since the cost function is usually used in conjunction with |
1340 |
|
automatic differentiation, the CPP option |
1341 |
|
{\bf ALLOW\_ADJOINT\_RUN} should be defined |
1342 |
|
(file {\it CPP\_OPTIONS.h}). |
1343 |
|
|
1344 |
\subsubsection{Initialisation} |
\subsubsection{Initialization} |
1345 |
% |
% |
1346 |
The initialisation of the {\it cost} package is readily enabled |
The initialization of the {\it cost} package is readily enabled |
1347 |
as soon as the CPP option {\bf ALLOW\_ADJOINT\_RUN} is defined. |
as soon as the CPP option {\bf ALLOW\_ADJOINT\_RUN} is defined. |
1348 |
% |
% |
1349 |
\begin{itemize} |
\begin{itemize} |
1374 |
} |
} |
1375 |
\\ |
\\ |
1376 |
This S/R |
This S/R |
1377 |
initialises the different cost function contributions. |
initializes the different cost function contributions. |
1378 |
The contribtion for the present example is {\bf objf\_tracer} |
The contribution for the present example is {\bf objf\_tracer} |
1379 |
which is defined on each tile (bi,bj). |
which is defined on each tile (bi,bj). |
1380 |
% |
% |
1381 |
\end{itemize} |
\end{itemize} |
1382 |
% |
% |
1383 |
\subsubsection{Incrementation} |
\subsubsection{Accumulation} |
1384 |
% |
% |
1385 |
\begin{itemize} |
\begin{itemize} |
1386 |
% |
% |
1428 |
tamc -output 'fc' ... |
tamc -output 'fc' ... |
1429 |
\end{verbatim} |
\end{verbatim} |
1430 |
|
|
1431 |
\begin{figure}[t!] |
%%%% \end{document} |
|
\input{part5/doc_ad_the_model} |
|
|
\label{fig:adthemodel} |
|
|
\caption{~} |
|
|
\end{figure} |
|
1432 |
|
|
1433 |
\begin{figure} |
\begin{figure} |
1434 |
\input{part5/doc_ad_the_main} |
\input{part5/doc_ad_the_main} |
|
\label{fig:adthemain} |
|
1435 |
\caption{~} |
\caption{~} |
1436 |
|
\label{fig:adthemain} |
1437 |
\end{figure} |
\end{figure} |
1438 |
|
|
1439 |
\subsection{The control variables (independent variables)} |
\subsection{The control variables (independent variables) |
1440 |
|
\label{section_ctrl}} |
1441 |
|
|
1442 |
The control variables are a subset of the model input |
The control variables are a subset of the model input |
1443 |
(initial conditions, boundary conditions, model parameters). |
(initial conditions, boundary conditions, model parameters). |
1444 |
Here we identify them with the variable $ \vec{u} $. |
Here we identify them with the variable $ \vec{u} $. |
1445 |
All intermediate variables whose derivative w.r.t. control |
All intermediate variables whose derivative w.r.t. control |
1446 |
variables don't vanish are called {\sf active variables}. |
variables do not vanish are called {\sf active variables}. |
1447 |
All subroutines whose derivative w.r.t. the control variables |
All subroutines whose derivative w.r.t. the control variables |
1448 |
don't vanish are called {\sf active routines}. |
don't vanish are called {\sf active routines}. |
1449 |
Read and write operations from and to file can be viewed |
Read and write operations from and to file can be viewed |
1451 |
active variables are written and from which active variables |
active variables are written and from which active variables |
1452 |
are read are called {\sf active files}. |
are read are called {\sf active files}. |
1453 |
All aspects relevant to the treatment of the control variables |
All aspects relevant to the treatment of the control variables |
1454 |
(parameter setting, initialisation, perturbation) |
(parameter setting, initialization, perturbation) |
1455 |
are controled by the package {\it pkg/ctrl}. |
are controlled by the package {\it pkg/ctrl}. |
1456 |
|
|
1457 |
|
\begin{figure}[h!] |
1458 |
|
\input{part5/doc_ctrl_flow} |
1459 |
|
\caption{~} |
1460 |
|
\label{fig:ctrlflow} |
1461 |
|
\end{figure} |
1462 |
|
|
1463 |
\subsubsection{genmake and CPP options} |
\subsubsection{genmake and CPP options} |
1464 |
% |
% |
1478 |
Each control variable is enabled via its own CPP option |
Each control variable is enabled via its own CPP option |
1479 |
in {\it ECCO\_CPPOPTIONS.h}. |
in {\it ECCO\_CPPOPTIONS.h}. |
1480 |
|
|
1481 |
\subsubsection{Initialisation} |
\subsubsection{Initialization} |
1482 |
% |
% |
1483 |
\begin{itemize} |
\begin{itemize} |
1484 |
% |
% |
1518 |
variables in the MITGCM need to be addressed. |
variables in the MITGCM need to be addressed. |
1519 |
First, in order to save memory, the control variable arrays |
First, in order to save memory, the control variable arrays |
1520 |
are not kept in memory, but rather read from file and added |
are not kept in memory, but rather read from file and added |
1521 |
to the initial (or first guess) fields. |
to the initial fields during the model initialization phase. |
1522 |
Similarly, the corresponding adjoint fields which represent |
Similarly, the corresponding adjoint fields which represent |
1523 |
the gradient of the cost function w.r.t. the control variables |
the gradient of the cost function w.r.t. the control variables |
1524 |
are written to to file. |
are written to file at the end of the adjoint integration. |
1525 |
Second, in addition to the files holding the 2-dim. and 3-dim. |
Second, in addition to the files holding the 2-dim. and 3-dim. |
1526 |
control variables and the gradient, a 1-dim. {\sf control vector} |
control variables and the corresponding cost gradients, |
1527 |
|
a 1-dim. {\sf control vector} |
1528 |
and {\sf gradient vector} are written to file. They contain |
and {\sf gradient vector} are written to file. They contain |
1529 |
only the wet points of the control variables and the corresponding |
only the wet points of the control variables and the corresponding |
1530 |
gradient. |
gradient. |
1531 |
This leads to a significant data compression. |
This leads to a significant data compression. |
1532 |
Furthermore, the control and the gradient vector can be passed to a |
Furthermore, an option is available |
1533 |
|
({\tt ALLOW\_NONDIMENSIONAL\_CONTROL\_IO}) to |
1534 |
|
non-dimensionalise the control and gradient vector, |
1535 |
|
which otherwise would contain different pieces of different |
1536 |
|
magnitudes and units. |
1537 |
|
Finally, the control and gradient vector can be passed to a |
1538 |
minimization routine if an update of the control variables |
minimization routine if an update of the control variables |
1539 |
is sought as part of a minimization exercise. |
is sought as part of a minimization exercise. |
1540 |
|
|
1545 |
|
|
1546 |
\subsubsection{Perturbation of the independent variables} |
\subsubsection{Perturbation of the independent variables} |
1547 |
% |
% |
1548 |
The dependency chain for differentiation starts |
The dependency flow for differentiation w.r.t. the controls |
1549 |
with adding a perturbation onto the the input variable, |
starts with adding a perturbation onto the input variable, |
1550 |
thus defining the independent or control variables for TAMC. |
thus defining the independent or control variables for TAMC. |
1551 |
Three classes of controls may be considered: |
Three types of controls may be considered: |
1552 |
% |
% |
1553 |
\begin{itemize} |
\begin{itemize} |
1554 |
% |
% |
1563 |
Consider as an example the initial tracer distribution |
Consider as an example the initial tracer distribution |
1564 |
{\bf tr1} as control variable. |
{\bf tr1} as control variable. |
1565 |
After {\bf tr1} has been initialised in |
After {\bf tr1} has been initialised in |
1566 |
{\it ini\_tr1} (dynamical variables including |
{\it ini\_tr1} (dynamical variables such as |
1567 |
temperature and salinity are initialised in {\it ini\_fields}), |
temperature and salinity are initialised in {\it ini\_fields}), |
1568 |
a perturbation anomaly is added to the field in S/R |
a perturbation anomaly is added to the field in S/R |
1569 |
{\it ctrl\_map\_ini} |
{\it ctrl\_map\_ini} |
1576 |
\end{split} |
\end{split} |
1577 |
\end{equation} |
\end{equation} |
1578 |
% |
% |
1579 |
In principle {\bf xx\_tr1} is a 3-dim. global array |
{\bf xx\_tr1} is a 3-dim. global array |
1580 |
holding the perturbation. In the case of a simple |
holding the perturbation. In the case of a simple |
1581 |
sensitivity study this array is identical to zero. |
sensitivity study this array is identical to zero. |
1582 |
However, it's specification is essential since TAMC |
However, it's specification is essential in the context |
1583 |
|
of automatic differentiation since TAMC |
1584 |
treats the corresponding line in the code symbolically |
treats the corresponding line in the code symbolically |
1585 |
when determining the differentiation chain and its origin. |
when determining the differentiation chain and its origin. |
1586 |
Thus, the variable names are part of the argument list |
Thus, the variable names are part of the argument list |
1598 |
and an 'active read' routine of the adjoint support |
and an 'active read' routine of the adjoint support |
1599 |
package {\it pkg/autodiff} is invoked. |
package {\it pkg/autodiff} is invoked. |
1600 |
The read-procedure is tagged with the variable |
The read-procedure is tagged with the variable |
1601 |
{\bf xx\_tr1\_dummy} enabbling TAMC to recognize the |
{\bf xx\_tr1\_dummy} enabling TAMC to recognize the |
1602 |
initialisation of the perturbation. |
initialization of the perturbation. |
1603 |
The modified call of TAMC thus reads |
The modified call of TAMC thus reads |
1604 |
% |
% |
1605 |
\begin{verbatim} |
\begin{verbatim} |
1620 |
to a variable assignment. Its derivative corresponds |
to a variable assignment. Its derivative corresponds |
1621 |
to a write statement of the adjoint variable. |
to a write statement of the adjoint variable. |
1622 |
The 'active file' routines have been designed |
The 'active file' routines have been designed |
1623 |
to support active read and corresponding active write |
to support active read and corresponding adjoint active write |
1624 |
operations. |
operations (and vice versa). |
1625 |
% |
% |
1626 |
\item |
\item |
1627 |
\fbox{ |
\fbox{ |
1638 |
Note however an important difference: |
Note however an important difference: |
1639 |
Since the boundary values are time dependent with a new |
Since the boundary values are time dependent with a new |
1640 |
forcing field applied at each time steps, |
forcing field applied at each time steps, |
1641 |
the general problem may be be thought of as |
the general problem may be thought of as |
1642 |
a new control variable at each time step, i.e. |
a new control variable at each time step |
1643 |
|
(or, if the perturbation is averaged over a certain period, |
1644 |
|
at each $ N $ timesteps), i.e. |
1645 |
\[ |
\[ |
1646 |
u_{\rm forcing} \, = \, |
u_{\rm forcing} \, = \, |
1647 |
\{ \, u_{\rm forcing} ( t_n ) \, \}_{ |
\{ \, u_{\rm forcing} ( t_n ) \, \}_{ |
1666 |
% |
% |
1667 |
This routine is not yet implemented, but would proceed |
This routine is not yet implemented, but would proceed |
1668 |
proceed along the same lines as the initial value sensitivity. |
proceed along the same lines as the initial value sensitivity. |
1669 |
|
The mixing parameters {\bf diffkr} and {\bf kapgm} |
1670 |
|
are currently added as controls in {\it ctrl\_map\_ini.F}. |
1671 |
% |
% |
1672 |
\end{itemize} |
\end{itemize} |
1673 |
% |
% |
1674 |
|
|
1675 |
\subsubsection{Output of adjoint variables and gradient} |
\subsubsection{Output of adjoint variables and gradient} |
1676 |
% |
% |
1677 |
Two ways exist to generate output of adjoint fields. |
Several ways exist to generate output of adjoint fields. |
1678 |
% |
% |
1679 |
\begin{itemize} |
\begin{itemize} |
1680 |
% |
% |
1681 |
\item |
\item |
1682 |
\fbox{ |
\fbox{ |
1683 |
\begin{minipage}{12cm} |
\begin{minipage}{12cm} |
1684 |
{\it ctrl\_pack}: |
{\it ctrl\_map\_ini, ctrl\_map\_forcing}: |
1685 |
\end{minipage} |
\end{minipage} |
1686 |
} |
} |
1687 |
\\ |
\\ |
|
At the end of the forward/adjoint integration, the S/R |
|
|
{\it ctrl\_pack} is called which mirrors S/R {\it ctrl\_unpack}. |
|
|
It writes the following files: |
|
|
% |
|
1688 |
\begin{itemize} |
\begin{itemize} |
1689 |
% |
% |
1690 |
\item {\bf xx\_...}: the control variable fields |
\item {\bf xx\_...}: the control variable fields \\ |
1691 |
|
Before the forward integration, the control |
1692 |
|
variables are read from file {\bf xx\_ ...} and added to |
1693 |
|
the model field. |
1694 |
% |
% |
1695 |
\item {\bf adxx\_...}: the adjoint variable fields, i.e. the gradient |
\item {\bf adxx\_...}: the adjoint variable fields, i.e. the gradient |
1696 |
$ \nabla _{u}{\cal J} $ for each control variable, |
$ \nabla _{u}{\cal J} $ for each control variable \\ |
1697 |
|
After the adjoint integration the corresponding adjoint |
1698 |
|
variables are written to {\bf adxx\_ ...}. |
1699 |
|
% |
1700 |
|
\end{itemize} |
1701 |
% |
% |
1702 |
\item {\bf vector\_ctrl}: the control vector |
\item |
1703 |
|
\fbox{ |
1704 |
|
\begin{minipage}{12cm} |
1705 |
|
{\it ctrl\_unpack, ctrl\_pack}: |
1706 |
|
\end{minipage} |
1707 |
|
} |
1708 |
|
\\ |
1709 |
|
% |
1710 |
|
\begin{itemize} |
1711 |
% |
% |
1712 |
\item {\bf vector\_grad}: the gradient vector |
\item {\bf vector\_ctrl}: the control vector \\ |
1713 |
|
At the very beginning of the model initialization, |
1714 |
|
the updated compressed control vector is read (or initialised) |
1715 |
|
and distributed to 2-dim. and 3-dim. control variable fields. |
1716 |
|
% |
1717 |
|
\item {\bf vector\_grad}: the gradient vector \\ |
1718 |
|
At the very end of the adjoint integration, |
1719 |
|
the 2-dim. and 3-dim. adjoint variables are read, |
1720 |
|
compressed to a single vector and written to file. |
1721 |
% |
% |
1722 |
\end{itemize} |
\end{itemize} |
1723 |
% |
% |
1729 |
} |
} |
1730 |
\\ |
\\ |
1731 |
In addition to writing the gradient at the end of the |
In addition to writing the gradient at the end of the |
1732 |
forward/adjoint integration, many more adjoint variables, |
forward/adjoint integration, many more adjoint variables |
1733 |
representing the Lagrange multipliers of the model state |
of the model state |
1734 |
w.r.t. the model state |
at intermediate times can be written using S/R |
|
at different times can be written using S/R |
|
1735 |
{\it addummy\_in\_stepping}. |
{\it addummy\_in\_stepping}. |
1736 |
This routine is part of the adjoint support package |
This routine is part of the adjoint support package |
1737 |
{\it pkg/autodiff} (cf.f. below). |
{\it pkg/autodiff} (cf.f. below). |
1745 |
Appropriate flow directives ({\it dummy\_in\_stepping.flow}) |
Appropriate flow directives ({\it dummy\_in\_stepping.flow}) |
1746 |
ensure that TAMC does not automatically |
ensure that TAMC does not automatically |
1747 |
generate {\it addummy\_in\_stepping} by trying to differentiate |
generate {\it addummy\_in\_stepping} by trying to differentiate |
1748 |
{\it dummy\_in\_stepping}, but rather takes the hand-written routine. |
{\it dummy\_in\_stepping}, but instead refers to |
1749 |
|
the hand-written routine. |
1750 |
|
|
1751 |
{\it dummy\_in\_stepping} is called in the forward code |
{\it dummy\_in\_stepping} is called in the forward code |
1752 |
at the beginning of each |
at the beginning of each |
1756 |
{\it addynamics}. |
{\it addynamics}. |
1757 |
|
|
1758 |
{\it addummy\_in\_stepping} includes the header files |
{\it addummy\_in\_stepping} includes the header files |
1759 |
{\it adffields.h, addynamics.h, adtr1.h}. |
{\it adcommon.h}. |
1760 |
These header files are also hand-written. They contain |
This header file is also hand-written. It contains |
1761 |
the common blocks {\bf /addynvars\_r/}, {\bf /addynvars\_cd/}, |
the common blocks |
1762 |
|
{\bf /addynvars\_r/}, {\bf /addynvars\_cd/}, |
1763 |
|
{\bf /addynvars\_diffkr/}, {\bf /addynvars\_kapgm/}, |
1764 |
{\bf /adtr1\_r/}, {\bf /adffields/}, |
{\bf /adtr1\_r/}, {\bf /adffields/}, |
1765 |
which have been extracted from the adjoint code to enable |
which have been extracted from the adjoint code to enable |
1766 |
access to the adjoint variables. |
access to the adjoint variables. |
1778 |
with the value of the cost function itself $ {\cal J}(u_{[k]}) $ |
with the value of the cost function itself $ {\cal J}(u_{[k]}) $ |
1779 |
at iteration step $ k $ serve |
at iteration step $ k $ serve |
1780 |
as input to a minimization routine (e.g. quasi-Newton method, |
as input to a minimization routine (e.g. quasi-Newton method, |
1781 |
conjugate gradient, ...) to compute an update in the |
conjugate gradient, ... \cite{gil-lem:89}) |
1782 |
|
to compute an update in the |
1783 |
control variable for iteration step $k+1$ |
control variable for iteration step $k+1$ |
1784 |
\[ |
\[ |
1785 |
u_{[k+1]} \, = \, u_{[0]} \, + \, \Delta u_{[k+1]} |
u_{[k+1]} \, = \, u_{[0]} \, + \, \Delta u_{[k+1]} |
1793 |
and the minimization routine. |
and the minimization routine. |
1794 |
|
|
1795 |
\begin{eqnarray*} |
\begin{eqnarray*} |
1796 |
\footnotesize |
\scriptsize |
1797 |
\begin{array}{ccccc} |
\begin{array}{ccccc} |
1798 |
u_{[0]} \,\, , \,\, \Delta u_{[k]} & ~ & ~ & ~ & ~ \\ |
u_{[0]} \,\, , \,\, \Delta u_{[k]} & ~ & ~ & ~ & ~ \\ |
1799 |
{\Big\downarrow} |
{\Big\downarrow} |
1810 |
{\cal J}_{[k]} = {\cal J} \left( M \left( u_{[k]} \right) \right)} \\ |
{\cal J}_{[k]} = {\cal J} \left( M \left( u_{[k]} \right) \right)} \\ |
1811 |
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
1812 |
\hline |
\hline |
1813 |
|
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
1814 |
|
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{{\Big\downarrow}} \\ |
1815 |
|
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
1816 |
\hline |
\hline |
1817 |
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
1818 |
\multicolumn{1}{|c}{ |
\multicolumn{1}{|c}{ |
1819 |
\nabla_u {\cal J}_{[k]} (\delta {\cal J}) = |
\nabla_u {\cal J}_{[k]} (\delta {\cal J}) = |
1820 |
T\!\!^{\ast} \cdot \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J})} & |
T^{\ast} \cdot \nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J})} & |
1821 |
\stackrel{\bf adjoint}{\mathbf \longleftarrow} & |
\stackrel{\bf adjoint}{\mathbf \longleftarrow} & |
1822 |
ad \, v_{[k]} (\delta {\cal J}) = |
ad \, v_{[k]} (\delta {\cal J}) = |
1823 |
\nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J}) & |
\nabla_v {\cal J} |_{v_{[k]}} (\delta {\cal J}) & |
1826 |
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
1827 |
\hline |
\hline |
1828 |
~ & ~ & ~ & ~ & ~ \\ |
~ & ~ & ~ & ~ & ~ \\ |
1829 |
~ & ~ & |
\hspace*{15ex}{\Bigg\downarrow} |
1830 |
{\cal J}_{[k]} \qquad {\Bigg\downarrow} \qquad \nabla_u {\cal J}_{[k]} |
\quad {\cal J}_{[k]}, \quad \nabla_u {\cal J}_{[k]} |
1831 |
& ~ & ~ \\ |
& ~ & ~ & ~ & ~ \\ |
1832 |
~ & ~ & ~ & ~ & ~ \\ |
~ & ~ & ~ & ~ & ~ \\ |
1833 |
\hline |
\hline |
1834 |
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
\multicolumn{1}{|c}{~} & ~ & ~ & ~ & \multicolumn{1}{c|}{~} \\ |
1856 |
|
|
1857 |
\vspace*{0.5cm} |
\vspace*{0.5cm} |
1858 |
|
|
1859 |
|
{\scriptsize |
1860 |
\begin{tabular}{ccccc} |
\begin{tabular}{ccccc} |
1861 |
{\bf vector\_ctrl\_$<$k$>$ } & ~ & ~ & ~ & ~ \\ |
{\bf vector\_ctrl\_$<$k$>$ } & ~ & ~ & ~ & ~ \\ |
1862 |
{\big\downarrow} & ~ & ~ & ~ & ~ \\ |
{\big\downarrow} & ~ & ~ & ~ & ~ \\ |
1867 |
\cline{3-3} |
\cline{3-3} |
1868 |
\multicolumn{1}{l}{\bf xx\_theta0...$<$k$>$} & ~ & |
\multicolumn{1}{l}{\bf xx\_theta0...$<$k$>$} & ~ & |
1869 |
\multicolumn{1}{|c|}{~} & ~ & ~ \\ |
\multicolumn{1}{|c|}{~} & ~ & ~ \\ |
1870 |
\multicolumn{1}{l}{\bf xx\_salt0...$<$k$>$} & $\longrightarrow$ & |
\multicolumn{1}{l}{\bf xx\_salt0...$<$k$>$} & |
1871 |
|
$\stackrel{\mbox{read}}{\longrightarrow}$ & |
1872 |
\multicolumn{1}{|c|}{forward integration} & ~ & ~ \\ |
\multicolumn{1}{|c|}{forward integration} & ~ & ~ \\ |
1873 |
\multicolumn{1}{l}{\bf \vdots} & ~ & \multicolumn{1}{|c|}{~} |
\multicolumn{1}{l}{\bf \vdots} & ~ & \multicolumn{1}{|c|}{~} |
1874 |
& ~ & ~ \\ |
& ~ & ~ \\ |
1875 |
\cline{3-3} |
\cline{3-3} |
1876 |
~ & ~ & ~ & ~ & ~ \\ |
~ & ~ & $\downarrow$ & ~ & ~ \\ |
1877 |
\cline{3-3} |
\cline{3-3} |
1878 |
~ & ~ & |
~ & ~ & |
1879 |
\multicolumn{1}{|c|}{~} & ~ & |
\multicolumn{1}{|c|}{~} & ~ & |
1880 |
\multicolumn{1}{l}{\bf adxx\_theta0...$<$k$>$} \\ |
\multicolumn{1}{l}{\bf adxx\_theta0...$<$k$>$} \\ |
1881 |
~ & ~ & \multicolumn{1}{|c|}{adjoint integration} & |
~ & ~ & \multicolumn{1}{|c|}{adjoint integration} & |
1882 |
$\longrightarrow$ & |
$\stackrel{\mbox{write}}{\longrightarrow}$ & |
1883 |
\multicolumn{1}{l}{\bf adxx\_salt0...$<$k$>$} \\ |
\multicolumn{1}{l}{\bf adxx\_salt0...$<$k$>$} \\ |
1884 |
~ & ~ & \multicolumn{1}{|c|}{~} |
~ & ~ & \multicolumn{1}{|c|}{~} |
1885 |
& ~ & \multicolumn{1}{l}{\bf \vdots} \\ |
& ~ & \multicolumn{1}{l}{\bf \vdots} \\ |
1891 |
~ & ~ & ~ & ~ & {\big\downarrow} \\ |
~ & ~ & ~ & ~ & {\big\downarrow} \\ |
1892 |
~ & ~ & ~ & ~ & {\bf vector\_grad\_$<$k$>$ } \\ |
~ & ~ & ~ & ~ & {\bf vector\_grad\_$<$k$>$ } \\ |
1893 |
\end{tabular} |
\end{tabular} |
1894 |
|
} |
1895 |
|
|
1896 |
\vspace*{0.5cm} |
\vspace*{0.5cm} |
1897 |
|
|
1898 |
|
|
1899 |
{\it ctrl\_unpack} reads in the updated control vector |
{\it ctrl\_unpack} reads the updated control vector |
1900 |
{\bf vector\_ctrl\_$<$k$>$}. |
{\bf vector\_ctrl\_$<$k$>$}. |
1901 |
It distributes the different control variables to |
It distributes the different control variables to |
1902 |
2-dim. and 3-dim. files {\it xx\_...$<$k$>$}. |
2-dim. and 3-dim. files {\it xx\_...$<$k$>$}. |
1903 |
During the forward integration the control variables |
At the start of the forward integration the control variables |
1904 |
are read from {\it xx\_...$<$k$>$}. |
are read from {\it xx\_...$<$k$>$} and added to the |
1905 |
Correspondingly, the adjoint fields are written |
field. |
1906 |
|
Correspondingly, at the end of the adjoint integration |
1907 |
|
the adjoint fields are written |
1908 |
to {\it adxx\_...$<$k$>$}, again via the active file routines. |
to {\it adxx\_...$<$k$>$}, again via the active file routines. |
1909 |
Finally, {\it ctrl\_pack} collects all adjoint field files |
Finally, {\it ctrl\_pack} collects all adjoint files |
1910 |
and writes them to the compressed vector file |
and writes them to the compressed vector file |
1911 |
{\bf vector\_grad\_$<$k$>$}. |
{\bf vector\_grad\_$<$k$>$}. |
1912 |
|
|
1914 |
|
|
1915 |
|
|
1916 |
|
|
1917 |
\subsection{Flow directives and adjoint support routines} |
\subsection{Flow directives and adjoint support routines \label{section_flowdir}} |
1918 |
|
|
1919 |
\subsection{Store directives and checkpointing} |
\subsection{Store directives and checkpointing \label{section_checkpointing}} |
1920 |
|
|
1921 |
\subsection{Gradient checks} |
\subsection{Gradient checks \label{section_grdchk}} |
1922 |
|
|
1923 |
\subsection{Second derivative generation via TAMC} |
\subsection{Second derivative generation via TAMC} |
1924 |
|
|